diff --git a/Shower/QTilde/Base/PartnerFinder.cc b/Shower/QTilde/Base/PartnerFinder.cc
--- a/Shower/QTilde/Base/PartnerFinder.cc
+++ b/Shower/QTilde/Base/PartnerFinder.cc
@@ -1,713 +1,710 @@
 // -*- C++ -*-
 //
 // PartnerFinder.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 PartnerFinder class.
 //
 
 #include "PartnerFinder.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
 #include "ThePEG/Repository/UseRandom.h" 
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Utilities/Debug.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 
 using namespace Herwig;
 
 DescribeClass<PartnerFinder,Interfaced>
 describePartnerFinder ("Herwig::PartnerFinder","HwShower.so");
 
 // some useful functions to avoid using #define
 namespace {
 
   // return bool if final-state particle
   inline bool FS(const tShowerParticlePtr a) {
     return a->isFinalState();
   }
 
   // return colour line pointer 
   inline Ptr<ThePEG::ColourLine>::transient_pointer  
   CL(const tShowerParticlePtr a, unsigned int index=0) {
     return a->colourInfo()->colourLines().empty() ? ThePEG::tColinePtr() :
       const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->colourLines()[index]);
   }
 
   // return colour line size
   inline size_t
   CLSIZE(const tShowerParticlePtr a) {
     return a->colourInfo()->colourLines().size();
   }
   
   inline Ptr<ThePEG::ColourLine>::transient_pointer
   ACL(const tShowerParticlePtr a, unsigned int index=0) {
     return a->colourInfo()->antiColourLines().empty() ?  ThePEG::tColinePtr() :
       const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->antiColourLines()[index]);
   }
 
   inline size_t
   ACLSIZE(const tShowerParticlePtr a) {
     return a->colourInfo()->antiColourLines().size();
   }
 }
 
 void PartnerFinder::persistentOutput(PersistentOStream & os) const {
   os << partnerMethod_ << QEDPartner_ << scaleChoice_;
 }
 
 void PartnerFinder::persistentInput(PersistentIStream & is, int) {
   is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_;
 }
 
 void PartnerFinder::Init() {
 
   static ClassDocumentation<PartnerFinder> documentation
     ("This class is responsible for finding the partners for each interaction types ",
      "and within the evolution scale range specified by the ShowerVariables ",
      "then to determine the initial evolution scales for each pair of partners.");
 
   static Switch<PartnerFinder,int> interfacePartnerMethod
     ("PartnerMethod",
      "Choice of partner finding method for gluon evolution.",
      &PartnerFinder::partnerMethod_, 0, false, false);
   static SwitchOption interfacePartnerMethodRandom
     (interfacePartnerMethod,
      "Random",
      "Choose partners of a gluon randomly.",
      0);
   static SwitchOption interfacePartnerMethodMaximum
     (interfacePartnerMethod,
      "Maximum",
      "Choose partner of gluon with largest angle.",
      1);
 
   static Switch<PartnerFinder,int> interfaceQEDPartner
     ("QEDPartner",
      "Control of which particles to use as the partner for QED radiation",
      &PartnerFinder::QEDPartner_, 0, false, false);
   static SwitchOption interfaceQEDPartnerAll
     (interfaceQEDPartner,
      "All",
      "Consider all possible choices which give a positive contribution"
      " in the soft limit.",
      0);
   static SwitchOption interfaceQEDPartnerIIandFF
     (interfaceQEDPartner,
      "IIandFF",
      "Only allow initial-initial or final-final combinations",
      1);
   static SwitchOption interfaceQEDPartnerIF
     (interfaceQEDPartner,
      "IF",
      "Only allow initial-final combinations",
      2);
 
   static Switch<PartnerFinder,int> interfaceScaleChoice
     ("ScaleChoice",
      "The choice of the evolution scales",
      &PartnerFinder::scaleChoice_, 0, false, false);
   static SwitchOption interfaceScaleChoicePartner
     (interfaceScaleChoice,
      "Partner",
      "Scale of all interactions is that of the evolution partner",
      0);
   static SwitchOption interfaceScaleChoiceDifferent
     (interfaceScaleChoice,
      "Different",
      "Allow each interaction to have different scales",
      1);
 }
 
 void PartnerFinder::setInitialEvolutionScales(const ShowerParticleVector &particles,
 					      const bool isDecayCase,
 					      ShowerInteraction type,
 					      const bool setPartners) {
   // clear the existing partners
   for(ShowerParticleVector::const_iterator cit = particles.begin();
       cit != particles.end(); ++cit) (*cit)->clearPartners();
   // set them
   if(type==ShowerInteraction::QCD) {
     setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
   }
   else if(type==ShowerInteraction::QED) {
     setInitialQEDEvolutionScales(particles,isDecayCase,setPartners);
   }
   else if(type==ShowerInteraction::EW) {
     setInitialEWEvolutionScales(particles,isDecayCase,false);
   }
   else if(type==ShowerInteraction::QEDQCD) {
     setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
     setInitialQEDEvolutionScales(particles,isDecayCase,false);
   }
   else if(type==ShowerInteraction::ALL) {
     setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
     setInitialQEDEvolutionScales(particles,isDecayCase,false);
     setInitialEWEvolutionScales(particles,isDecayCase,false);
   }
   else
     assert(false);
   // \todo EW scales here
   // print out for debugging
   if(Debug::level>=10) {
     for(ShowerParticleVector::const_iterator cit = particles.begin();
 	cit != particles.end(); ++cit) {
       generator()->log() << "Particle: " << **cit << "\n";
       if(!(**cit).partner()) continue;
       generator()->log() << "Primary partner: " << *(**cit).partner() << "\n";
       for(vector<ShowerParticle::EvolutionPartner>::const_iterator it= (**cit).partners().begin();
 	  it!=(**cit).partners().end();++it) {
 	generator()->log() << static_cast<long>(it->type) << " "
 			   << it->weight << " " 
 			   << it->scale/GeV << " " 
  			   << *(it->partner) 
 			   << "\n";
       }
     }
     generator()->log() << flush;
   }
 }
 
 void PartnerFinder::setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
                                                  const bool isDecayCase,
                                                  const bool setPartners) {
   // Loop over  particles  and consider only coloured particles which don't
   // have already their colour partner fixed and that don't have children
   // (the latter requirement is relaxed in the case isDecayCase is true). 
   // Build a map which has as key one of these particles (i.e. a pointer 
   // to a ShowerParticle object) and as a corresponding value the vector
   // of all its possible *normal* candidate colour partners, defined as follows:
   //  --- have colour, and no children (this is not required in the case 
   //      isDecayCase is true);
   //  --- if both are initial (incoming) state particles, then the (non-null) colourLine() 
   //      of one of them must match the (non-null) antiColourLine() of the other.
   //  --- if one is an initial (incoming) state particle and the other is
   //      a final (outgoing) state particle,  then both must have the
   //      same (non-null) colourLine() or the same (non-null) antiColourLine();
   // Notice that this definition exclude the special case of baryon-violating
   // processes (as in R-parity Susy), which will show up as particles
   // without candidate colour partners, and that we will be treated a part later
   // (this means that no modifications of the following loop is needed!
 
   for ( const auto & sp : particles ) {
     // Skip colourless particles
     if(!sp->data().coloured()) continue;
     // find the partners
     auto partners = findQCDPartners(sp,particles);
     // must have a partner
     if(partners.empty()) {
       throw Exception() << "`Failed to make colour connections in " 
 			<< "PartnerFinder::setQCDInitialEvolutionScales"
 			<< *sp
 			<< Exception::eventerror;
     }
     // Calculate the evolution scales for all possible pairs of of particles
     vector<pair<Energy,Energy> > scales;
     int position = -1;
     for(size_t ix=0; ix< partners.size(); ++ix) {
       scales.push_back(calculateInitialEvolutionScales(ShowerPPair(sp, partners[ix].second),isDecayCase));
       if (!setPartners && partners[ix].second) position = ix; 
     }
     assert(setPartners || position >= 0);
     
     // set partners if required
     if (setPartners) {
       // In the case of more than one candidate colour partners,
       //	       there are now two approaches to choosing the partner. The
       //	       first method is based on two assumptions:
       //               1) the choice of which is the colour partner is done
       //                  *randomly* between the available candidates.
       //               2) the choice of which is the colour partner is done
       //                  *independently* from each particle: in other words,
       //                  if for a particle "i" its selected colour partner is  
       //                  the particle "j", then the colour partner of "j" 
       //                  does not have to be necessarily "i".
       // 	      The second method always chooses the furthest partner
       //	      for hard gluons and gluinos.
       // random choice
       if( partnerMethod_ == 0 ) {
 	// random choice of partner
 	position = UseRandom::irnd(partners.size());
       }
       // take the one with largest angle
       else if (partnerMethod_ == 1 ) {
 	if (sp->perturbative() == 1 && 
 	    sp->dataPtr()->iColour()==PDT::Colour8 ) {
 	  assert(partners.size()==2);
 	  // Determine largest angle
 	  double maxAngle(0.);
 	  for(unsigned int ix=0;ix<partners.size();++ix) {
 	    double angle = sp->momentum().vect().
 	      angle(partners[ix].second->momentum().vect());
 	    if(angle>maxAngle) {
 	      maxAngle = angle;
 	      position = ix;
 	    }
 	  }
 	}
 	else position = UseRandom::irnd(partners.size());
       }
       else assert(false);
       // set the evolution partner
       sp->partner(partners[position].second);
     }
     
     // primary partner set, set the others and do the scale
     for(size_t ix=0; ix<partners.size(); ++ix) {
       sp->addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,1.,partners[ix].first,
 						      scales[ix].first));
     }
     
     // set scales for all interactions to that of the partner, default
     Energy scale = scales[position].first;
     for(unsigned int ix=0;ix<partners.size();++ix) {
       if(partners[ix].first==ShowerPartnerType::QCDColourLine) {
 	sp->scales().QCD_c = 
 	  sp->scales().QCD_c_noAO = 
 	  (scaleChoice_==0 ? scale : scales[ix].first);
       }
       else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) {
 	sp->scales().QCD_ac = 
 	  sp->scales().QCD_ac_noAO =
 	  (scaleChoice_==0 ? scale : scales[ix].first);
       }
       else assert(false);
     }
   }
 }
 
 void PartnerFinder::setInitialQEDEvolutionScales(const ShowerParticleVector &particles,
 						 const bool isDecayCase,
 						 const bool setPartners) {
   // loop over all the particles
   for(const auto & sp : particles) {
     // not charged or photon continue
     if(!sp->dataPtr()->charged()) continue;
     // find the potential partners
     vector<pair<double,tShowerParticlePtr> > partners = findQEDPartners(sp,particles,isDecayCase);
     if(partners.empty()) {
       throw Exception() << "Failed to find partner in "
 			<< "PartnerFinder::setQEDInitialEvolutionScales"
 			<< *sp << Exception::eventerror;
     }
     // calculate the probabilities
     double prob(0.);
     for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
     // normalise
     for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
     // set the partner if required
     int position(-1);
     // use QCD partner if set
     if(!setPartners&&sp->partner()) {
       for(unsigned int ix=0;ix<partners.size();++ix) {
 	if(sp->partner()==partners[ix].second) {
 	  position = ix;
 	  break;
 	}
       }
     }
     // set the partner
     if(setPartners||!sp->partner()||position<0) {
       prob = UseRandom::rnd();
       for(unsigned int ix=0;ix<partners.size();++ix) {
  	if(partners[ix].first>prob) {
 	  position = ix;
 	  break;
 	}
 	prob -= partners[ix].first;
       }
       if(position>=0&&(setPartners||!sp->partner())) {
 	sp->partner(partners[position].second);
       }
     }
     // must have a partner
     if(position<0) throw Exception() << "Failed to find partner in "
 				 << "PartnerFinder::setQEDInitialEvolutionScales"
 				 << *sp << Exception::eventerror; 
     // Calculate the evolution scales for all possible pairs of of particles
     vector<pair<Energy,Energy> > scales;
     for(unsigned int ix=0;ix< partners.size();++ix) {
       scales.push_back(calculateInitialEvolutionScales(ShowerPPair(sp,partners[ix].second),
 						       isDecayCase));
     }
     // store all the possible partners
     for(unsigned int ix=0;ix<partners.size();++ix) {
       sp->addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
 							  partners[ix].first,
 							  ShowerPartnerType::QED,
 							  scales[ix].first));
     }
     // set scales
     sp->scales().QED      = scales[position].first;
     sp->scales().QED_noAO = scales[position].first;
   }
 }
 
 pair<Energy,Energy> PartnerFinder::
 calculateInitialEvolutionScales(const ShowerPPair &particlePair,
                                 const bool isDecayCase, int key) {
   bool FS1=FS(particlePair.first),FS2= FS(particlePair.second);
   if(FS1 && FS2){
     return calculateFinalFinalScales(particlePair.first->momentum(),particlePair.second->momentum(), key);
   }
   else if(FS1 && !FS2) {
     pair<Energy,Energy> rval = calculateInitialFinalScales(particlePair.second->momentum(),
                                                            particlePair.first->momentum(),
                                                            isDecayCase);
     return { rval.second, rval.first };
   }
   else if(!FS1 &&FS2)
     return calculateInitialFinalScales(particlePair.first->momentum(),particlePair.second->momentum(),isDecayCase);
   else
     return calculateInitialInitialScales(particlePair.first->momentum(),particlePair.second->momentum());
 }
 
 vector< pair<ShowerPartnerType, tShowerParticlePtr> > 
 PartnerFinder::findQCDPartners(tShowerParticlePtr particle,
 			       const ShowerParticleVector &particles) {
   vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners;
   for(const auto & sp : particles) {
     if(!sp->data().coloured() || particle==sp) continue;
     // one initial-state and one final-state particle
     if(FS(particle) != FS(sp)) {
       // loop over all the colours of both particles
       for(size_t ix=0; ix<CLSIZE(particle); ++ix) {
 	for(size_t jx=0; jx<CLSIZE(sp); ++jx) {
 	  if((CL(particle,ix) && CL(particle,ix)==CL(sp,jx))) {
 	    partners.push_back({ ShowerPartnerType::    QCDColourLine, sp });
 	  }
 	}
       }
       //loop over all the anti-colours of both particles
       for(size_t ix=0; ix<ACLSIZE(particle); ++ix) {
 	for(size_t jx=0; jx<ACLSIZE(sp); ++jx) {
 	  if((ACL(particle,ix) && ACL(particle,ix)==ACL(sp,jx))) {
 	    partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
 	  }
 	}
       }
     }
     // two initial-state or two final-state particles
     else {
       //loop over the colours of the first particle and the anti-colours of the other
       for(size_t ix=0; ix<CLSIZE(particle); ++ix){
 	for(size_t jx=0; jx<ACLSIZE(sp); ++jx){
 	  if(CL(particle,ix) && CL(particle,ix)==ACL(sp,jx)) {
 	    partners.push_back({ ShowerPartnerType::    QCDColourLine, sp });
 	  }
 	}
       }
       //loop over the anti-colours of the first particle and the colours of the other
       for(size_t ix=0; ix<ACLSIZE(particle); ++ix){
 	for(size_t jx=0; jx<CLSIZE(sp); jx++){
 	  if(ACL(particle,ix) && ACL(particle,ix)==CL(sp,jx)) {
 	    partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
 	  }
 	}
       }
     }
   }
   // if we haven't found any partners look for RPV
   if (partners.empty()) {
     // special for RPV
     tColinePtr col = CL(particle); 
     if(FS(particle)&&col&&col->sourceNeighbours().first) {
       tColinePair cpair = col->sourceNeighbours();
       for(const auto & sp : particles) {
 	if(( FS(sp) && ( CL(sp) == cpair.first || CL(sp)  == cpair.second))||
 	   (!FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second ))) {
 	  partners.push_back({ ShowerPartnerType::    QCDColourLine, sp });
 	}
       }
     }
     else if(col&&col->sinkNeighbours().first) {
       tColinePair cpair = col->sinkNeighbours();
       for(const auto & sp : particles) {
 	if(( FS(sp) && (ACL(sp) == cpair.first || ACL(sp)  == cpair.second))||
 	   (!FS(sp) && ( CL(sp) == cpair.first ||  CL(sp) == cpair.second))) {
 	  partners.push_back({ ShowerPartnerType::    QCDColourLine, sp });
 	}
       }
     }
     col = ACL(particle);
     if(FS(particle)&&col&&col->sinkNeighbours().first) {
       tColinePair cpair = col->sinkNeighbours();
       for(const auto & sp : particles) {
 	if(( FS(sp) && (ACL(sp) == cpair.first || ACL(sp)  == cpair.second))||
 	   (!FS(sp) && ( CL(sp) == cpair.first ||  CL(sp) == cpair.second ))) {
 	  partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
 	}
       }
     }
     else if(col&&col->sourceNeighbours().first) {
       tColinePair cpair = col->sourceNeighbours();
       for(const auto & sp : particles) {
 	if(( FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))||
 	   (!FS(sp) && (ACL(sp) == cpair.first ||ACL(sp) == cpair.second))) {
 	  partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp });
 	}
       }
     }
   }
   // return the partners
   return partners;
 }
 
 vector< pair<double, tShowerParticlePtr> >
 PartnerFinder::findQEDPartners(tShowerParticlePtr particle,
 			       const ShowerParticleVector & particles,
 			       const bool isDecayCase) {
   vector< pair<double, tShowerParticlePtr> > partners;
   const double pcharge = 
     particle->id()==ParticleID::gamma ? 1 : double(particle->data().iCharge());
   vector< pair<double, tShowerParticlePtr> > photons;
   for (const auto & sp : particles) {
     if (particle == sp) continue;
     if (sp->id()==ParticleID::gamma) photons.push_back(make_pair(1.,sp));
     if (!sp->data().charged() ) continue;
     double charge = pcharge*double((sp)->data().iCharge());
     if ( FS(particle) != FS(sp) ) charge *=-1.;
     if ( QEDPartner_ != 0 && !isDecayCase ) {
       // only include II and FF as requested
       if ( QEDPartner_ == 1 && FS(particle) != FS(sp) )
 	continue;
       // only include IF is requested
       else if (QEDPartner_ == 2 && FS(particle) == FS(sp) )
 	continue;
     }
     if (particle->id()==ParticleID::gamma) charge = -abs(charge);
     // only keep positive dipoles
     if (charge<0.) partners.push_back({ -charge, sp });
   }
   if (particle->id()==ParticleID::gamma && partners.empty())
     return photons;
   return partners;
 }
 
 void PartnerFinder::setInitialEWEvolutionScales(const ShowerParticleVector &particles,
 						const bool isDecayCase,
 						const bool setPartners) {
   // loop over all the particles
   for(ShowerParticleVector::const_iterator cit = particles.begin();
       cit != particles.end(); ++cit) {
     // if not weakly interacting continue
     if( !weaklyInteracting( (**cit).dataPtr()))
       continue;
     // find the potential partners
     vector<pair<double,tShowerParticlePtr> > partners = findEWPartners(*cit,particles,isDecayCase);
     if(partners.empty()) {
       throw Exception() << "Failed to find partner in "
   			<< "PartnerFinder::setEWInitialEvolutionScales"
   			<< (**cit) << Exception::eventerror;
     }
     // calculate the probabilities
     double prob(0.);
     for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
     // normalise
     for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
     // set the partner if required
     int position(-1);
     // use QCD partner if set
     if(!setPartners&&(*cit)->partner()) {
       for(unsigned int ix=0;ix<partners.size();++ix) {
   	if((*cit)->partner()==partners[ix].second) {
   	  position = ix;
   	  break;
   	}
       }
     }
     // set the partner
     if(setPartners||!(*cit)->partner()||position<0) {
       prob = UseRandom::rnd();
       for(unsigned int ix=0;ix<partners.size();++ix) {
   	if(partners[ix].first>prob) {
   	  position = ix;
   	  break;
   	}
   	prob -= partners[ix].first;
       }
       if(position>=0&&(setPartners||!(*cit)->partner())) {
   	(*cit)->partner(partners[position].second);
       }
     }
     // must have a partner
     if(position<0) throw Exception() << "Failed to find partner in "
   				 << "PartnerFinder::setEWInitialEvolutionScales"
   				 << (**cit) << Exception::eventerror;
     // Calculate the evolution scales for all possible pairs of of particles
-    vector<pair<Energy,Energy> > scalesW,scalesZ;
+    vector<pair<Energy,Energy> > scales;
     for(unsigned int ix=0;ix< partners.size();++ix) {
-      scalesZ.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
+      scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
 							isDecayCase, 0));
-      scalesW.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
-							isDecayCase, 1));
     }
     // store all the possible partners
     for(unsigned int ix=0;ix<partners.size();++ix) {
       (**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
   							  partners[ix].first,
   							  ShowerPartnerType::EW,
-  							  scalesW[ix].first));
+  							  scales[ix].first));
     }
     // set scales
-    (**cit).scales().EW_Z    = scalesZ[position].first;
-    (**cit).scales().EW_W    = scalesW[position].first;
+    (**cit).scales().EW      = scales[position].first;
   }
 }
 
 vector< pair<double, tShowerParticlePtr> >
 PartnerFinder::findEWPartners(tShowerParticlePtr particle,
 			       const ShowerParticleVector &particles,
 			       const bool isDecayCase ) {
   vector< pair<double, tShowerParticlePtr> > partners;
   ShowerParticleVector::const_iterator cjt;
   for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
     if(!weaklyInteracting((*cjt)->dataPtr()) ||
        particle == *cjt) continue;
     if( QEDPartner_ != 0 && !isDecayCase ) {
       // only include II and FF as requested
       if( QEDPartner_ == 1 && FS(particle) != FS(*cjt) )
 	continue;
       // only include IF is requested
       else if(QEDPartner_ == 2 && FS(particle) == FS(*cjt) )
 	continue;
     }
     double charge = 1.;
     // only keep positive dipoles
     if(charge>0.) partners.push_back(make_pair(charge,*cjt));
   }
   return partners;
 }
 
 pair<Energy,Energy> 
 PartnerFinder::calculateFinalFinalScales(
         const Lorentz5Momentum & p1,
         const Lorentz5Momentum & p2, int key) 
 {
   static const double eps=1e-7;
   // Using JHEP 12(2003)045 we find that we need ktilde = 1/2(1+b-c+lambda)
   // ktilde = qtilde^2/Q^2 therefore qtilde = sqrt(ktilde*Q^2)
   // find momenta in rest frame of system
   // calculate quantities for the scales
   Energy2 Q2 = (p1+p2).m2();
   double b = p1.mass2()/Q2;
   double c = p2.mass2()/Q2;
   if(b<0.) {
     if(b<-eps) {
       throw Exception() << "Negative Mass squared b = " << b
 			<< "in PartnerFinder::calculateFinalFinalScales()"
 			<< Exception::eventerror;
     }
     b = 0.;
   }
   if(c<0.) {
     if(c<-eps) {
       throw Exception() << "Negative Mass squared c = " << c
 			<< "in PartnerFinder::calculateFinalFinalScales()"
 			<< Exception::eventerror;
     }
     c = 0.;
   }
   // KMH & PR - 16 May 2008 - swapped lambda calculation from 
   // double lam=2.*p1.vect().mag()/Q; to sqrt(kallen(1,b,c)), 
   // which should be identical for p1 & p2 onshell in their COM
   // but in the inverse construction for the Nason method, this
   // was not the case, leading to misuse. 
   const double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c))
                    *(sqrt(b)-1.-sqrt(c))*(sqrt(c)-1.-sqrt(b)));
   
   Energy firstQ, secondQ;
   double kappab(0.), kappac(0.);
   //key = 0; // symmetric case pre-selection 
   switch(key) {
   case 0: // symmetric case
     firstQ  = sqrt(0.5*Q2*(1.+b-c+lam));
     secondQ = sqrt(0.5*Q2*(1.-b+c+lam));
     break;
   case 1: // maximum emission from both legs
     kappab=4.*(1.-2.*sqrt(c)-b+c);
     kappac=4.*(1.-2.*sqrt(b)-c+b);
     firstQ  = sqrt(Q2*kappab);
     secondQ = sqrt(Q2*kappac);
     break;
   default:
     assert(false);
   }     
   // calculate the scales
   return pair<Energy,Energy>(firstQ, secondQ);
 }
 
 pair<Energy,Energy>
 PartnerFinder::calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
 					   const bool isDecayCase) {
   if(!isDecayCase) { 
     // In this case from JHEP 12(2003)045 we find the conditions
     // ktilde_b = (1+c) and ktilde_c = (1+2c)
     // We also find that c = m_c^2/Q^2. The process is a+b->c where
     // particle a is not colour connected (considered as a colour singlet).
     // Therefore we simply find that q_b = sqrt(Q^2+m_c^2) and 
     // q_c = sqrt(Q^2+2 m_c^2)
     // We also assume that the first particle in the pair is the initial
     // state particle and the second is the final state one c 
     const Energy2  mc2 = sqr(pc.mass());
     const Energy2  Q2  = -(pb-pc).m2();
     return { sqrt(Q2+mc2), sqrt(Q2+2*mc2) };
   }
   else {    
     // In this case from JHEP 12(2003)045 we find, for the decay
     // process b->c+a(neutral), the condition
     // (ktilde_b-1)*(ktilde_c-c)=(1/4)*sqr(1-a+c+lambda). 
     // We also assume that the first particle in the pair is the initial
     // state particle (b) and the second is the final state one (c).
     //  - We find maximal phase space coverage through emissions from 
     //    c if we set ktilde_c = 4.*(sqr(1.-sqrt(a))-c)
     //  - We find the most 'symmetric' way to populate the phase space
     //    occurs for (ktilde_b-1)=(ktilde_c-c)=(1/2)*(1-a+c+lambda) 
     //  - We find the most 'smooth' way to populate the phase space
     //    occurs for...
     Energy2 mb2(sqr(pb.mass()));
     double a=(pb-pc).m2()/mb2;
     double c=sqr(pc.mass())/mb2;
     double lambda   = 1. + a*a + c*c - 2.*a - 2.*c - 2.*a*c;
     lambda = sqrt(max(lambda,0.));
     const double PROD     = 0.25*sqr(1. - a + c + lambda);
     int key = 0;
     double ktilde_b, ktilde_c, cosi = 0.;
     switch (key) {
     case 0: // the 'symmetric' choice
       ktilde_c = 0.5*(1-a+c+lambda) + c ;
       ktilde_b = 1.+PROD/(ktilde_c-c)   ;
       break;
     case 1: // the 'maximal' choice
       ktilde_c = 4.0*(sqr(1.-sqrt(a))-c);
       ktilde_b = 1.+PROD/(ktilde_c-c)   ;
       break;
     case 2: // the 'smooth' choice
       c = max(c,1.*GeV2/mb2);
       cosi = (sqr(1-sqrt(c))-a)/lambda;
       ktilde_b = 2.0/(1.0-cosi);
       ktilde_c = (1.0-a+c+lambda)*(1.0+c-a-lambda*cosi)/(2.0*(1.0+cosi));
       break;
     }
     return { sqrt(mb2*ktilde_b), sqrt(mb2*ktilde_c) };
   }
 }
 
 pair<Energy,Energy>
 PartnerFinder::calculateInitialInitialScales(const Lorentz5Momentum& p1, const Lorentz5Momentum& p2) {
   // This case is quite simple. From JHEP 12(2003)045 we find the condition
   // that ktilde_b = ktilde_c = 1. In this case we have the process
   // b+c->a so we need merely boost to the CM frame of the two incoming
   // particles and then qtilde is equal to the energy in that frame
   const Energy Q = sqrt((p1+p2).m2());
   return {Q,Q};
 }
diff --git a/Shower/QTilde/Base/ShowerParticle.cc b/Shower/QTilde/Base/ShowerParticle.cc
--- a/Shower/QTilde/Base/ShowerParticle.cc
+++ b/Shower/QTilde/Base/ShowerParticle.cc
@@ -1,585 +1,584 @@
 // -*- C++ -*-
 //
 // ShowerParticle.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 ShowerParticle class.
 //
 
 #include "ShowerParticle.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
 #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
 #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
 #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
 #include <ostream>
 
 using namespace Herwig;
 
 PPtr ShowerParticle::clone() const {
   return new_ptr(*this);
 }
 
 PPtr ShowerParticle::fullclone() const {
   return new_ptr(*this);
 }
 
 ClassDescription<ShowerParticle> ShowerParticle::initShowerParticle;
 // Definition of the static class description member.
 
 void ShowerParticle::vetoEmission(ShowerPartnerType, Energy scale) {
   scales_.QED         = min(scale,scales_.QED        );
   scales_.QED_noAO    = min(scale,scales_.QED_noAO   );
   scales_.QCD_c       = min(scale,scales_.QCD_c      );
   scales_.QCD_c_noAO  = min(scale,scales_.QCD_c_noAO );
   scales_.QCD_ac      = min(scale,scales_.QCD_ac     );
   scales_.QCD_ac_noAO = min(scale,scales_.QCD_ac_noAO);
-  scales_.EW_Z        = min(scale,scales_.EW_Z       );
-  scales_.EW_W        = min(scale,scales_.EW_W       );
+  scales_.EW          = min(scale,scales_.EW         );
   if(spinInfo()) spinInfo()->undecay();
 }
 
 void ShowerParticle::addPartner(EvolutionPartner in ) {
   partners_.push_back(in); 
 }
 
 namespace {
 
 LorentzRotation boostToShower(Lorentz5Momentum & porig, tShowerBasisPtr basis) {
   LorentzRotation output; 
   assert(basis);
   if(basis->frame()==ShowerBasis::BackToBack) {
     // we are doing the evolution in the back-to-back frame
     // work out the boostvector
     Boost boostv(-(basis->pVector()+basis->nVector()).boostVector());
     // momentum of the parton
     Lorentz5Momentum ptest(basis->pVector());
     // construct the Lorentz boost
     output = LorentzRotation(boostv);
     ptest *= output;
     Axis axis(ptest.vect().unit());
     // now rotate so along the z axis as needed for the splitting functions
     if(axis.perp2()>1e-10) {
       double sinth(sqrt(1.-sqr(axis.z())));
       output.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
     }
     else if(axis.z()<0.) {
       output.rotate(Constants::pi,Axis(1.,0.,0.));
     }
     porig = output*basis->pVector();
     porig.setX(ZERO);
     porig.setY(ZERO);
   }
   else {
     output = LorentzRotation(-basis->pVector().boostVector());
     porig = output*basis->pVector();
     porig.setX(ZERO);
     porig.setY(ZERO);
     porig.setZ(ZERO);
   }
   return output;
 }
 
 RhoDMatrix bosonMapping(ShowerParticle & particle,
 			const Lorentz5Momentum & porig,
 			VectorSpinPtr vspin,
 			const LorentzRotation & rot,
 			Helicity::Direction dir) {
   // rotate the original basis
   vector<LorentzPolarizationVector> sbasis;
   for(unsigned int ix=0;ix<3;++ix) {
     sbasis.push_back(vspin->getProductionBasisState(ix));
     sbasis.back().transform(rot);
   }
   // splitting basis
   vector<LorentzPolarizationVector> fbasis;
   bool massless(particle.id()==ParticleID::g||particle.id()==ParticleID::gamma);
   VectorWaveFunction wave(porig,particle.dataPtr(),dir==outgoing ? incoming : outgoing);
   for(unsigned int ix=0;ix<3;++ix) {
     if(massless&&ix==1) {
       fbasis.push_back(LorentzPolarizationVector());
     }
     else {
       wave.reset(ix,vector_phase);
       fbasis.push_back(wave.wave());
     }
   }
   // work out the mapping
   RhoDMatrix mapping=RhoDMatrix(PDT::Spin1,false);
   for(unsigned int ix=0;ix<3;++ix) {
     for(unsigned int iy=0;iy<3;++iy) {
       mapping(ix,iy)= -conj(sbasis[ix].dot(fbasis[iy]));
       if(particle.id()<0) mapping(ix,iy)=conj(mapping(ix,iy));
     }
   }
   return mapping;
 }
 
 RhoDMatrix fermionMapping(ShowerParticle & particle,
 			  const Lorentz5Momentum & porig,
 			  FermionSpinPtr fspin,
 			  const LorentzRotation & rot) {
   // extract the original basis states
   vector<LorentzSpinor<SqrtEnergy> > sbasis;
   for(unsigned int ix=0;ix<2;++ix) {
     sbasis.push_back(fspin->getProductionBasisState(ix));
     sbasis.back().transform(rot);
   }
   // calculate the states in the splitting basis
   vector<LorentzSpinor<SqrtEnergy> > fbasis;
   SpinorWaveFunction wave(porig,particle.dataPtr(),
 			  particle.id()>0 ? incoming : outgoing);
   for(unsigned int ix=0;ix<2;++ix) {
     wave.reset(ix);
     fbasis.push_back(wave.dimensionedWave());
   }
   RhoDMatrix mapping=RhoDMatrix(PDT::Spin1Half,false);
   for(unsigned int ix=0;ix<2;++ix) {
     if(fbasis[0].s2()==complex<SqrtEnergy>()) {
       mapping(ix,0) = sbasis[ix].s3()/fbasis[0].s3();
       mapping(ix,1) = sbasis[ix].s2()/fbasis[1].s2();
     }
     else {
       mapping(ix,0) = sbasis[ix].s2()/fbasis[0].s2();
       mapping(ix,1) = sbasis[ix].s3()/fbasis[1].s3();
     }
   }
   return mapping;
 }
 
 FermionSpinPtr createFermionSpinInfo(ShowerParticle & particle,
 				     const Lorentz5Momentum & porig,
 				     const LorentzRotation & rot,
 				     Helicity::Direction dir) {
   // calculate the splitting basis for the branching
   // and rotate back to construct the basis states
   LorentzRotation rinv = rot.inverse();
   SpinorWaveFunction wave;
   if(particle.id()>0)
     wave=SpinorWaveFunction(porig,particle.dataPtr(),incoming);
   else
     wave=SpinorWaveFunction(porig,particle.dataPtr(),outgoing);
   FermionSpinPtr fspin = new_ptr(FermionSpinInfo(particle.momentum(),dir==outgoing));
   for(unsigned int ix=0;ix<2;++ix) {
     wave.reset(ix);
     LorentzSpinor<SqrtEnergy> basis = wave.dimensionedWave();
     basis.transform(rinv);
     fspin->setBasisState(ix,basis);
   }
   particle.spinInfo(fspin);
   return fspin;
 }
 
 VectorSpinPtr createVectorSpinInfo(ShowerParticle & particle,
 				   const Lorentz5Momentum & porig,
 				   const LorentzRotation & rot,
 				   Helicity::Direction dir) {
   // calculate the splitting basis for the branching
   // and rotate back to construct the basis states
   LorentzRotation rinv = rot.inverse();
   bool massless(particle.id()==ParticleID::g||particle.id()==ParticleID::gamma);
   VectorWaveFunction wave(porig,particle.dataPtr(),dir);
   VectorSpinPtr vspin = new_ptr(VectorSpinInfo(particle.momentum(),dir==outgoing));
   for(unsigned int ix=0;ix<3;++ix) {
     LorentzPolarizationVector basis;
     if(massless&&ix==1) {
       basis = LorentzPolarizationVector();
     }
     else {
       wave.reset(ix,vector_phase);
       basis = wave.wave();
     }
     basis *= rinv;
     vspin->setBasisState(ix,basis);
   }
   particle.spinInfo(vspin);
   vspin->  DMatrix() = RhoDMatrix(PDT::Spin1);
   vspin->rhoMatrix() = RhoDMatrix(PDT::Spin1);
   if(massless) {
     vspin->  DMatrix()(0,0) = 0.5;
     vspin->rhoMatrix()(0,0) = 0.5;
     vspin->  DMatrix()(2,2) = 0.5;
     vspin->rhoMatrix()(2,2) = 0.5;
   }
   return vspin;
 }
 
 }
 
 RhoDMatrix ShowerParticle::extractRhoMatrix(bool forward) {
   // get the spin density matrix and the mapping
   RhoDMatrix mapping;
   SpinPtr inspin;
   bool needMapping = getMapping(inspin,mapping);
   // set the decayed flag
   inspin->decay();
   // get the spin density matrix
   RhoDMatrix rho = forward ? inspin->rhoMatrix() : inspin->DMatrix();
   // map to the shower basis if needed  
   if(needMapping) {
     RhoDMatrix rhop(rho.iSpin(),false);
     for(int ixb=0;ixb<rho.iSpin();++ixb) {
   	for(int iyb=0;iyb<rho.iSpin();++iyb) {
   	  if(mapping(iyb,ixb)==0.)continue;
 	  for(int iya=0;iya<rho.iSpin();++iya) {
             if(rho(iya,iyb)==0.)continue;
 	    for(int ixa=0;ixa<rho.iSpin();++ixa) {
   	      rhop(ixa,ixb) += rho(iya,iyb)*mapping(iya,ixa)*conj(mapping(iyb,ixb));
   	    }
   	  }
   	}
     }
     rhop.normalize();
     rho = rhop;
   }
   return rho;
 }
 
 bool ShowerParticle::getMapping(SpinPtr & output, RhoDMatrix & mapping) {
   // if the particle is not from the hard process
   if(!this->perturbative()) {
     // mapping is the identity
     output=this->spinInfo();
     mapping=RhoDMatrix(this->dataPtr()->iSpin());
     if(output) {
       return false;
     }
     else {
       Lorentz5Momentum porig;
       LorentzRotation rot = boostToShower(porig,showerBasis());
       Helicity::Direction dir = this->isFinalState() ? outgoing : incoming;
       if(this->dataPtr()->iSpin()==PDT::Spin0) {
 	assert(false);
       }
       else if(this->dataPtr()->iSpin()==PDT::Spin1Half) {
 	output = createFermionSpinInfo(*this,porig,rot,dir);
       }
       else if(this->dataPtr()->iSpin()==PDT::Spin1) {
 	output = createVectorSpinInfo(*this,porig,rot,dir);
       }
       else {
 	assert(false);
       }
       return false;
     }
   }
   // if particle is final-state and is from the hard process
   else if(this->isFinalState()) {
     assert(this->perturbative()==1 || this->perturbative()==2);
     // get transform to shower frame
     Lorentz5Momentum porig;
     LorentzRotation rot = boostToShower(porig,showerBasis());
     // the rest depends on the spin of the particle
     PDT::Spin spin(this->dataPtr()->iSpin());
     mapping=RhoDMatrix(spin,false);
     // do the spin dependent bit
     if(spin==PDT::Spin0) {
       ScalarSpinPtr sspin=dynamic_ptr_cast<ScalarSpinPtr>(this->spinInfo());
       if(!sspin) {
 	ScalarWaveFunction::constructSpinInfo(this,outgoing,true);
       }
       output=this->spinInfo();
       return false;
     }
     else if(spin==PDT::Spin1Half) {
       FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
       // spin info exists get information from it
       if(fspin) {
 	output=fspin;
 	mapping = fermionMapping(*this,porig,fspin,rot);
 	return true;
       }
       // spin info does not exist create it
       else {
 	output = createFermionSpinInfo(*this,porig,rot,outgoing);
 	return false;
       }
     }
     else if(spin==PDT::Spin1) {
       VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
       // spin info exists get information from it
       if(vspin) {
 	output=vspin;
 	mapping = bosonMapping(*this,porig,vspin,rot,outgoing);
 	return true; 
       }
       else {
 	output = createVectorSpinInfo(*this,porig,rot,outgoing);
 	return false;
       }
     }
     // not scalar/fermion/vector
     else
       assert(false);
   }
   // incoming to hard process
   else if(this->perturbative()==1 && !this->isFinalState()) {
     // get the basis vectors
     // get transform to shower frame
     Lorentz5Momentum porig;
     LorentzRotation rot = boostToShower(porig,showerBasis());
     porig *= this->x();
     // the rest depends on the spin of the particle
     PDT::Spin spin(this->dataPtr()->iSpin());
     mapping=RhoDMatrix(spin);
     // do the spin dependent bit
     if(spin==PDT::Spin0) {
       cerr << "testing spin 0 not yet implemented " << endl;
       assert(false);
     }
     // spin-1/2
     else if(spin==PDT::Spin1Half) {
       FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
       // spin info exists get information from it
       if(fspin) {
 	output=fspin;
 	mapping = fermionMapping(*this,porig,fspin,rot);
 	return true;
       }
       // spin info does not exist create it
       else {
 	output = createFermionSpinInfo(*this,porig,rot,incoming);
 	return false;
       }
     }
     // spin-1
     else if(spin==PDT::Spin1) {
       VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
       // spinInfo exists map it
       if(vspin) {
 	output=vspin;
 	mapping = bosonMapping(*this,porig,vspin,rot,incoming);
 	return true;
       }
       // create the spininfo
       else {
 	output = createVectorSpinInfo(*this,porig,rot,incoming);
 	return false;
       }
     }
     assert(false);
   }
   // incoming to decay
   else if(this->perturbative() == 2 && !this->isFinalState()) { 
     // get the basis vectors
     Lorentz5Momentum porig;
     LorentzRotation rot=boostToShower(porig,showerBasis());
     // the rest depends on the spin of the particle
     PDT::Spin spin(this->dataPtr()->iSpin());
     mapping=RhoDMatrix(spin);
     // do the spin dependent bit
     if(spin==PDT::Spin0) {
       cerr << "testing spin 0 not yet implemented " << endl;
       assert(false);
     }
     // spin-1/2
     else if(spin==PDT::Spin1Half) {
     //   FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
     //   // spin info exists get information from it
     //   if(fspin) {
     // 	output=fspin;
     // 	mapping = fermionMapping(*this,porig,fspin,rot);
     // 	return true;
     //   // spin info does not exist create it
     //   else {
     // 	output = createFermionSpinInfo(*this,porig,rot,incoming);
     // 	return false;
     //   }
     // }
       assert(false);
     }
     // // spin-1
     // else if(spin==PDT::Spin1) {
     //   VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
     //   // spinInfo exists map it
     //   if(vspin) {
     // 	output=vspin;
     // 	mapping = bosonMapping(*this,porig,vspin,rot);
     // 	return true;
     //   }
     //   // create the spininfo
     //   else {
     // 	output = createVectorSpinInfo(*this,porig,rot,incoming);
     // 	return false;
     //   }
     // }
     // assert(false);
     assert(false);
   }
   else
     assert(false);
   return true;
 }
 
 void ShowerParticle::constructSpinInfo(bool timeLike) {
   // now construct the required spininfo and calculate the basis states
   PDT::Spin spin(dataPtr()->iSpin());
   if(spin==PDT::Spin0) {
     ScalarWaveFunction::constructSpinInfo(this,outgoing,timeLike);
   }
   // calculate the basis states and construct the SpinInfo for a spin-1/2 particle
   else if(spin==PDT::Spin1Half) {
     // outgoing particle
     if(id()>0) {
       vector<LorentzSpinorBar<SqrtEnergy> > stemp;
       SpinorBarWaveFunction::calculateWaveFunctions(stemp,this,outgoing);
       SpinorBarWaveFunction::constructSpinInfo(stemp,this,outgoing,timeLike);
     }
     // outgoing antiparticle
     else {
       vector<LorentzSpinor<SqrtEnergy> > stemp;
       SpinorWaveFunction::calculateWaveFunctions(stemp,this,outgoing);
       SpinorWaveFunction::constructSpinInfo(stemp,this,outgoing,timeLike);
     }
   }
   // calculate the basis states and construct the SpinInfo for a spin-1 particle
   else if(spin==PDT::Spin1) {
     bool massless(id()==ParticleID::g||id()==ParticleID::gamma);
     vector<Helicity::LorentzPolarizationVector> vtemp;
     VectorWaveFunction::calculateWaveFunctions(vtemp,this,outgoing,massless,vector_phase);
     VectorWaveFunction::constructSpinInfo(vtemp,this,outgoing,timeLike,massless);
   }
   else {
     throw Exception() << "Spins higher than 1 are not yet implemented in " 
 		      << "FS_QtildaShowerKinematics1to2::constructVertex() "
 		      << Exception::runerror;
   }
 }
 
 void ShowerParticle::initializeDecay() {
   Lorentz5Momentum p, n, ppartner, pcm;
   assert(perturbative()!=1);
   // this is for the initial decaying particle
   if(perturbative()==2) {
     ShowerBasisPtr newBasis;
     p = momentum();
     Lorentz5Momentum ppartner(partner()->momentum());
     // removed to make inverse recon work properly
     //if(partner->thePEGBase()) ppartner=partner->thePEGBase()->momentum();
     pcm=ppartner;
     Boost boost(p.findBoostToCM());
     pcm.boost(boost);
     n = Lorentz5Momentum( ZERO,0.5*p.mass()*pcm.vect().unit()); 
     n.boost( -boost);
     newBasis = new_ptr(ShowerBasis());
     newBasis->setBasis(p,n,ShowerBasis::Rest);
     showerBasis(newBasis,false);
   }
   else {
     showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0])->showerBasis(),true);
   }
 }
 
 void ShowerParticle::initializeInitialState(PPtr parent) {
   // For the time being we are considering only 1->2 branching
   Lorentz5Momentum p, n, pthis, pcm;
   assert(perturbative()!=2);
   if(perturbative()==1) {
     ShowerBasisPtr newBasis;
     // find the partner and its momentum
     if(!partner()) return;
     if(partner()->isFinalState()) {
       Lorentz5Momentum pa = -momentum()+partner()->momentum();
       Lorentz5Momentum pb =  momentum();
       Energy scale=parent->momentum().t();
       Lorentz5Momentum pbasis(ZERO,parent->momentum().vect().unit()*scale);
       Axis axis(pa.vect().unit());
       LorentzRotation rot;
       double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
       if(axis.perp2()>1e-20) {
 	rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
 	rot.rotateX(Constants::pi);
       }
       if(abs(1.-pa.e()/pa.vect().mag())>1e-6) rot.boostZ( pa.e()/pa.vect().mag());
       pb *= rot;
       if(pb.perp2()/GeV2>1e-20) {
 	Boost trans = -1./pb.e()*pb.vect();
 	trans.setZ(0.);
 	rot.boost(trans);
       }
       pbasis *=rot;
       rot.invert();
       n = rot*Lorentz5Momentum(ZERO,-pbasis.vect());
       p = rot*Lorentz5Momentum(ZERO, pbasis.vect());
     }
     else {
       pcm = parent->momentum();
       p = Lorentz5Momentum(ZERO, pcm.vect());
       n = Lorentz5Momentum(ZERO, -pcm.vect());
     }
     newBasis = new_ptr(ShowerBasis());
     newBasis->setBasis(p,n,ShowerBasis::BackToBack);
     showerBasis(newBasis,false);
   } 
   else {
     showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(children()[0])->showerBasis(),true);
   }
 }
 
 void ShowerParticle::initializeFinalState() {
   // set the basis vectors
   Lorentz5Momentum p,n;
   if(perturbative()!=0) {
     ShowerBasisPtr newBasis;
     // find the partner() and its momentum
     if(!partner()) return;
     Lorentz5Momentum ppartner(partner()->momentum());
     // momentum of the emitting particle
     p = momentum();
     Lorentz5Momentum pcm;
     // if the partner is a final-state particle then the reference
     // vector is along the partner in the rest frame of the pair
     if(partner()->isFinalState()) {
       Boost boost=(p + ppartner).findBoostToCM();
       pcm = ppartner;
       pcm.boost(boost);
       n = Lorentz5Momentum(ZERO,pcm.vect());
       n.boost( -boost);
     }
     else if(!partner()->isFinalState()) {
       // if the partner is an initial-state particle then the reference
       // vector is along the partner which should be massless
       if(perturbative()==1) {
 	n = Lorentz5Momentum(ZERO,ppartner.vect());
       }
       // if the partner is an initial-state decaying particle then the reference
       // vector is along the backwards direction in rest frame of decaying particle
       else {
 	Boost boost=ppartner.findBoostToCM();
 	pcm = p;
 	pcm.boost(boost);
 	n = Lorentz5Momentum( ZERO, -pcm.vect()); 
 	n.boost( -boost);
       } 
     }
     newBasis = new_ptr(ShowerBasis());
     newBasis->setBasis(p,n,ShowerBasis::BackToBack);
     showerBasis(newBasis,false);
   }
   else if(initiatesTLS()) {
     showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0]->children()[0])->showerBasis(),true);
   }
   else  {
     showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0])->showerBasis(),true);
   }
 }
 
 void ShowerParticle::setShowerMomentum(bool timeLike) {
   Energy m = this->mass() > ZERO ? this->mass() : this->data().mass();
   // calculate the momentum of the assuming on-shell
   Energy2 pt2 = sqr(this->showerParameters().pt);
   double alpha = timeLike ? this->showerParameters().alpha : this->x();
   double beta = 0.5*(sqr(m) + pt2 - sqr(alpha)*showerBasis()->pVector().m2())/(alpha*showerBasis()->p_dot_n());
   Lorentz5Momentum porig=showerBasis()->sudakov2Momentum(alpha,beta,
 							 this->showerParameters().ptx,
 							 this->showerParameters().pty);
   porig.setMass(m);
   this->set5Momentum(porig);
 }
diff --git a/Shower/QTilde/Base/ShowerParticle.h b/Shower/QTilde/Base/ShowerParticle.h
--- a/Shower/QTilde/Base/ShowerParticle.h
+++ b/Shower/QTilde/Base/ShowerParticle.h
@@ -1,531 +1,529 @@
 // -*- C++ -*-
 //
 // ShowerParticle.h 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.
 //
 #ifndef HERWIG_ShowerParticle_H
 #define HERWIG_ShowerParticle_H
 //
 // This is the declaration of the ShowerParticle class.
 //
 
 #include "ThePEG/EventRecord/Particle.h"
 #include "Herwig/Shower/QTilde/SplittingFunctions/SplittingFunction.fh"
 #include "Herwig/Shower/QTilde/ShowerConfig.h"
 #include "Herwig/Shower/QTilde/Kinematics/ShowerBasis.h"
 #include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.h"
 #include "ShowerParticle.fh"
 #include <iosfwd>
 
 namespace Herwig {
 
 using namespace ThePEG;
 
   /** \ingroup Shower
    *  This class represents a particle in the showering process.
    *  It inherits from the Particle class of ThePEG and has some  
    *  specifics information useful only during the showering process.
    * 
    *  Notice that: 
    *    - for forward evolution, it is clear what is meant by parent/child; 
    *      for backward evolution, however, it depends whether we want 
    *      to keep a physical picture or a Monte-Carlo effective one. 
    *      In the former case, an incoming particle (emitting particle)  
    *      splits into an emitted particle and the emitting particle after 
    *      the emission: the latter two are then children of the 
    *      emitting particle, the parent. In the Monte-Carlo effective  
    *      picture, we have that the particle close to the hard subprocess, 
    *      with higher (space-like) virtuality, splits into an emitted particle 
    *      and the emitting particle at lower virtuality: the latter two are, 
    *      in this case, the children of the first one, the parent. However we
    *      choose a more physical picture where the new emitting particle is the
    *      parented of the emitted final-state particle and the original emitting
    *      particle.
    *    - the pointer to a SplitFun object is set only in the case 
    *      that the particle has undergone a shower emission. This is similar to
    *      the case of the decay of a normal Particle where 
    *      the pointer to a Decayer object is set only in the case 
    *      that the particle has undergone to a decay. 
    *      In the case of particle connected directly to the hard subprocess, 
    *      there is no pointer to the hard subprocess, but there is a method 
    *      isFromHardSubprocess() which returns true only in this case.
    *
    *  @see Particle
    *  @see ShowerConfig
    *  @see ShowerKinematics
    */
 class ShowerParticle: public Particle {
 
 public:
 
   /**
    *  Struct for all the info on an evolution partner
    */
   struct EvolutionPartner {
 
     /**
      *  Constructor
      */
     EvolutionPartner(tShowerParticlePtr p,double w, ShowerPartnerType t,
 		     Energy s) : partner(p), weight(w), type(t), scale(s)
     {}
 
     /**
      * The partner
      */
     tShowerParticlePtr partner;
     
     /**
      *  Weight
      */
     double weight;
 
     /**
      *  Type
      */
     ShowerPartnerType type;
 
     /**
      *  The assoicated evolution scale
      */
     Energy scale;
   };
 
   /**
    *  Struct to store the evolution scales
    */
   struct EvolutionScales {
 
     /**
      *  Constructor
      */
     EvolutionScales() : QED(),QCD_c(),QCD_ac(),
 			QED_noAO(),QCD_c_noAO(),QCD_ac_noAO()
     {}
 
     /**
      *  QED scale
      */
     Energy QED;
 
     /**
      * QCD colour scale
      */
     Energy QCD_c;
 
     /**
      *  QCD anticolour scale
      */
     Energy QCD_ac;
 
     /**
      *  QED scale
      */
     Energy QED_noAO;
 
     /**
      * QCD colour scale
      */
     Energy QCD_c_noAO;
 
     /**
      *  QCD anticolour scale
      */
     Energy QCD_ac_noAO;
 
     /**
      *   EW scales
      */
-    Energy EW_Z;
-    Energy EW_W;
+    Energy EW;
     
   };
 
 
   /** @name Construction and descruction functions. */
   //@{
 
   /**
    * Standard Constructor. Note that the default constructor is
    * private - there is no particle without a pointer to a
    * ParticleData object.
    * @param x the ParticleData object
    * @param fs  Whether or not the particle is an inital or final-state particle
    * @param tls Whether or not the particle initiates a time-like shower
    */
   ShowerParticle(tcEventPDPtr x, bool fs, bool tls=false) 
     : Particle(x), _isFinalState(fs),
       _perturbative(0), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
       _vMass(ZERO), _thePEGBase() {}
 
   /**
    * Copy constructor from a ThePEG Particle
    * @param x ThePEG particle
    * @param pert Where the particle came from
    * @param fs Whether or not the particle is an inital or final-state particle
    * @param tls Whether or not the particle initiates a time-like shower
    */
   ShowerParticle(const Particle & x, unsigned int pert, bool fs, bool tls=false)
     : Particle(x), _isFinalState(fs),
     _perturbative(pert), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
     _vMass(ZERO), _thePEGBase(&x) {}
   //@}
 
 public:
 
   /**
    *  Set a preliminary momentum for the particle
    */
   void setShowerMomentum(bool timelike);
 
   /**
    *  Construct the spin info object for a shower particle
    */
   void constructSpinInfo(bool timelike);
 
   /**
    * Perform any initial calculations needed after the branching has been selected
    */
   void initializeDecay();
 
   /**
    * Perform any initial calculations needed after the branching has been selected
    * @param parent The beam particle
    */
   void initializeInitialState(PPtr parent);
 
   /**
    * Perform any initial calculations needed after the branching has been selected
    */
   void initializeFinalState();
 
   /**
    *   Access/Set various flags about the state of the particle
    */
   //@{
   /**
    * Access the flag that tells if the particle is final state
    * or initial state.
    */
   bool isFinalState() const { return _isFinalState; }
 
   /**
    * Access the flag that tells if the particle is initiating a
    * time like shower when it has been emitted in an initial state shower.
    */
   bool initiatesTLS() const { return _initiatesTLS; }
 
   /**
    * Access the flag which tells us where the particle came from
    * This is 0 for a particle produced in the shower, 1 if the particle came
    * from the hard sub-process and 2 is it came from a decay.
    */
   unsigned int perturbative() const { return _perturbative; }
   //@}
 
   /**
    * Set/Get the momentum fraction for initial-state particles
    */
   //@{
   /**
    *  For an initial state particle get the fraction of the beam momentum
    */
   void x(double x) { _x = x; }
 
   /**
    *  For an initial state particle set the fraction of the beam momentum
    */
   double x() const { return _x; }
   //@}
 
   /**
    * Set/Get methods for the ShowerKinematics objects
    */
   //@{
   /**
    * Access/ the ShowerKinematics object.
    */
   const ShoKinPtr & showerKinematics() const { return _showerKinematics; }
 
 
   /**
    * Set the ShowerKinematics object.
    */
   void showerKinematics(const ShoKinPtr in) { _showerKinematics = in; }
   //@}
 
   /**
    * Set/Get methods for the ShowerBasis objects
    */
   //@{
   /**
    * Access/ the ShowerBasis object.
    */
   const ShowerBasisPtr & showerBasis() const { return _showerBasis; }
 
 
   /**
    * Set the ShowerBasis object.
    */
   void showerBasis(const ShowerBasisPtr in, bool copy) {
     if(!copy) 
       _showerBasis = in;
     else {
       _showerBasis = new_ptr(ShowerBasis());
       _showerBasis->setBasis(in->pVector(),in->nVector(),in->frame());
     } 
   }
   //@}
 
   /**
    *  Members relating to the initial evolution scale and partner for the particle
    */
   //@{
   /**
    *  Veto emission at a given scale 
    */
   void vetoEmission(ShowerPartnerType type, Energy scale);
 
   /**
    *  Access to the evolution scales
    */
   const EvolutionScales & scales() const {return scales_;} 
 
   /**
    *  Access to the evolution scales
    */
   EvolutionScales & scales() {return scales_;} 
 
   /**
    * Return the virtual mass\f$
    */
   Energy virtualMass() const { return _vMass; }
 
   /**
    *  Set the virtual mass
    */
   void virtualMass(Energy mass) { _vMass = mass; }
 
   /** 
    * Return the partner
    */
   tShowerParticlePtr partner() const { return _partner; }
 
   /**
    * Set the partner
    */
   void partner(const tShowerParticlePtr partner) { _partner = partner; } 
 
   /**
    *  Get the possible partners 
    */
   vector<EvolutionPartner> & partners() { return partners_; }
 
   /**
    *  Add a possible partners 
    */
   void addPartner(EvolutionPartner in );
 
   /**
    *  Clear the evolution partners
    */
   void clearPartners() { partners_.clear(); }
     
   /** 
    * Return the progenitor of the shower
    */
   tShowerParticlePtr progenitor() const { return _progenitor; }
 
   /**
    * Set the progenitor of the shower
    */
   void progenitor(const tShowerParticlePtr progenitor) { _progenitor = progenitor; } 
   //@}
 
 
   /**
    *  Members to store and provide access to variables for a specific
    *  shower evolution scheme
    */
   //@{
   struct Parameters {
     Parameters() : alpha(1.), beta(), ptx(), pty(), pt() {}
     double alpha;
     double beta;
     Energy ptx;
     Energy pty;
     Energy pt;
   };
 
 
   /**
    *  Set the vector containing dimensionless variables
    */
   Parameters & showerParameters() { return _parameters; }
   //@}
 
   /**
    *  If this particle came from the hard process get a pointer to ThePEG particle
    *  it came from
    */
   const tcPPtr thePEGBase() const { return _thePEGBase; }
  
 public:
 
   /**
    *  Extract the rho matrix including mapping needed in the shower
    */
   RhoDMatrix extractRhoMatrix(bool forward);
 
   /**
    * For a particle which came from the hard process get the spin density and
    * the mapping required to the basis used in the Shower
    * @param rho The \f$\rho\f$ matrix
    * @param mapping The mapping
    * @param showerkin The ShowerKinematics object
    */
   bool getMapping(SpinPtr &, RhoDMatrix & map);
 
 protected:
 
   /**
    * Standard clone function.
    */
   virtual PPtr clone() const;
 
   /**
    * Standard clone function.
    */
   virtual PPtr fullclone() const;
 
 private:
 
   /**
    * The static object used to initialize the description of this class.
    * Indicates that this is a concrete class with persistent data.
    */
   static ClassDescription<ShowerParticle> initShowerParticle;
 
   /**
    * The assignment operator is private and must never be called.
    * In fact, it should not even be implemented.
    */
   ShowerParticle & operator=(const ShowerParticle &) = delete;
 
 private:
 
   /**
    *  Whether the particle is in the final or initial state
    */
   bool _isFinalState;
 
   /**
    *  Whether the particle came from 
    */
   unsigned int _perturbative;
 
   /**
    *  Does a particle produced in the backward shower initiate a time-like shower 
    */
   bool _initiatesTLS;
 
   /**
    * Dimensionless parameters
    */
   Parameters _parameters;
 
   /**
    *  The beam energy fraction for particle's in the initial state
    */
   double _x;
 
   /**
    *  The shower kinematics for the particle
    */
   ShoKinPtr _showerKinematics;
 
   /**
    *  The shower basis for the particle
    */
   ShowerBasisPtr _showerBasis;
 
   /**
    *  Storage of the evolution scales
    */
   EvolutionScales scales_;
 
   /**
    *  Virtual mass
    */
   Energy _vMass;
 
   /**
    *  Partners
    */
   tShowerParticlePtr _partner;
 
   /**
    *  Pointer to ThePEG Particle this ShowerParticle was created from
    */
   const tcPPtr _thePEGBase;
   
   /**
    *  Progenitor
    */   
   tShowerParticlePtr _progenitor;
 
   /**
    *  Partners
    */
   vector<EvolutionPartner> partners_;
     
 };
 
 inline ostream & operator<<(ostream & os, const ShowerParticle::EvolutionScales & es) {
   os << "Scales: QED=" << es.QED / GeV
      << " QCD_c=" << es.QCD_c / GeV
      << " QCD_ac=" << es.QCD_ac / GeV
-     << " EW_Z=" << es.EW_Z / GeV
-     << " EW_W=" << es.EW_W / GeV
+     << " EW=" << es.EW / GeV
      << " QED_noAO=" << es.QED_noAO / GeV
      << " QCD_c_noAO=" << es.QCD_c_noAO / GeV
      << " QCD_ac_noAO=" << es.QCD_ac_noAO / GeV
      << '\n';
   return os;
 }
 
 }
 
 #include "ThePEG/Utilities/ClassTraits.h"
 
 namespace ThePEG {
 
 /** @cond TRAITSPECIALIZATIONS */
 
 /** This template specialization informs ThePEG about the
  *  base classes of ShowerParticle. */
 template <>
 struct BaseClassTrait<Herwig::ShowerParticle,1> {
   /** Typedef of the first base class of ShowerParticle. */
   typedef Particle NthBase;
 };
 
 /** This template specialization informs ThePEG about the name of
  *  the ShowerParticle class and the shared object where it is defined. */
 template <>
 struct ClassTraits<Herwig::ShowerParticle>
   : public ClassTraitsBase<Herwig::ShowerParticle> {
   /** Return a platform-independent class name */
   static string className() { return "Herwig::ShowerParticle"; }
   /** Create a Event object. */
   static TPtr create() { return TPtr::Create(Herwig::ShowerParticle(tcEventPDPtr(),true)); }
 };
 
 /** @endcond */
 
 }
 
 #endif /* HERWIG_ShowerParticle_H */
diff --git a/Shower/QTilde/QTildeShowerHandler.cc b/Shower/QTilde/QTildeShowerHandler.cc
--- a/Shower/QTilde/QTildeShowerHandler.cc
+++ b/Shower/QTilde/QTildeShowerHandler.cc
@@ -1,3237 +1,3230 @@
 // -*- C++ -*-
 //
 // QTildeShowerHandler.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 QTildeShowerHandler class.
 //
 
 #include "QTildeShowerHandler.h"
 #include "ThePEG/Interface/Deleted.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/Interface/RefVector.h"
 #include "ThePEG/Interface/Parameter.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 "ThePEG/Utilities/EnumIO.h"
 #include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
 #include "Herwig/PDF/MPIPDF.h"
 #include "Herwig/PDF/MinBiasPDF.h"
 #include "Herwig/Shower/QTilde/Base/ShowerTree.h"
 #include "Herwig/Shower/QTilde/Base/HardTree.h"
 #include "Herwig/Shower/QTilde/Kinematics/KinematicsReconstructor.h"
 #include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
 #include "Herwig/PDF/HwRemDecayer.h"
 #include "Herwig/Shower/QTilde/Base/ShowerVertex.h"
 #include "ThePEG/Repository/CurrentGenerator.h"
 #include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h"
 #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
 #include "ThePEG/PDF/PartonExtractor.h"
 #include "Herwig/Shower/RealEmissionProcess.h"
 #include "Herwig/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h"
 #include "Herwig/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h"
 #include "Herwig/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h"
 
 
 using namespace Herwig;
 
 bool QTildeShowerHandler::_hardEmissionWarn = true;
 bool QTildeShowerHandler::_missingTruncWarn = true;
 
 QTildeShowerHandler::QTildeShowerHandler() :
   _maxtry(100), _meCorrMode(1), _evolutionScheme(1),
   _hardVetoReadOption(false),
   _iptrms(ZERO), _beta(0.), _gamma(ZERO), _iptmax(),
   _limitEmissions(0), _initialenhance(1.), _finalenhance(1.),
   _nReWeight(100), _reWeight(false),
   interaction_(ShowerInteraction::QEDQCD),
   _trunc_Mode(true), _hardEmission(1),
   _softOpt(2), _hardPOWHEG(false), muPt(ZERO)
 {}
 
 QTildeShowerHandler::~QTildeShowerHandler() {}
 
 IBPtr QTildeShowerHandler::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr QTildeShowerHandler::fullclone() const {
   return new_ptr(*this);
 }
 
 void QTildeShowerHandler::persistentOutput(PersistentOStream & os) const {
   os << _splittingGenerator << _maxtry 
      << _meCorrMode << _hardVetoReadOption
      << _limitEmissions << _softOpt << _hardPOWHEG
      << ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV) 
      << _vetoes << _fullShowerVetoes << _nReWeight << _reWeight
      << _trunc_Mode << _hardEmission << _evolutionScheme 
      << ounit(muPt,GeV) << oenum(interaction_) 
      << _reconstructor << _partnerfinder;
 }
 
 void QTildeShowerHandler::persistentInput(PersistentIStream & is, int) {
   is >> _splittingGenerator >> _maxtry 
      >> _meCorrMode >> _hardVetoReadOption
      >> _limitEmissions >> _softOpt >> _hardPOWHEG
      >> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
      >> _vetoes >> _fullShowerVetoes >> _nReWeight >> _reWeight
      >> _trunc_Mode >> _hardEmission >> _evolutionScheme
      >> iunit(muPt,GeV) >> ienum(interaction_)
      >> _reconstructor >> _partnerfinder;
 }
 
 
 // The following static variable is needed for the type
 // description system in ThePEG. 
 DescribeClass<QTildeShowerHandler,ShowerHandler>
 describeHerwigQTildeShowerHandler("Herwig::QTildeShowerHandler", "HwShower.so");
 
 void QTildeShowerHandler::Init() {
 
   static ClassDocumentation<QTildeShowerHandler> documentation
     ("TheQTildeShowerHandler class is the main class"
      " for the angular-ordered parton shower",
      "The Shower evolution was performed using an algorithm described in "
      "\\cite{Marchesini:1983bm,Marchesini:1987cf,Gieseke:2003rz,Bahr:2008pv}.",
      "%\\cite{Marchesini:1983bm}\n"
      "\\bibitem{Marchesini:1983bm}\n"
      "  G.~Marchesini and B.~R.~Webber,\n"
      "  ``Simulation Of QCD Jets Including Soft Gluon Interference,''\n"
      "  Nucl.\\ Phys.\\  B {\\bf 238}, 1 (1984).\n"
      "  %%CITATION = NUPHA,B238,1;%%\n"
      "%\\cite{Marchesini:1987cf}\n"
      "\\bibitem{Marchesini:1987cf}\n"
      "  G.~Marchesini and B.~R.~Webber,\n"
      "   ``Monte Carlo Simulation of General Hard Processes with Coherent QCD\n"
      "  Radiation,''\n"
      "  Nucl.\\ Phys.\\  B {\\bf 310}, 461 (1988).\n"
      "  %%CITATION = NUPHA,B310,461;%%\n"
      "%\\cite{Gieseke:2003rz}\n"
      "\\bibitem{Gieseke:2003rz}\n"
      "  S.~Gieseke, P.~Stephens and B.~Webber,\n"
      "  ``New formalism for QCD parton showers,''\n"
      "  JHEP {\\bf 0312}, 045 (2003)\n"
      "  [arXiv:hep-ph/0310083].\n"
      "  %%CITATION = JHEPA,0312,045;%%\n"
      );
 
   static Reference<QTildeShowerHandler,SplittingGenerator> 
     interfaceSplitGen("SplittingGenerator", 
 		      "A reference to the SplittingGenerator object", 
 		      &Herwig::QTildeShowerHandler::_splittingGenerator,
 		      false, false, true, false);
 
   static Parameter<QTildeShowerHandler,unsigned int> interfaceMaxTry
     ("MaxTry",
      "The maximum number of attempts to generate the shower from a"
      " particular ShowerTree",
      &QTildeShowerHandler::_maxtry, 100, 1, 100000,
      false, false, Interface::limited);
 
   static Parameter<QTildeShowerHandler,unsigned int> interfaceNReWeight
     ("NReWeight",
      "The number of attempts for the shower when reweighting",
      &QTildeShowerHandler::_nReWeight, 100, 10, 10000,
      false, false, Interface::limited);
 
   static Switch<QTildeShowerHandler, unsigned int> ifaceMECorrMode
     ("MECorrMode",
      "Choice of the ME Correction Mode",
      &QTildeShowerHandler::_meCorrMode, 1, false, false);
   static SwitchOption on
     (ifaceMECorrMode,"HardPlusSoft","hard+soft on", 1);
   static SwitchOption hard
     (ifaceMECorrMode,"Hard","only hard on", 2);
   static SwitchOption soft
     (ifaceMECorrMode,"Soft","only soft on", 3);
 
   static Switch<QTildeShowerHandler, bool> ifaceHardVetoReadOption
     ("HardVetoReadOption",
      "Apply read-in scale veto to all collisions or just the primary one?",
      &QTildeShowerHandler::_hardVetoReadOption, false, false, false);
   static SwitchOption AllCollisions
     (ifaceHardVetoReadOption,
      "AllCollisions",
      "Read-in pT veto applied to primary and secondary collisions.",
      false);
   static SwitchOption PrimaryCollision
     (ifaceHardVetoReadOption,
      "PrimaryCollision",
      "Read-in pT veto applied to primary but not secondary collisions.",
      true);
 
   static Parameter<QTildeShowerHandler, Energy> ifaceiptrms
     ("IntrinsicPtGaussian",
      "RMS of intrinsic pT of Gaussian distribution:\n"
      "2*(1-Beta)*exp(-sqr(intrinsicpT/RMS))/sqr(RMS)",
      &QTildeShowerHandler::_iptrms, GeV, ZERO, ZERO, 1000000.0*GeV,
      false, false, Interface::limited);
 
   static Parameter<QTildeShowerHandler, double> ifacebeta
     ("IntrinsicPtBeta",
      "Proportion of inverse quadratic distribution in generating intrinsic pT.\n"
      "(1-Beta) is the proportion of Gaussian distribution",
      &QTildeShowerHandler::_beta, 0, 0, 1,
      false, false, Interface::limited);
 
   static Parameter<QTildeShowerHandler, Energy> ifacegamma
     ("IntrinsicPtGamma",
      "Parameter for inverse quadratic:\n"
      "2*Beta*Gamma/(sqr(Gamma)+sqr(intrinsicpT))",
      &QTildeShowerHandler::_gamma,GeV, ZERO, ZERO, 100000.0*GeV,
      false, false, Interface::limited);
 
   static Parameter<QTildeShowerHandler, Energy> ifaceiptmax
     ("IntrinsicPtIptmax",
      "Upper bound on intrinsic pT for inverse quadratic",
      &QTildeShowerHandler::_iptmax,GeV, ZERO, ZERO, 100000.0*GeV,
      false, false, Interface::limited);
 
   static RefVector<QTildeShowerHandler,ShowerVeto> ifaceVetoes
     ("Vetoes",
      "The vetoes to be checked during showering",
      &QTildeShowerHandler::_vetoes, -1,
      false,false,true,true,false);
  
   static RefVector<QTildeShowerHandler,FullShowerVeto> interfaceFullShowerVetoes
     ("FullShowerVetoes",
      "The vetos to be appliede on the full final state of the shower",
      &QTildeShowerHandler::_fullShowerVetoes, -1, false, false, true, false, false);
 
   static Switch<QTildeShowerHandler,unsigned int> interfaceLimitEmissions
     ("LimitEmissions",
      "Limit the number and type of emissions for testing",
      &QTildeShowerHandler::_limitEmissions, 0, false, false);
   static SwitchOption interfaceLimitEmissionsNoLimit
     (interfaceLimitEmissions,
      "NoLimit",
      "Allow an arbitrary number of emissions",
      0);
   static SwitchOption interfaceLimitEmissionsOneInitialStateEmission
     (interfaceLimitEmissions,
      "OneInitialStateEmission",
      "Allow one emission in the initial state and none in the final state",
      1);
   static SwitchOption interfaceLimitEmissionsOneFinalStateEmission
     (interfaceLimitEmissions,
      "OneFinalStateEmission",
      "Allow one emission in the final state and none in the initial state",
      2);
   static SwitchOption interfaceLimitEmissionsHardOnly
     (interfaceLimitEmissions,
      "HardOnly",
      "Only allow radiation from the hard ME correction",
      3);
   static SwitchOption interfaceLimitEmissionsOneEmission
     (interfaceLimitEmissions,
      "OneEmission",
      "Allow one emission in either the final state or initial state, but not both",
      4);
 
   static Switch<QTildeShowerHandler,bool> interfaceTruncMode
     ("TruncatedShower", "Include the truncated shower?", 
      &QTildeShowerHandler::_trunc_Mode, 1, false, false);
   static SwitchOption interfaceTruncMode0
     (interfaceTruncMode,"No","Truncated Shower is OFF", 0);
   static SwitchOption interfaceTruncMode1
     (interfaceTruncMode,"Yes","Truncated Shower is ON", 1);
 
   static Switch<QTildeShowerHandler,int> interfaceHardEmission
     ("HardEmission",
      "Whether to use ME corrections or POWHEG for the hardest emission",
      &QTildeShowerHandler::_hardEmission, 0, false, false);
   static SwitchOption interfaceHardEmissionNone
     (interfaceHardEmission,
      "None",
      "No Corrections",
      0);
   static SwitchOption interfaceHardEmissionMECorrection
     (interfaceHardEmission,
      "MECorrection",
      "Old fashioned ME correction",
      1);
   static SwitchOption interfaceHardEmissionPOWHEG
     (interfaceHardEmission,
      "POWHEG",
      "Powheg style hard emission",
      2);
 
   static Switch<QTildeShowerHandler,ShowerInteraction> interfaceInteractions
     ("Interactions",
      "The interactions to be used in the shower",
      &QTildeShowerHandler::interaction_, ShowerInteraction::QEDQCD, false, false);
   static SwitchOption interfaceInteractionsQCD
     (interfaceInteractions,
      "QCD",
      "Only QCD radiation",
      ShowerInteraction::QCD);
   static SwitchOption interfaceInteractionsQED
     (interfaceInteractions,
      "QED",
      "Only QEd radiation",
      ShowerInteraction::QED);
   static SwitchOption interfaceInteractionEWOnly
     (interfaceInteractions,
      "EWOnly",
      "Only EW",
      ShowerInteraction::EW);
   static SwitchOption interfaceInteractionQEDQCD
     (interfaceInteractions,
      "QEDQCD",
      "QED and QCD",
      ShowerInteraction::QEDQCD);
   static SwitchOption interfaceInteractionALL
     (interfaceInteractions,
      "ALL",
      "QED, QCD and EW",
      ShowerInteraction::ALL);
   
   static Deleted<QTildeShowerHandler> delReconstructionOption
     ("ReconstructionOption", "The old reconstruction option switch has been replaced with"
      " the new EvolutionScheme switch, see  arXiv:1904.11866 for details");
 
   static Switch<QTildeShowerHandler,unsigned int> interfaceEvolutionScheme
     ("EvolutionScheme",
      "The scheme to interpret the evolution variable in the case of multple emission.",
      &QTildeShowerHandler::_evolutionScheme, 1, false, false);
   static SwitchOption interfaceEvolutionSchemepT
     (interfaceEvolutionScheme,
      "pT",
      "pT scheme",
      0);
   static SwitchOption interfaceEvolutionSchemeDotProduct
     (interfaceEvolutionScheme,
      "DotProduct",
      "Dot-product scheme",
      1);
   static SwitchOption interfaceEvolutionSchemeQ2
     (interfaceEvolutionScheme,
      "Q2",
      "Q2 scheme",
      2);
 
   static Switch<QTildeShowerHandler,unsigned int> interfaceSoftCorrelations
     ("SoftCorrelations",
      "Option for the treatment of soft correlations in the parton shower",
      &QTildeShowerHandler::_softOpt, 2, false, false);
   static SwitchOption interfaceSoftCorrelationsNone
     (interfaceSoftCorrelations,
      "No",
      "No soft correlations",
      0);
   static SwitchOption interfaceSoftCorrelationsFull
     (interfaceSoftCorrelations,
      "Full",
      "Use the full eikonal",
      1);
   static SwitchOption interfaceSoftCorrelationsSingular
     (interfaceSoftCorrelations,
      "Singular",
      "Use original Webber-Marchisini form",
      2);
 
   static Switch<QTildeShowerHandler,bool> interfaceHardPOWHEG
     ("HardPOWHEG",
      "Treatment of powheg emissions which are too hard to have a shower interpretation",
      &QTildeShowerHandler::_hardPOWHEG, false, false, false);
   static SwitchOption interfaceHardPOWHEGAsShower
     (interfaceHardPOWHEG,
      "AsShower",
      "Still interpret as shower emissions",
      false);
   static SwitchOption interfaceHardPOWHEGRealEmission
     (interfaceHardPOWHEG,
      "RealEmission",
      "Generate shower from the real emmission configuration",
      true);
 
   static Reference<QTildeShowerHandler,KinematicsReconstructor> interfaceKinematicsReconstructor
     ("KinematicsReconstructor",
      "Reference to the KinematicsReconstructor object",
      &QTildeShowerHandler::_reconstructor, false, false, true, false, false);
 
   static Reference<QTildeShowerHandler,PartnerFinder> interfacePartnerFinder
     ("PartnerFinder",
      "Reference to the PartnerFinder object",
      &QTildeShowerHandler::_partnerfinder, false, false, true, false, false);
 
 }
 
 tPPair QTildeShowerHandler::cascade(tSubProPtr sub,
 				    XCPtr xcomb) {
   // use me for reference in tex file etc
   useMe();
   prepareCascade(sub);
   // set things up in the base class
   resetWeights();
   hard_=ShowerTreePtr();
   decay_.clear();
   done_.clear();
   // start of the try block for the whole showering process
   unsigned int countFailures=0;
   while (countFailures<maxtry()) {
     try {
       decay_.clear();
       done_.clear();
       PerturbativeProcessPtr hard;
       DecayProcessMap decay;
       splitHardProcess(firstInteraction() ? tagged() :
       		       tPVector(currentSubProcess()->outgoing().begin(),
       				currentSubProcess()->outgoing().end()),
       		       hard,decay);
       ShowerTree::constructTrees(hard_,decay_,hard,decay);
       // if no hard process
       if(!hard_)  throw Exception() << "Shower starting with a decay"
 				    << "is not implemented" 
 				    << Exception::runerror;
       // perform the shower for the hard process
       showerHardProcess(hard_,xcomb);
       done_.push_back(hard_);
       hard_->updateAfterShower(decay_);
       // if no decaying particles to shower break out of the loop
       if(decay_.empty()) break;
       // shower the decay products
       while(!decay_.empty()) {
 	// find particle whose production process has been showered
 	ShowerDecayMap::iterator dit = decay_.begin();
 	while(!dit->second->parent()->hasShowered() && dit!=decay_.end()) ++dit;
 	assert(dit!=decay_.end());
 	// get the particle
 	ShowerTreePtr decayingTree = dit->second;
 	// remove it from the multimap
 	decay_.erase(dit);
 	// make sure the particle has been decayed
 	QTildeShowerHandler::decay(decayingTree,decay_);
 	// now shower the decay
 	showerDecay(decayingTree);
 	done_.push_back(decayingTree);
 	decayingTree->updateAfterShower(decay_);
       }
       // suceeded break out of the loop
       break;
     }
     catch (KinematicsReconstructionVeto) {
       resetWeights();
       ++countFailures;
     }
     catch ( ... ) {
       hard_=ShowerTreePtr();
       decay_.clear();
       done_.clear();
       throw;
     }
   }
   // if loop exited because of too many tries, throw event away
   if (countFailures >= maxtry()) {
     resetWeights();
     hard_=ShowerTreePtr();
     decay_.clear();
     done_.clear();
     throw Exception() << "Too many tries for main while loop "
 		      << "in QTildeShowerHandler::cascade()." 
 		      << Exception::eventerror; 	
   }
   //enter the particles in the event record
   fillEventRecord();
   // clear storage
   hard_=ShowerTreePtr();
   decay_.clear();
   done_.clear();
   // non hadronic case return
   if (!isResolvedHadron(incomingBeams().first ) && 
       !isResolvedHadron(incomingBeams().second) )
     return incomingBeams();
   // remake the remnants (needs to be after the colours are sorted
   //                       out in the insertion into the event record)
   if ( firstInteraction() ) return remakeRemnant(sub->incoming());
   //Return the new pair of incoming partons. remakeRemnant is not
   //necessary here, because the secondary interactions are not yet
   //connected to the remnants.
   return make_pair(findFirstParton(sub->incoming().first ),
 		   findFirstParton(sub->incoming().second));
 }
 
 void QTildeShowerHandler::fillEventRecord() {
   // create a new step 
   StepPtr pstep = newStep();
   assert(!done_.empty());
   assert(done_[0]->isHard());
   // insert the steps
   for(unsigned int ix=0;ix<done_.size();++ix) {
     done_[ix]->fillEventRecord(pstep,doISR(),doFSR());
   }
 }
 
 HardTreePtr QTildeShowerHandler::generateCKKW(ShowerTreePtr ) const {
   return HardTreePtr();
 }
 
 void QTildeShowerHandler::doinit() {
   ShowerHandler::doinit();
   // check on the reweighting
   for(unsigned int ix=0;ix<_fullShowerVetoes.size();++ix) {
     if(_fullShowerVetoes[ix]->behaviour()==1) {
       _reWeight = true;
       break;
     }
   }
   if(_reWeight && maximumTries()<_nReWeight) {
     throw Exception() << "Reweight being performed in the shower but the number of attempts for the"
 		      << "shower is less than that for the reweighting.\n"
 		      << "Maximum number of attempt for the shower "
 		      << fullName() << ":MaxTry is " << maximumTries() << "\nand for reweighting is "
 		      << fullName() << ":NReWeight is " << _nReWeight << "\n"
 		      << "we recommend the number of attempts is 10 times the number for reweighting\n"
 		      << Exception::runerror;
   }
   ShowerTree::_vmin2 = vMin();
   ShowerTree::_spaceTime = includeSpaceTime();
 }
 
 void QTildeShowerHandler::doinitrun() {
   ShowerHandler::doinitrun();
   ShowerTree::_vmin2 = vMin();
   ShowerTree::_spaceTime = includeSpaceTime();
 }
   
 void QTildeShowerHandler::generateIntrinsicpT(vector<ShowerProgenitorPtr> particlesToShower) {
   _intrinsic.clear();
   if ( !ipTon() || !doISR() ) return;
   // don't do anything for the moment for secondary scatters
   if( !firstInteraction() ) return;
   // generate intrinsic pT
   for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
     // only consider initial-state particles
     if(particlesToShower[ix]->progenitor()->isFinalState()) continue;
     if(!particlesToShower[ix]->progenitor()->dataPtr()->coloured()) continue;
     Energy ipt;
     if(UseRandom::rnd() > _beta) {
       ipt=_iptrms*sqrt(-log(UseRandom::rnd()));
     }
     else {
       ipt=_gamma*sqrt(pow(1.+sqr(_iptmax/_gamma), UseRandom::rnd())-1.);
     }
     pair<Energy,double> pt = make_pair(ipt,UseRandom::rnd(Constants::twopi));
     _intrinsic[particlesToShower[ix]] = pt;
   }
 }
 
 void QTildeShowerHandler::setupMaximumScales(const vector<ShowerProgenitorPtr> & p,
 				 XCPtr xcomb) {
   // let POWHEG events radiate freely
   if(_hardEmission==2&&hardTree()) {
     vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
     for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(Constants::MaxEnergy);
     return;
   }
   // return if no vetos
   if (!restrictPhasespace()) return; 
   // find out if hard partonic subprocess.
   bool isPartonic(false); 
   map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator 
     cit = _currenttree->incomingLines().begin();
   Lorentz5Momentum pcm;
   for(; cit!=currentTree()->incomingLines().end(); ++cit) {
     pcm += cit->first->progenitor()->momentum();
     isPartonic |= cit->first->progenitor()->coloured();
   }
   // find minimum pt from hard process, the maximum pt from all outgoing
   // coloured lines (this is simpler and more general than
   // 2stu/(s^2+t^2+u^2)).  Maximum scale for scattering processes will
   // be transverse mass.
   Energy ptmax = generator()->maximumCMEnergy();
   // general case calculate the scale  
   if ( !hardScaleIsMuF() || (hardVetoReadOption()&&!firstInteraction()) ) {
     // scattering process
     if(currentTree()->isHard()) {
       assert(xcomb);
       // coloured incoming particles
       if (isPartonic) {
 	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator 
 	  cjt = currentTree()->outgoingLines().begin();
 	for(; cjt!=currentTree()->outgoingLines().end(); ++cjt) {
 	  if (cjt->first->progenitor()->coloured())
 	    ptmax = min(ptmax,cjt->first->progenitor()->momentum().mt());
 	}
       }
       if (ptmax == generator()->maximumCMEnergy() ) ptmax = pcm.m();
       if(hardScaleIsMuF()&&hardVetoReadOption()&&
 	 !firstInteraction()) {
 	ptmax=min(ptmax,sqrt(xcomb->lastShowerScale()));
       }
     } 
     // decay, incoming() is the decaying particle.
     else { 
       ptmax = currentTree()->incomingLines().begin()->first
 	->progenitor()->momentum().mass(); 
     }
   }
   // hepeup.SCALUP is written into the lastXComb by the
   // LesHouchesReader itself - use this by user's choice. 
   // Can be more general than this. 
   else {
     if(currentTree()->isHard()) {
       assert(xcomb);
       ptmax = sqrt( xcomb->lastShowerScale() );
     }
     else {
       ptmax = currentTree()->incomingLines().begin()->first
 	->progenitor()->momentum().mass(); 
     }
   }
   ptmax *= hardScaleFactor();
   // set maxHardPt for all progenitors.  For partonic processes this
   // is now the max pt in the FS, for non-partonic processes or
   // processes with no coloured FS the invariant mass of the IS
   vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
   for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(ptmax);
 }
 
 void QTildeShowerHandler::setupHardScales(const vector<ShowerProgenitorPtr> & p,
 					  XCPtr xcomb) {
   if ( hardScaleIsMuF() &&
        (!hardVetoReadOption() || firstInteraction()) ) {
     Energy hardScale = ZERO;
     if(currentTree()->isHard()) {
       assert(xcomb);
       hardScale = sqrt( xcomb->lastShowerScale() );
     }
     else {
       hardScale = currentTree()->incomingLines().begin()->first
 	->progenitor()->momentum().mass(); 
     }
     hardScale *= hardScaleFactor();
     vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
     for (; ckt != p.end(); ckt++) (*ckt)->hardScale(hardScale);
     muPt = hardScale;
   }
 }
 
 void QTildeShowerHandler::showerHardProcess(ShowerTreePtr hard, XCPtr xcomb) {
   _hardme = HwMEBasePtr();
   // extract the matrix element
   tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(xcomb);
   if(lastXC) {
     _hardme = dynamic_ptr_cast<HwMEBasePtr>(lastXC->matrixElement());
   }
   _decayme = HwDecayerBasePtr();
   // set the current tree
   currentTree(hard);
   hardTree(HardTreePtr());
   // work out the type of event
   currentTree()->xcombPtr(dynamic_ptr_cast<StdXCombPtr>(xcomb));
   currentTree()->identifyEventType();
   checkFlags();
   // generate the showering
   doShowering(true,xcomb);
 }
 
 RealEmissionProcessPtr QTildeShowerHandler::hardMatrixElementCorrection(bool hard) {
   // set the initial enhancement factors for the soft correction
   _initialenhance = 1.;
   _finalenhance   = 1.;
   // see if we can get the correction from the matrix element
   // or decayer
   RealEmissionProcessPtr real;
   if(hard) {
     if(_hardme&&_hardme->hasMECorrection()) {
       _hardme->initializeMECorrection(_currenttree->perturbativeProcess(),
      				      _initialenhance,_finalenhance);
       if(hardMEC())
      	real = 
 	  _hardme->applyHardMatrixElementCorrection(_currenttree->perturbativeProcess());
     }
   }
   else {
     if(_decayme&&_decayme->hasMECorrection()) {
       _decayme->initializeMECorrection(_currenttree->perturbativeProcess(),
    				       _initialenhance,_finalenhance);
       if(hardMEC())
    	real = _decayme->applyHardMatrixElementCorrection(_currenttree->perturbativeProcess());
     }
   }
   return real;
 }
 
 ShowerParticleVector QTildeShowerHandler::createTimeLikeChildren(tShowerParticlePtr, IdList ids) {
   // Create the ShowerParticle objects for the two children of
   // the emitting particle; set the parent/child relationship
   // if same as definition create particles, otherwise create cc
   ShowerParticleVector children;
   for(unsigned int ix=0;ix<2;++ix) {
     children.push_back(new_ptr(ShowerParticle(ids[ix+1],true)));
     if(children[ix]->id()==_progenitor->id()&&!ids[ix+1]->stable()&&abs(ids[ix+1]->id())!=ParticleID::tauminus)
       children[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
     else
       children[ix]->set5Momentum(Lorentz5Momentum(ids[ix+1]->mass()));
   }
   return children;
 }
 
 bool QTildeShowerHandler::timeLikeShower(tShowerParticlePtr particle, 
 					 ShowerInteraction type,
 					 Branching fb, bool first) {
   // don't do anything if not needed
   if(_limitEmissions == 1 || hardOnly() || 
      ( _limitEmissions == 2 && _nfs != 0) ||
      ( _limitEmissions == 4 && _nfs + _nis != 0) ) {
     if(particle->spinInfo()) particle->spinInfo()->develop();
     return false;
   }
   // generate the emission
   ShowerParticleVector children;
   // generate the emission
   if(!fb.kinematics) 
     fb = selectTimeLikeBranching(particle,type,HardBranchingPtr());
   // no emission, return
   if(!fb.kinematics) {
     if(particle->spinInfo()) particle->spinInfo()->develop();
     return false;
   }
   Branching fc[2] = {Branching(),Branching()};
 
   assert(fb.kinematics);
   // has emitted
   // Assign the shower kinematics to the emitting particle.
   particle->showerKinematics(fb.kinematics);
   // check highest pT
   if(fb.kinematics->pT()>progenitor()->highestpT())
     progenitor()->highestpT(fb.kinematics->pT());
   // create the children
   children = createTimeLikeChildren(particle,fb.ids);
   // update the children
   particle->showerKinematics()->
     updateChildren(particle, children,fb.type);
   // update number of emissions
   ++_nfs;
   if(_limitEmissions!=0) {
     if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
     if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
     if(particle->spinInfo()) particle->spinInfo()->develop();
     return true;
   }
   // select branchings for children
   fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
   fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
   // shower the first  particle
   if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
   if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
   // shower the second particle
   if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
   if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
   particle->showerKinematics()->updateParent(particle, children,_evolutionScheme,fb.type);
   // branching has happened
   if(first&&!children.empty())
     particle->showerKinematics()->resetChildren(particle,children);
   if(particle->spinInfo()) particle->spinInfo()->develop();
   return true;
 }
 
 bool 
 QTildeShowerHandler::spaceLikeShower(tShowerParticlePtr particle, PPtr beam,
 				     ShowerInteraction type) {
   //using the pdf's associated with the ShowerHandler assures, that
   //modified pdf's are used for the secondary interactions via 
   //CascadeHandler::resetPDFs(...)
   tcPDFPtr pdf;
   if(firstPDF().particle() == _beam)
     pdf = firstPDF().pdf();
   if(secondPDF().particle() == _beam)
     pdf = secondPDF().pdf();
   Energy freeze = pdfFreezingScale();
   // don't do anything if not needed
   if(_limitEmissions == 2  || hardOnly() ||
      ( _limitEmissions == 1 && _nis != 0 ) ||
      ( _limitEmissions == 4 && _nis + _nfs != 0 ) ) {
     if(particle->spinInfo()) particle->spinInfo()->develop();
     return false;
   }
   Branching bb;
   // generate branching
   while (true) {
     bb=_splittingGenerator->chooseBackwardBranching(*particle,beam,
 						    _initialenhance,
 						    _beam,type,
 						    pdf,freeze);
     // return if no emission
     if(!bb.kinematics) {
       if(particle->spinInfo()) particle->spinInfo()->develop();
       return false;
     }
     // if not vetoed break
     if(!spaceLikeVetoed(bb,particle)) break;
     // otherwise reset scale and continue
     particle->vetoEmission(bb.type,bb.kinematics->scale());
     if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
   }
   // assign the splitting function and shower kinematics
   particle->showerKinematics(bb.kinematics);
   if(bb.kinematics->pT()>progenitor()->highestpT())
     progenitor()->highestpT(bb.kinematics->pT());
   // For the time being we are considering only 1->2 branching
   // particles as in Sudakov form factor
   tcPDPtr part[2]={bb.ids[0],bb.ids[2]};
   // Now create the actual particles, make the otherChild a final state
   // particle, while the newParent is not
   ShowerParticlePtr newParent  = new_ptr(ShowerParticle(part[0],false));
   ShowerParticlePtr otherChild = new_ptr(ShowerParticle(part[1],true,true));
   ShowerParticleVector theChildren;
   theChildren.push_back(particle); 
   theChildren.push_back(otherChild);
   //this updates the evolution scale
   particle->showerKinematics()->
     updateParent(newParent, theChildren,_evolutionScheme,bb.type);
   // update the history if needed
   _currenttree->updateInitialStateShowerProduct(_progenitor,newParent);
   _currenttree->addInitialStateBranching(particle,newParent,otherChild);
   // for the reconstruction of kinematics, parent/child
   // relationships are according to the branching process:
   // now continue the shower
   ++_nis;
   bool emitted = _limitEmissions==0 ? 
     spaceLikeShower(newParent,beam,type) : false;
   if(newParent->spinInfo()) newParent->spinInfo()->develop();
   // now reconstruct the momentum
   if(!emitted) {
     if(_intrinsic.find(_progenitor)==_intrinsic.end()) {
       bb.kinematics->updateLast(newParent,ZERO,ZERO);
     }
     else {
       pair<Energy,double> kt=_intrinsic[_progenitor];
       bb.kinematics->updateLast(newParent,
 				kt.first*cos(kt.second),
 				kt.first*sin(kt.second));
     }
   }
   particle->showerKinematics()->
     updateChildren(newParent, theChildren,bb.type);
   if(_limitEmissions!=0) {
     if(particle->spinInfo()) particle->spinInfo()->develop();
     return true;
   }
   // perform the shower of the final-state particle
   timeLikeShower(otherChild,type,Branching(),true);
   updateHistory(otherChild);
   if(theChildren[1]->spinInfo()) theChildren[1]->spinInfo()->develop();
   // return the emitted
   if(particle->spinInfo()) particle->spinInfo()->develop();
   return true;
 }
 
 void QTildeShowerHandler::showerDecay(ShowerTreePtr decay) {
   // work out the type of event
   currentTree()->xcombPtr(StdXCombPtr());
   currentTree()->identifyEventType();
   _decayme = HwDecayerBasePtr();
   _hardme  = HwMEBasePtr();
   // find the decayer
   // try the normal way if possible
   tDMPtr dm = decay->incomingLines().begin()->first->original()   ->decayMode();
   if(!dm) dm = decay->incomingLines().begin()->first->copy()      ->decayMode();
   if(!dm) dm = decay->incomingLines().begin()->first->progenitor()->decayMode();
   // otherwise make a string and look it up
   if(!dm) {
     string tag = decay->incomingLines().begin()->first->original()->dataPtr()->name() 
       + "->";
     OrderedParticles outgoing;
     for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator 
 	  it=decay->outgoingLines().begin();it!=decay->outgoingLines().end();++it) {
       if(abs(decay->incomingLines().begin()->first->original()->id()) == ParticleID::t &&
 	 abs(it->first->original()->id())==ParticleID::Wplus &&
 	 decay->treelinks().size() == 1) {
 	ShowerTreePtr Wtree = decay->treelinks().begin()->first;
 	for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
 	      it2=Wtree->outgoingLines().begin();it2!=Wtree->outgoingLines().end();++it2) {
 	  outgoing.insert(it2->first->original()->dataPtr());
 	}
       }
       else {
 	outgoing.insert(it->first->original()->dataPtr());
       }
     }
     for(OrderedParticles::const_iterator it=outgoing.begin(); it!=outgoing.end();++it) {
       if(it!=outgoing.begin()) tag += ",";
       tag +=(**it).name();
     }
     tag += ";";
     dm = findDecayMode(tag);
   }
   if(dm) _decayme = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
   // set the ShowerTree to be showered
   currentTree(decay);
   decay->applyTransforms();
   hardTree(HardTreePtr());
   // generate the showering
   doShowering(false,XCPtr());
   // if no vetos 
   // force calculation of spin correlations
   SpinPtr spInfo = decay->incomingLines().begin()->first->progenitor()->spinInfo();
   if(spInfo) {
     if(!spInfo->developed()) spInfo->needsUpdate();
     spInfo->develop();
   }
 }
 
 bool QTildeShowerHandler::spaceLikeDecayShower(tShowerParticlePtr particle,
 				   const ShowerParticle::EvolutionScales & maxScales,
 				   Energy minmass,ShowerInteraction type,
 				   Branching fb) {
   // don't do anything if not needed
   if(_limitEmissions == 1 || hardOnly() || 
      ( _limitEmissions == 3 && _nis != 0) ||
      ( _limitEmissions == 4 && _nfs + _nis != 0) ) {
     return false;
   }
   // generate the emission
   ShowerParticleVector children;
   // generate the emission
   if(!fb.kinematics) 
     fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,
 				       HardBranchingPtr());
   // no emission, return
   if(!fb.kinematics) return false;
   Branching fc[2];
   if(particle->virtualMass()==ZERO) 
     particle->virtualMass(_progenitor->progenitor()->mass());
   fc[0] = Branching();
   fc[1] = Branching();
   assert(fb.kinematics);
   // has emitted
   // Assign the shower kinematics to the emitting particle.
   particle->showerKinematics(fb.kinematics);
   if(fb.kinematics->pT()>progenitor()->highestpT())
     progenitor()->highestpT(fb.kinematics->pT());
   // create the ShowerParticle objects for the two children
   children = createTimeLikeChildren(particle,fb.ids);
   // updateChildren the children
   particle->showerKinematics()->
     updateChildren(particle, children, fb.type);
   // select branchings for children
   fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,
 					type,HardBranchingPtr());
   fc[1] = selectTimeLikeBranching      (children[1],type,HardBranchingPtr());
   // old default
   ++_nis;
   // shower the first  particle
   _currenttree->updateInitialStateShowerProduct(_progenitor,children[0]);
   _currenttree->addInitialStateBranching(particle,children[0],children[1]);
   if(fc[0].kinematics) spaceLikeDecayShower(children[0],maxScales,minmass,type,Branching());
   // shower the second particle
   if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],true);
   updateHistory(children[1]);
   // TODO NEED AN UPDATE HERE FOR RECONOPT!=0
   // branching has happened
   return true;
 }
 
 vector<ShowerProgenitorPtr> QTildeShowerHandler::setupShower(bool hard) {
   RealEmissionProcessPtr real;
   // generate hard me if needed
   if(_hardEmission==1) {
     real = hardMatrixElementCorrection(hard);
     if(real&&!real->outgoing().empty()) setupMECorrection(real);
   }
   // generate POWHEG hard emission if needed
   else if(_hardEmission==2) 
     hardestEmission(hard);
   // set the initial colour partners
   setEvolutionPartners(hard,interaction_,false);
   // get the particles to be showered
   vector<ShowerProgenitorPtr> particlesToShower = 
     currentTree()->extractProgenitors();
   // return the answer
   return particlesToShower;
 }
 
 void QTildeShowerHandler::setEvolutionPartners(bool hard,ShowerInteraction type,
 					       bool clear) {
   // match the particles in the ShowerTree and hardTree
   if(hardTree() && !hardTree()->connect(currentTree()))
     throw Exception() << "Can't match trees in "
 		      << "QTildeShowerHandler::setEvolutionPartners()"
 		      << Exception::eventerror;
   // extract the progenitors
   vector<ShowerParticlePtr> particles = 
     currentTree()->extractProgenitorParticles();
   // clear the partners if needed
   if(clear) {
     for(unsigned int ix=0;ix<particles.size();++ix) {
       particles[ix]->partner(ShowerParticlePtr());
       particles[ix]->clearPartners();
     }
   }
   // sort out the colour partners
   if(hardTree()) {
     // find the partner
     for(unsigned int ix=0;ix<particles.size();++ix) {
       tShowerParticlePtr partner = hardTree()->particles()[particles[ix]]->branchingParticle()->partner();
       if(!partner) continue;
       for(map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
 	    it=hardTree()->particles().begin();
 	  it!=hardTree()->particles().end();++it) {
 	if(it->second->branchingParticle()==partner) {
 	  particles[ix]->partner(it->first);
 	  break;
 	}
       }
       if(!particles[ix]->partner()) 
 	throw Exception() << "Can't match partners in "
 			  << "QTildeShowerHandler::setEvolutionPartners()"
 			  << Exception::eventerror;
     }
   }
 
   // Set the initial evolution scales
   partnerFinder()->
     setInitialEvolutionScales(particles,!hard,interaction_,!_hardtree);
   if(hardTree() && _hardPOWHEG) {
     bool tooHard=false;
     map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator 
       eit=hardTree()->particles().end();
     for(unsigned int ix=0;ix<particles.size();++ix) {
       map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator 
 	mit = hardTree()->particles().find(particles[ix]);
       Energy hardScale(ZERO);
       ShowerPartnerType type(ShowerPartnerType::Undefined);
       // final-state
       if(particles[ix]->isFinalState()) {
 	if(mit!= eit && !mit->second->children().empty()) {
 	  hardScale = mit->second->scale();
 	  type = mit->second->type();
 	}
       }
       // initial-state
       else {
 	if(mit!= eit && mit->second->parent()) {
 	  hardScale = mit->second->parent()->scale();
 	  type = mit->second->parent()->type();
 	}
       }
       if(type!=ShowerPartnerType::Undefined) {
 	if(type==ShowerPartnerType::QED) {
 	  tooHard |= particles[ix]->scales().QED_noAO<hardScale;
 	}
 	else if(type==ShowerPartnerType::QCDColourLine) {
 	  tooHard |= particles[ix]->scales().QCD_c_noAO<hardScale;
 	}
 	else if(type==ShowerPartnerType::QCDAntiColourLine) {
 	  tooHard |= particles[ix]->scales().QCD_ac_noAO<hardScale;
 	}
-	else if(type==ShowerPartnerType::EW && 
-	       (particles[1]->id()==ParticleID::Z0 
-		   || particles[2]->id()==ParticleID::Z0)) {
-	  tooHard |= particles[ix]->scales().EW_Z<hardScale;
-	}
-	else if(type==ShowerPartnerType::EW && 
-	       (abs(particles[1]->id())==ParticleID::Wplus 
-	       || abs(particles[2]->id())==ParticleID::Wplus)) {
-	  tooHard |= particles[ix]->scales().EW_W<hardScale;
+	else if(type==ShowerPartnerType::EW) {
+	  tooHard |= particles[ix]->scales().EW<hardScale;
 	}
       }
     }
     if(tooHard) convertHardTree(hard,type);
   }
 }
 
 void QTildeShowerHandler::updateHistory(tShowerParticlePtr particle) {
   if(!particle->children().empty()) {
     ShowerParticleVector theChildren;
     for(unsigned int ix=0;ix<particle->children().size();++ix) {
       ShowerParticlePtr part = dynamic_ptr_cast<ShowerParticlePtr>
 	(particle->children()[ix]);
       theChildren.push_back(part);
     }
     // update the history if needed
     if(particle==_currenttree->getFinalStateShowerProduct(_progenitor))
       _currenttree->updateFinalStateShowerProduct(_progenitor,
 						  particle,theChildren);
     _currenttree->addFinalStateBranching(particle,theChildren);
     for(unsigned int ix=0;ix<theChildren.size();++ix)
       updateHistory(theChildren[ix]);
   }
 }
 
 bool QTildeShowerHandler::startTimeLikeShower(ShowerInteraction type) {
   // initialize basis vectors etc
   if(!progenitor()->progenitor()->partner()) return false;
   progenitor()->progenitor()->initializeFinalState();
   if(hardTree()) {
     map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator 
       eit=hardTree()->particles().end(),
       mit = hardTree()->particles().find(progenitor()->progenitor());
     if( mit != eit && !mit->second->children().empty() ) {
       bool output=truncatedTimeLikeShower(progenitor()->progenitor(),
 					  mit->second ,type,Branching(),true);
       if(output) updateHistory(progenitor()->progenitor());
       return output;
     }
   }
   // do the shower
   bool output = hardOnly() ? false :
     timeLikeShower(progenitor()->progenitor() ,type,Branching(),true) ;
   if(output) updateHistory(progenitor()->progenitor());
   return output;
 }
 
 bool QTildeShowerHandler::startSpaceLikeShower(PPtr parent, ShowerInteraction type) {
   // initialise the basis vectors
   if(!progenitor()->progenitor()->partner()) return false;
   progenitor()->progenitor()->initializeInitialState(parent);
   if(hardTree()) {
     map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator 
       eit =hardTree()->particles().end(),
       mit = hardTree()->particles().find(progenitor()->progenitor());
     if( mit != eit && mit->second->parent() ) {
       return truncatedSpaceLikeShower( progenitor()->progenitor(),
 				       parent, mit->second->parent(), type );
     } 
   }
   // perform the shower
   return  hardOnly() ? false :
     spaceLikeShower(progenitor()->progenitor(),parent,type);
 }
 
 bool QTildeShowerHandler::
 startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
 			  Energy minimumMass,ShowerInteraction type) {
   // set up the particle basis vectors
   if(!progenitor()->progenitor()->partner()) return false;
   progenitor()->progenitor()->initializeDecay();
   if(hardTree()) {
     map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator 
       eit =hardTree()->particles().end(),
       mit = hardTree()->particles().find(progenitor()->progenitor());
     if( mit != eit && mit->second->parent() ) {
       HardBranchingPtr branch=mit->second;
       while(branch->parent()) branch=branch->parent();
       return truncatedSpaceLikeDecayShower(progenitor()->progenitor(),maxScales,
 					   minimumMass, branch ,type, Branching());
     }
   }
   // perform the shower
   return  hardOnly() ? false :
     spaceLikeDecayShower(progenitor()->progenitor(),maxScales,minimumMass,type,Branching());
 }
 
 bool QTildeShowerHandler::timeLikeVetoed(const Branching & fb,
 			     ShowerParticlePtr particle) {
   // work out type of interaction
   ShowerInteraction type = convertInteraction(fb.type);
   // check whether emission was harder than largest pt of hard subprocess
   if ( restrictPhasespace() && fb.kinematics->pT() > _progenitor->maxHardPt() ) 
     return true;
   // soft matrix element correction veto
   if( softMEC()) {
     if(_hardme && _hardme->hasMECorrection()) {
       if(_hardme->softMatrixElementVeto(particle,
 					_progenitor->progenitor(),
 					particle->isFinalState(),
 					_progenitor->highestpT(),
 					fb.ids, fb.kinematics->z(),
 					fb.kinematics->scale(),
 					fb.kinematics->pT()))
 	return true;
     }
     else if(_decayme && _decayme->hasMECorrection()) {
       if(_decayme->softMatrixElementVeto(particle,
 					 _progenitor->progenitor(),
 					 particle->isFinalState(),
 					 _progenitor->highestpT(),
 					 fb.ids, fb.kinematics->z(),
 					 fb.kinematics->scale(),
 					 fb.kinematics->pT()))
 	return true;
     }
   }
   // veto on maximum pt
   if(fb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
   // general vetos
   if (fb.kinematics && !_vetoes.empty()) {
     bool vetoed=false;
     for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
 	 v != _vetoes.end(); ++v) {
       bool test = (**v).vetoTimeLike(_progenitor,particle,fb,currentTree());
       switch((**v).vetoType()) {
       case ShowerVeto::Emission:
 	vetoed |= test;
 	break;
       case ShowerVeto::Shower:
 	if(test) throw VetoShower();
 	break;
       case ShowerVeto::Event:
 	if(test) throw Veto();
 	break;
       }
     }
     if(vetoed) return true;
   }
   if ( firstInteraction() &&
        profileScales() ) {
     double weight = 
       profileScales()->
       hardScaleProfile(_progenitor->hardScale(),fb.kinematics->pT());
     if ( UseRandom::rnd() > weight )
       return true;
   }
   return false;
 }
 
 bool QTildeShowerHandler::spaceLikeVetoed(const Branching & bb,
 			      ShowerParticlePtr particle) {
   // work out type of interaction
   ShowerInteraction type = convertInteraction(bb.type);
   // check whether emission was harder than largest pt of hard subprocess
   if (restrictPhasespace() && bb.kinematics->pT() > _progenitor->maxHardPt())
     return true;
   // apply the soft correction
   if( softMEC() && _hardme && _hardme->hasMECorrection() ) {
     if(_hardme->softMatrixElementVeto(particle,
 				      _progenitor->progenitor(),
 				      particle->isFinalState(),
 				      _progenitor->highestpT(),
 				      bb.ids, bb.kinematics->z(),
 				      bb.kinematics->scale(),
 				      bb.kinematics->pT()))
       return true;
   }
   // the more general vetos
 
   // check vs max pt for the shower
   if(bb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
 
   if (!_vetoes.empty()) {
     bool vetoed=false;
     for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
 	 v != _vetoes.end(); ++v) {
       bool test = (**v).vetoSpaceLike(_progenitor,particle,bb,currentTree());
       switch ((**v).vetoType()) {
       case ShowerVeto::Emission:
 	vetoed |= test;
 	break;
       case ShowerVeto::Shower:
 	if(test) throw VetoShower();
 	break;
       case ShowerVeto::Event:
 	if(test) throw Veto();
 	break;
       }
     }
     if (vetoed) return true;
   }
   if ( firstInteraction() &&
        profileScales() ) {
     double weight = 
       profileScales()->
       hardScaleProfile(_progenitor->hardScale(),bb.kinematics->pT());
     if ( UseRandom::rnd() > weight )
       return true;
   }
   return false;
 }
 
 bool QTildeShowerHandler::spaceLikeDecayVetoed( const Branching & fb,
 				    ShowerParticlePtr particle) {
   // work out type of interaction
   ShowerInteraction type = convertInteraction(fb.type);
   // apply the soft correction
   if( softMEC() && _decayme && _decayme->hasMECorrection() ) {
     if(_decayme->softMatrixElementVeto(particle,
 				       _progenitor->progenitor(),
 				       particle->isFinalState(),
 				       _progenitor->highestpT(),
 				       fb.ids, fb.kinematics->z(),
 				       fb.kinematics->scale(),
 				       fb.kinematics->pT()))
       return true;
   }
   // veto on hardest pt in the shower
   if(fb.kinematics->pT()> _progenitor->maximumpT(type)) return true;
   // general vetos
   if (!_vetoes.empty()) {
     bool vetoed=false;
     for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
 	 v != _vetoes.end(); ++v) {
       bool test = (**v).vetoSpaceLike(_progenitor,particle,fb,currentTree());
       switch((**v).vetoType()) {
       case ShowerVeto::Emission:
 	vetoed |= test;
 	break;
       case ShowerVeto::Shower:
 	if(test) throw VetoShower();
 	break;
       case ShowerVeto::Event:
 	if(test) throw Veto();
 	break;
       }
       if (vetoed) return true;
     }
   }
   return false;
 }
 
 void QTildeShowerHandler::hardestEmission(bool hard) {
   HardTreePtr ISRTree;
   // internal POWHEG in production or decay
   if( (( _hardme  &&  _hardme->hasPOWHEGCorrection()!=0 ) ||
        ( _decayme && _decayme->hasPOWHEGCorrection()!=0 ) ) ) {
     RealEmissionProcessPtr real;
     unsigned int type(0);
     // production
     if(_hardme) {
       assert(hard);
       real = _hardme->generateHardest( currentTree()->perturbativeProcess(),
 				       interaction_);
       type = _hardme->hasPOWHEGCorrection();
     }
     // decay
     else {
       assert(!hard);
       real = _decayme->generateHardest( currentTree()->perturbativeProcess() );
       type = _decayme->hasPOWHEGCorrection();
     }
     if(real) {
       // set up ther hard tree
       if(!real->outgoing().empty()) _hardtree = new_ptr(HardTree(real));
       // set up the vetos
       currentTree()->setVetoes(real->pT(),type);
     }
     // store initial state POWHEG radiation
     if(_hardtree && _hardme && _hardme->hasPOWHEGCorrection()==1) 
       ISRTree = _hardtree;
   }
   else if (hard) {
     // Get minimum pT cutoff used in shower approximation
     Energy maxpt = 1.*GeV;
 
     if ( currentTree()->showerApproximation() ) {
       int colouredIn  = 0;
       int colouredOut = 0;
       for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it
 	     = currentTree()->outgoingLines().begin();
 	   it != currentTree()->outgoingLines().end(); ++it ) {
 	if( it->second->coloured() ) ++colouredOut;
       }
       for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
 	     = currentTree()->incomingLines().begin();
 	   it != currentTree()->incomingLines().end(); ++it ) {
 	if( it->second->coloured() ) ++colouredIn;
       }
       if ( currentTree()->showerApproximation()->ffPtCut() == currentTree()->showerApproximation()->fiPtCut() &&
 	   currentTree()->showerApproximation()->ffPtCut() == currentTree()->showerApproximation()->iiPtCut() )
 	maxpt = currentTree()->showerApproximation()->ffPtCut();
       else if ( colouredIn == 2 && colouredOut == 0 )
 	maxpt = currentTree()->showerApproximation()->iiPtCut();
       else if ( colouredIn == 0 && colouredOut > 1 )
 	maxpt = currentTree()->showerApproximation()->ffPtCut();
       else if ( colouredIn == 2 && colouredOut == 1 )
 	maxpt = min(currentTree()->showerApproximation()->iiPtCut(), currentTree()->showerApproximation()->fiPtCut());
       else if ( colouredIn == 1 && colouredOut > 1 )
 	maxpt = min(currentTree()->showerApproximation()->ffPtCut(), currentTree()->showerApproximation()->fiPtCut());
       else
 	maxpt = min(min(currentTree()->showerApproximation()->iiPtCut(), currentTree()->showerApproximation()->fiPtCut()), 
 		    currentTree()->showerApproximation()->ffPtCut());
     }
 
     // Generate hardtree from born and real emission subprocesses
     _hardtree = generateCKKW(currentTree());
 
     // Find transverse momentum of hardest emission
     if (_hardtree){
       for(set<HardBranchingPtr>::iterator it=_hardtree->branchings().begin();
      	  it!=_hardtree->branchings().end();++it) {
       	if ((*it)->parent() && (*it)->status()==HardBranching::Incoming)
       	  maxpt=(*it)->branchingParticle()->momentum().perp();
       	if ((*it)->children().size()==2 && (*it)->status()==HardBranching::Outgoing){
 	  if ((*it)->branchingParticle()->id()!=21 &&
 	      abs((*it)->branchingParticle()->id())>5 ){
 	    if ((*it)->children()[0]->branchingParticle()->id()==21 ||
 		abs((*it)->children()[0]->branchingParticle()->id())<6)
 	      maxpt=(*it)->children()[0]->branchingParticle()->momentum().perp();
 	    else if ((*it)->children()[1]->branchingParticle()->id()==21 ||
 		     abs((*it)->children()[1]->branchingParticle()->id())<6)
 	      maxpt=(*it)->children()[1]->branchingParticle()->momentum().perp();
 	  }
 	  else {
 	    if ( abs((*it)->branchingParticle()->id())<6){
 	      if (abs((*it)->children()[0]->branchingParticle()->id())<6)
 		maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp();
 	      else 
 		maxpt = (*it)->children()[0]->branchingParticle()->momentum().perp();
 	    }
 	    else maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp();
 	  }
       	}
       } 
     }
      
     
     // Hardest (pt) emission should be the first powheg emission.
     maxpt=min(sqrt(lastXCombPtr()->lastShowerScale()),maxpt);
 
     // set maximum pT for subsequent emissions from S events
     if ( currentTree()->isPowhegSEvent() ) {
       for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it
 	     = currentTree()->outgoingLines().begin(); 
 	   it != currentTree()->outgoingLines().end(); ++it ) {
 	if( ! it->second->coloured() ) continue;
 	it->first->maximumpT(maxpt, ShowerInteraction::QCD  );
       }  
       for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
 	     = currentTree()->incomingLines().begin(); 
 	   it != currentTree()->incomingLines().end(); ++it ) {
 	if( ! it->second->coloured() ) continue;
 	it->first->maximumpT(maxpt, ShowerInteraction::QCD );
       }
     }
   }
   else 
     _hardtree = generateCKKW(currentTree());
 
   // if hard me doesn't have a FSR powheg 
   // correction use decay powheg correction
   if (_hardme && _hardme->hasPOWHEGCorrection()<2) {
     addFSRUsingDecayPOWHEG(ISRTree);
   }
   // connect the trees
   if(_hardtree) {
     connectTrees(currentTree(),_hardtree,hard); 
   }
 }
 
 void QTildeShowerHandler::addFSRUsingDecayPOWHEG(HardTreePtr ISRTree) {
   // check for intermediate colour singlet resonance
   const ParticleVector inter =  _hardme->subProcess()->intermediates();
   if (inter.size()!=1 || inter[0]->momentum().m2()/GeV2 < 0 || 
       inter[0]->dataPtr()->iColour()!=PDT::Colour0) {
     return;
   }
    
   // ignore cases where outgoing particles are not coloured
   map<ShowerProgenitorPtr, tShowerParticlePtr > out = currentTree()->outgoingLines();
   if (out.size() != 2 ||
       out. begin()->second->dataPtr()->iColour()==PDT::Colour0 ||
       out.rbegin()->second->dataPtr()->iColour()==PDT::Colour0) {
     return;
   }
 
   // look up decay mode
   tDMPtr dm;
   string tag;
   string inParticle = inter[0]->dataPtr()->name() + "->";
   vector<string> outParticles;
   outParticles.push_back(out.begin ()->first->progenitor()->dataPtr()->name());
   outParticles.push_back(out.rbegin()->first->progenitor()->dataPtr()->name());
   for (int it=0; it<2; ++it){
     tag = inParticle + outParticles[it] + "," + outParticles[(it+1)%2] + ";";
     dm = generator()->findDecayMode(tag);
     if(dm) break;
   }
 
   // get the decayer
   HwDecayerBasePtr decayer;   
   if(dm) decayer = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
   // check if decayer has a FSR POWHEG correction 
   if (!decayer || decayer->hasPOWHEGCorrection()<2) {
     return;
   }
   // generate the hardest emission
   // create RealEmissionProcess
   PPtr in = new_ptr(*inter[0]);
   RealEmissionProcessPtr newProcess(new_ptr(RealEmissionProcess()));
   newProcess->bornIncoming().push_back(in);
   newProcess->bornOutgoing().push_back(out.begin ()->first->progenitor());
   newProcess->bornOutgoing().push_back(out.rbegin()->first->progenitor());
   // generate the FSR
   newProcess = decayer->generateHardest(newProcess);
   HardTreePtr FSRTree;
   if(newProcess) {
     // set up ther hard tree
     if(!newProcess->outgoing().empty()) FSRTree = new_ptr(HardTree(newProcess));
     // set up the vetos
     currentTree()->setVetoes(newProcess->pT(),2);
   }
   if(!FSRTree) return;
   
   // if there is no ISRTree make _hardtree from FSRTree
   if (!ISRTree){
     vector<HardBranchingPtr> inBranch,hardBranch;
     for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
 	  cit =currentTree()->incomingLines().begin();
 	cit!=currentTree()->incomingLines().end();++cit ) {
       inBranch.push_back(new_ptr(HardBranching(cit->second,SudakovPtr(),
 					       HardBranchingPtr(),
 					       HardBranching::Incoming)));
       inBranch.back()->beam(cit->first->original()->parents()[0]);
       hardBranch.push_back(inBranch.back());
     }
     if(inBranch[0]->branchingParticle()->dataPtr()->coloured()) {
       inBranch[0]->colourPartner(inBranch[1]);
       inBranch[1]->colourPartner(inBranch[0]);
     }
     for(set<HardBranchingPtr>::iterator it=FSRTree->branchings().begin();
 	it!=FSRTree->branchings().end();++it) {
       if((**it).branchingParticle()->id()!=in->id()) 
 	hardBranch.push_back(*it);
     } 
     hardBranch[2]->colourPartner(hardBranch[3]);
     hardBranch[3]->colourPartner(hardBranch[2]);
     HardTreePtr newTree = new_ptr(HardTree(hardBranch,inBranch,
 					   ShowerInteraction::QCD));            
     _hardtree = newTree;    
   }
   
   // Otherwise modify the ISRTree to include the emission in FSRTree
   else {
     vector<tShowerParticlePtr> FSROut, ISROut;   
     set<HardBranchingPtr>::iterator itFSR, itISR;
     // get outgoing particles 
     for(itFSR =FSRTree->branchings().begin();
 	itFSR!=FSRTree->branchings().end();++itFSR){
       if ((**itFSR).status()==HardBranching::Outgoing) 
 	FSROut.push_back((*itFSR)->branchingParticle());
     }     
     for(itISR =ISRTree->branchings().begin();
 	itISR!=ISRTree->branchings().end();++itISR){
       if ((**itISR).status()==HardBranching::Outgoing) 
 	ISROut.push_back((*itISR)->branchingParticle());
     }
     
     // find COM frame formed by outgoing particles
     LorentzRotation eventFrameFSR, eventFrameISR;
     eventFrameFSR = ((FSROut[0]->momentum()+FSROut[1]->momentum()).findBoostToCM());  
     eventFrameISR = ((ISROut[0]->momentum()+ISROut[1]->momentum()).findBoostToCM());
 
     // find rotation between ISR and FSR frames
     int j=0;
     if (ISROut[0]->id()!=FSROut[0]->id()) j=1;
     eventFrameISR.rotateZ( (eventFrameFSR*FSROut[0]->momentum()).phi()-
 			   (eventFrameISR*ISROut[j]->momentum()).phi() );
     eventFrameISR.rotateY( (eventFrameFSR*FSROut[0]->momentum()).theta()-
 			   (eventFrameISR*ISROut[j]->momentum()).theta() );
     eventFrameISR.invert();
     
     for (itFSR=FSRTree->branchings().begin();
 	 itFSR!=FSRTree->branchings().end();++itFSR){
       if ((**itFSR).branchingParticle()->id()==in->id()) continue;
       for (itISR =ISRTree->branchings().begin();
 	   itISR!=ISRTree->branchings().end();++itISR){
 	if ((**itISR).status()==HardBranching::Incoming) continue;
 	if ((**itFSR).branchingParticle()->id()==
 	    (**itISR).branchingParticle()->id()){
 	  // rotate FSRTree particle to ISRTree event frame
 	  (**itISR).branchingParticle()->setMomentum(eventFrameISR*
 						     eventFrameFSR*
 						     (**itFSR).branchingParticle()->momentum());
 	  (**itISR).branchingParticle()->rescaleMass();
 	  // add the children of the FSRTree particles to the ISRTree
 	  if(!(**itFSR).children().empty()){
 	    (**itISR).addChild((**itFSR).children()[0]);
 	    (**itISR).addChild((**itFSR).children()[1]);
 	    // rotate momenta to ISRTree event frame
 	    (**itISR).children()[0]->branchingParticle()->setMomentum(eventFrameISR*
 								      eventFrameFSR*
 								      (**itFSR).children()[0]->branchingParticle()->momentum());
 	    (**itISR).children()[1]->branchingParticle()->setMomentum(eventFrameISR*
 								      eventFrameFSR*
 								      (**itFSR).children()[1]->branchingParticle()->momentum());
 	  }
 	}
       }	
     }
     _hardtree = ISRTree;
   }
 }
 
 bool QTildeShowerHandler::truncatedTimeLikeShower(tShowerParticlePtr particle,
 				      HardBranchingPtr branch,
 				      ShowerInteraction type,
 				      Branching fb, bool first) {
   // select a branching if we don't have one
   if(!fb.kinematics)
     fb = selectTimeLikeBranching(particle,type,branch);
   // must be an emission, the forced one it not a truncated one
   assert(fb.kinematics);
   ShowerParticleVector children;
   Branching fc[2] = {Branching(),Branching()};
   // Assign the shower kinematics to the emitting particle.
   particle->showerKinematics(fb.kinematics);
   if(fb.kinematics->pT()>progenitor()->highestpT())
     progenitor()->highestpT(fb.kinematics->pT());
   // create the children
   children = createTimeLikeChildren(particle,fb.ids);
   // update the children
   particle->showerKinematics()->
     updateChildren(particle, children,fb.type);
   // select branchings for children
   if(!fc[0].kinematics) {
     // select branching for first particle
     if(!fb.hard && fb.iout ==1 )
       fc[0] = selectTimeLikeBranching(children[0],type,branch);
     else if(fb.hard && !branch->children()[0]->children().empty() )
       fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]);
     else
       fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
   }
   // select branching for the second particle
   if(!fc[1].kinematics) {
     // select branching for first particle
     if(!fb.hard && fb.iout ==2 )
       fc[1] = selectTimeLikeBranching(children[1],type,branch);
     else if(fb.hard && !branch->children()[1]->children().empty() )
       fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]);
     else
       fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
   }
   // shower the first  particle
   if(fc[0].kinematics) {
     // the parent has truncated emission and following line
     if(!fb.hard && fb.iout == 1)
       truncatedTimeLikeShower(children[0],branch,type,fc[0],false);
     // hard emission and subsquent hard emissions
     else if(fb.hard && !branch->children()[0]->children().empty() )
       truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
     // normal shower
     else
       timeLikeShower(children[0],type,fc[0],false);
   }
   if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
   // shower the second particle
   if(fc[1].kinematics) {
     // the parent has truncated emission and following line
     if(!fb.hard && fb.iout == 2)
       truncatedTimeLikeShower(children[1],branch,type,fc[1],false);
     // hard emission and subsquent hard emissions
     else if(fb.hard && !branch->children()[1]->children().empty() )
       truncatedTimeLikeShower(children[1],branch->children()[1],type,fc[1],false);
     else
       timeLikeShower(children[1],type,fc[1],false);
   }
   if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
   // branching has happened
   particle->showerKinematics()->updateParent(particle, children,_evolutionScheme,fb.type); 
   if(first&&!children.empty())
     particle->showerKinematics()->resetChildren(particle,children);
   if(particle->spinInfo()) particle->spinInfo()->develop();
   // TODO NEED AN UPDATE HERE FOR RECONOPT!=0 ?
   return true;
 }
 
 bool QTildeShowerHandler::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam,
 				       HardBranchingPtr branch,
 				       ShowerInteraction type) {
   tcPDFPtr pdf;
   if(firstPDF().particle()  == beamParticle())
     pdf = firstPDF().pdf();
   if(secondPDF().particle() == beamParticle())
     pdf = secondPDF().pdf();
   Energy freeze = pdfFreezingScale();
   Branching bb;
   // parameters of the force branching
   double z(0.);
   HardBranchingPtr timelike;
   for( unsigned int ix = 0; ix < branch->children().size(); ++ix ) {
     if( branch->children()[ix]->status() ==HardBranching::Outgoing) {
       timelike = branch->children()[ix];
     }
     if( branch->children()[ix]->status() ==HardBranching::Incoming )
       z = branch->children()[ix]->z();
   }
   // generate truncated branching
   tcPDPtr part[2];
   if(z>=0.&&z<=1.) {
     while (true) {
       if( !isTruncatedShowerON() || hardOnly() ) break;
       bb = splittingGenerator()->chooseBackwardBranching( *particle, 
 							  beam, 1., beamParticle(), 
 							  type , pdf,freeze);
       if( !bb.kinematics || bb.kinematics->scale() < branch->scale() ) {
 	bb = Branching();
 	break;
       }
       // particles as in Sudakov form factor
       part[0] = bb.ids[0];
       part[1] = bb.ids[2];
       double zsplit = bb.kinematics->z();
       // apply the vetos for the truncated shower
       // if doesn't carry most of momentum
       ShowerInteraction type2 = convertInteraction(bb.type);
       if(type2==branch->sudakov()->interactionType() &&
 	 zsplit < 0.5) {
 	particle->vetoEmission(bb.type,bb.kinematics->scale());
 	continue;
       }
       // others
       if( part[0]->id() != particle->id() || // if particle changes type
 	  bb.kinematics->pT() > progenitor()->maximumpT(type2) ||   // pt veto
 	  bb.kinematics->scale() < branch->scale()) { // angular ordering veto
 	particle->vetoEmission(bb.type,bb.kinematics->scale());
 	continue;
       }
       // and those from the base class
       if(spaceLikeVetoed(bb,particle)) {
 	particle->vetoEmission(bb.type,bb.kinematics->scale());
 	continue;
       }
       break;
     }
   }
   if( !bb.kinematics ) {
     //do the hard emission
     ShoKinPtr kinematics = new_ptr(IS_QTildeShowerKinematics1to2(
     	branch->scale(), z, branch->phi(),
     	branch->children()[0]->pT(), branch->sudakov() 
     	));
     // assign the splitting function and shower kinematics
     particle->showerKinematics( kinematics );
     if(kinematics->pT()>progenitor()->highestpT())
       progenitor()->highestpT(kinematics->pT());
     // For the time being we are considering only 1->2 branching
     // Now create the actual particles, make the otherChild a final state
     // particle, while the newParent is not
     ShowerParticlePtr newParent = 
       new_ptr( ShowerParticle( branch->branchingParticle()->dataPtr(), false ) );
     ShowerParticlePtr otherChild = 
       new_ptr( ShowerParticle( timelike->branchingParticle()->dataPtr(),
     			       true, true ) );
     ShowerParticleVector theChildren;
     theChildren.push_back( particle ); 
     theChildren.push_back( otherChild );
     particle->showerKinematics()->
       updateParent( newParent, theChildren,_evolutionScheme, branch->type());
     // update the history if needed
     currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
     currentTree()->addInitialStateBranching( particle, newParent, otherChild );
     // for the reconstruction of kinematics, parent/child
     // relationships are according to the branching process:
     // now continue the shower
     bool emitted=false;
     if(!hardOnly()) {
       if( branch->parent() ) {
     	emitted = truncatedSpaceLikeShower( newParent, beam, branch->parent() , type);
       }
       else {
     	emitted = spaceLikeShower( newParent, beam , type);
       }
     }
     if( !emitted ) {
       if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
     	kinematics->updateLast( newParent, ZERO, ZERO );
       }
       else {
     	pair<Energy,double> kt = intrinsicpT()[progenitor()];
     	kinematics->updateLast( newParent,
     				kt.first*cos( kt.second ),
     				kt.first*sin( kt.second ) );
       }
     }
     particle->showerKinematics()->
       updateChildren( newParent, theChildren,bb.type);
     if(hardOnly()) return true;
     // perform the shower of the final-state particle
     if( timelike->children().empty() ) {
       timeLikeShower( otherChild , type,Branching(),true);
     }
     else {
       truncatedTimeLikeShower( otherChild, timelike , type,Branching(), true);
     }
     updateHistory(otherChild);
     // return the emitted
     return true;
   }
   // assign the splitting function and shower kinematics
   particle->showerKinematics( bb.kinematics );
   if(bb.kinematics->pT()>progenitor()->highestpT())
     progenitor()->highestpT(bb.kinematics->pT());
   // For the time being we are considering only 1->2 branching
   // Now create the actual particles, make the otherChild a final state
   // particle, while the newParent is not
   ShowerParticlePtr newParent = new_ptr( ShowerParticle( part[0], false ) );
   ShowerParticlePtr otherChild = new_ptr( ShowerParticle( part[1], true, true ) );
   ShowerParticleVector theChildren; 
   theChildren.push_back( particle ); 
   theChildren.push_back( otherChild );
   particle->showerKinematics()->
     updateParent( newParent, theChildren,_evolutionScheme, bb.type);
   // update the history if needed
   currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
   currentTree()->addInitialStateBranching( particle, newParent, otherChild );
   // for the reconstruction of kinematics, parent/child
   // relationships are according to the branching process:
   // now continue the shower
   bool emitted = truncatedSpaceLikeShower( newParent, beam, branch,type);
   // now reconstruct the momentum
   if( !emitted ) {
     if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
       bb.kinematics->updateLast( newParent, ZERO, ZERO );
     }
     else {
       pair<Energy,double> kt = intrinsicpT()[ progenitor() ];
       bb.kinematics->updateLast( newParent,
 				 kt.first*cos( kt.second ),
 				 kt.first*sin( kt.second ) );
     }
   }
   particle->showerKinematics()->
     updateChildren( newParent, theChildren, bb.type);
   // perform the shower of the final-state particle
   timeLikeShower( otherChild , type,Branching(),true);
   updateHistory(otherChild);
   // return the emitted
   return true;
 }
 
 bool QTildeShowerHandler::
 truncatedSpaceLikeDecayShower(tShowerParticlePtr particle, 
 			      const ShowerParticle::EvolutionScales & maxScales,
 			      Energy minmass, HardBranchingPtr branch,
 			      ShowerInteraction type, Branching fb) {
   // select a branching if we don't have one
   if(!fb.kinematics)
     fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,branch);
   // must be an emission, the forced one it not a truncated one
   assert(fb.kinematics);
   ShowerParticleVector children;
   Branching fc[2]={Branching(),Branching()};
   // Assign the shower kinematics to the emitting particle.
   particle->showerKinematics(fb.kinematics);
   if(fb.kinematics->pT()>progenitor()->highestpT())
     progenitor()->highestpT(fb.kinematics->pT());
   // create the ShowerParticle objects for the two children
   children = createTimeLikeChildren(particle,fb.ids);
   // updateChildren the children
   particle->showerKinematics()->
     updateChildren(particle, children, fb.type);
   // select branchings for children
   if(!fc[0].kinematics) {
     if(children[0]->id()==particle->id()) {
       // select branching for first particle
       if(!fb.hard)
 	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,branch);
       else if(fb.hard && ! branch->children()[0]->children().empty() )
 	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,
 					      branch->children()[0]);
       else
 	fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,
 					      HardBranchingPtr());
     }
     else {
       // select branching for first particle
       if(fb.hard && !branch->children()[0]->children().empty() )
 	fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]);
       else
 	fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
     }
   }
   // select branching for the second particle
   if(!fc[1].kinematics) {
     if(children[1]->id()==particle->id()) {
       // select branching for first particle
       if(!fb.hard)
 	   fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,branch);
       else if(fb.hard && ! branch->children()[1]->children().empty() )
 	fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,
 					      branch->children()[1]);
       else
 	fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,
 					      HardBranchingPtr());
     }
     else {
       if(fb.hard && !branch->children()[1]->children().empty() )
 	fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]);
       else
 	fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
     }
   }
   // old default
   // update the history if needed
   currentTree()->updateInitialStateShowerProduct(progenitor(),children[0]);
   currentTree()->addInitialStateBranching(particle,children[0],children[1]);
   // shower the first  particle
   if(fc[0].kinematics) {
     if(children[0]->id()==particle->id()) {
       if(!fb.hard)
 	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
 				       branch,type,fc[0]);
       else if(fb.hard && ! branch->children()[0]->children().empty() )
 	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
 				       branch->children()[0],type,fc[0]);
       else
 	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[0]);
     }
     else {
       if(fb.hard && !branch->children()[0]->children().empty() )
 	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
       // normal shower
       else
 	timeLikeShower(children[0],type,fc[0],false);
     }
   }
   // shower the second  particle
   if(fc[1].kinematics) {
     if(children[0]->id()==particle->id()) {
       if(!fb.hard)
 	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
 				       branch,type,fc[1]);
       else if(fb.hard && ! branch->children()[0]->children().empty() )
 	truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
 				       branch->children()[0],type,fc[1]);
       else
 	spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[1]);
     }
     else {
       if(fb.hard && !branch->children()[0]->children().empty() )
 	truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[1],false);
       // normal shower
       else
 	timeLikeShower(children[0],type,fc[1],false);
     }
   }
   updateHistory(children[1]);
   // TODO DO WE NEED A CHECK IF RECONPT !=0
   return true;
 }
 
 void QTildeShowerHandler::connectTrees(ShowerTreePtr showerTree, 
 				       HardTreePtr hardTree, bool hard ) {
   ShowerParticleVector particles;
   // find the Sudakovs
   for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
       cit!=hardTree->branchings().end();++cit) {
     // Sudakovs for ISR
     if((**cit).parent()&&(**cit).status()==HardBranching::Incoming) {
       ++_nis;
       array<long,3> br;
       br[0] = (**cit).parent()->branchingParticle()->id();
       br[1] = (**cit).          branchingParticle()->id();
       br[2] = (**cit).parent()->children()[0]==*cit ?
 	(**cit).parent()->children()[1]->branchingParticle()->id() :
 	(**cit).parent()->children()[0]->branchingParticle()->id();
       BranchingList branchings = splittingGenerator()->initialStateBranchings();
       if(br[1]<0&&br[0]==br[1]) {
 	br[0] = abs(br[0]);
 	br[1] = abs(br[1]);
       }
       else if(br[1]<0) {
 	br[1] = -br[1];
 	br[2] = -br[2];
       }
       long index = abs(br[1]);
       SudakovPtr sudakov;
       for(BranchingList::const_iterator cjt = branchings.lower_bound(index); 
 	  cjt != branchings.upper_bound(index); ++cjt ) {
 	IdList ids = cjt->second.particles;
 	if(ids[0]->id()==br[0]&&ids[1]->id()==br[1]&&ids[2]->id()==br[2]) {
 	  sudakov=cjt->second.sudakov;
 	  break;
 	}
       }
       if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
 				     << "QTildeShowerHandler::connectTrees() for ISR" 
 				     << Exception::runerror;
       (**cit).parent()->sudakov(sudakov);
     }
     // Sudakovs for FSR
     else if(!(**cit).children().empty()) {
       ++_nfs;
       array<long,3> br;
       br[0] = (**cit)               .branchingParticle()->id();
       br[1] = (**cit).children()[0]->branchingParticle()->id();
       br[2] = (**cit).children()[1]->branchingParticle()->id();
       BranchingList branchings = splittingGenerator()->finalStateBranchings();
       if(br[0]<0) {
 	br[0] = abs(br[0]);
 	br[1] = abs(br[1]);
 	br[2] = abs(br[2]);
       }
       long index = br[0];
       SudakovPtr sudakov;
       for(BranchingList::const_iterator cjt = branchings.lower_bound(index); 
 	  cjt != branchings.upper_bound(index); ++cjt ) {
 	IdList ids = cjt->second.particles;
 	if(ids[0]->id()==br[0]&&ids[1]->id()==br[1]&&ids[2]->id()==br[2]) {
 	  sudakov=cjt->second.sudakov;
 	  break;
 	}
       }
       if(!sudakov) {
 	throw Exception() << "Can't find Sudakov for the hard emission in "
 			  << "QTildeShowerHandler::connectTrees()" 
 			  << Exception::runerror;
       }
       (**cit).sudakov(sudakov);
     }
   }
   // calculate the evolution scale
   for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
       cit!=hardTree->branchings().end();++cit) {
     particles.push_back((*cit)->branchingParticle());
   }
   partnerFinder()->
     setInitialEvolutionScales(particles,!hard,interaction_,true);
   hardTree->partnersSet(true);
   // inverse reconstruction
   if(hard) {
     kinematicsReconstructor()->
       deconstructHardJets(hardTree,interaction_);
   }
   else
     kinematicsReconstructor()->
       deconstructDecayJets(hardTree,interaction_);
   // now reset the momenta of the showering particles
   vector<ShowerProgenitorPtr> particlesToShower=showerTree->extractProgenitors();
   // match them
   map<ShowerProgenitorPtr,HardBranchingPtr> partners;
   for(set<HardBranchingPtr>::const_iterator bit=hardTree->branchings().begin();
       bit!=hardTree->branchings().end();++bit) {
     Energy2 dmin( 1e30*GeV2 );
     ShowerProgenitorPtr partner;
     for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
 	pit!=particlesToShower.end();++pit) {
       if(partners.find(*pit)!=partners.end()) continue;
       if( (**bit).branchingParticle()->id() !=  (**pit).progenitor()->id() ) continue;
       if( (**bit).branchingParticle()->isFinalState() !=
 	  (**pit).progenitor()->isFinalState() ) continue;
       if( (**pit).progenitor()->isFinalState() ) {
 	Energy2 dtest =
 	  sqr( (**pit).progenitor()->momentum().x() - (**bit).showerMomentum().x() ) +
 	  sqr( (**pit).progenitor()->momentum().y() - (**bit).showerMomentum().y() ) +
 	  sqr( (**pit).progenitor()->momentum().z() - (**bit).showerMomentum().z() ) +
 	  sqr( (**pit).progenitor()->momentum().t() - (**bit).showerMomentum().t() );
 	// add mass difference for identical particles (e.g. Z0 Z0 production)
 	dtest += 1e10*sqr((**pit).progenitor()->momentum().m()-(**bit).showerMomentum().m());
 	if( dtest < dmin ) {
 	  partner = *pit;
 	  dmin = dtest;
 	}
       }
       else {
 	// ensure directions are right
 	if((**pit).progenitor()->momentum().z()/(**bit).showerMomentum().z()>ZERO) {
 	  partner = *pit;
 	  break;
 	}
       }
     }
     if(!partner) throw Exception() << "Failed to match shower and hard trees in QTildeShowerHandler::hardestEmission"
 				   << Exception::eventerror;
     partners[partner] = *bit;
   }
   for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
       pit!=particlesToShower.end();++pit) {
     HardBranchingPtr partner = partners[*pit];
     if((**pit).progenitor()->dataPtr()->stable()) {
       (**pit).progenitor()->set5Momentum(partner->showerMomentum());
       (**pit).copy()->set5Momentum(partner->showerMomentum());
     }
     else {
       Lorentz5Momentum oldMomentum = (**pit).progenitor()->momentum();
       Lorentz5Momentum newMomentum = partner->showerMomentum();
       LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
       (**pit).progenitor()->transform(boost);
       (**pit).copy()      ->transform(boost);
       boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
       (**pit).progenitor()->transform(boost);
       (**pit).copy()      ->transform(boost);
     }
   }
   // correction boosts for daughter trees
   for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator 
 	tit  = showerTree->treelinks().begin();
       tit != showerTree->treelinks().end();++tit) {
     ShowerTreePtr decayTree = tit->first;
     map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator 
       cit = decayTree->incomingLines().begin();
     // reset the momentum of the decay particle
     Lorentz5Momentum oldMomentum = cit->first->progenitor()->momentum();
     Lorentz5Momentum newMomentum = tit->second.second->momentum();
     LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
     decayTree->transform(boost,true);
     boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
     decayTree->transform(boost,true);
   }
 }
 
 void QTildeShowerHandler::doShowering(bool hard,XCPtr xcomb) {
   // zero number of emissions
   _nis = _nfs = 0;
   // if MC@NLO H event and limited emissions
   // indicate both final and initial state emission
   if ( currentTree()->isMCatNLOHEvent() && _limitEmissions != 0 ) {
     _nis = _nfs = 1;
   }
   // extract particles to shower
   vector<ShowerProgenitorPtr> particlesToShower(setupShower(hard));
   // check if we should shower
   bool colCharge = false;
   for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
     if(particlesToShower[ix]->progenitor()->dataPtr()->coloured() ||
        particlesToShower[ix]->progenitor()->dataPtr()->charged()) {
       colCharge = true;
       break;
     }
   }
   if(!colCharge) {
     _currenttree->hasShowered(true);
     return;
   }
   // setup the maximum scales for the shower
   if (restrictPhasespace()) setupMaximumScales(particlesToShower,xcomb);
   // set the hard scales for the profiles
   setupHardScales(particlesToShower,xcomb);
   // specific stuff for hard processes and decays
   Energy minmass(ZERO), mIn(ZERO);
   // hard process generate the intrinsic p_T once and for all
   if(hard) {
     generateIntrinsicpT(particlesToShower);
   }
   // decay compute the minimum mass of the final-state
   else {
     for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
       if(particlesToShower[ix]->progenitor()->isFinalState()) {
         if(particlesToShower[ix]->progenitor()->dataPtr()->stable()) {
           auto dm=  ShowerHandler::currentHandler()->retConstituentMasses()?
 	    particlesToShower[ix]->progenitor()->dataPtr()->constituentMass():
 	    particlesToShower[ix]->progenitor()->dataPtr()->mass();
           minmass += dm;
 	}
         else
 	  minmass += particlesToShower[ix]->progenitor()->mass();
       }
       else {
 	mIn = particlesToShower[ix]->progenitor()->mass();
       }
     }
     // throw exception if decay can't happen
     if ( minmass > mIn ) {
       throw Exception() << "QTildeShowerHandler.cc: Mass of decaying particle is "
 			<< "below constituent masses of decay products."
 			<< Exception::eventerror;
     }
   }
   // setup for reweighted
   bool reWeighting = _reWeight && hard && ShowerHandler::currentHandler()->firstInteraction();
   double eventWeight=0.;
   unsigned int nTryReWeight(0);
   // create random particle vector (only need to do once)
   vector<ShowerProgenitorPtr> tmp;
   unsigned int nColouredIncoming = 0;
   while(particlesToShower.size()>0){
     unsigned int xx=UseRandom::irnd(particlesToShower.size());
     tmp.push_back(particlesToShower[xx]);
     particlesToShower.erase(particlesToShower.begin()+xx);
   }
   particlesToShower=tmp;
   for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
     if(!particlesToShower[ix]->progenitor()->isFinalState() &&
        particlesToShower[ix]->progenitor()->coloured()) ++nColouredIncoming;
   }
   bool switchRecon = hard && nColouredIncoming !=1;
   // main shower loop
   unsigned int ntry(0);
   bool reconstructed = false;
   do {
     // clear results of last attempt if needed
     if(ntry!=0) {
       currentTree()->clear();
       setEvolutionPartners(hard,interaction_,true);
       _nis = _nfs = 0;
       // if MC@NLO H event and limited emissions
       // indicate both final and initial state emission
       if ( currentTree()->isMCatNLOHEvent() && _limitEmissions != 0 ) {
 	_nis = _nfs = 1;
       }
       for(unsigned int ix=0; ix<particlesToShower.size();++ix) {
 	SpinPtr spin = particlesToShower[ix]->progenitor()->spinInfo();
 	if(spin && spin->decayVertex() &&
 	   dynamic_ptr_cast<tcSVertexPtr>(spin->decayVertex())) {
 	  spin->decayVertex(VertexPtr());
 	}
       }
       for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
 	if(particlesToShower[ix]->progenitor()->isFinalState() ||
 	   (hard && !particlesToShower[ix]->progenitor()->isFinalState())) {
 	  if(particlesToShower[ix]->progenitor()->spinInfo())
 	    particlesToShower[ix]->progenitor()->spinInfo()->reset();
 	}
       }
     }
     // loop over particles
     for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
       // extract the progenitor
       progenitor(particlesToShower[ix]);
       // final-state radiation
       if(progenitor()->progenitor()->isFinalState()) {
 	if(!doFSR()) continue;
 	// perform shower
 	progenitor()->hasEmitted(startTimeLikeShower(interaction_));
       }
       // initial-state radiation
       else {
 	if(!doISR()) continue;
 	// hard process
 	if(hard) {
 	  // get the PDF
 	  setBeamParticle(_progenitor->beam());
 	  if(!beamParticle()) {
 	    throw Exception() << "Incorrect type of beam particle in "
 			      << "QTildeShowerHandler::doShowering(). "
 			      << "This should not happen for conventional choices but may happen if you have used a"
 			      << " non-default choice and have not changed the create ParticleData line in the input files"
 			      << " for this particle to create BeamParticleData."
 			      << Exception::runerror;
 	  }
 	  // perform the shower
 	  // set the beam particle
 	  tPPtr beamparticle=progenitor()->original();
 	  if(!beamparticle->parents().empty()) 
 	    beamparticle=beamparticle->parents()[0];
 	  // generate the shower
 	  progenitor()->hasEmitted(startSpaceLikeShower(beamparticle,
 							interaction_));
 	}
 	// decay
 	else {
 	  // skip colour and electrically neutral particles
 	  if(!progenitor()->progenitor()->dataPtr()->coloured() &&
 	     !progenitor()->progenitor()->dataPtr()->charged()) {
 	    progenitor()->hasEmitted(false);
 	    continue;
 	  }
 	  // perform shower
 	  // set the scales correctly. The current scale is the maximum scale for
 	  // emission not the starting scale
 	  ShowerParticle::EvolutionScales maxScales(progenitor()->progenitor()->scales());
 	  progenitor()->progenitor()->scales() = ShowerParticle::EvolutionScales();
 	  if(progenitor()->progenitor()->dataPtr()->charged()) {
 	    progenitor()->progenitor()->scales().QED      = progenitor()->progenitor()->mass();
 	    progenitor()->progenitor()->scales().QED_noAO = progenitor()->progenitor()->mass();
 	  }
 	  if(progenitor()->progenitor()->hasColour()) {
 	    progenitor()->progenitor()->scales().QCD_c       = progenitor()->progenitor()->mass();
 	    progenitor()->progenitor()->scales().QCD_c_noAO  = progenitor()->progenitor()->mass();
 	  }
 	  if(progenitor()->progenitor()->hasAntiColour()) {
 	    progenitor()->progenitor()->scales().QCD_ac      = progenitor()->progenitor()->mass();
 	    progenitor()->progenitor()->scales().QCD_ac_noAO = progenitor()->progenitor()->mass();
 	  }
 	  // perform the shower
 	  progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass,
 							     interaction_));
 	}
       }
     }
     // do the kinematic reconstruction, checking if it worked
     reconstructed = hard ?
       kinematicsReconstructor()->
       reconstructHardJets (currentTree(),intrinsicpT(),interaction_,
 			   switchRecon && ntry>maximumTries()/2) :
       kinematicsReconstructor()->
       reconstructDecayJets(currentTree(),interaction_);
     if(!reconstructed) continue;
     // apply vetos on the full shower
     for(vector<FullShowerVetoPtr>::const_iterator it=_fullShowerVetoes.begin();
 	it!=_fullShowerVetoes.end();++it) {
       int veto = (**it).applyVeto(currentTree());
       if(veto<0) continue;
       // veto the shower
       if(veto==0) {
 	reconstructed = false;
 	break;
       }
       // veto the shower and reweight
       else if(veto==1) {
 	reconstructed = false;
 	break;
       }
       // veto the event
       else if(veto==2) {
 	throw Veto();
       }
     }
     if(reWeighting) {
       if(reconstructed) eventWeight += 1.;
       reconstructed=false;
       ++nTryReWeight;
       if(nTryReWeight==_nReWeight) {
 	reWeighting = false;
 	if(eventWeight==0.) throw Veto();
       }
     }
   }
   while(!reconstructed&&maximumTries()>++ntry);
   // check if failed to generate the shower
   if(ntry==maximumTries()) {
     if(hard)
       throw ShowerHandler::ShowerTriesVeto(ntry);
     else
       throw Exception() << "Failed to generate the shower after "
 			<< ntry << " attempts in QTildeShowerHandler::showerDecay()"
 			<< Exception::eventerror;
   }
   // handle the weights and apply any reweighting required
   if(nTryReWeight>0) {
     tStdEHPtr seh = dynamic_ptr_cast<tStdEHPtr>(generator()->currentEventHandler());
     static bool first = true;
     if(seh) {
       seh->reweight(eventWeight/double(nTryReWeight));
     }
     else if(first) {
       generator()->log() << "Reweighting the shower only works with internal Herwig7 processes"
     			 << "Presumably you are showering Les Houches Events. These will not be"
     			 << "reweighted\n";
       first = false;
     }
   }
   // tree has now showered
   _currenttree->hasShowered(true);
   hardTree(HardTreePtr());
 }
 
 void QTildeShowerHandler:: convertHardTree(bool hard,ShowerInteraction type) {
   map<ColinePtr,ColinePtr> cmap;
   // incoming particles
   for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator 
 	cit=currentTree()->incomingLines().begin();cit!=currentTree()->incomingLines().end();++cit) {
     map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator 
       mit = hardTree()->particles().find(cit->first->progenitor());
     // put the colour lines in the map
     ShowerParticlePtr oldParticle = cit->first->progenitor();
     ShowerParticlePtr newParticle = mit->second->branchingParticle();
     ColinePtr cLine = oldParticle->    colourLine();
     ColinePtr aLine = oldParticle->antiColourLine();
     if(newParticle->colourLine() &&
        cmap.find(newParticle->    colourLine())==cmap.end())
       cmap[newParticle->    colourLine()] = cLine;
     if(newParticle->antiColourLine() &&
        cmap.find(newParticle->antiColourLine())==cmap.end())
       cmap[newParticle->antiColourLine()] = aLine;
     // check whether or not particle emits
     bool emission = mit->second->parent();
     if(emission) {
       if(newParticle->colourLine()) {
 	ColinePtr ctemp = newParticle->    colourLine();
 	ctemp->removeColoured(newParticle);
       }
       if(newParticle->antiColourLine()) {
 	ColinePtr ctemp = newParticle->antiColourLine();
 	ctemp->removeAntiColoured(newParticle);
       }
       newParticle = mit->second->parent()->branchingParticle();
     }
     // get the new colour lines
     ColinePtr newCLine,newALine;
     // sort out colour lines
     if(newParticle->colourLine()) {
       ColinePtr ctemp = newParticle->    colourLine();
       ctemp->removeColoured(newParticle);
       if(cmap.find(ctemp)!=cmap.end()) {
 	newCLine = cmap[ctemp];
       }
       else {
 	newCLine = new_ptr(ColourLine());
 	cmap[ctemp] = newCLine;
       }
     }
     // and anticolour lines
     if(newParticle->antiColourLine()) {
       ColinePtr ctemp = newParticle->antiColourLine();
       ctemp->removeAntiColoured(newParticle);
       if(cmap.find(ctemp)!=cmap.end()) {
 	newALine = cmap[ctemp];
       }
       else {
 	newALine = new_ptr(ColourLine());
 	cmap[ctemp] = newALine;
       }
     }
     // remove colour lines from old particle
     if(aLine) {
       aLine->removeAntiColoured(cit->first->copy());
       aLine->removeAntiColoured(cit->first->progenitor());
     }
     if(cLine) {
       cLine->removeColoured(cit->first->copy());
       cLine->removeColoured(cit->first->progenitor());
     }
     // add particle to colour lines
     if(newCLine) newCLine->addColoured    (newParticle);
     if(newALine) newALine->addAntiColoured(newParticle);
     // insert new particles
     cit->first->copy(newParticle);
     ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,false)));
     cit->first->progenitor(sp);
     currentTree()->incomingLines()[cit->first]=sp;
     cit->first->perturbative(!emission);
     // and the emitted particle if needed
     if(emission) {
       ShowerParticlePtr newOut = mit->second->parent()->children()[1]->branchingParticle();
       if(newOut->colourLine()) {
 	ColinePtr ctemp = newOut->    colourLine();
 	ctemp->removeColoured(newOut);
 	assert(cmap.find(ctemp)!=cmap.end());
 	cmap[ctemp]->addColoured    (newOut);
       }
       if(newOut->antiColourLine()) {
 	ColinePtr ctemp = newOut->antiColourLine();
 	ctemp->removeAntiColoured(newOut);
 	assert(cmap.find(ctemp)!=cmap.end());
 	cmap[ctemp]->addAntiColoured(newOut);
       }
       ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true));
       ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout));
       out->perturbative(false);
       currentTree()->outgoingLines().insert(make_pair(out,sout));
     }
     if(hard) {
       // sort out the value of x
       if(mit->second->beam()->momentum().z()>ZERO) {
 	sp->x(newParticle->momentum(). plus()/mit->second->beam()->momentum(). plus());
       }
       else {
 	sp->x(newParticle->momentum().minus()/mit->second->beam()->momentum().minus());
       }
     }
   }
   // outgoing particles
   for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator 
 	cit=currentTree()->outgoingLines().begin();cit!=currentTree()->outgoingLines().end();++cit) {
     map<tShowerTreePtr,pair<tShowerProgenitorPtr,
 			    tShowerParticlePtr> >::const_iterator tit;
     for(tit  = currentTree()->treelinks().begin();
 	tit != currentTree()->treelinks().end();++tit) {
       if(tit->second.first && tit->second.second==cit->first->progenitor())
 	break;
     }
     map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator 
       mit = hardTree()->particles().find(cit->first->progenitor());
     if(mit==hardTree()->particles().end()) continue;
     // put the colour lines in the map
     ShowerParticlePtr oldParticle = cit->first->progenitor();
     ShowerParticlePtr newParticle = mit->second->branchingParticle();
     ShowerParticlePtr newOut;
     ColinePtr cLine = oldParticle->    colourLine();
     ColinePtr aLine = oldParticle->antiColourLine();
     if(newParticle->colourLine() &&
        cmap.find(newParticle->    colourLine())==cmap.end())
       cmap[newParticle->    colourLine()] = cLine;
     if(newParticle->antiColourLine() &&
        cmap.find(newParticle->antiColourLine())==cmap.end())
       cmap[newParticle->antiColourLine()] = aLine;
     // check whether or not particle emits
     bool emission = !mit->second->children().empty();
     if(emission) {
       if(newParticle->colourLine()) {
 	ColinePtr ctemp = newParticle->    colourLine();
 	ctemp->removeColoured(newParticle);
       }
       if(newParticle->antiColourLine()) {
 	ColinePtr ctemp = newParticle->antiColourLine();
 	ctemp->removeAntiColoured(newParticle);
       }
       newParticle = mit->second->children()[0]->branchingParticle();
       newOut      = mit->second->children()[1]->branchingParticle();
       if(newParticle->id()!=oldParticle->id()&&newParticle->id()==newOut->id())
 	swap(newParticle,newOut);
     }
     // get the new colour lines
     ColinePtr newCLine,newALine;
     // sort out colour lines
     if(newParticle->colourLine()) {
       ColinePtr ctemp = newParticle->    colourLine();
       ctemp->removeColoured(newParticle);
       if(cmap.find(ctemp)!=cmap.end()) {
 	newCLine = cmap[ctemp];
       }
       else {
 	newCLine = new_ptr(ColourLine());
 	cmap[ctemp] = newCLine;
       }
     }
     // and anticolour lines
     if(newParticle->antiColourLine()) {
       ColinePtr ctemp = newParticle->antiColourLine();
       ctemp->removeAntiColoured(newParticle);
       if(cmap.find(ctemp)!=cmap.end()) {
 	newALine = cmap[ctemp];
       }
       else {
 	newALine = new_ptr(ColourLine());
 	cmap[ctemp] = newALine;
       }
     }
     // remove colour lines from old particle
     if(aLine) {
       aLine->removeAntiColoured(cit->first->copy());
       aLine->removeAntiColoured(cit->first->progenitor());
     }
     if(cLine) {
       cLine->removeColoured(cit->first->copy());
       cLine->removeColoured(cit->first->progenitor());
     }
     // special for unstable particles
     if(newParticle->id()==oldParticle->id() &&
        (tit!=currentTree()->treelinks().end()||!oldParticle->dataPtr()->stable())) {
       Lorentz5Momentum oldMomentum = oldParticle->momentum();
       Lorentz5Momentum newMomentum = newParticle->momentum();
       LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
       if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false);
       oldParticle->transform(boost);
       boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
       oldParticle->transform(boost);
       if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false);
       newParticle=oldParticle;
     }
     // add particle to colour lines
     if(newCLine) newCLine->addColoured    (newParticle);
     if(newALine) newALine->addAntiColoured(newParticle);
     // insert new particles
     cit->first->copy(newParticle);
     ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,true)));
     cit->first->progenitor(sp);
     currentTree()->outgoingLines()[cit->first]=sp;
     cit->first->perturbative(!emission);
     // and the emitted particle if needed
     if(emission) {
       if(newOut->colourLine()) {
 	ColinePtr ctemp = newOut->    colourLine();
 	ctemp->removeColoured(newOut);
 	assert(cmap.find(ctemp)!=cmap.end());
 	cmap[ctemp]->addColoured    (newOut);
       }
       if(newOut->antiColourLine()) {
 	ColinePtr ctemp = newOut->antiColourLine();
 	ctemp->removeAntiColoured(newOut);
 	assert(cmap.find(ctemp)!=cmap.end());
 	cmap[ctemp]->addAntiColoured(newOut);
       }
       ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true));
       ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout));
       out->perturbative(false);
       currentTree()->outgoingLines().insert(make_pair(out,sout));
     }
     // update any decay products
     if(tit!=currentTree()->treelinks().end())
       currentTree()->updateLink(tit->first,make_pair(cit->first,sp));
   }
   // reset the tree
   currentTree()->resetShowerProducts();
   // reextract the particles and set the colour partners
   vector<ShowerParticlePtr> particles = 
     currentTree()->extractProgenitorParticles();
   // clear the partners
   for(unsigned int ix=0;ix<particles.size();++ix) {
     particles[ix]->partner(ShowerParticlePtr());
     particles[ix]->clearPartners();
   }
   // clear the tree
   hardTree(HardTreePtr());
   // Set the initial evolution scales
   partnerFinder()->
     setInitialEvolutionScales(particles,!hard,type,!_hardtree);
 }
 
 Branching QTildeShowerHandler::selectTimeLikeBranching(tShowerParticlePtr particle,
 						       ShowerInteraction type,
 						       HardBranchingPtr branch) {
   Branching fb;
   unsigned int iout=0;
   while (true) {
     // break if doing truncated shower and no truncated shower needed
     if(branch && (!isTruncatedShowerON()||hardOnly())) break;
     fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type);
     // no emission break
     if(!fb.kinematics) break;
     // special for truncated shower
     if(branch) {
       // check haven't evolved too far
       if(fb.kinematics->scale() < branch->scale()) {
 	fb=Branching();
 	break;
       }
       // find the truncated line
       iout=0;
       if(fb.ids[1]->id()!=fb.ids[2]->id()) {
 	if(fb.ids[1]->id()==particle->id())       iout=1;
 	else if (fb.ids[2]->id()==particle->id()) iout=2;
       }
       else if(fb.ids[1]->id()==particle->id()) {
 	if(fb.kinematics->z()>0.5) iout=1;
 	else                       iout=2;
       }
       // apply the vetos for the truncated shower
       // no flavour changing branchings
       if(iout==0) {
 	particle->vetoEmission(fb.type,fb.kinematics->scale());
 	continue;
       }
       double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
       // only if same interaction for forced branching
       ShowerInteraction type2 = convertInteraction(fb.type);
       // and evolution
       if(type2==branch->sudakov()->interactionType()) {
 	if(zsplit < 0.5 || // hardest line veto
 	   fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
 	  particle->vetoEmission(fb.type,fb.kinematics->scale());
 	  continue;
 	}
       }
       // pt veto
       if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
 	particle->vetoEmission(fb.type,fb.kinematics->scale());
 	continue;
       }
     }
     // standard vetos for all emissions
     if(timeLikeVetoed(fb,particle)) {
       particle->vetoEmission(fb.type,fb.kinematics->scale());
       if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
       continue;
     }
     // special for already decayed particles
     // don't allow flavour changing branchings
     bool vetoDecay = false;
     for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,
 	  tShowerParticlePtr> >::const_iterator tit  = currentTree()->treelinks().begin();
 	tit != currentTree()->treelinks().end();++tit) {
       if(tit->second.first == progenitor()) {
 	map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator 
 	  it = currentTree()->outgoingLines().find(progenitor());
 	if(it!=currentTree()->outgoingLines().end() && particle == it->second &&
 	   fb.ids[0]!=fb.ids[1] && fb.ids[1]!=fb.ids[2]) {
 	  vetoDecay = true;
 	  break;
 	}
       }
     }
     if(vetoDecay) {
       particle->vetoEmission(fb.type,fb.kinematics->scale());
       if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
       continue;
     }
     break;
   }
   // normal case
   if(!branch) {
     if(fb.kinematics) fb.hard = false;
     return fb;
   }
   // truncated emission
   if(fb.kinematics) {
     fb.hard = false;
     fb.iout = iout;
     return fb;
   }
   // otherwise need to return the hard emission
   // construct the kinematics for the hard emission
   ShoKinPtr showerKin = new_ptr(FS_QTildeShowerKinematics1to2(
     branch->scale(),
     branch->children()[0]->z(),
     branch->phi(),
     branch->children()[0]->pT(),
     branch->sudakov()
     ));						 
 
   IdList idlist(3);
   idlist[0] = particle->dataPtr();
   idlist[1] = branch->children()[0]->branchingParticle()->dataPtr();
   idlist[2] = branch->children()[1]->branchingParticle()->dataPtr();
   fb = Branching( showerKin, idlist, branch->sudakov(),branch->type() );
   fb.hard = true;
   fb.iout=0;
   // return it
   return fb;
 }
 
 Branching QTildeShowerHandler::selectSpaceLikeDecayBranching(tShowerParticlePtr particle,
 						 const ShowerParticle::EvolutionScales & maxScales,
 						 Energy minmass,ShowerInteraction type,
 						 HardBranchingPtr branch) {
   Branching fb;
   unsigned int iout=0;
   while (true) {
     // break if doing truncated shower and no truncated shower needed
     if(branch && (!isTruncatedShowerON()||hardOnly())) break;
     // select branching
     fb=_splittingGenerator->chooseDecayBranching(*particle,maxScales,minmass,
 						 _initialenhance,type);
     // return if no radiation
     if(!fb.kinematics) break;
     // special for truncated shower
     if(branch) {
       // check haven't evolved too far
       if(fb.kinematics->scale() < branch->scale()) {
 	fb=Branching();
 	break;
       }
       // find the truncated line
       iout=0;
       if(fb.ids[1]->id()!=fb.ids[2]->id()) {
 	if(fb.ids[1]->id()==particle->id())       iout=1;
 	else if (fb.ids[2]->id()==particle->id()) iout=2;
       }
       else if(fb.ids[1]->id()==particle->id()) {
 	if(fb.kinematics->z()>0.5) iout=1;
 	else                       iout=2;
       }
       // apply the vetos for the truncated shower
       // no flavour changing branchings
       if(iout==0) {
 	particle->vetoEmission(fb.type,fb.kinematics->scale());
 	continue;
       }
       ShowerInteraction type2 = convertInteraction(fb.type);
       double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
       if(type2==branch->sudakov()->interactionType()) {
 	if(zsplit < 0.5 || // hardest line veto
 	   fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
 	  particle->vetoEmission(fb.type,fb.kinematics->scale());
 	  continue;
 	}
       }
       // pt veto
       if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
 	particle->vetoEmission(fb.type,fb.kinematics->scale());
 	continue;
       }
     }
     // if not vetoed break
     if(spaceLikeDecayVetoed(fb,particle)) {
       // otherwise reset scale and continue
       particle->vetoEmission(fb.type,fb.kinematics->scale());
       continue;
     }
     break;
   }
   // normal case
   if(!branch) {
     if(fb.kinematics) fb.hard = false;
     return fb;
   }
   // truncated emission
   if(fb.kinematics) {
     fb.hard = false;
     fb.iout = iout;
     return fb;
   }
   // otherwise need to return the hard emission
   // construct the kinematics for the hard emission
   ShoKinPtr showerKin = new_ptr(Decay_QTildeShowerKinematics1to2(
                                 branch->scale(),
 				branch->children()[0]->z(),
 				branch->phi(),
 				branch->children()[0]->pT(),
 				branch->sudakov()));
    IdList idlist(3);
    idlist[0] = particle->dataPtr();
    idlist[1] = branch->children()[0]->branchingParticle()->dataPtr();
    idlist[2] = branch->children()[1]->branchingParticle()->dataPtr();
    // create the branching
    fb = Branching( showerKin, idlist, branch->sudakov(),ShowerPartnerType::QCDColourLine  );
    fb.hard=true;
    fb.iout=0;
    // return it
   return fb;
 }
 
 void QTildeShowerHandler::checkFlags() {
   string error = "Inconsistent hard emission set-up in QTildeShowerHandler::showerHardProcess(). "; 
   if ( ( currentTree()->isMCatNLOSEvent() || currentTree()->isMCatNLOHEvent() ) ) {
     if (_hardEmission ==2 )
       throw Exception() << error
 			<< "Cannot generate POWHEG matching with MC@NLO shower "
 			<< "approximation.  Add 'set QTildeShowerHandler:HardEmission 0' to input file."
 			<< Exception::runerror;
     if ( canHandleMatchboxTrunc() )
       throw Exception() << error
 			<< "Cannot use truncated qtilde shower with MC@NLO shower "
 			<< "approximation.  Set LHCGenerator:EventHandler"
 			<< ":CascadeHandler to '/Herwig/Shower/ShowerHandler' or "
 			<< "'/Herwig/Shower/Dipole/DipoleShowerHandler'."
 			<< Exception::runerror;
   }
   else if ( ((currentTree()->isPowhegSEvent() || currentTree()->isPowhegHEvent()) ) &&
 	    _hardEmission != 2){
     if ( canHandleMatchboxTrunc())
       throw Exception() << error
 			<< "Unmatched events requested for POWHEG shower "
 			<< "approximation.  Set QTildeShowerHandler:HardEmission to "
 			<< "'POWHEG'."
 			<< Exception::runerror;
     else if (_hardEmissionWarn) {
       _hardEmissionWarn = false;
       _hardEmission=2;
       throw Exception() << error
 			<< "Unmatched events requested for POWHEG shower "
 			<< "approximation. Changing QTildeShowerHandler:HardEmission from "
 			<< _hardEmission << " to 2"
 			<< Exception::warning;
     }
   }
 
   if ( currentTree()->isPowhegSEvent() || currentTree()->isPowhegHEvent()) {
     if (currentTree()->showerApproximation()->needsTruncatedShower() &&
 	!canHandleMatchboxTrunc() )
       throw Exception() << error
 			<< "Current shower handler cannot generate truncated shower.  "
 			<< "Set Generator:EventHandler:CascadeHandler to "
 			<< "'/Herwig/Shower/PowhegShowerHandler'."
 			<< Exception::runerror;
   }
   else if ( currentTree()->truncatedShower() && _missingTruncWarn) {
     _missingTruncWarn=false;
     throw Exception() << "Warning: POWHEG shower approximation used without "
 		      << "truncated shower.  Set Generator:EventHandler:"
 		      << "CascadeHandler to '/Herwig/Shower/PowhegShowerHandler' and "
 		      << "'MEMatching:TruncatedShower Yes'."
 		      << Exception::warning;   
   }
   // else if ( !dipme && _hardEmissionMode > 1 && 
   // 	    firstInteraction())
   //   throw Exception() << error
   // 		      << "POWHEG matching requested for LO events.  Include "
   // 		      << "'set Factory:ShowerApproximation MEMatching' in input file."
   // 		      << Exception::runerror;
 }
 
 
 tPPair QTildeShowerHandler::remakeRemnant(tPPair oldp){
   // get the parton extractor
   PartonExtractor & pex = *lastExtractor();
   // get the new partons
   tPPair newp = make_pair(findFirstParton(oldp.first ), 
 			  findFirstParton(oldp.second));
   // if the same do nothing
   if(newp == oldp) return oldp;
   // Creates the new remnants and returns the new PartonBinInstances
   // ATTENTION Broken here for very strange configuration
   PBIPair newbins = pex.newRemnants(oldp, newp, newStep());
   newStep()->addIntermediate(newp.first);
   newStep()->addIntermediate(newp.second);
   // return the new partons
   return newp;
 }
 
 PPtr QTildeShowerHandler::findFirstParton(tPPtr seed) const{
   if(seed->parents().empty()) return seed;
   tPPtr parent = seed->parents()[0];
   //if no parent there this is a loose end which will 
   //be connected to the remnant soon.
   if(!parent || parent == incomingBeams().first || 
      parent == incomingBeams().second ) return seed;
   else return findFirstParton(parent);
 }
 
 void QTildeShowerHandler::decay(ShowerTreePtr tree, ShowerDecayMap & decay) {
   // must be one incoming particle
   assert(tree->incomingLines().size()==1);
   // apply any transforms
   tree->applyTransforms();
   // if already decayed return
   if(!tree->outgoingLines().empty()) return;
   // now we need to replace the particle with a new copy after the shower
   // find particle after the shower
   map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
     tit = tree->parent()->treelinks().find(tree);
   assert(tit!=tree->parent()->treelinks().end());
   ShowerParticlePtr newparent=tit->second.second;
   PerturbativeProcessPtr newProcess =  new_ptr(PerturbativeProcess());
   newProcess->incoming().push_back(make_pair(newparent,PerturbativeProcessPtr()));
   DecayProcessMap decayMap;
   ShowerHandler::decay(newProcess,decayMap);
   ShowerTree::constructTrees(tree,decay,newProcess,decayMap);
 }
 
 namespace {
 
   ShowerProgenitorPtr
   findFinalStateLine(ShowerTreePtr tree, long id, Lorentz5Momentum momentum) { 
     map<ShowerProgenitorPtr,tShowerParticlePtr>::iterator partner;
     Energy2 dmin(1e30*GeV2);
     for(map<ShowerProgenitorPtr,tShowerParticlePtr>::iterator 
 	  cit =tree->outgoingLines().begin(); cit!=tree->outgoingLines().end(); ++cit) {
       if(cit->second->id()!=id) continue;
       Energy2 test = 
 	sqr(cit->second->momentum().x()-momentum.x())+
 	sqr(cit->second->momentum().y()-momentum.y())+
 	sqr(cit->second->momentum().z()-momentum.z())+
 	sqr(cit->second->momentum().t()-momentum.t());
       if(test<dmin) {
 	dmin    = test;
 	partner = cit;
       }
     }
     return partner->first;
   }
 
   ShowerProgenitorPtr 
   findInitialStateLine(ShowerTreePtr tree, long id, Lorentz5Momentum momentum) { 
     map<ShowerProgenitorPtr,ShowerParticlePtr>::iterator partner;
     Energy2 dmin(1e30*GeV2);
     for(map<ShowerProgenitorPtr,ShowerParticlePtr>::iterator 
 	  cit =tree->incomingLines().begin(); cit!=tree->incomingLines().end(); ++cit) {
       if(cit->second->id()!=id) continue;
       Energy2 test = 
 	sqr(cit->second->momentum().x()-momentum.x())+
 	sqr(cit->second->momentum().y()-momentum.y())+
 	sqr(cit->second->momentum().z()-momentum.z())+
 	sqr(cit->second->momentum().t()-momentum.t());
       if(test<dmin) {
 	dmin    = test;
 	partner = cit;
       }
     }
     return partner->first;
   }
 
   void fixSpectatorColours(PPtr newSpect,ShowerProgenitorPtr oldSpect,
 			   ColinePair & cline,ColinePair & aline, bool reconnect) {
     cline.first  = oldSpect->progenitor()->colourLine();
     cline.second = newSpect->colourLine();
     aline.first  = oldSpect->progenitor()->antiColourLine();
     aline.second = newSpect->antiColourLine();
     if(!reconnect) return;
     if(cline.first) {
       cline.first ->removeColoured(oldSpect->copy());
       cline.first ->removeColoured(oldSpect->progenitor());
       cline.second->removeColoured(newSpect);
       cline.first ->addColoured(newSpect);
     }
     if(aline.first) {
       aline.first ->removeAntiColoured(oldSpect->copy());
       aline.first ->removeAntiColoured(oldSpect->progenitor());
       aline.second->removeAntiColoured(newSpect);
       aline.first ->addAntiColoured(newSpect);
     }
   }
 
   void fixInitialStateEmitter(ShowerTreePtr tree, PPtr newEmit,PPtr emitted, ShowerProgenitorPtr emitter,
 			      ColinePair cline,ColinePair aline,double x) {
     // sort out the colours
     if(emitted->dataPtr()->iColour()==PDT::Colour8) {
       // emitter
       if(cline.first && cline.first == emitter->progenitor()->antiColourLine() &&
 	 cline.second !=newEmit->antiColourLine()) {
 	// sort out not radiating line
 	ColinePtr col = emitter->progenitor()->colourLine();
 	if(col) {
 	  col->removeColoured(emitter->copy());
 	  col->removeColoured(emitter->progenitor());
 	  newEmit->colourLine()->removeColoured(newEmit);
 	  col->addColoured(newEmit);
 	}
       }
       else if(aline.first && aline.first == emitter->progenitor()->colourLine() &&
 	      aline.second !=newEmit->colourLine()) {
 	// sort out not radiating line
 	ColinePtr anti = emitter->progenitor()->antiColourLine();
 	if(anti) {
 	  anti->removeAntiColoured(emitter->copy());
 	  anti->removeAntiColoured(emitter->progenitor());
 	  newEmit->colourLine()->removeAntiColoured(newEmit);
 	  anti->addAntiColoured(newEmit);
 	}
       }
       else
 	assert(false);
       // emitted
       if(cline.first && cline.second==emitted->colourLine()) {
 	cline.second->removeColoured(emitted);
 	cline.first->addColoured(emitted);
       }
       else if(aline.first && aline.second==emitted->antiColourLine()) {
 	aline.second->removeAntiColoured(emitted);
 	aline.first->addAntiColoured(emitted);
       }
       else
 	assert(false);
     }
     else {
       if(emitter->progenitor()->antiColourLine() ) {
 	ColinePtr col = emitter->progenitor()->antiColourLine();
 	col->removeAntiColoured(emitter->copy());
 	col->removeAntiColoured(emitter->progenitor());
  	if(newEmit->antiColourLine()) {
 	  newEmit->antiColourLine()->removeAntiColoured(newEmit);
 	  col->addAntiColoured(newEmit);
 	}
 	else if (emitted->colourLine()) {
 	  emitted->colourLine()->removeColoured(emitted);
 	  col->addColoured(emitted);
 	}
 	else
 	  assert(false);
       }
       if(emitter->progenitor()->colourLine() ) {
 	ColinePtr col = emitter->progenitor()->colourLine();
 	col->removeColoured(emitter->copy());
 	col->removeColoured(emitter->progenitor());
  	if(newEmit->colourLine()) {
 	  newEmit->colourLine()->removeColoured(newEmit);
 	  col->addColoured(newEmit);
 	}
 	else if (emitted->antiColourLine()) {
 	  emitted->antiColourLine()->removeAntiColoured(emitted);
 	  col->addAntiColoured(emitted);
 	}
 	else
 	  assert(false);
       }
     }
     // update the emitter
     emitter->copy(newEmit);
     ShowerParticlePtr sp = new_ptr(ShowerParticle(*newEmit,1,false));
     sp->x(x);
     emitter->progenitor(sp);
     tree->incomingLines()[emitter]=sp;
     emitter->perturbative(false);
     // add emitted
     sp=new_ptr(ShowerParticle(*emitted,1,true));
     ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(emitter->original(),emitted,sp));
     gluon->perturbative(false);
     tree->outgoingLines().insert(make_pair(gluon,sp));
   }
 
   void fixFinalStateEmitter(ShowerTreePtr tree, PPtr newEmit,PPtr emitted, ShowerProgenitorPtr emitter,
 			    ColinePair cline,ColinePair aline) {
     map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
     // special case if decayed
     for(tit  = tree->treelinks().begin(); tit != tree->treelinks().end();++tit) {
       if(tit->second.first && tit->second.second==emitter->progenitor())
 	break;
     }
     // sort out the colour lines
     if(cline.first && cline.first == emitter->progenitor()->antiColourLine() &&
        cline.second !=newEmit->antiColourLine()) {
       // sort out not radiating line
       ColinePtr col = emitter->progenitor()->colourLine();
       if(col) {
 	col->removeColoured(emitter->copy());
 	col->removeColoured(emitter->progenitor());
 	newEmit->colourLine()->removeColoured(newEmit);
 	col->addColoured(newEmit);
       }
     }
     else if(aline.first && aline.first == emitter->progenitor()->colourLine() &&
 	    aline.second !=newEmit->colourLine()) {
       // sort out not radiating line
       ColinePtr anti = emitter->progenitor()->antiColourLine();
       if(anti) {
 	anti->removeAntiColoured(emitter->copy());
 	anti->removeAntiColoured(emitter->progenitor());
 	newEmit->colourLine()->removeAntiColoured(newEmit);
 	anti->addAntiColoured(newEmit);
       }
     }
     else
       assert(false);
     // update the emitter
     emitter->copy(newEmit);
     ShowerParticlePtr sp = new_ptr(ShowerParticle(*newEmit,1,true));
     emitter->progenitor(sp);
     tree->outgoingLines()[emitter]=sp;
     emitter->perturbative(false);
     // update for decaying particles
     if(tit!=tree->treelinks().end())
       tree->updateLink(tit->first,make_pair(emitter,sp));
     // add the emitted particle
     // sort out the colour
     if(cline.first && cline.second==emitted->antiColourLine()) {
       cline.second->removeAntiColoured(emitted);
       cline.first->addAntiColoured(emitted);
     }
     else if(aline.first && aline.second==emitted->colourLine()) {
       aline.second->removeColoured(emitted);
       aline.first->addColoured(emitted);
     }
     else
       assert(false);
     sp=new_ptr(ShowerParticle(*emitted,1,true));
     ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(emitter->original(),
 						       emitted,sp));
     gluon->perturbative(false);
     tree->outgoingLines().insert(make_pair(gluon,sp));
   }
   
 }
 
 void QTildeShowerHandler::setupMECorrection(RealEmissionProcessPtr real) {
   assert(real);
   // II emission
   if(real->emitter()   < real->incoming().size() &&
      real->spectator() < real->incoming().size()) {
     // recoiling system
     for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator 
 	   cjt= currentTree()->outgoingLines().begin();
 	 cjt != currentTree()->outgoingLines().end();++cjt ) {
       cjt->first->progenitor()->transform(real->transformation());
       cjt->first->copy()->transform(real->transformation());
     }
     // the the radiating system
     ShowerProgenitorPtr emitter,spectator;
     unsigned int iemit  = real->emitter();
     unsigned int ispect = real->spectator();
     int ig = int(real->emitted())-int(real->incoming().size());
     emitter = findInitialStateLine(currentTree(),
 				   real->bornIncoming()[iemit]->id(),
 				   real->bornIncoming()[iemit]->momentum());
     spectator = findInitialStateLine(currentTree(),
 				     real->bornIncoming()[ispect]->id(),
 				     real->bornIncoming()[ispect]->momentum());
     // sort out the colours
     ColinePair cline,aline;
     fixSpectatorColours(real->incoming()[ispect],spectator,cline,aline,true);
     // update the spectator
     spectator->copy(real->incoming()[ispect]);
     ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->incoming()[ispect],1,false)));
     sp->x(ispect ==0 ? real->x().first :real->x().second);
     spectator->progenitor(sp);
     currentTree()->incomingLines()[spectator]=sp;
     spectator->perturbative(true);
     // now for the emitter
     fixInitialStateEmitter(currentTree(),real->incoming()[iemit],real->outgoing()[ig],
 			   emitter,cline,aline,iemit ==0 ? real->x().first :real->x().second);
   }
   // FF emission
   else if(real->emitter()   >= real->incoming().size() &&
 	  real->spectator() >= real->incoming().size()) {
     assert(real->outgoing()[real->emitted()-real->incoming().size()]->id()==ParticleID::g);
     // find the emitter and spectator in the shower tree
     ShowerProgenitorPtr emitter,spectator;
     int iemit = int(real->emitter())-int(real->incoming().size());
     emitter = findFinalStateLine(currentTree(),
 				 real->bornOutgoing()[iemit]->id(),
 				 real->bornOutgoing()[iemit]->momentum());
     int ispect = int(real->spectator())-int(real->incoming().size());
     spectator = findFinalStateLine(currentTree(),
 				   real->bornOutgoing()[ispect]->id(),
 				   real->bornOutgoing()[ispect]->momentum());
     map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
     // first the spectator
     // special case if decayed
     for(tit  = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) {
       if(tit->second.first && tit->second.second==spectator->progenitor())
     	break;
     }
     // sort out the colours
     ColinePair cline,aline;
     fixSpectatorColours(real->outgoing()[ispect],spectator,cline,aline,true);
     // update the spectator
     spectator->copy(real->outgoing()[ispect]);
     ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->outgoing()[ispect],1,true)));
     spectator->progenitor(sp);
     currentTree()->outgoingLines()[spectator]=sp;
     spectator->perturbative(true);
     // update for decaying particles
     if(tit!=currentTree()->treelinks().end())
       currentTree()->updateLink(tit->first,make_pair(spectator,sp));
     // now the emitting particle
     int ig = int(real->emitted())-int(real->incoming().size());
     fixFinalStateEmitter(currentTree(),real->outgoing()[iemit],
 			 real->outgoing()[ig],
 			 emitter,cline,aline);
   }
   // IF emission
   else {
     // scattering process
     if(real->incoming().size()==2) {
       ShowerProgenitorPtr emitter,spectator;
       unsigned int iemit  = real->emitter();
       unsigned int ispect = real->spectator();
       int ig = int(real->emitted())-int(real->incoming().size());
       ColinePair cline,aline;
       // incoming spectator
       if(ispect<2) {
 	spectator = findInitialStateLine(currentTree(),
 					 real->bornIncoming()[ispect]->id(),
 					 real->bornIncoming()[ispect]->momentum());
 	fixSpectatorColours(real->incoming()[ispect],spectator,cline,aline,true);
 	// update the spectator
 	spectator->copy(real->incoming()[ispect]);
 	ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->incoming()[ispect],1,false)));
 	sp->x(ispect ==0 ? real->x().first :real->x().second);
 	spectator->progenitor(sp);
 	currentTree()->incomingLines()[spectator]=sp;
 	spectator->perturbative(true);
       }
       // outgoing spectator
       else {
 	spectator = findFinalStateLine(currentTree(),
 				       real->bornOutgoing()[ispect-real->incoming().size()]->id(),
 				       real->bornOutgoing()[ispect-real->incoming().size()]->momentum());
 	// special case if decayed
 	map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
 	for(tit  = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) {
 	  if(tit->second.first && tit->second.second==spectator->progenitor())
 	    break;
 	}
 	fixSpectatorColours(real->outgoing()[ispect-real->incoming().size()],spectator,cline,aline,true);
 	// update the spectator
 	spectator->copy(real->outgoing()[ispect-real->incoming().size()]);
 	ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->outgoing()[ispect-real->incoming().size()],1,true)));
 	spectator->progenitor(sp);
 	currentTree()->outgoingLines()[spectator]=sp;
 	spectator->perturbative(true);
 	// update for decaying particles
 	if(tit!=currentTree()->treelinks().end())
 	  currentTree()->updateLink(tit->first,make_pair(spectator,sp));
       }
       // incoming emitter
       if(iemit<2) {
 	emitter = findInitialStateLine(currentTree(),
 				       real->bornIncoming()[iemit]->id(),
 				       real->bornIncoming()[iemit]->momentum());
 	fixInitialStateEmitter(currentTree(),real->incoming()[iemit],real->outgoing()[ig],
 			       emitter,aline,cline,iemit ==0 ? real->x().first :real->x().second);
       }
       // outgoing emitter
       else {
 	emitter = findFinalStateLine(currentTree(),
 				     real->bornOutgoing()[iemit-real->incoming().size()]->id(),
 				     real->bornOutgoing()[iemit-real->incoming().size()]->momentum());
 	fixFinalStateEmitter(currentTree(),real->outgoing()[iemit-real->incoming().size()],
 			     real->outgoing()[ig],emitter,aline,cline);
       }
     }
     // decay process
     else {
       assert(real->spectator()==0);
       unsigned int iemit  = real->emitter()-real->incoming().size();
       int ig = int(real->emitted())-int(real->incoming().size());
       ColinePair cline,aline;
       // incoming spectator
       ShowerProgenitorPtr spectator = findInitialStateLine(currentTree(),
 							   real->bornIncoming()[0]->id(),
 							   real->bornIncoming()[0]->momentum());
       fixSpectatorColours(real->incoming()[0],spectator,cline,aline,false);
       // find the emitter
       ShowerProgenitorPtr emitter = 
 	findFinalStateLine(currentTree(),
 			   real->bornOutgoing()[iemit]->id(),
 			   real->bornOutgoing()[iemit]->momentum());
       // recoiling system
       for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator 
 	     cjt= currentTree()->outgoingLines().begin();
 	   cjt != currentTree()->outgoingLines().end();++cjt ) {
 	if(cjt->first==emitter) continue;
 	cjt->first->progenitor()->transform(real->transformation());
 	cjt->first->copy()->transform(real->transformation());
       }
       // sort out the emitter
       fixFinalStateEmitter(currentTree(),real->outgoing()[iemit],
 			   real->outgoing()[ig],emitter,aline,cline);
     }
   }
   // clean up the shower tree
   _currenttree->resetShowerProducts();
 }
diff --git a/Shower/QTilde/SplittingFunctions/SplittingFunction.cc b/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
--- a/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
+++ b/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
@@ -1,1071 +1,1048 @@
 // -*- C++ -*-
 //
 // SplittingFunction.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 SplittingFunction class.
 //
 
 #include "SplittingFunction.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Repository/UseRandom.h"
 #include "ThePEG/Utilities/EnumIO.h"
 #include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 
 using namespace Herwig;
 
 DescribeAbstractClass<SplittingFunction,Interfaced>
 describeSplittingFunction ("Herwig::SplittingFunction","");
 
 void SplittingFunction::Init() {
 
   static ClassDocumentation<SplittingFunction> documentation
     ("The SplittingFunction class is the based class for 1->2 splitting functions"
      " in Herwig");
 
   static Switch<SplittingFunction,ColourStructure> interfaceColourStructure
     ("ColourStructure",
      "The colour structure for the splitting function",
      &SplittingFunction::_colourStructure, Undefined, false, false);
   static SwitchOption interfaceColourStructureTripletTripletOctet
     (interfaceColourStructure,
      "TripletTripletOctet",
      "3 -> 3 8",
      TripletTripletOctet);
   static SwitchOption interfaceColourStructureOctetOctetOctet
     (interfaceColourStructure,
      "OctetOctetOctet",
      "8 -> 8 8",
      OctetOctetOctet);
   static SwitchOption interfaceColourStructureOctetTripletTriplet
     (interfaceColourStructure,
      "OctetTripletTriplet",
      "8 -> 3 3bar",
      OctetTripletTriplet);
   static SwitchOption interfaceColourStructureTripletOctetTriplet
     (interfaceColourStructure,
      "TripletOctetTriplet",
      "3 -> 8 3",
      TripletOctetTriplet); 
   static SwitchOption interfaceColourStructureSextetSextetOctet
     (interfaceColourStructure,
      "SextetSextetOctet",
      "6 -> 6 8",
      SextetSextetOctet);
   
   static SwitchOption interfaceColourStructureChargedChargedNeutral
     (interfaceColourStructure,
      "ChargedChargedNeutral",
      "q -> q 0",
      ChargedChargedNeutral);
   static SwitchOption interfaceColourStructureNeutralChargedCharged
     (interfaceColourStructure,
      "NeutralChargedCharged",
      "0 -> q qbar",
      NeutralChargedCharged);
   static SwitchOption interfaceColourStructureChargedNeutralCharged
     (interfaceColourStructure,
      "ChargedNeutralCharged",
      "q -> 0 q",
      ChargedNeutralCharged);
   static SwitchOption interfaceColourStructureEW
     (interfaceColourStructure,
      "EW",
      "q -> q W/Z",
      EW);
 
   static Switch<SplittingFunction,ShowerInteraction> 
     interfaceInteractionType
     ("InteractionType",
      "Type of the interaction",
      &SplittingFunction::_interactionType, 
      ShowerInteraction::UNDEFINED, false, false);
   static SwitchOption interfaceInteractionTypeQCD
     (interfaceInteractionType,
      "QCD","QCD",ShowerInteraction::QCD);
   static SwitchOption interfaceInteractionTypeQED
     (interfaceInteractionType,
      "QED","QED",ShowerInteraction::QED);
   static SwitchOption interfaceInteractionTypeEW
     (interfaceInteractionType,
      "EW","EW",ShowerInteraction::EW);
 
   static Switch<SplittingFunction,bool> interfaceAngularOrdered
     ("AngularOrdered",
      "Whether or not this interaction is angular ordered, "
      "normally only g->q qbar and gamma-> f fbar are the only ones which aren't.",
      &SplittingFunction::angularOrdered_, true, false, false);
   static SwitchOption interfaceAngularOrderedYes
     (interfaceAngularOrdered,
      "Yes",
      "Interaction is angular ordered",
      true);
   static SwitchOption interfaceAngularOrderedNo
     (interfaceAngularOrdered,
      "No",
      "Interaction isn't angular ordered",
      false);
 
   static Switch<SplittingFunction,unsigned int> interfaceScaleChoice
     ("ScaleChoice",
      "The scale choice to be used",
      &SplittingFunction::scaleChoice_, 2, false, false);
   static SwitchOption interfaceScaleChoicepT
     (interfaceScaleChoice,
      "pT",
      "pT of the branching",
      0);
   static SwitchOption interfaceScaleChoiceQ2
     (interfaceScaleChoice,
      "Q2",
      "Q2 of the branching",
      1);
   static SwitchOption interfaceScaleChoiceFromAngularOrdering
     (interfaceScaleChoice,
      "FromAngularOrdering",
      "If angular order use pT, otherwise Q2",
      2);
 
   static Switch<SplittingFunction,bool> interfaceStrictAO
     ("StrictAO",
      "Whether or not to apply strict angular-ordering,"
      " i.e. for QED even in QCD emission, and vice versa",
      &SplittingFunction::strictAO_, true, false, false);
   static SwitchOption interfaceStrictAOYes
     (interfaceStrictAO,
      "Yes",
      "Apply strict ordering",
      true);
   static SwitchOption interfaceStrictAONo
     (interfaceStrictAO,
      "No",
      "Don't apply strict ordering",
      false);
 }
 
 void SplittingFunction::persistentOutput(PersistentOStream & os) const {
    os << oenum(_interactionType)
       << oenum(_colourStructure) << _colourFactor
       << angularOrdered_ << scaleChoice_ << strictAO_;
 }
 
 void SplittingFunction::persistentInput(PersistentIStream & is, int) {
   is >> ienum(_interactionType)
      >>	ienum(_colourStructure) >> _colourFactor
      >> angularOrdered_ >> scaleChoice_ >> strictAO_;
 }
 
 void SplittingFunction::colourConnection(tShowerParticlePtr parent,
                                          tShowerParticlePtr first,
                                          tShowerParticlePtr second,
 					 ShowerPartnerType partnerType, 
                                          const bool back) const {
   if(_colourStructure==TripletTripletOctet) {
     if(!back) {
       ColinePair cparent = ColinePair(parent->colourLine(), 
                                       parent->antiColourLine());
       // ensure input consistency
       assert((  cparent.first && !cparent.second && 
 		partnerType==ShowerPartnerType::QCDColourLine) || 
              ( !cparent.first &&  cparent.second && 
 		partnerType==ShowerPartnerType::QCDAntiColourLine));
       // q -> q g
       if(cparent.first) {
         ColinePtr newline=new_ptr(ColourLine());
         cparent.first->addColoured(second);
         newline->addColoured     ( first);
         newline->addAntiColoured (second);
       }
       // qbar -> qbar g
       else {
         ColinePtr newline=new_ptr(ColourLine());
         cparent.second->addAntiColoured(second);
         newline->addColoured(second);
         newline->addAntiColoured(first);
       }
       // Set progenitor
       first->progenitor(parent->progenitor());
       second->progenitor(parent->progenitor());
     }
     else {
       ColinePair cfirst = ColinePair(first->colourLine(), 
                                      first->antiColourLine());
       // ensure input consistency
       assert((  cfirst.first && !cfirst.second && 
 		partnerType==ShowerPartnerType::QCDColourLine) || 
              ( !cfirst.first &&  cfirst.second && 
 		partnerType==ShowerPartnerType::QCDAntiColourLine));
       // q -> q g
       if(cfirst.first) {
         ColinePtr newline=new_ptr(ColourLine());
         cfirst.first->addAntiColoured(second);
         newline->addColoured(second);
         newline->addColoured(parent);
       }
       // qbar -> qbar g
       else {
         ColinePtr newline=new_ptr(ColourLine());
         cfirst.second->addColoured(second);
         newline->addAntiColoured(second);
         newline->addAntiColoured(parent);
       }
       // Set progenitor
       parent->progenitor(first->progenitor());
       second->progenitor(first->progenitor());
     }
   }
   else if(_colourStructure==OctetOctetOctet) {
     if(!back) {
       ColinePair cparent = ColinePair(parent->colourLine(), 
                                       parent->antiColourLine());
       // ensure input consistency
       assert(cparent.first&&cparent.second);
       // ensure first gluon is hardest
       if( first->id()==second->id() && parent->showerKinematics()->z()<0.5 ) 
 	swap(first,second);
       // colour line radiates
       if(partnerType==ShowerPartnerType::QCDColourLine) {
 	// The colour line is radiating
 	ColinePtr newline=new_ptr(ColourLine());
 	cparent.first->addColoured(second);
 	cparent.second->addAntiColoured(first);
 	newline->addColoured(first);
 	newline->addAntiColoured(second);
       }
       // anti colour line radiates
       else if(partnerType==ShowerPartnerType::QCDAntiColourLine) {
 	ColinePtr newline=new_ptr(ColourLine());
 	cparent.first->addColoured(first);
 	cparent.second->addAntiColoured(second);
 	newline->addColoured(second);
 	newline->addAntiColoured(first);
       }
       else
 	assert(false);
     }
     else {
       ColinePair cfirst = ColinePair(first->colourLine(), 
                                      first->antiColourLine());
       // ensure input consistency
       assert(cfirst.first&&cfirst.second);
       // The colour line is radiating
       if(partnerType==ShowerPartnerType::QCDColourLine) {
 	ColinePtr newline=new_ptr(ColourLine());
 	cfirst.first->addAntiColoured(second);
 	cfirst.second->addAntiColoured(parent);
 	newline->addColoured(parent);
 	newline->addColoured(second);
       }
       // anti colour line radiates
       else if(partnerType==ShowerPartnerType::QCDAntiColourLine) {
 	ColinePtr newline=new_ptr(ColourLine());
 	cfirst.first->addColoured(parent);
 	cfirst.second->addColoured(second);
 	newline->addAntiColoured(second);
 	newline->addAntiColoured(parent);
       }
       else
 	assert(false);
     }    
   }
   else if(_colourStructure == OctetTripletTriplet) {
     if(!back) {
       ColinePair cparent = ColinePair(parent->colourLine(), 
                                       parent->antiColourLine());
       // ensure input consistency
       assert(cparent.first&&cparent.second);
       cparent.first ->addColoured    ( first);
       cparent.second->addAntiColoured(second);
       // Set progenitor
       first->progenitor(parent->progenitor());
       second->progenitor(parent->progenitor());
     }
     else {
       ColinePair cfirst = ColinePair(first->colourLine(), 
                                      first->antiColourLine());
       // ensure input consistency
       assert(( cfirst.first && !cfirst.second) ||
              (!cfirst.first &&  cfirst.second));
       // g -> q qbar
       if(cfirst.first) {
 	ColinePtr newline=new_ptr(ColourLine());
 	cfirst.first->addColoured(parent);
 	newline->addAntiColoured(second);
 	newline->addAntiColoured(parent);	
       }
       // g -> qbar q
       else {
         ColinePtr newline=new_ptr(ColourLine());
         cfirst.second->addAntiColoured(parent);
         newline->addColoured(second);
         newline->addColoured(parent);
       }
       // Set progenitor
       parent->progenitor(first->progenitor());
       second->progenitor(first->progenitor());
     }
   }
   else if(_colourStructure == TripletOctetTriplet) {
     if(!back) {
       ColinePair cparent = ColinePair(parent->colourLine(), 
                                       parent->antiColourLine());
       // ensure input consistency
       assert(( cparent.first && !cparent.second) || 
              (!cparent.first &&  cparent.second));
       // q -> g q
       if(cparent.first) {
         ColinePtr newline=new_ptr(ColourLine());
         cparent.first->addColoured(first);
         newline->addColoured    (second);
         newline->addAntiColoured( first);
       }
       // qbar -> g qbar
       else {
 	ColinePtr newline=new_ptr(ColourLine());
 	cparent.second->addAntiColoured(first);
 	newline->addColoured    ( first);
 	newline->addAntiColoured(second);	
       }
       // Set progenitor
       first->progenitor(parent->progenitor());
       second->progenitor(parent->progenitor());
     }
     else {
       ColinePair cfirst = ColinePair(first->colourLine(), 
                                      first->antiColourLine());
       // ensure input consistency
       assert(cfirst.first&&cfirst.second);
       // q -> g q
       if(parent->id()>0) {
         cfirst.first ->addColoured(parent);
         cfirst.second->addColoured(second);
       }
       else {
         cfirst.first ->addAntiColoured(second);
         cfirst.second->addAntiColoured(parent);
       }
       // Set progenitor
       parent->progenitor(first->progenitor());
       second->progenitor(first->progenitor());
     }
   }
   else if(_colourStructure==SextetSextetOctet) {
     //make sure we're not doing backward evolution
     assert(!back);
 
     //make sure something sensible
     assert(parent->colourLine() || parent->antiColourLine());
    
     //get the colour lines or anti-colour lines
     bool isAntiColour=true;
     ColinePair cparent;
     if(parent->colourLine()) {
       cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->colourLines()[0]), 
 			   const_ptr_cast<tColinePtr>(parent->colourInfo()->colourLines()[1]));
       isAntiColour=false;
     }
     else {
       cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->antiColourLines()[0]), 
 			   const_ptr_cast<tColinePtr>(parent->colourInfo()->antiColourLines()[1]));
     }
     
     //check for sensible input
     //    assert(cparent.first && cparent.second);
 
     // sextet has 2 colour lines
     if(!isAntiColour) {
       //pick at random which of the colour topolgies to take
       double topology = UseRandom::rnd();
       if(topology < 0.25) {
         ColinePtr newline=new_ptr(ColourLine());
         cparent.first->addColoured(second);
         cparent.second->addColoured(first);
         newline->addColoured(first);
         newline->addAntiColoured(second);
       }
       else if(topology >=0.25 && topology < 0.5) {
         ColinePtr newline=new_ptr(ColourLine());
         cparent.first->addColoured(first);
         cparent.second->addColoured(second);
         newline->addColoured(first);
         newline->addAntiColoured(second); 
       }
       else if(topology >= 0.5 && topology < 0.75) {
         ColinePtr newline=new_ptr(ColourLine());
         cparent.first->addColoured(second);
         cparent.second->addColoured(first); 
         newline->addColoured(first); 
         newline->addAntiColoured(second); 
       }
       else {
         ColinePtr newline=new_ptr(ColourLine());
         cparent.first->addColoured(first);
         cparent.second->addColoured(second);
         newline->addColoured(first);
         newline->addAntiColoured(second);
       }
     }
     // sextet has 2 anti-colour lines
     else {
       double topology = UseRandom::rnd();
       if(topology < 0.25){
         ColinePtr newline=new_ptr(ColourLine());
         cparent.first->addAntiColoured(second);
         cparent.second->addAntiColoured(first);
         newline->addAntiColoured(first);
         newline->addColoured(second);
       }
       else if(topology >=0.25 && topology < 0.5) {
         ColinePtr newline=new_ptr(ColourLine());
         cparent.first->addAntiColoured(first);
         cparent.second->addAntiColoured(second);
         newline->addAntiColoured(first);
         newline->addColoured(second); 
       }
       else if(topology >= 0.5 && topology < 0.75) {
         ColinePtr newline=new_ptr(ColourLine());
         cparent.first->addAntiColoured(second);
         cparent.second->addAntiColoured(first);
         newline->addAntiColoured(first);
         newline->addColoured(second); 
       }
       else {
         ColinePtr newline=new_ptr(ColourLine());
         cparent.first->addAntiColoured(first);
         cparent.second->addAntiColoured(second);
         newline->addAntiColoured(first);
         newline->addColoured(second);
       }
     }   
   }
   else if(_colourStructure == ChargedChargedNeutral) {
     if(!parent->data().coloured()) return;
     if(!back) {
       ColinePair cparent = ColinePair(parent->colourLine(), 
 				      parent->antiColourLine());
       // q -> q g
       if(cparent.first) {
 	cparent.first->addColoured(first);
       }
       // qbar -> qbar g
       if(cparent.second) {
 	cparent.second->addAntiColoured(first);
       }
     }
     else {
       ColinePair cfirst = ColinePair(first->colourLine(), 
 				     first->antiColourLine());
       // q -> q g
       if(cfirst.first) {
 	cfirst.first->addColoured(parent);
       }
       // qbar -> qbar g
       if(cfirst.second) {
 	cfirst.second->addAntiColoured(parent);
       }
     }
   }
   else if(_colourStructure == ChargedNeutralCharged) {
     if(!parent->data().coloured()) return;
     if(!back) {
       ColinePair cparent = ColinePair(parent->colourLine(), 
 				      parent->antiColourLine());
       // q -> q g
       if(cparent.first) {
 	cparent.first->addColoured(second);
       }
       // qbar -> qbar g
       if(cparent.second) {
 	cparent.second->addAntiColoured(second);
       }
     }
     else {
       if (second->dataPtr()->iColour()==PDT::Colour3 ) {
 	ColinePtr newline=new_ptr(ColourLine());
 	newline->addColoured(second);
 	newline->addColoured(parent);
       }
       else if (second->dataPtr()->iColour()==PDT::Colour3bar ) {
 	ColinePtr newline=new_ptr(ColourLine());
 	newline->addAntiColoured(second);
 	newline->addAntiColoured(parent);
       }
     }
   }
   else if(_colourStructure == NeutralChargedCharged ) {
     if(!back) {
       if(first->dataPtr()->coloured()) {
 	ColinePtr newline=new_ptr(ColourLine());
 	if(first->dataPtr()->iColour()==PDT::Colour3) {
 	  newline->addColoured    (first );
 	  newline->addAntiColoured(second);
 	}
 	else if (first->dataPtr()->iColour()==PDT::Colour3bar) {
 	  newline->addColoured    (second);
 	  newline->addAntiColoured(first );
 	}
 	else
 	  assert(false);
       }
     }
     else {   
       ColinePair cfirst = ColinePair(first->colourLine(), 
 				     first->antiColourLine());
       // gamma -> q qbar
       if(cfirst.first) {
 	cfirst.first->addAntiColoured(second);
       }
       // gamma -> qbar q
       else if(cfirst.second) {
 	cfirst.second->addColoured(second);
       }
       else 
 	assert(false);
     }
   }
   else if(_colourStructure == EW) {
     if(!parent->data().coloured()) return;
     if(!back) {
       ColinePair cparent = ColinePair(parent->colourLine(), 
 				      parent->antiColourLine());
       // q -> q g
       if(cparent.first) {
 	cparent.first->addColoured(first);
       }
       // qbar -> qbar g
       if(cparent.second) {
 	cparent.second->addAntiColoured(first);
       }
     }
     else {
       ColinePair cfirst = ColinePair(first->colourLine(), 
 				     first->antiColourLine());
       // q -> q g
       if(cfirst.first) {
 	cfirst.first->addColoured(parent);
       }
       // qbar -> qbar g
       if(cfirst.second) {
 	cfirst.second->addAntiColoured(parent);
       }
     }
   }
   else {
     assert(false);
   }
 }
 
 void SplittingFunction::doinit() {
   Interfaced::doinit();
   assert(_interactionType!=ShowerInteraction::UNDEFINED);
   assert((_colourStructure>0&&_interactionType==ShowerInteraction::QCD) ||
 	 (_colourStructure<0&&(_interactionType==ShowerInteraction::QED ||
 			       _interactionType==ShowerInteraction::EW)) );
   if(_colourFactor>0.) return;
   // compute the colour factors if need
   if(_colourStructure==TripletTripletOctet) {
     _colourFactor = 4./3.;
   }
   else if(_colourStructure==OctetOctetOctet) {
     _colourFactor = 3.;
   }
   else if(_colourStructure==OctetTripletTriplet) {
     _colourFactor = 0.5;
   }
   else if(_colourStructure==TripletOctetTriplet) {
     _colourFactor = 4./3.;
   }
   else if(_colourStructure==SextetSextetOctet) {
     _colourFactor = 10./3.;
   }
   else if(_colourStructure<0) {
     _colourFactor = 1.;
   }
   else {
     assert(false);
   }
 }
 
 bool SplittingFunction::checkColours(const IdList & ids) const {
   if(_colourStructure==TripletTripletOctet) {
     if(ids[0]!=ids[1]) return false;
     if((ids[0]->iColour()==PDT::Colour3||ids[0]->iColour()==PDT::Colour3bar) &&
        ids[2]->iColour()==PDT::Colour8) return true;
     return false;
   }
   else if(_colourStructure==OctetOctetOctet) {
     for(unsigned int ix=0;ix<3;++ix) {
       if(ids[ix]->iColour()!=PDT::Colour8) return false;
     }
     return true;
   }
   else if(_colourStructure==OctetTripletTriplet) {
     if(ids[0]->iColour()!=PDT::Colour8) return false;
     if(ids[1]->iColour()==PDT::Colour3&&ids[2]->iColour()==PDT::Colour3bar)
       return true;
     if(ids[1]->iColour()==PDT::Colour3bar&&ids[2]->iColour()==PDT::Colour3)
       return true;
     return false;
   }
   else if(_colourStructure==TripletOctetTriplet) {
     if(ids[0]!=ids[2]) return false;
     if((ids[0]->iColour()==PDT::Colour3||ids[0]->iColour()==PDT::Colour3bar) &&
        ids[1]->iColour()==PDT::Colour8) return true;
     return false;
   }
   else if(_colourStructure==SextetSextetOctet) {
     if(ids[0]!=ids[1]) return false;
     if((ids[0]->iColour()==PDT::Colour6 || ids[0]->iColour()==PDT::Colour6bar) &&
        ids[2]->iColour()==PDT::Colour8) return true;
     return false;
   }
   else if(_colourStructure==ChargedChargedNeutral) {
     if(ids[0]!=ids[1]) return false;
     if(ids[2]->iCharge()!=0) return false;
     if(ids[0]->iCharge()==ids[1]->iCharge()) return true;
     return false;
   }
   else if(_colourStructure==ChargedNeutralCharged) {
     if(ids[0]!=ids[2]) return false;
     if(ids[1]->iCharge()!=0) return false;
     if(ids[0]->iCharge()==ids[2]->iCharge()) return true;
     return false;
   }
   else if(_colourStructure==NeutralChargedCharged) {
     if(ids[1]->id()!=-ids[2]->id()) return false;
     if(ids[0]->iCharge()!=0) return false;
     if(ids[1]->iCharge()==-ids[2]->iCharge()) return true;
     return false;
   }
   else {
     assert(false);
   }
   return false;
 }
 
 namespace {
 
   bool hasColour(tPPtr p) {
     PDT::Colour colour = p->dataPtr()->iColour();
     return colour==PDT::Colour3 || colour==PDT::Colour8 || colour == PDT::Colour6;
   }
 
   bool hasAntiColour(tPPtr p) {
     PDT::Colour colour = p->dataPtr()->iColour();
     return colour==PDT::Colour3bar || colour==PDT::Colour8 || colour == PDT::Colour6bar;
   }
   
 }
 
 void SplittingFunction::evaluateFinalStateScales(ShowerPartnerType partnerType,
 						 Energy scale, double z,
 						 tShowerParticlePtr parent,
 						 tShowerParticlePtr emitter,
 						 tShowerParticlePtr emitted) {
   // identify emitter and emitted
   double zEmitter = z, zEmitted = 1.-z;
   bool bosonSplitting(false);
   // special for g -> gg, particle highest z is emitter
   if(emitter->id() == emitted->id() && emitter->id() == parent->id() &&
      zEmitted > zEmitter) {
     swap(zEmitted,zEmitter);
     swap( emitted, emitter);
   }
   // otherwise if particle ID same
   else if(emitted->id()==parent->id()) {
     swap(zEmitted,zEmitter);
     swap( emitted, emitter);
   }
   // no real emitter/emitted
   else if(emitter->id()!=parent->id()) {
     bosonSplitting = true;
   }
   // may need to add angularOrder flag here
   // now the various scales
   // QED
   if(partnerType==ShowerPartnerType::QED) {
     assert(colourStructure()==ChargedChargedNeutral ||
 	   colourStructure()==ChargedNeutralCharged ||
 	   colourStructure()==NeutralChargedCharged );
     // normal case
     if(!bosonSplitting) {
       assert(colourStructure()==ChargedChargedNeutral ||
 	     colourStructure()==ChargedNeutralCharged );
       // set the scales
       // emitter
       emitter->scales().QED         = zEmitter*scale;
       emitter->scales().QED_noAO    =          scale;
       if(strictAO_)
 	emitter->scales().QCD_c       = min(zEmitter*scale,parent->scales().QCD_c      );
       else
 	emitter->scales().QCD_c       = min(         scale,parent->scales().QCD_c      );
       emitter->scales().QCD_c_noAO  = min(scale,parent->scales().QCD_c_noAO );
       if(strictAO_)
 	emitter->scales().QCD_ac      = min(zEmitter*scale,parent->scales().QCD_ac     );
       else
 	emitter->scales().QCD_ac      = min(         scale,parent->scales().QCD_ac     );
       emitter->scales().QCD_ac_noAO = min(scale,parent->scales().QCD_ac_noAO);
-      emitter->scales().EW_Z        = min(scale,parent->scales().EW_Z       );
-      emitter->scales().EW_W        = min(scale,parent->scales().EW_W       );
+      emitter->scales().EW          = min(scale,parent->scales().EW         );
       // emitted 
       emitted->scales().QED         = zEmitted*scale;
       emitted->scales().QED_noAO    =          scale;
       emitted->scales().QCD_c       = ZERO;
       emitted->scales().QCD_c_noAO  = ZERO;
       emitted->scales().QCD_ac      = ZERO;
       emitted->scales().QCD_ac_noAO = ZERO;
-      emitted->scales().EW_Z        = min(scale,parent->scales().EW_Z       );
-      emitted->scales().EW_W        = min(scale,parent->scales().EW_W       );
+      emitted->scales().EW          = min(scale,parent->scales().EW         );
     }
     // gamma -> f fbar
     else {
       assert(colourStructure()==NeutralChargedCharged );
       // emitter
       emitter->scales().QED           = zEmitter*scale;
       emitter->scales().QED_noAO      =          scale;
       if(hasColour(emitter)) {
 	emitter->scales().QCD_c       = zEmitter*scale;
 	emitter->scales().QCD_c_noAO  =          scale;
       }
       if(hasAntiColour(emitter)) {
 	emitter->scales().QCD_ac      = zEmitter*scale;
 	emitter->scales().QCD_ac_noAO =          scale;
       }
-      emitter->scales().EW_Z          = zEmitter*scale;
-      emitter->scales().EW_W          = zEmitter*scale;
+      emitter->scales().EW            = zEmitter*scale;
       // emitted 
       emitted->scales().QED           = zEmitted*scale;
       emitted->scales().QED_noAO      =          scale;
       if(hasColour(emitted)) {
 	emitted->scales().QCD_c       = zEmitted*scale;
 	emitted->scales().QCD_c_noAO  =          scale;
       }
       if(hasAntiColour(emitted)) {
 	emitted->scales().QCD_ac      = zEmitted*scale;
 	emitted->scales().QCD_ac_noAO =          scale;
       }
-      emitted->scales().EW_Z          = zEmitted*scale;
-      emitted->scales().EW_W          = zEmitted*scale;
+      emitted->scales().EW            = zEmitted*scale;
     }
   }
   // QCD
   else if (partnerType==ShowerPartnerType::QCDColourLine ||
 	   partnerType==ShowerPartnerType::QCDAntiColourLine) {
    // normal case eg q -> q g and g -> g g
     if(!bosonSplitting) {
       if(strictAO_)
 	emitter->scales().QED         = min(zEmitter*scale,parent->scales().QED     );
       else
 	emitter->scales().QED         = min(         scale,parent->scales().QED     );
-      emitter->scales().QED_noAO    = min(scale,parent->scales().QED_noAO);
-      emitter->scales().EW_Z        = min(scale,parent->scales().EW_Z   );
-      emitter->scales().EW_W        = min(scale,parent->scales().EW_W   );
+      emitter->scales().QED_noAO      = min(scale,parent->scales().QED_noAO);
+      emitter->scales().EW            = min(scale,parent->scales().EW     );
       if(partnerType==ShowerPartnerType::QCDColourLine) {
 	emitter->scales().QCD_c       = zEmitter*scale;
 	emitter->scales().QCD_c_noAO  =          scale;
 	emitter->scales().QCD_ac      = min(zEmitter*scale,parent->scales().QCD_ac     );
 	emitter->scales().QCD_ac_noAO = min(         scale,parent->scales().QCD_ac_noAO);
       }
       else {
 	emitter->scales().QCD_c       = min(zEmitter*scale,parent->scales().QCD_c      );
 	emitter->scales().QCD_c_noAO  = min(         scale,parent->scales().QCD_c_noAO );
 	emitter->scales().QCD_ac      = zEmitter*scale;
 	emitter->scales().QCD_ac_noAO =          scale;
       }
       // emitted 
       emitted->scales().QED         = ZERO;
       emitted->scales().QED_noAO    = ZERO;
       emitted->scales().QCD_c       = zEmitted*scale;
       emitted->scales().QCD_c_noAO  =          scale;
       emitted->scales().QCD_ac      = zEmitted*scale;
       emitted->scales().QCD_ac_noAO =          scale;
-      emitted->scales().EW_Z        = min(scale,parent->scales().EW_Z   );
-      emitted->scales().EW_W        = min(scale,parent->scales().EW_W   );
+      emitted->scales().EW          = min(scale,parent->scales().EW     );
     }
     // g -> q qbar
     else {
       // emitter
       if(emitter->dataPtr()->charged()) {
 	emitter->scales().QED         = zEmitter*scale;
 	emitter->scales().QED_noAO    =          scale;
       }
-      emitter->scales().EW_Z          = zEmitter*scale;
-      emitter->scales().EW_W          = zEmitter*scale;
+      emitter->scales().EW            = zEmitter*scale;
       emitter->scales().QCD_c         = zEmitter*scale;
       emitter->scales().QCD_c_noAO    =          scale;
       emitter->scales().QCD_ac        = zEmitter*scale;
       emitter->scales().QCD_ac_noAO   =          scale;
       // emitted 
       if(emitted->dataPtr()->charged()) {
 	emitted->scales().QED         = zEmitted*scale;
 	emitted->scales().QED_noAO    =          scale;
       }
-      emitted->scales().EW_Z          = zEmitted*scale;
-      emitted->scales().EW_W          = zEmitted*scale;
+      emitted->scales().EW            = zEmitted*scale;
       emitted->scales().QCD_c         = zEmitted*scale;
       emitted->scales().QCD_c_noAO    =          scale;
       emitted->scales().QCD_ac        = zEmitted*scale;
       emitted->scales().QCD_ac_noAO   =          scale;
     }
   }
   else if(partnerType==ShowerPartnerType::EW) {
     // EW
-    emitter->scales().EW_Z        = zEmitter*scale;
-    emitter->scales().EW_W        = zEmitter*scale;
-    emitted->scales().EW_Z        = zEmitted*scale;
-    emitted->scales().EW_W        = zEmitted*scale;
+    emitter->scales().EW          = zEmitter*scale;
+    emitted->scales().EW          = zEmitted*scale;
     // QED
     // W radiation AO
     if(emitted->dataPtr()->charged()) {
       emitter->scales().QED       = zEmitter*scale;
       emitter->scales().QED_noAO  =          scale;
       emitted->scales().QED       = zEmitted*scale;
       emitted->scales().QED_noAO  =          scale;
     }
     // Z don't
     else {
       emitter->scales().QED         = min(scale,parent->scales().QED     );
       emitter->scales().QED_noAO    = min(scale,parent->scales().QED_noAO);
       emitted->scales().QED         = ZERO;
       emitted->scales().QED_noAO    = ZERO;
     }
     // QCD
     emitter->scales().QCD_c       = min(scale,parent->scales().QCD_c      );
     emitter->scales().QCD_c_noAO  = min(scale,parent->scales().QCD_c_noAO );
     emitter->scales().QCD_ac      = min(scale,parent->scales().QCD_ac     );
     emitter->scales().QCD_ac_noAO = min(scale,parent->scales().QCD_ac_noAO);
     emitted->scales().QCD_c       = ZERO;
     emitted->scales().QCD_c_noAO  = ZERO;
     emitted->scales().QCD_ac      = ZERO;
     emitted->scales().QCD_ac_noAO = ZERO;
   }
   else
     assert(false);
 }
 
 void SplittingFunction::evaluateInitialStateScales(ShowerPartnerType partnerType,
 						   Energy scale, double z,
 						   tShowerParticlePtr parent,
 						   tShowerParticlePtr spacelike,
 						   tShowerParticlePtr timelike) {
   // scale for time-like child
   Energy AOScale = (1.-z)*scale;
   // QED
   if(partnerType==ShowerPartnerType::QED) {
     if(parent->id()==spacelike->id()) {
       // parent
       parent   ->scales().QED         =   scale;
       parent   ->scales().QED_noAO    =   scale;
       parent   ->scales().QCD_c       = min(scale,spacelike->scales().QCD_c      );
       parent   ->scales().QCD_c_noAO  = min(scale,spacelike->scales().QCD_c_noAO );
       parent   ->scales().QCD_ac      = min(scale,spacelike->scales().QCD_ac     );
       parent   ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
       // timelike
       timelike->scales().QED         = AOScale;
       timelike->scales().QED_noAO    =   scale;
       timelike->scales().QCD_c       =    ZERO;
       timelike->scales().QCD_c_noAO  =    ZERO;
       timelike->scales().QCD_ac      =    ZERO;
       timelike->scales().QCD_ac_noAO =    ZERO;
     }
     else if(parent->id()==timelike->id()) {
       parent   ->scales().QED         =   scale;
       parent   ->scales().QED_noAO    =   scale;
       if(hasColour(parent)) {
 	parent   ->scales().QCD_c       = scale;
 	parent   ->scales().QCD_c_noAO  = scale;
       }
       if(hasAntiColour(parent)) {
 	parent   ->scales().QCD_ac      = scale;
 	parent   ->scales().QCD_ac_noAO = scale;
       }
       // timelike 
       timelike->scales().QED         = AOScale;
       timelike->scales().QED_noAO    =   scale;
       if(hasColour(timelike)) {
 	timelike->scales().QCD_c       = AOScale;
 	timelike->scales().QCD_c_noAO  =   scale;
       }
       if(hasAntiColour(timelike)) {
 	timelike->scales().QCD_ac      = AOScale;
 	timelike->scales().QCD_ac_noAO =   scale;
       }
     }
     else {
       parent   ->scales().QED         = scale;
       parent   ->scales().QED_noAO    = scale;
       parent   ->scales().QCD_c       = ZERO ;
       parent   ->scales().QCD_c_noAO  = ZERO ;
       parent   ->scales().QCD_ac      = ZERO ;
       parent   ->scales().QCD_ac_noAO = ZERO ;
       // timelike 
       timelike->scales().QED         = AOScale;
       timelike->scales().QED_noAO    =   scale;
       if(hasColour(timelike)) {
 	timelike->scales().QCD_c       = min(AOScale,spacelike->scales().QCD_ac     );
 	timelike->scales().QCD_c_noAO  = min(  scale,spacelike->scales().QCD_ac_noAO);
       }
       if(hasAntiColour(timelike)) {
 	timelike->scales().QCD_ac      = min(AOScale,spacelike->scales().QCD_c      );
 	timelike->scales().QCD_ac_noAO = min(  scale,spacelike->scales().QCD_c_noAO );
       }
     }
   }
   // QCD
   else if (partnerType==ShowerPartnerType::QCDColourLine ||
 	   partnerType==ShowerPartnerType::QCDAntiColourLine) {
     // timelike 
     if(timelike->dataPtr()->charged()) {
       timelike->scales().QED         = AOScale;
       timelike->scales().QED_noAO    =   scale;
     }
     if(hasColour(timelike)) {
       timelike->scales().QCD_c       = AOScale;
       timelike->scales().QCD_c_noAO  =   scale;
     }
     if(hasAntiColour(timelike)) {
       timelike->scales().QCD_ac      = AOScale;
       timelike->scales().QCD_ac_noAO =   scale;
     }
     if(parent->id()==spacelike->id()) {
       parent   ->scales().QED         = min(scale,spacelike->scales().QED        );
       parent   ->scales().QED_noAO    = min(scale,spacelike->scales().QED_noAO   );
       parent   ->scales().QCD_c       = min(scale,spacelike->scales().QCD_c      );
       parent   ->scales().QCD_c_noAO  = min(scale,spacelike->scales().QCD_c_noAO );
       parent   ->scales().QCD_ac      = min(scale,spacelike->scales().QCD_ac     );
       parent   ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
     }
     else {
       if(parent->dataPtr()->charged()) {
 	parent   ->scales().QED         = scale;
 	parent   ->scales().QED_noAO    = scale;
       }
       if(hasColour(parent)) {
 	parent   ->scales().QCD_c      = scale;
 	parent   ->scales().QCD_c_noAO  = scale;
       }
       if(hasAntiColour(parent)) {
 	parent   ->scales().QCD_ac      = scale;
 	parent   ->scales().QCD_ac_noAO = scale;
       }
     }
   }
   else if(partnerType==ShowerPartnerType::EW) {
     if(abs(spacelike->id())!=ParticleID::Wplus && 
        spacelike->id() !=ParticleID::Z0 ) {
       // QCD scales
       parent   ->scales().QCD_c       = min(scale,spacelike->scales().QCD_c      );
       parent   ->scales().QCD_c_noAO  = min(scale,spacelike->scales().QCD_c_noAO );
       parent   ->scales().QCD_ac      = min(scale,spacelike->scales().QCD_ac     );
       parent   ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
       timelike->scales().QCD_c       =    ZERO;
       timelike->scales().QCD_c_noAO  =    ZERO;
       timelike->scales().QCD_ac      =    ZERO;
       timelike->scales().QCD_ac_noAO =    ZERO;
       // QED scales
       if(timelike->id()==ParticleID::Z0) {
 	parent   ->scales().QED       = min(scale,spacelike->scales().QED      );
 	parent   ->scales().QED_noAO  = min(scale,spacelike->scales().QED_noAO );
 	timelike->scales().QED       =    ZERO;
 	timelike->scales().QED_noAO  =    ZERO;
       }
       else {
 	parent   ->scales().QED         =   scale;
 	parent   ->scales().QED_noAO    =   scale;
 	timelike->scales().QED          = AOScale;
 	timelike->scales().QED_noAO     =   scale;
       }
       // EW scales
-      if(timelike->id()==ParticleID::Z0) {
-	parent   ->scales().EW_Z       =   scale;
-	timelike->scales().EW_Z        = AOScale;
-      }
-      if(timelike->id()==ParticleID::Wplus){
-	parent   ->scales().EW_W       =   scale;
-	timelike->scales().EW_W        = AOScale;
-      }
+	parent   ->scales().EW         =   scale;
+	timelike->scales().EW          = AOScale;
     }
     else assert(false);
   }
   else
     assert(false);
 }
 
 void SplittingFunction::evaluateDecayScales(ShowerPartnerType partnerType,
 					    Energy scale, double z,
 					    tShowerParticlePtr parent,
 					    tShowerParticlePtr spacelike,
 					    tShowerParticlePtr timelike) {
   assert(parent->id()==spacelike->id());
   // angular-ordered scale for 2nd child
   Energy AOScale = (1.-z)*scale;
   // QED
   if(partnerType==ShowerPartnerType::QED) {
     // timelike
     timelike->scales().QED         = AOScale;
     timelike->scales().QED_noAO    =   scale;
     timelike->scales().QCD_c       =    ZERO;
     timelike->scales().QCD_c_noAO  =    ZERO;
     timelike->scales().QCD_ac      =    ZERO;
     timelike->scales().QCD_ac_noAO =    ZERO;
-    timelike->scales().EW_Z        = ZERO;
-    timelike->scales().EW_W        = ZERO;
+    timelike->scales().EW          = ZERO;
     // spacelike
     spacelike->scales().QED         =   scale;
     spacelike->scales().QED_noAO    =   scale;
-    spacelike->scales().EW_Z        = max(scale,parent->scales().EW_Z       );
-    spacelike->scales().EW_W        = max(scale,parent->scales().EW_W       );
+    spacelike->scales().EW          = max(scale,parent->scales().EW         );
   }
   // QCD
   else if(partnerType==ShowerPartnerType::QCDColourLine ||
 	  partnerType==ShowerPartnerType::QCDAntiColourLine) {
     // timelike 
     timelike->scales().QED         = ZERO;
     timelike->scales().QED_noAO    = ZERO;
     timelike->scales().QCD_c       = AOScale;
     timelike->scales().QCD_c_noAO  =   scale;
     timelike->scales().QCD_ac      = AOScale;
     timelike->scales().QCD_ac_noAO =   scale;
-    timelike->scales().EW_Z        = ZERO;
-    timelike->scales().EW_W        = ZERO;
+    timelike->scales().EW          = ZERO;
     // spacelike
     spacelike->scales().QED         = max(scale,parent->scales().QED        );
     spacelike->scales().QED_noAO    = max(scale,parent->scales().QED_noAO   );
-    spacelike->scales().EW_Z        = max(scale,parent->scales().EW_Z       );
-    spacelike->scales().EW_W        = max(scale,parent->scales().EW_W       );
+    spacelike->scales().EW          = max(scale,parent->scales().EW         );
   }
   else if(partnerType==ShowerPartnerType::EW) {
     // EW
-    timelike->scales().EW_Z         = AOScale;
-    timelike->scales().EW_W         = AOScale;
-    spacelike->scales().EW_Z        = max(scale,parent->scales().EW_Z       );
-    spacelike->scales().EW_W        = max(scale,parent->scales().EW_W       );
+    timelike->scales().EW           = AOScale;
+    spacelike->scales().EW          = max(scale,parent->scales().EW         );
     // QCD
     timelike->scales().QCD_c       =    ZERO;
     timelike->scales().QCD_c_noAO  =    ZERO;
     timelike->scales().QCD_ac      =    ZERO;
     timelike->scales().QCD_ac_noAO =    ZERO;
-    timelike->scales().EW_Z        = ZERO;
-    timelike->scales().EW_W        = ZERO;
+    timelike->scales().EW          = ZERO;
     // QED
     timelike->scales().QED         = ZERO;
     timelike->scales().QED_noAO    = ZERO;
     spacelike->scales().QED         = max(scale,parent->scales().QED        );
     spacelike->scales().QED_noAO    = max(scale,parent->scales().QED_noAO   );
   }
   else 
     assert(false);
   spacelike->scales().QCD_c       = max(scale,parent->scales().QCD_c      );
   spacelike->scales().QCD_c_noAO  = max(scale,parent->scales().QCD_c_noAO );
   spacelike->scales().QCD_ac      = max(scale,parent->scales().QCD_ac     );
   spacelike->scales().QCD_ac_noAO = max(scale,parent->scales().QCD_ac_noAO);
 }
diff --git a/Shower/QTilde/SplittingFunctions/SplittingGenerator.cc b/Shower/QTilde/SplittingFunctions/SplittingGenerator.cc
--- a/Shower/QTilde/SplittingFunctions/SplittingGenerator.cc
+++ b/Shower/QTilde/SplittingFunctions/SplittingGenerator.cc
@@ -1,656 +1,648 @@
 // -*- C++ -*-
 //
 // SplittingGenerator.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 SplittingGenerator class.
 //
 
 #include "SplittingGenerator.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/Command.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Utilities/StringUtils.h"
 #include "ThePEG/Repository/Repository.h"
 #include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
 #include "Herwig/Shower/ShowerHandler.h"
 #include "ThePEG/Utilities/Rebinder.h"
 #include <cassert>
 #include "ThePEG/Utilities/DescribeClass.h"
 
 using namespace Herwig;
 
 namespace {
 
 bool checkInteraction(ShowerInteraction allowed,
 		      ShowerInteraction splitting) {
   if(allowed==ShowerInteraction::ALL)
     return true;
   else if(allowed==ShowerInteraction::QEDQCD &&
 	  (splitting==ShowerInteraction::QED ||
 	   splitting==ShowerInteraction::QCD ))
     return true;
   else if(allowed == splitting)
     return true;
   else
     return false;
 }
 
 }
 
 DescribeClass<SplittingGenerator,Interfaced>
 describeSplittingGenerator ("Herwig::SplittingGenerator","");
 
 IBPtr SplittingGenerator::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr SplittingGenerator::fullclone() const {
   return new_ptr(*this);
 }
 
 void SplittingGenerator::persistentOutput(PersistentOStream & os) const {
   os << _bbranchings << _fbranchings << _deTuning;
 }
 
 void SplittingGenerator::persistentInput(PersistentIStream & is, int) {
   is >> _bbranchings >> _fbranchings >> _deTuning;
 }
 
 void SplittingGenerator::Init() {
 
   static ClassDocumentation<SplittingGenerator> documentation
     ("There class is responsible for initializing the Sudakov form factors ",
      "and generating splittings.");
 
   static Command<SplittingGenerator> interfaceAddSplitting
     ("AddFinalSplitting",
      "Adds another splitting to the list of splittings considered "
      "in the shower. Command is a->b,c; Sudakov",
      &SplittingGenerator::addFinalSplitting);
 
   static Command<SplittingGenerator> interfaceAddInitialSplitting
     ("AddInitialSplitting",
      "Adds another splitting to the list of initial splittings to consider "
      "in the shower. Command is a->b,c; Sudakov. Here the particle a is the "
      "particle that is PRODUCED by the splitting. b is the initial state "
      "particle that is splitting in the shower.",
      &SplittingGenerator::addInitialSplitting);
 
   static Command<SplittingGenerator> interfaceDeleteSplitting
     ("DeleteFinalSplitting",
      "Deletes a splitting from the list of splittings considered "
      "in the shower. Command is a->b,c; Sudakov",
      &SplittingGenerator::deleteFinalSplitting);
 
   static Command<SplittingGenerator> interfaceDeleteInitialSplitting
     ("DeleteInitialSplitting",
      "Deletes a splitting from the list of initial splittings to consider "
      "in the shower. Command is a->b,c; Sudakov. Here the particle a is the "
      "particle that is PRODUCED by the splitting. b is the initial state "
      "particle that is splitting in the shower.",
      &SplittingGenerator::deleteInitialSplitting);
 
   static Parameter<SplittingGenerator,double> interfaceDetuning
     ("Detuning",
      "The Detuning parameter to make the veto algorithm less efficient to improve the weight variations",
      &SplittingGenerator::_deTuning, 1.0, 1.0, 10.0,
      false, false, Interface::limited);
 
 }
 
 string SplittingGenerator::addSplitting(string arg, bool final) {
   string partons = StringUtils::car(arg);
   string sudakov = StringUtils::cdr(arg);
   vector<tPDPtr> products;
   string::size_type next = partons.find("->");
   if(next == string::npos) 
     return "Error: Invalid string for splitting " + arg;
   if(partons.find(';') == string::npos) 
     return "Error: Invalid string for splitting " + arg;
   tPDPtr parent = Repository::findParticle(partons.substr(0,next));
   partons = partons.substr(next+2);
   do {
     next = min(partons.find(','), partons.find(';'));
     tPDPtr pdp = Repository::findParticle(partons.substr(0,next));
     partons = partons.substr(next+1);
     if(pdp) products.push_back(pdp);
     else return "Error: Could not create splitting from " + arg;
   } while(partons[0] != ';' && partons.size());
   SudakovPtr s;
   s = dynamic_ptr_cast<SudakovPtr>(Repository::TraceObject(sudakov));
   if(!s) return "Error: Could not load Sudakov " + sudakov + '\n';
   IdList ids;
   ids.push_back(parent);
   for(vector<tPDPtr>::iterator it = products.begin(); it!=products.end(); ++it)
     ids.push_back(*it);
   // check splitting can handle this
   if(!s->splittingFn()->accept(ids)) 
     return "Error: Sudakov " + sudakov + " SplittingFunction can't handle particles\n";
   // add to map
   addToMap(ids,s,final);
   return "";
 }
 
 string SplittingGenerator::deleteSplitting(string arg, bool final) {
   string partons = StringUtils::car(arg);
   string sudakov = StringUtils::cdr(arg);
   vector<tPDPtr> products;
   string::size_type next = partons.find("->");
   if(next == string::npos) 
     return "Error: Invalid string for splitting " + arg;
   if(partons.find(';') == string::npos) 
     return "Error: Invalid string for splitting " + arg;
   tPDPtr parent = Repository::findParticle(partons.substr(0,next));
   partons = partons.substr(next+2);
   do {
     next = min(partons.find(','), partons.find(';'));
     tPDPtr pdp = Repository::findParticle(partons.substr(0,next));
     partons = partons.substr(next+1);
     if(pdp) products.push_back(pdp);
     else return "Error: Could not create splitting from " + arg;
   } while(partons[0] != ';' && partons.size());
   SudakovPtr s;
   s = dynamic_ptr_cast<SudakovPtr>(Repository::TraceObject(sudakov));
   if(!s) return "Error: Could not load Sudakov " + sudakov + '\n';
   IdList ids;
   ids.push_back(parent);
   for(vector<tPDPtr>::iterator it = products.begin(); it!=products.end(); ++it)
     ids.push_back(*it);
   // check splitting can handle this
   if(!s->splittingFn()->accept(ids)) 
     return "Error: Sudakov " + sudakov + " SplittingFunction can't handle particles\n";
   // delete from map
   deleteFromMap(ids,s,final);
   return "";
 }
 
 void SplittingGenerator::addToMap(const IdList &ids, const SudakovPtr &s, bool final) {
   if(!final) {
       // search if the branching was already included.
     auto binsert =BranchingInsert(abs(ids[1]->id()),BranchingElement(s,ids));
       // get the range of already inserted splittings.
     auto eqrange=_bbranchings.equal_range(binsert.first);
     for(auto it = eqrange.first; it != eqrange.second; ++it){
       if((*it).second == binsert.second)
         throw Exception()<<"SplittingGenerator: Trying to insert existing splitting.\n"
         << Exception::setuperror;
     }
     _bbranchings.insert(binsert);
     s->addSplitting(ids);
   }
   else {
       // search if the branching was already included.
     auto binsert =BranchingInsert(abs(ids[0]->id()),BranchingElement(s,ids));
       // get the range of already inserted splittings.
     auto eqrange=_fbranchings.equal_range(binsert.first);
     for(auto it = eqrange.first; it != eqrange.second; ++it){
       if((*it).second ==binsert.second)
         throw Exception()<<"SplittingGenerator: Trying to insert existing splitting.\n"
         << Exception::setuperror;
     }
     
     _fbranchings.insert(binsert);
     s->addSplitting(ids);
   }
 }
 
 void SplittingGenerator::deleteFromMap(const IdList &ids, 
 				       const SudakovPtr &s, bool final) {
   bool didRemove=false;
   if(!final) {
     pair<BranchingList::iterator,BranchingList::iterator> 
       range = _bbranchings.equal_range(abs(ids[1]->id()));
     for(BranchingList::iterator it=range.first;
 	it!=range.second&&it!=_bbranchings.end();++it) {
       if(it->second.sudakov==s&&it->second.particles==ids) {
 	BranchingList::iterator it2=it;
 	--it;
 	_bbranchings.erase(it2);
     didRemove=true;
       }
     }
     s->removeSplitting(ids);
   }
   else {
     pair<BranchingList::iterator,BranchingList::iterator> 
       range = _fbranchings.equal_range(abs(ids[0]->id()));
     for(BranchingList::iterator it=range.first;
 	it!=range.second&&it!=_fbranchings.end();++it) {
       if(it->second.sudakov==s&&it->second.particles==ids) {
 	BranchingList::iterator it2 = it;
 	--it;
 	_fbranchings.erase(it2);
     didRemove=true;
       }
     }
     s->removeSplitting(ids);
   }
   if (!didRemove)
     throw Exception()<<"SplittingGenerator: Try to remove non existing splitting.\n"
                       << Exception::setuperror;
 }
 
 Branching SplittingGenerator::chooseForwardBranching(ShowerParticle &particle,
 						     double enhance,
 						     ShowerInteraction type) const {
   RhoDMatrix rho;
   bool rhoCalc(false);
   Energy newQ = ZERO;
   ShoKinPtr kinematics = ShoKinPtr();
   ShowerPartnerType partnerType(ShowerPartnerType::Undefined);
   SudakovPtr sudakov   = SudakovPtr();
   IdList ids;
   // First, find the eventual branching, corresponding to the highest scale.
   long index = abs(particle.data().id());
   // if no branchings return empty branching struct
   if( _fbranchings.find(index) == _fbranchings.end() ) 
     return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined);
   // otherwise select branching
   for(BranchingList::const_iterator cit = _fbranchings.lower_bound(index); 
       cit != _fbranchings.upper_bound(index); ++cit) {
     // check either right interaction or doing both
     if(!checkInteraction(type,cit->second.sudakov->interactionType())) continue;
     if(!rhoCalc) {
       rho = particle.extractRhoMatrix(true);
       rhoCalc = true;
     }
     // whether or not this interaction should be angular ordered
     bool angularOrdered = cit->second.sudakov->splittingFn()->angularOrdered();
     ShoKinPtr newKin;
     ShowerPartnerType type;
     IdList particles = particle.id()!=cit->first ? cit->second.conjugateParticles : cit->second.particles;
     // work out which starting scale we need
     if(cit->second.sudakov->interactionType()==ShowerInteraction::QED) {
       type = ShowerPartnerType::QED;
       Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO;
       newKin = cit->second.sudakov->
     	generateNextTimeBranching(startingScale,particles,rho,enhance,
 				  _deTuning);
     }
     else if(cit->second.sudakov->interactionType()==ShowerInteraction::QCD) {
       // special for octets
       if(particle.dataPtr()->iColour()==PDT::Colour8) {
 	// octet -> octet octet
 	if(cit->second.sudakov->splittingFn()->colourStructure()==OctetOctetOctet) {
     	  type = ShowerPartnerType::QCDColourLine;
 	  Energy startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
     	  newKin= cit->second.sudakov->
     	    generateNextTimeBranching(startingScale,particles,rho,0.5*enhance,
 				      _deTuning);
 	  startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
     	  ShoKinPtr newKin2 = cit->second.sudakov->
 	    generateNextTimeBranching(startingScale,particles,rho,0.5*enhance,
 				      _deTuning);
     	  // pick the one with the highest scale
     	  if( ( newKin && newKin2 && newKin2->scale() > newKin->scale()) ||
     	      (!newKin && newKin2) ) {
     	    newKin = newKin2;
     	    type = ShowerPartnerType::QCDAntiColourLine;
     	  }
     	}
      	// other g -> q qbar
      	else {
      	  Energy startingScale = angularOrdered ? 
 	    max(particle.scales().QCD_c     , particle.scales().QCD_ac     ) : 
 	    max(particle.scales().QCD_c_noAO, particle.scales().QCD_ac_noAO);
     	  newKin= cit->second.sudakov->
     	    generateNextTimeBranching(startingScale,particles,rho,enhance,
 				      _deTuning);
     	  type = UseRandom::rndbool() ? 
     	    ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
 	}
       }
       // everything else q-> qg etc
       else {
 	Energy startingScale;
 	if(particle.hasColour()) {
 	  type = ShowerPartnerType::QCDColourLine;
 	  startingScale = angularOrdered ? particle.scales().QCD_c  : particle.scales().QCD_c_noAO;
 	}
 	else {
 	  type = ShowerPartnerType::QCDAntiColourLine;
 	  startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
 	}
 	newKin= cit->second.sudakov->
 	  generateNextTimeBranching(startingScale,particles,rho,enhance,
 				    _deTuning);
       }
     }
     else if(cit->second.sudakov->interactionType()==ShowerInteraction::EW) {
       type = ShowerPartnerType::EW;
-      Energy startingScale = (abs(particles[1]->id())==ParticleID::Wplus|| 
-			      abs(particles[2]->id())==ParticleID::Wplus) ? particle.scales().EW_W :  particle.scales().EW_Z;
+      Energy startingScale = particle.scales().EW;
       newKin = cit->second.sudakov->
     	generateNextTimeBranching(startingScale,particles,rho,enhance,_deTuning);
     }
     // shouldn't be anything else
     else
       assert(false);
     // if no kinematics contine
     if(!newKin) continue;
     // select highest scale 
     if( newKin->scale() > newQ ) {
       kinematics  = newKin;
       newQ        = newKin->scale();
       ids         = particles;
       sudakov     = cit->second.sudakov;
       partnerType = type;
     }
   }
   // return empty branching if nothing happened
   if(!kinematics) {
     if ( particle.spinInfo() ) particle.spinInfo()->undecay();
     return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
 		     ShowerPartnerType::Undefined);
   }
   // if not hard generate phi
   kinematics->phi(sudakov->generatePhiForward(particle,ids,kinematics,rho));
   // and return it
   return Branching(kinematics, ids,sudakov,partnerType);
 }
 
 Branching SplittingGenerator::
 chooseDecayBranching(ShowerParticle &particle,
 		     const ShowerParticle::EvolutionScales & stoppingScales,
 		     Energy minmass, double enhance,
 		     ShowerInteraction interaction) const {
   RhoDMatrix rho(particle.dataPtr()->iSpin());
   Energy newQ = Constants::MaxEnergy;
   ShoKinPtr kinematics;
   SudakovPtr sudakov;
   ShowerPartnerType partnerType(ShowerPartnerType::Undefined);
   IdList ids;
   // First, find the eventual branching, corresponding to the lowest scale.
   long index = abs(particle.data().id());
   // if no branchings return empty branching struct
   if(_fbranchings.find(index) == _fbranchings.end()) 
     return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined);
   // otherwise select branching
   for(BranchingList::const_iterator cit = _fbranchings.lower_bound(index); 
       cit != _fbranchings.upper_bound(index); ++cit)  {
     // check interaction doesn't change flavour
     if(cit->second.particles[1]->id()!=index&&cit->second.particles[2]->id()!=index) continue;
     // check either right interaction or doing both
     if(!checkInteraction(interaction,cit->second.sudakov->interactionType())) continue;
     // whether or not this interaction should be angular ordered
     bool angularOrdered = cit->second.sudakov->splittingFn()->angularOrdered();
     ShoKinPtr newKin;
     IdList particles = particle.id()!=cit->first ? cit->second.conjugateParticles : cit->second.particles;
     ShowerPartnerType type;
     // work out which starting scale we need
     if(cit->second.sudakov->interactionType()==ShowerInteraction::QED) {
       type = ShowerPartnerType::QED;
       Energy stoppingScale = angularOrdered ? stoppingScales.QED    : stoppingScales.QED_noAO;
       Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO;
       if(startingScale < stoppingScale ) { 
     	newKin = cit->second.sudakov->
     	  generateNextDecayBranching(startingScale,stoppingScale,minmass,particles,rho,
 				     enhance,_deTuning);
       }
     }
     else if(cit->second.sudakov->interactionType()==ShowerInteraction::QCD) {
       // special for octets
       if(particle.dataPtr()->iColour()==PDT::Colour8) {
 	// octet -> octet octet
 	if(cit->second.sudakov->splittingFn()->colourStructure()==OctetOctetOctet) {
 	  Energy stoppingColour = angularOrdered ? stoppingScales.QCD_c     : stoppingScales.QCD_c_noAO;
 	  Energy stoppingAnti   = angularOrdered ? stoppingScales.QCD_ac    : stoppingScales.QCD_ac_noAO;
 	  Energy startingColour = angularOrdered ? particle.scales().QCD_c  : particle.scales().QCD_c_noAO;
 	  Energy startingAnti   = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
 	  type = ShowerPartnerType::QCDColourLine;
 	  if(startingColour<stoppingColour) {
 	    newKin= cit->second.sudakov->	
 	      generateNextDecayBranching(startingColour,stoppingColour,minmass,
 					 particles,rho,0.5*enhance,_deTuning);
 	  }
 	  ShoKinPtr newKin2; 
 	  if(startingAnti<stoppingAnti) {
 	    newKin2 = cit->second.sudakov->
 	      generateNextDecayBranching(startingAnti,stoppingAnti,minmass,
 					 particles,rho,0.5*enhance,_deTuning);
 	  }
 	  // pick the one with the lowest scale
 	  if( (newKin&&newKin2&&newKin2->scale()<newKin->scale()) ||
 	      (!newKin&&newKin2) ) {
 	    newKin = newKin2;
 	    type = ShowerPartnerType::QCDAntiColourLine;
 	  }
 	}
 	// other
 	else {
 	  assert(false);
 	}
       }
       // everything else
       else {
 	Energy startingScale,stoppingScale;
 	if(particle.hasColour()) {
 	  type = ShowerPartnerType::QCDColourLine;
 	  stoppingScale = angularOrdered ? stoppingScales.QCD_c     : stoppingScales.QCD_c_noAO;
 	  startingScale = angularOrdered ? particle.scales().QCD_c  : particle.scales().QCD_c_noAO;
 	}
 	else {
 	  type = ShowerPartnerType::QCDAntiColourLine;
 	  stoppingScale = angularOrdered ? stoppingScales.QCD_ac    : stoppingScales.QCD_ac_noAO;
 	  startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
 	}
 	if(startingScale < stoppingScale ) {
 	  newKin = cit->second.sudakov->
 	    generateNextDecayBranching(startingScale,stoppingScale,minmass,particles,rho,
 				       enhance,_deTuning);
 	}
       }
     }
     else if(cit->second.sudakov->interactionType()==ShowerInteraction::EW) {
       type = ShowerPartnerType::EW;
       Energy stoppingScale, startingScale;
-      if(abs(particles[1]->id())==ParticleID::Wplus|| abs(particles[2]->id())==ParticleID::Wplus) {
-         stoppingScale = stoppingScales.EW_W;
-	 startingScale = particle.scales().EW_W;
-      }
-      else {
-	stoppingScale = stoppingScales.EW_Z;
-	startingScale = particle.scales().EW_Z;
-      }
+	stoppingScale = stoppingScales.EW;
+	startingScale = particle.scales().EW;
       if(startingScale < stoppingScale ) { 
     	newKin = cit->second.sudakov->
     	  generateNextDecayBranching(startingScale,stoppingScale,minmass,particles,rho,enhance,_deTuning);
       }
     }
     // shouldn't be anything else
     else
       assert(false);
     if(!newKin) continue;
     // select highest scale
     if(newKin->scale() < newQ ) {
       newQ = newKin->scale();
       ids = particles;
       kinematics=newKin;
       sudakov=cit->second.sudakov;
       partnerType = type;
     }
   }
   // return empty branching if nothing happened
   if(!kinematics)  return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
 				    ShowerPartnerType::Undefined);
   // and generate phi
   kinematics->phi(sudakov->generatePhiDecay(particle,ids,kinematics,rho));
   // and return it
   return Branching(kinematics, ids,sudakov,partnerType);
 }
 
 Branching SplittingGenerator::
 chooseBackwardBranching(ShowerParticle &particle,PPtr,
 			double enhance,
 			Ptr<BeamParticleData>::transient_const_pointer beam,
 			ShowerInteraction type,
 			tcPDFPtr pdf, Energy freeze) const {
   RhoDMatrix rho;
   bool rhoCalc(false);
   Energy newQ=ZERO;
   ShoKinPtr kinematics=ShoKinPtr();
   ShowerPartnerType partnerType(ShowerPartnerType::Undefined);
   SudakovPtr sudakov;
   IdList ids;
   // First, find the eventual branching, corresponding to the highest scale.
   long index = abs(particle.id());
   // if no possible branching return
   if(_bbranchings.find(index) == _bbranchings.end())
     return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined);
   // otherwise select branching
   for(BranchingList::const_iterator cit = _bbranchings.lower_bound(index); 
       cit != _bbranchings.upper_bound(index); ++cit ) {
     // check either right interaction or doing both
     if(!checkInteraction(type,cit->second.sudakov->interactionType())) continue;
     // setup the PDF
     cit->second.sudakov->setPDF(pdf,freeze);
     //calc rho as needed
     if(!rhoCalc) {
       rho = particle.extractRhoMatrix(false);
       rhoCalc = true;
     }
     // whether or not this interaction should be angular ordered
     bool angularOrdered = cit->second.sudakov->splittingFn()->angularOrdered();
     ShoKinPtr newKin;
     IdList particles = particle.id()!=cit->first ? cit->second.conjugateParticles : cit->second.particles;
     ShowerPartnerType type;
     if(cit->second.sudakov->interactionType()==ShowerInteraction::QED) {
       type = ShowerPartnerType::QED;
       Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO;
       newKin=cit->second.sudakov->
     	generateNextSpaceBranching(startingScale,particles, particle.x(),rho,enhance,
 				   beam,_deTuning);
     }
     else if(cit->second.sudakov->interactionType()==ShowerInteraction::QCD) { 
       // special for octets
       if(particle.dataPtr()->iColour()==PDT::Colour8) {
 	// octet -> octet octet
 	if(cit->second.sudakov->splittingFn()->colourStructure()==OctetOctetOctet) {
     	  type = ShowerPartnerType::QCDColourLine;
 	  Energy startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
 	  newKin = cit->second.sudakov->
 	    generateNextSpaceBranching(startingScale,particles, particle.x(),rho,0.5*enhance,
 				       beam,_deTuning);
 	  startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
 	  ShoKinPtr newKin2 = cit->second.sudakov->
 	    generateNextSpaceBranching(startingScale,particles, particle.x(),rho,
 				       0.5*enhance,beam,_deTuning);
 	  // pick the one with the highest scale
 	  if( (newKin&&newKin2&&newKin2->scale()>newKin->scale()) ||
 	      (!newKin&&newKin2) ) {
 	    newKin = newKin2;
 	    type = ShowerPartnerType::QCDAntiColourLine;
 	  }
 	}
 	else {
      	  Energy startingScale = angularOrdered ? 
 	    max(particle.scales().QCD_c     , particle.scales().QCD_ac    ) : 
 	    max(particle.scales().QCD_c_noAO, particle.scales().QCD_ac_noAO);
 	  type = UseRandom::rndbool() ? 
 	    ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
 	  newKin=cit->second.sudakov->
 	    generateNextSpaceBranching(startingScale,particles, particle.x(),rho,enhance,beam,_deTuning);
 	}
       }
       // everything else
       else {
 	Energy startingScale;
 	if(particle.hasColour()) {
 	  type = ShowerPartnerType::QCDColourLine;
 	  startingScale = angularOrdered ? particle.scales().QCD_c  : particle.scales().QCD_c_noAO;
 	}
 	else {
 	  type = ShowerPartnerType::QCDAntiColourLine;
 	  startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
 	}
     	newKin=cit->second.sudakov->
     	  generateNextSpaceBranching(startingScale,particles,particle.x(),rho,enhance,beam,_deTuning);
       }
     }
     else if(cit->second.sudakov->interactionType()==ShowerInteraction::EW) {
       type = ShowerPartnerType::EW;
-      Energy startingScale = (abs(particles[1]->id())==ParticleID::Wplus|| 
-			      abs(particles[2]->id())==ParticleID::Wplus) ? particle.scales().EW_W :  particle.scales().EW_Z;
+      Energy startingScale = particle.scales().EW;
       newKin=cit->second.sudakov->
     	generateNextSpaceBranching(startingScale,particles,particle.x(),rho,enhance,beam,_deTuning);
     }
     // shouldn't be anything else
     else
       assert(false);
     // if no kinematics contine
     if(!newKin) continue;
     // select highest scale
     if(newKin->scale() > newQ) {
       newQ = newKin->scale();
       kinematics=newKin;
       ids = particles;
       sudakov=cit->second.sudakov;
       partnerType = type;
     }
   }
   // return empty branching if nothing happened
   if(!kinematics) {
     if ( particle.spinInfo() ) particle.spinInfo()->undecay();
     return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
 		     ShowerPartnerType::Undefined);
   }
   // initialize the ShowerKinematics 
   // and generate phi
   kinematics->phi(sudakov->generatePhiBackward(particle,ids,kinematics,rho));
   // return the answer
   return Branching(kinematics, ids,sudakov,partnerType);
 }
 
 void SplittingGenerator::rebind(const TranslationMap & trans) {
   BranchingList::iterator cit;
   for(cit=_fbranchings.begin();cit!=_fbranchings.end();++cit) {
     (cit->second).sudakov=trans.translate((cit->second).sudakov);
     for(unsigned int ix=0;ix<(cit->second).particles.size();++ix) {
       (cit->second).particles[ix]=trans.translate((cit->second).particles[ix]);
     }
     for(unsigned int ix=0;ix<(cit->second).conjugateParticles.size();++ix) {
       (cit->second).conjugateParticles[ix]=trans.translate((cit->second).conjugateParticles[ix]);
     }
   }
   for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit) {
     (cit->second).sudakov=trans.translate((cit->second).sudakov);
     for(unsigned int ix=0;ix<(cit->second).particles.size();++ix) {
       (cit->second).particles[ix]=trans.translate((cit->second).particles[ix]);
     }
     for(unsigned int ix=0;ix<(cit->second).conjugateParticles.size();++ix) {
       (cit->second).conjugateParticles[ix]=trans.translate((cit->second).conjugateParticles[ix]);
     }
   }
   Interfaced::rebind(trans);
 }
 
 IVector SplittingGenerator::getReferences() {
   IVector ret = Interfaced::getReferences();
   BranchingList::iterator cit;
   for(cit=_fbranchings.begin();cit!=_fbranchings.end();++cit) {
     ret.push_back((cit->second).sudakov);
     for(unsigned int ix=0;ix<(cit->second).particles.size();++ix) 
       ret.push_back(const_ptr_cast<tPDPtr>((cit->second).particles[ix]));
     for(unsigned int ix=0;ix<(cit->second).conjugateParticles.size();++ix) 
       ret.push_back(const_ptr_cast<tPDPtr>((cit->second).conjugateParticles[ix]));
   }
   for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit) {
     ret.push_back((cit->second).sudakov);
     for(unsigned int ix=0;ix<(cit->second).particles.size();++ix) 
       ret.push_back(const_ptr_cast<tPDPtr>((cit->second).particles[ix]));
     for(unsigned int ix=0;ix<(cit->second).conjugateParticles.size();++ix) 
       ret.push_back(const_ptr_cast<tPDPtr>((cit->second).conjugateParticles[ix]));
   }
   return ret;
 }