diff --git a/MatrixElement/Matchbox/Matching/QTildeMatching.cc b/MatrixElement/Matchbox/Matching/QTildeMatching.cc
--- a/MatrixElement/Matchbox/Matching/QTildeMatching.cc
+++ b/MatrixElement/Matchbox/Matching/QTildeMatching.cc
@@ -1,525 +1,525 @@
 // -*- C++ -*-
 //
 // QTildeMatching.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 QTildeMatching class.
 //
 
 #include "QTildeMatching.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/EventRecord/Particle.h"
 #include "ThePEG/Repository/UseRandom.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 
 
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
 #include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
 
 using namespace Herwig;
 
 QTildeMatching::QTildeMatching() 
   : theCorrectForXZMismatch(true) {}
 
 QTildeMatching::~QTildeMatching() {}
 
 IBPtr QTildeMatching::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr QTildeMatching::fullclone() const {
   return new_ptr(*this);
 }
 
 void QTildeMatching::checkCutoff() {
   if ( showerTildeKinematics() ) {
     showerTildeKinematics()->
       prepare(realCXComb(),bornCXComb());
     showerTildeKinematics()->dipole(dipole());
     showerTildeKinematics()->getShowerVariables();
   }
 }
 
 void QTildeMatching::getShowerVariables() {
 
   // already filled from checkCutoff in this case
   if ( showerTildeKinematics() )
     return;
 
   // get the shower variables
   calculateShowerVariables();
 
   // check for the cutoff
   dipole()->isAboveCutoff(isAboveCutoff());
 
   // get the hard scale
   dipole()->showerHardScale(hardScale());
 
   // check for phase space
   dipole()->isInShowerPhasespace(isInShowerPhasespace());
 
 }
 
 bool QTildeMatching::isInShowerPhasespace() const {
 
   assert((theQTildeSudakov->cutOffOption() == 0 || theQTildeSudakov->cutOffOption() == 2) && 
 	 "implementation only provided for default and pt cutoff");
 
   Energy qtildeHard = ZERO;
 
   Energy qtilde = dipole()->showerScale();
   assert(!dipole()->showerParameters().empty());
   double z = dipole()->showerParameters()[0];
 
   // FF
   if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) {
     qtildeHard = 
       theQTildeFinder->
       calculateFinalFinalScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
 				bornCXComb()->meMomenta()[dipole()->bornSpectator()],
 				bornCXComb()->mePartonData()[dipole()->bornEmitter()]->iColour() == PDT::Colour3).first;
   }
 
   // FI
   if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) {
     qtildeHard = 
       theQTildeFinder->
       calculateInitialFinalScales(bornCXComb()->meMomenta()[dipole()->bornSpectator()],
 				  bornCXComb()->meMomenta()[dipole()->bornEmitter()],false).second;
   }
 
   // IF
   if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) {
     qtildeHard = 
       theQTildeFinder->
       calculateInitialFinalScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
 				  bornCXComb()->meMomenta()[dipole()->bornSpectator()],false).first;
     if ( z < (dipole()->bornEmitter() == 0 ? bornCXComb()->lastX1() : bornCXComb()->lastX2()) )
       return false;
   }
 
   // II
   if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2 ) {
     qtildeHard = 
       theQTildeFinder->
       calculateInitialInitialScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
 				    bornCXComb()->meMomenta()[dipole()->bornSpectator()]).first;
     if ( z < (dipole()->bornEmitter() == 0 ? bornCXComb()->lastX1() : bornCXComb()->lastX2()) )
       return false;
   }
 
 
 
   Energy Qg = theQTildeSudakov->kinScale();
   Energy2 pt2 = ZERO;
   if ( dipole()->bornEmitter() > 1 ) {
     Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
     if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
       pt2 = sqr(z*(1.-z)*qtilde) - sqr(mu);
     else
       pt2 = sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg);
   }
   if ( dipole()->bornEmitter() < 2 ) {
     pt2 = sqr((1.-z)*qtilde) - z*sqr(Qg);
   }
 
   if ( pt2 < max(theQTildeSudakov->pT2min(),sqr(safeCut()) ))
     return false;
 
   bool hardVeto = restrictPhasespace() && sqrt(pt2) >= dipole()->showerHardScale();
   return qtilde <= qtildeHard && !hardVeto;
 
 }
 
 bool QTildeMatching::isAboveCutoff() const {
 
   assert((theQTildeSudakov->cutOffOption() == 0 || theQTildeSudakov->cutOffOption() == 2) && 
 	 "implementation only provided for default and pt cutoff");
   Energy qtilde = dipole()->showerScale();
   assert(!dipole()->showerParameters().empty());
   double z = dipole()->showerParameters()[0];
   Energy Qg = theQTildeSudakov->kinScale();
   if ( dipole()->bornEmitter() > 1 ) {
     Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
     if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
       return sqr(z*(1.-z)*qtilde) - sqr(mu) >= 
              max(theQTildeSudakov->pT2min(),sqr(safeCut()));
     else
       return sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg) >= 
              max(theQTildeSudakov->pT2min(),sqr(safeCut()));
   }
   if ( dipole()->bornEmitter() < 2 ) {
     return
       sqr((1.-z)*qtilde) - z*sqr(Qg) >= 
       max(theQTildeSudakov->pT2min(),sqr(safeCut()));
   }
   return false;
 }
 
 CrossSection QTildeMatching::dSigHatDR() const {
 
   assert(!dipole()->showerParameters().empty());
   pair<Energy2,double> vars = 
     make_pair(sqr(dipole()->showerScale()),
 	      dipole()->showerParameters()[0]);
 
   pair<int,int> ij(dipole()->bornEmitter(),
 		   dipole()->bornSpectator());
   double ccme2 =
      dipole()->underlyingBornME()->largeNColourCorrelatedME2(ij,theLargeNBasis);
    
    if(ccme2==0.)return 0.*nanobarn;
       
    double lnme2=dipole()->underlyingBornME()->largeNME2(theLargeNBasis);
    if(lnme2==0){
      generator()->log() <<"\nQTildeMatching: ";
      generator()->log() <<"\n  largeNME2 is ZERO, while largeNColourCorrelatedME2 is not ZERO." ;
      generator()->log() <<"\n  This is too seriuos.\n" ;
      generator()->log() << Exception::runerror;
    }
   
    ccme2 *=
      dipole()->underlyingBornME()->me2() /lnme2;
      
 
   Energy2 prop = ZERO;
   if ( dipole()->bornEmitter() > 1 ) {
     prop =
       (realCXComb()->meMomenta()[dipole()->realEmitter()] +
        realCXComb()->meMomenta()[dipole()->realEmission()]).m2()
       - bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2();
   } else {
     prop = 
       2.*vars.second*(realCXComb()->meMomenta()[dipole()->realEmitter()]*
 		      realCXComb()->meMomenta()[dipole()->realEmission()]);
   }
 
   // note alphas included downstream from subtractionScaleWeight()
   double xme2 = -8.*Constants::pi*ccme2*splitFn(vars)*realXComb()->lastSHat()/prop;
   xme2 *= 
     pow(realCXComb()->lastSHat() / bornCXComb()->lastSHat(),
 	bornCXComb()->mePartonData().size()-4.);
 
   double bornPDF = bornPDFWeight(dipole()->underlyingBornME()->lastScale());
   if ( bornPDF == 0.0 )
     return ZERO;
 
   xme2 *= bornPDF;
   
   xme2 *= dipole()->realEmissionME()->finalStateSymmetry() /
     dipole()->underlyingBornME()->finalStateSymmetry(); 
 
   // take care of mismatch between z and x as we are approaching the
   // hard phase space boundary
   // TODO get rid of this useless scale option business and simplify PDF handling in here
   if ( dipole()->bornEmitter() < 2 && theCorrectForXZMismatch ) {
     Energy2 emissionScale = ZERO;
     if ( emissionScaleInSubtraction() == showerScale ) {
       emissionScale = showerFactorizationScale();
     } else if ( emissionScaleInSubtraction() == realScale ) {
       emissionScale = dipole()->realEmissionME()->lastScale();
     } else if ( emissionScaleInSubtraction() == bornScale ) {
       emissionScale = dipole()->underlyingBornME()->lastScale();
     }
     double xzMismatch = 
       dipole()->subtractionParameters()[0] / dipole()->showerParameters()[0];
     double realCorrectedPDF = 
       dipole()->bornEmitter() == 0 ?
       dipole()->realEmissionME()->pdf1(emissionScale,theExtrapolationX,
 				       xzMismatch) :
       dipole()->realEmissionME()->pdf2(emissionScale,theExtrapolationX,
 				       xzMismatch);
     double realPDF = 
       dipole()->bornEmitter() == 0 ?
       dipole()->realEmissionME()->pdf1(emissionScale,theExtrapolationX,1.0) :
       dipole()->realEmissionME()->pdf2(emissionScale,theExtrapolationX,1.0);
     if ( realPDF == 0.0 || realCorrectedPDF == 0.0 )
       return ZERO;
     xme2 *= realCorrectedPDF / realPDF;
   }
 
   Energy qtilde = sqrt(vars.first);
   double z = vars.second;
   Energy2 pt2 = ZERO;
   Energy Qg = theQTildeSudakov->kinScale();
   if ( dipole()->bornEmitter() > 1 ) {
     Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
     if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
       pt2 = sqr(z*(1.-z)*qtilde) - sqr(mu);
     else
       pt2 = sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg);
   }
   if ( dipole()->bornEmitter() < 2 ) {
     pt2 = sqr((1.-z)*qtilde) - z*sqr(Qg);
   }
   assert(pt2 >= ZERO);
 
   if ( profileScales() )
     xme2 *= profileScales()->hardScaleProfile(dipole()->showerHardScale(),sqrt(pt2));
 
   CrossSection res = 
     sqr(hbarc) * 
     realXComb()->jacobian() * 
     subtractionScaleWeight() *
     xme2 /
     (2. * realXComb()->lastSHat());
 
   return res;
 
 }
 
 double QTildeMatching::me2() const {
   throw Exception() << "QTildeMatching::me2(): Not intented to use. Disable the ShowerApproximationGenerator."
 		    << Exception::runerror;
   return 0.;
 }
 
 void QTildeMatching::calculateShowerVariables() const {
 
   Lorentz5Momentum n;
   Energy2 Q2 = ZERO;
 
   const Lorentz5Momentum& pb = bornCXComb()->meMomenta()[dipole()->bornEmitter()];
   const Lorentz5Momentum& pc = bornCXComb()->meMomenta()[dipole()->bornSpectator()];
 
   if ( dipole()->bornEmitter() > 1 ) {
     Q2 = (pb+pc).m2();
   } else {
     Q2 = -(pb-pc).m2();
   }
 
   if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) {
     double b = sqr(bornCXComb()->meMomenta()[dipole()->bornEmitter()].m())/Q2;
     double c = sqr(bornCXComb()->meMomenta()[dipole()->bornSpectator()].m())/Q2;
     double lambda = sqrt(1.+sqr(b)+sqr(c)-2.*b-2.*c-2.*b*c);
     n = (1.-0.5*(1.-b+c-lambda))*pc - 0.5*(1.-b+c-lambda)*pb;
   }
 
   if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) {
     n = bornCXComb()->meMomenta()[dipole()->bornSpectator()];
   }
 
   if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) {
     double c = sqr(bornCXComb()->meMomenta()[dipole()->bornSpectator()].m())/Q2;
     n = (1.+c)*pc - c*pb;
   }
 
   if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2 ) {
     n = bornCXComb()->meMomenta()[dipole()->bornSpectator()];
   }
 
   // the light-cone condition is numerically not very stable, so we
   // explicitly push it on the light-cone here
   n.setMass(ZERO);
   n.rescaleEnergy();
 
   double z = 0.0;
 
   if ( dipole()->bornEmitter() > 1 ) {
     z = 1. - 
     (n*realCXComb()->meMomenta()[dipole()->realEmission()])/
     (n*bornCXComb()->meMomenta()[dipole()->bornEmitter()]);
   } else {
     z = 1. - 
     (n*realCXComb()->meMomenta()[dipole()->realEmission()])/
     (n*realCXComb()->meMomenta()[dipole()->realEmitter()]);
   }
   // allow small violations (numerical inaccuracies)
   if ( z <= 0 && z >= -1e-6 ) {
     z = std::numeric_limits<double>::epsilon();
   } else if ( z >= 1 && z <= 1+1e-6 ) {
     z = 1-std::numeric_limits<double>::epsilon();
   }
 
   Energy2 qtilde2 = ZERO;
   Energy2 q2 = ZERO;
 
   if ( dipole()->bornEmitter() > 1 ) {
     q2 = 
       (realCXComb()->meMomenta()[dipole()->realEmitter()] + realCXComb()->meMomenta()[dipole()->realEmission()]).m2();
     qtilde2 = (q2 - bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2())/(z*(1.-z));
   } else {
     q2 = 
       -(realCXComb()->meMomenta()[dipole()->realEmitter()] - realCXComb()->meMomenta()[dipole()->realEmission()]).m2();
     qtilde2 = (q2 + bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2())/(1.-z);
   }
   if ( qtilde2 < ZERO ) {
     qtilde2 = ZERO;
   }
 
   assert(qtilde2 >= ZERO && z > 0.0 && z < 1.0);
 
   dipole()->showerScale(sqrt(qtilde2));
   dipole()->showerParameters().resize(1);
   dipole()->showerParameters()[0] = z;
 
 }
 
 double QTildeMatching::splitFn(const pair<Energy2,double>& vars) const {
   const Energy2& qtilde2 = vars.first;
   const double z = vars.second;
   double Nc = SM().Nc();
 
   // final state branching
   if ( dipole()->bornEmitter() > 1 ) {
     // final state quark quark branching
     if ( abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 7 ) {
       Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
       return
 	((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(z)-2.*sqr(m)/(z*qtilde2))/(1.-z);
     }
     // final state gluon branching
     if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g ) {
       if ( realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
 	// ATTENTION the factor 2 here is intentional as it cancels to the 1/2
 	// stemming from the large-N colour correlator
 	return 2.*Nc*(z/(1.-z)+(1.-z)/z+z*(1.-z));
       }
       if ( abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
 	Energy m = realCXComb()->mePartonData()[dipole()->realEmission()]->hardProcessMass();
 	return (1./2.)*(1.-2.*z*(1.-z)+2.*sqr(m)/(z*(1.-z)*qtilde2));
       }
     }
     // final state squark branching
     if ((abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) > 1000000 && 
 	 abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 1000007) ||
 	(abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) > 2000000 && 
 	 abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 2000007)){
       Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
       return ((sqr(Nc)-1.)/Nc)*(z-sqr(m)/(z*qtilde2))/(1.-z);
     }
     // final state gluino branching
     if (bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == 1000021){
       Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
       return Nc*(1.+sqr(z)-2.*sqr(m)/(z*qtilde2))/(1.-z);
     }
   }
 
   // initial state branching
   if ( dipole()->bornEmitter() < 2 ) {
     // g/g
     if ( realCXComb()->mePartonData()[dipole()->realEmitter()]->id() == ParticleID::g &&
 	 realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
       // see above for factor of 2
       return 2.*Nc*(z/(1.-z)+(1.-z)/z+z*(1.-z));
     }
     // q/q
     if ( abs(realCXComb()->mePartonData()[dipole()->realEmitter()]->id()) < 7 &&
 	 realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
       return
 	((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(z))/(1.-z);
     }
     // g/q
     if ( realCXComb()->mePartonData()[dipole()->realEmitter()]->id() == ParticleID::g &&
 	 abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
       return (1./2.)*(1.-2.*z*(1.-z));
     }
     // q/g
     if ( abs(realCXComb()->mePartonData()[dipole()->realEmitter()]->id()) < 7 &&
 	 abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
       return
 	((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(1.-z))/z;
     }
   }
 
   return 0.0;
 
 }
 
 // If needed, insert default implementations of virtual function defined
 // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
 
 void QTildeMatching::doinit() {
   assert(theShowerHandler && theQTildeFinder && theQTildeSudakov);
   theShowerHandler->init();
   theQTildeFinder->init();
   theQTildeSudakov->init();
   hardScaleFactor(theShowerHandler->hardScaleFactor());
   factorizationScaleFactor(theShowerHandler->factorizationScaleFactor());
   renormalizationScaleFactor(theShowerHandler->renormalizationScaleFactor());
   profileScales(theShowerHandler->profileScales());
   restrictPhasespace(theShowerHandler->restrictPhasespace());
   hardScaleIsMuF(theShowerHandler->hardScaleIsMuF());
   ShowerApproximation::doinit();
 }
 
 void QTildeMatching::doinitrun() {
   assert(theShowerHandler && theQTildeFinder && theQTildeSudakov);
   theShowerHandler->initrun();
   theQTildeFinder->initrun();
   theQTildeSudakov->initrun();
   ShowerApproximation::doinitrun();
 }
 
 
 void QTildeMatching::persistentOutput(PersistentOStream & os) const {
   os << theQTildeFinder << theQTildeSudakov
      << theShowerHandler << theCorrectForXZMismatch;
 }
 
 void QTildeMatching::persistentInput(PersistentIStream & is, int) {
   is >> theQTildeFinder >> theQTildeSudakov
      >> theShowerHandler >> theCorrectForXZMismatch;
 }
 
 
 // *** Attention *** The following static variable is needed for the type
 // description system in ThePEG. Please check that the template arguments
 // are correct (the class and its base class), and that the constructor
 // arguments are correct (the class name and the name of the dynamically
 // loadable library where the class implementation can be found).
 DescribeClass<QTildeMatching,Herwig::ShowerApproximation>
   describeHerwigQTildeMatching("Herwig::QTildeMatching", "HwShower.so HwQTildeMatching.so");
 
 void QTildeMatching::Init() {
 
   static ClassDocumentation<QTildeMatching> documentation
     ("QTildeMatching implements NLO matching with the default shower.");
 
-  static Reference<QTildeMatching,QTildeFinder> interfaceQTildeFinder
+  static Reference<QTildeMatching,PartnerFinder> interfaceQTildeFinder
     ("QTildeFinder",
      "Set the partner finder to calculate hard scales.",
      &QTildeMatching::theQTildeFinder, false, false, true, false, false);
   interfaceQTildeFinder.rank(-1);
 
   static Reference<QTildeMatching,QTildeSudakov> interfaceQTildeSudakov
     ("QTildeSudakov",
      "Set the partner finder to calculate hard scales.",
      &QTildeMatching::theQTildeSudakov, false, false, true, false, false);
   interfaceQTildeSudakov.rank(-1);
 
   static Reference<QTildeMatching,ShowerHandler> interfaceShowerHandler
     ("ShowerHandler",
      "The QTilde shower handler to use.",
      &QTildeMatching::theShowerHandler, false, false, true, true, false);
   interfaceShowerHandler.rank(-1);
 
   static Switch<QTildeMatching,bool> interfaceCorrectForXZMismatch
     ("CorrectForXZMismatch",
      "Correct for x/z mismatch near hard phase space boundary.",
      &QTildeMatching::theCorrectForXZMismatch, true, false, false);
   static SwitchOption interfaceCorrectForXZMismatchYes
     (interfaceCorrectForXZMismatch,
      "Yes",
      "Include the correction factor.",
      true);
   static SwitchOption interfaceCorrectForXZMismatchNo
     (interfaceCorrectForXZMismatch,
      "No",
      "Do not include the correction factor.",
      false);
   interfaceCorrectForXZMismatch.rank(-1);
 
 }
 
diff --git a/MatrixElement/Matchbox/Matching/QTildeMatching.h b/MatrixElement/Matchbox/Matching/QTildeMatching.h
--- a/MatrixElement/Matchbox/Matching/QTildeMatching.h
+++ b/MatrixElement/Matchbox/Matching/QTildeMatching.h
@@ -1,201 +1,201 @@
 // -*- C++ -*-
 //
 // QTildeMatching.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_QTildeMatching_H
 #define Herwig_QTildeMatching_H
 //
 // This is the declaration of the QTildeMatching class.
 //
 
 #include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
 #include "Herwig/Shower/ShowerHandler.h"
-#include "Herwig/Shower/QTilde/Default/QTildeFinder.h"
+#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
 #include "Herwig/Shower/QTilde/Default/QTildeSudakov.h"
 
 namespace Herwig {
 
 using namespace ThePEG;
 
 /**
  * \ingroup Matchbox
  * \author Simon Platzer
  *
  * \brief QTildeMatching implements NLO matching with the default shower.
  *
  */
 class QTildeMatching: public Herwig::ShowerApproximation {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * The default constructor.
    */
   QTildeMatching();
 
   /**
    * The destructor.
    */
   virtual ~QTildeMatching();
   //@}
 
 public:
 
   /**
    * Return the shower approximation to the real emission cross
    * section for the given pair of Born and real emission
    * configurations.
    */
   virtual CrossSection dSigHatDR() const;
 
   /**
    * Return the shower approximation splitting kernel for the given
    * pair of Born and real emission configurations in units of the
    * Born center of mass energy squared, and including a weight to
    * project onto the splitting given by the dipole used.
    */
   virtual double me2() const;
 
   /**
    * Determine if the configuration is below or above the cutoff.
    */
   virtual void checkCutoff();
 
   /**
    * Determine all kinematic variables which are not provided by the
    * dipole kinematics; store all shower variables in the respective
    * dipole object for later use.
    */
   virtual void getShowerVariables();
 
 protected:
 
   /**
    * Return true, if the shower was able to generate an emission
    * leading from the given Born to the given real emission process.
    */
   virtual bool isInShowerPhasespace() const;
 
   /**
    * Return true, if the shower emission leading from the given Born
    * to the given real emission process would have been generated
    * above the shower's infrared cutoff.
    */
   virtual bool isAboveCutoff() const;
 
   /**
    * Calculate qtilde^2 and z for the splitting considered
    */
   void calculateShowerVariables() const;
 
   /**
    * Return the splitting function as a function of the kinematic
    * variables
    */
   double splitFn(const pair<Energy2,double>&) const; 
 
 public:
 
   /** @name Functions used by the persistent I/O system. */
   //@{
   /**
    * Function used to write out object persistently.
    * @param os the persistent output stream written to.
    */
   void persistentOutput(PersistentOStream & os) const;
 
   /**
    * Function used to read in object persistently.
    * @param is the persistent input stream read from.
    * @param version the version number of the object when written.
    */
   void persistentInput(PersistentIStream & is, int version);
   //@}
 
   /**
    * The standard Init function used to initialize the interfaces.
    * Called exactly once for each class by the class description system
    * before the main function starts or
    * when this class is dynamically loaded.
    */
   static void Init();
 
 protected:
 
   /** @name Clone Methods. */
   //@{
   /**
    * Make a simple clone of this object.
    * @return a pointer to the new object.
    */
   virtual IBPtr clone() const;
 
   /** Make a clone of this object, possibly modifying the cloned object
    * to make it sane.
    * @return a pointer to the new object.
    */
   virtual IBPtr fullclone() const;
   //@}
 
 
 // If needed, insert declarations of virtual function defined in the
 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
 
 protected:
 
   /** @name Standard Interfaced functions. */
   //@{
   /**
    * Initialize this object after the setup phase before saving an
    * EventGenerator to disk.
    * @throws InitException if object could not be initialized properly.
    */
   virtual void doinit();
 
   /**
    * Initialize this object. Called in the run phase just before
    * a run begins.
    */
   virtual void doinitrun();
   //@}
 
 private:
 
   /**
    * The assignment operator is private and must never be called.
    * In fact, it should not even be implemented.
    */
   QTildeMatching & operator=(const QTildeMatching &);
 
   /**
    * The shower handler to be used
    */
   Ptr<ShowerHandler>::ptr theShowerHandler;
 
   /**
    * The qtilde partner finder for calculating the hard scales
    */
-  Ptr<QTildeFinder>::ptr theQTildeFinder;
+  Ptr<PartnerFinder>::ptr theQTildeFinder;
 
   /**
    * The qtilde Sudakov to access the cutoff
    */
   Ptr<QTildeSudakov>::ptr theQTildeSudakov;
 
   /**
    * True, if PDF weight should be corrected for z/x mismatch at the
    * hard phase space boundary
    */
   bool theCorrectForXZMismatch;
 
 };
 
 }
 
 #endif /* Herwig_QTildeMatching_H */
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,536 +1,757 @@
 // -*- 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;
 
-DescribeAbstractClass<PartnerFinder,Interfaced>
+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 unsigned int
   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 unsigned int
   ACLSIZE(const tShowerParticlePtr a) {
     return a->colourInfo()->antiColourLines().size();
   }
 }
 
 void PartnerFinder::persistentOutput(PersistentOStream & os) const {
-  os << partnerMethod_ << QEDPartner_ << scaleChoice_;
+  os << partnerMethod_ << QEDPartner_ << scaleChoice_
+     << _finalFinalConditions << _initialFinalDecayConditions
+     << _initialInitialConditions;
 }
 
 void PartnerFinder::persistentInput(PersistentIStream & is, int) {
-  is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_;
+  is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_
+     >> _finalFinalConditions >> _initialFinalDecayConditions
+     >>_initialInitialConditions;
 }
 
 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);
 
+
+  static Switch<PartnerFinder,unsigned int> interfaceFinalFinalConditions
+    ("FinalFinalConditions",
+     "The initial conditions for the shower of a final-final colour connection",
+     &PartnerFinder::_finalFinalConditions, 0, false, false);
+  static SwitchOption interfaceFinalFinalConditionsSymmetric
+    (interfaceFinalFinalConditions,
+     "Symmetric",
+     "The symmetric choice",
+     0);
+  static SwitchOption interfaceFinalFinalConditionsColoured
+    (interfaceFinalFinalConditions,
+     "Coloured",
+     "Maximal radiation from the coloured particle",
+     1);
+  static SwitchOption interfaceFinalFinalConditionsAntiColoured
+    (interfaceFinalFinalConditions,
+     "AntiColoured",
+     "Maximal emission from the anticoloured particle",
+     2);
+  static SwitchOption interfaceFinalFinalConditionsRandom
+    (interfaceFinalFinalConditions,
+     "Random",
+     "Randomly selected maximal emission from one of the particles",
+     3);
+
+  static Switch<PartnerFinder,unsigned int> interfaceInitialFinalDecayConditions
+    ("InitialFinalDecayConditions",
+     "The initial conditions for the shower of an initial-final"
+     " decay colour connection.",
+     &PartnerFinder::_initialFinalDecayConditions, 0, false, false);
+  static SwitchOption interfaceInitialFinalDecayConditionsSymmetric
+    (interfaceInitialFinalDecayConditions,
+     "Symmetric",
+     "The symmetric choice",
+     0);
+  static SwitchOption interfaceInitialFinalDecayConditionsMaximal
+    (interfaceInitialFinalDecayConditions,
+     "Maximal",
+     "Maximal radiation from the decay product",
+     1);
+  static SwitchOption interfaceInitialFinalDecayConditionsSmooth
+    (interfaceInitialFinalDecayConditions,
+     "Smooth",
+     "Smooth matching in the soft limit",
+     2);
+
+  static Switch<PartnerFinder,unsigned int> interfaceInitialInitialConditions
+    ("InitialInitialConditions",
+     "The initial conditions for the shower of an initial-initial"
+     " colour connection.",
+     &PartnerFinder::_initialInitialConditions, 0, false, false);
+  static SwitchOption interfaceInitialInitialConditionsSymmetric
+    (interfaceInitialInitialConditions,
+     "Symmetric",
+     "The symmetric choice",
+     0);
+  static SwitchOption interfaceInitialInitialConditionsMaximiseB
+    (interfaceInitialInitialConditions,
+     "MaximiseB",
+     "Maximal radiation from parton b",
+     1);
+  static SwitchOption interfaceInitialInitialConditionsMaximiseC
+    (interfaceInitialInitialConditions,
+     "MaximiseC",
+     "Maximal radiation from parton c",
+     2);
 }
 
 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 {
     setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
     setInitialQEDEvolutionScales(particles,isDecayCase,false);
   }
   // 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) {
   // set the partners and the scales
   if(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!)
     ShowerParticleVector::const_iterator cit, cjt;
     for(cit = particles.begin(); cit != particles.end(); ++cit) {
       // Skip colourless particles
       if(!(*cit)->data().coloured()) continue;
       // find the partners
       vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners = 
 	findQCDPartners(*cit,particles);
       // must have a partner
       if(partners.empty()) {
 	throw Exception() << "`Failed to make colour connections in " 
                           << "PartnerFinder::setQCDInitialEvolutionScales"
                           << (**cit)
                           << 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(*cit,partners[ix].second),
 							 isDecayCase));
       }
       // 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.
       // store the choice
       int position(-1);
       // 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 ((*cit)->perturbative() == 1 && 
 	    (*cit)->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 = (*cit)->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
       (*cit)->partner(partners[position].second);
       for(unsigned int ix=0;ix<partners.size();++ix) {
 	(**cit).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) {
 	  (**cit).scales().QCD_c = 
 	    (**cit).scales().QCD_c_noAO = 
 	    (scaleChoice_==0 ? scale : scales[ix].first);
 	}
 	else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) {
 	  (**cit).scales().QCD_ac = 
 	    (**cit).scales().QCD_ac_noAO =
 	    (scaleChoice_==0 ? scale : scales[ix].first);
 	}
 	else
 	  assert(false);
       }
     }
   }
   // primary partner set, set the others and do the scale
   else {
     for(ShowerParticleVector::const_iterator cit = particles.begin();
         cit != particles.end(); ++cit) {
       // Skip colourless particles
       if(!(*cit)->data().coloured()) continue;
       // find the partners
       vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners = 
 	findQCDPartners(*cit,particles);
       // must have a partner
       if(partners.empty()) {
         throw Exception() << "`Failed to make colour connections in " 
                           << "PartnerFinder::setQCDInitialEvolutionScales"
                           << (**cit)
                           << Exception::eventerror;
       }
       // Calculate the evolution scales for all possible pairs of of particles
       vector<pair<Energy,Energy> > scales;
       int position(-1);
       for(unsigned int ix=0;ix< partners.size();++ix) {
 	if(partners[ix].second) position = ix;
 	scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
 							 isDecayCase));
       }
       assert(position>=0);
       for(unsigned int ix=0;ix<partners.size();++ix) {
 	(**cit).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) {
 	  (**cit).scales().QCD_c = 
 	    (**cit).scales().QCD_c_noAO = 
 	    (scaleChoice_==0 ? scale : scales[ix].first);
 	}
 	else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) {
 	  (**cit).scales().QCD_ac = 
 	    (**cit).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(ShowerParticleVector::const_iterator cit = particles.begin();
       cit != particles.end(); ++cit) {
     // not charged or photon continue
     if(!(**cit).dataPtr()->charged()) continue;
     // find the potential partners
     vector<pair<double,tShowerParticlePtr> > partners = findQEDPartners(*cit,particles,isDecayCase);
     if(partners.empty()) {
       throw Exception() << "Failed to find partner in " 
 			<< "PartnerFinder::setQEDInitialEvolutionScales"
 			<< (**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::setQEDInitialEvolutionScales"
 				 << (**cit) << 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(*cit,partners[ix].second),
 						       isDecayCase));
     }
     // 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::QED,
 							  scales[ix].first));
     }
     // set scales
     (**cit).scales().QED      = scales[position].first;
     (**cit).scales().QED_noAO = scales[position].first;
   }
 }
 
 pair<Energy,Energy> PartnerFinder::
 calculateInitialEvolutionScales(const ShowerPPair &particlePair,
                                 const bool isDecayCase) {
   bool FS1=FS(particlePair.first),FS2= FS(particlePair.second);
 
-  if(FS1 && FS2)
-    return calculateFinalFinalScales(particlePair);
+  if(FS1 && FS2) {
+    bool colouredFirst =
+     particlePair.first->colourLine()&&
+     particlePair.first->colourLine()==particlePair.second->antiColourLine();
+   return calculateFinalFinalScales(particlePair.first->momentum(),
+                                    particlePair.second->momentum(),
+ 				   colouredFirst);
+  }
   else if(FS1 && !FS2) {
     ShowerPPair a(particlePair.second, particlePair.first);
-    pair<Energy,Energy> rval = calculateInitialFinalScales(a,isDecayCase);
+    pair<Energy,Energy> rval = calculateInitialFinalScales(a.first->momentum(),
+                                                           a.second->momentum(),
+                                                           isDecayCase);
     return pair<Energy,Energy>(rval.second,rval.first);
   }
   else if(!FS1 &&FS2)
-    return calculateInitialFinalScales(particlePair,isDecayCase);
+    return calculateInitialFinalScales(particlePair.first->momentum(),particlePair.second->momentum(),isDecayCase);
   else
-    return calculateInitialInitialScales(particlePair);
+    return calculateInitialInitialScales(particlePair.first->momentum(),particlePair.second->momentum());
 }
 
 vector< pair<ShowerPartnerType, tShowerParticlePtr> > 
 PartnerFinder::findQCDPartners(tShowerParticlePtr particle,
 			       const ShowerParticleVector &particles) {
   vector< pair<ShowerPartnerType, tShowerParticlePtr> > partners;
   ShowerParticleVector::const_iterator cjt;
   for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
     if(!(*cjt)->data().coloured() || particle==*cjt) continue;
     // one initial-state and one final-state particle
     if(FS(particle) != FS(*cjt)) {
       // loop over all the colours of both particles
       for(unsigned int ix=0; ix<CLSIZE(particle); ++ix) {
 	for(unsigned int jx=0; jx<CLSIZE(*cjt); ++jx) {
 	  if((CL(particle,ix) && CL(particle,ix)==CL(*cjt,jx))) {
 	    partners.push_back(make_pair(ShowerPartnerType::    QCDColourLine,*cjt));
 	  }
 	}
       }
       //loop over all the anti-colours of both particles
       for(unsigned int ix=0; ix<ACLSIZE(particle); ++ix) {
 	for(unsigned int jx=0; jx<ACLSIZE(*cjt); ++jx) {
 	  if((ACL(particle,ix) && ACL(particle,ix)==ACL(*cjt,jx))) {
 	    partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
 	  }
 	}
       }
     }
     // 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(unsigned int ix=0; ix<CLSIZE(particle); ++ix){
 	for(unsigned int jx=0; jx<ACLSIZE(*cjt); ++jx){
 	  if(CL(particle,ix) && CL(particle,ix)==ACL(*cjt,jx)) {
 	    partners.push_back(make_pair(ShowerPartnerType::    QCDColourLine,*cjt));
 	  }
 	}
       }
       //loop over the anti-colours of the first particle and the colours of the other
       for(unsigned int ix=0; ix<ACLSIZE(particle); ++ix){
 	for(unsigned int jx=0; jx<CLSIZE(*cjt); jx++){
 	  if(ACL(particle,ix) && ACL(particle,ix)==CL(*cjt,jx)) {
 	    partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
 	  }
 	}
       }
     }
   }
   // 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(cjt=particles.begin();cjt!=particles.end();++cjt) {
 	if(( FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt)  == cpair.second))||
 	   (!FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second ))) {
 	  partners.push_back(make_pair(ShowerPartnerType::    QCDColourLine,*cjt));
 	}
       }
     }
     else if(col&&col->sinkNeighbours().first) {
       tColinePair cpair = col->sinkNeighbours();
       for(cjt=particles.begin();cjt!=particles.end();++cjt) {
 	if(( FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt)  == cpair.second))||
 	   (!FS(*cjt) && ( CL(*cjt) == cpair.first ||  CL(*cjt) == cpair.second))) {
 	  partners.push_back(make_pair(ShowerPartnerType::    QCDColourLine,*cjt));    
 	}
       }
     }
     col = ACL(particle);
     if(FS(particle)&&col&&col->sinkNeighbours().first) {
       tColinePair cpair = col->sinkNeighbours();
       for(cjt=particles.begin();cjt!=particles.end();++cjt) {
 	if(( FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt)  == cpair.second))||
 	   (!FS(*cjt) && ( CL(*cjt) == cpair.first ||  CL(*cjt) == cpair.second ))) {
 	  partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));  
 	}
       }
     }
     else if(col&&col->sourceNeighbours().first) {
       tColinePair cpair = col->sourceNeighbours();
       for(cjt=particles.begin();cjt!=particles.end();++cjt) {
 	if(( FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))||
 	   (!FS(*cjt) && (ACL(*cjt) == cpair.first ||ACL(*cjt) == cpair.second))) {
 	  partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));     
 	}
       }
     }
   }
   // return the partners
   return partners;
 }
 
 vector< pair<double, tShowerParticlePtr> > 
 PartnerFinder::findQEDPartners(tShowerParticlePtr particle,
 			       const ShowerParticleVector &particles,
 			       const bool isDecayCase) {
   vector< pair<double, tShowerParticlePtr> > partners;
   ShowerParticleVector::const_iterator cjt;
   double pcharge = particle->id()==ParticleID::gamma ? 1 : double(particle->data().iCharge());
   vector< pair<double, tShowerParticlePtr> > photons;
   for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
     if(particle == *cjt) continue;
     if((**cjt).id()==ParticleID::gamma) photons.push_back(make_pair(1.,*cjt));
     if(!(*cjt)->data().charged() ) continue;
     double charge = pcharge*double((*cjt)->data().iCharge());
     if( FS(particle) != FS(*cjt) ) charge *=-1.;
     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;
     }
     if(particle->id()==ParticleID::gamma) charge = -abs(charge);
     // only keep positive dipoles
     if(charge<0.) partners.push_back(make_pair(-charge,*cjt));
   }
   if(particle->id()==ParticleID::gamma&& partners.empty()) {
     return photons;
   }
   return partners;
 }
+
+pair<Energy,Energy> 
+PartnerFinder::calculateFinalFinalScales(
+        const Lorentz5Momentum & p1,
+        const Lorentz5Momentum & p2,
+	bool colouredFirst) 
+{
+  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. 
+  double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c))
+                 *(sqrt(b)-1.-sqrt(c))*(sqrt(c)-1.-sqrt(b)));
+  // symmetric case
+  unsigned int iopt=finalFinalConditions();
+  Energy firstQ,secondQ;
+  if(iopt==0) {
+    firstQ  = sqrt(0.5*Q2*(1.+b-c+lam));
+    secondQ = sqrt(0.5*Q2*(1.-b+c+lam));
+  }
+  // assymetric choice
+  else {
+    double kappab,kappac;
+    // calculate kappa with coloured line getting maximum
+    if((iopt==1&&colouredFirst)|| // first particle coloured+maximal for coloured
+       (iopt==2&&!colouredFirst)|| // first particle anticoloured+maximal for acoloured
+       (iopt==3&&UseRandom::rndbool(0.5))) { // random choice
+      kappab=4.*(1.-2.*sqrt(c)-b+c);
+      kappac=c+0.25*sqr(1.-b-c+lam)/(kappab-b);
+    }
+    else {
+      kappac=4.*(1.-2.*sqrt(b)-c+b);
+      kappab=b+0.25*sqr(1.-b-c+lam)/(kappac-c);
+    }
+    // calculate the scales
+    firstQ  = sqrt(Q2*kappab);
+    secondQ = sqrt(Q2*kappac);
+  }
+  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 
+    Energy2  mc2 = sqr(pc.mass());
+    Energy2  Q2  = -(pb-pc).m2();
+    return pair<Energy,Energy>(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.));
+    double PROD     = 0.25*sqr(1. - a + c + lambda);
+    double ktilde_b, ktilde_c,cosi(0.);
+    switch(initialFinalDecayConditions()) {
+    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 is a problem if very small here use 1GeV as minimum
+      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;
+    default:
+      throw Exception() << "Invalid option for decay shower's phase space"
+			<< " PartnerFinder::calculateInitialFinalScales"
+			<< Exception::abortnow;
+    }
+    return pair<Energy,Energy>(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
+  Energy Q = sqrt((p1+p2).m2());
+  if(_initialInitialConditions==1) {
+    return pair<Energy,Energy>(sqrt(2.0)*Q,sqrt(0.5)*Q);
+  } else if(_initialInitialConditions==2) {
+    return pair<Energy,Energy>(sqrt(0.5)*Q,sqrt(2.0)*Q);
+  } else {
+    return pair<Energy,Energy>(Q,Q);
+  }
+}
diff --git a/Shower/QTilde/Base/PartnerFinder.h b/Shower/QTilde/Base/PartnerFinder.h
--- a/Shower/QTilde/Base/PartnerFinder.h
+++ b/Shower/QTilde/Base/PartnerFinder.h
@@ -1,208 +1,297 @@
 // -*- C++ -*-
 //
 // PartnerFinder.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_PartnerFinder_H
 #define HERWIG_PartnerFinder_H
 //
 // This is the declaration of the PartnerFinder class.
 //
 #include "Herwig/Shower/QTilde/ShowerConfig.h"
 #include "ThePEG/Interface/Interfaced.h"
 #include "PartnerFinder.fh"
 
 namespace Herwig {
 
 using namespace ThePEG;
 
 /**
  *  typedef of a pair of particle for calculating the evolution scales
  */
 typedef pair<tShowerParticlePtr,tShowerParticlePtr> ShowerPPair;
 
 /** \ingroup Shower
  *
  *  This class is responsible of two related tasks: 
  *  -  it finds the partners 
  *  -  for each pair of partners (and interaction therefore)
  *      it sets the initial evolution scales of both of them.
  *
  *  In general the finding of the partners is performed by this class but
  *  the calculation of the initial evolution scales should be implemented
  *  for different shower evolution models in classes inheriting from this one.
  *  Notice that a given particle has, in general, a different partner
  *  for each different interaction; however, given a partner, its 
  *  initial evolution scale, Q, is purely a kinematical relationship 
  *  between the pair, without dependence on the dynamics (i.e. type of interaction).
  *
  * @see \ref PartnerFinderInterfaces "The interfaces"
  * defined for PartnerFinder.
  */
 class PartnerFinder: public Interfaced {
 
 public:
 
   /**
    * The default constructor.
    */
-  PartnerFinder() : partnerMethod_(0), QEDPartner_(0), scaleChoice_(0) {}
+  PartnerFinder() : partnerMethod_(0), QEDPartner_(0), scaleChoice_(0),
+                    _finalFinalConditions(0),
+		    _initialFinalDecayConditions(0),
+		    _initialInitialConditions(0) {}
 
   /**
    * Given in input a collection of particles (ShowerParticle objects),
    * each of these methods set the initial evolution scales of those particles, 
    * between the ones given in input, that do not have yet their
    * evolution scale set. 
    * The input collection of particles can be either the full collection of 
    * showering particles (kept in the main class ShowerHandler,
    * in the case isDecayCase is false, or simply, in the case isDecayCase 
    * is true, the decaying particle and its decay products.    
    * The methods returns true, unless something wrong (inconsistencies,
    * or undefined values) happens.
    *
    * These methods are virtual but in most cases inheriting classes should not
    * need to overide them as they simply find the relevant partner and call
    * one of the calculateScale members to calculate the scale.
    */
   //@{
   /**
    * Set the initial scales
    * @param particles        The particles to be considered
    * @param isDecayCase      Whether or not this is a decay
    * @param setPartners Whether to set the colour partners or just the scales
    */
   virtual void setInitialEvolutionScales(const ShowerParticleVector &particles,
 					 const bool isDecayCase,
 					 ShowerInteraction,
 					 const bool setPartners=true);
   //@}
+protected:
+
+  /** @name Clone Methods. */
+  //@{
+  /**
+   * Make a simple clone of this object.
+   * @return a pointer to the new object.
+   */
+  virtual IBPtr clone() const {return new_ptr(*this);}
+
+  /** Make a clone of this object, possibly modifying the cloned object
+   * to make it sane.
+   * @return a pointer to the new object.
+   */
+  virtual IBPtr fullclone() const {return new_ptr(*this);}
+  //@}
+
 
 public:
 
   /** @name Functions used by the persistent I/O system. */
   //@{
   /**
    * Function used to write out object persistently.
    * @param os the persistent output stream written to.
    */
   void persistentOutput(PersistentOStream & os) const;
 
   /**
    * Function used to read in object persistently.
    * @param is the persistent input stream read from.
    * @param version the version number of the object when written.
    */
   void persistentInput(PersistentIStream & is, int version);
   //@}
 
   /**
    * The standard Init function used to initialize the interfaces.
    * Called exactly once for each class by the class description system
    * before the main function starts or
    * when this class is dynamically loaded.
    */
   static void Init();
 
 protected:
 
   /**
+   * Access function for the initial conditions for the shower
+   */
+  //@{
+  /**
+   * Initial conditions for the shower of a final-final colour connection
+   * - 0 is the symmetric choice
+   * - 1 is maximal emmision from the coloured particle
+   * - 2 is maximal emmision from the anticoloured particle
+   * - 3 is randomly selected maximal emmision
+   */
+  unsigned int finalFinalConditions() const 
+  {return _finalFinalConditions;}
+
+  /**
+   * Initial conditions for the shower of an initial-final decay colour connection
+   * - 0 is the symmetric choice
+   * - 1 is maximal emission from the decay product
+   * - 2 is the smooth choice
+   */
+  unsigned int initialFinalDecayConditions() const
+  {return _initialFinalDecayConditions;}
+  //@}
+
+protected:
+
+  /**
    *  Members to set the scales for different interactions
    */
   //@{
   /**
    *  Set initial scales for a QCD interaction
    */
   virtual void setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
 					    const bool isDecayCase,
 					    const bool setPartners=true);
 
   /**
    *  Set initial scales for a QED interaction
    */
   virtual void setInitialQEDEvolutionScales(const ShowerParticleVector &particles,
 					    const bool isDecayCase,
 					    const bool setPartners=true);
   //@}
 
   /**
    *  Find the QCD partners
    * @param particle The particle to find the partners for
    * @param particles The full set of particles to search
    */
   vector< pair<ShowerPartnerType, tShowerParticlePtr> > 
   findQCDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles);
 
   /**
    *  Find the QED partners
    * @param particle The particle to find the partners for
    * @param particles The full set of particles to search
    */
   vector< pair<double, tShowerParticlePtr> > 
   findQEDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles,
 		  const bool isDecayCase);
 
+public:
   /**
    * Given a pair of particles, supposedly partners w.r.t. an interaction,
    * this method returns their initial evolution scales as a pair.
    * If something wrong happens, it returns the null (ZERO,ZERO) pair. 
    * This method is used by the above setXXXInitialEvolutionScales 
    * methods.
    * These methods must be overiden in inheriting classes
    */
   //@{
   /**
    *  General method to calculate the initial evolution scales
    */
-  virtual pair<Energy,Energy> calculateInitialEvolutionScales(const ShowerPPair &,
-							      const bool isDecayCase);
+  pair<Energy,Energy> calculateInitialEvolutionScales(const ShowerPPair &,
+						      const bool isDecayCase);
 
   /**
-   *  Calculate the initial evolution scales for two final-state particles
+   *  Calculate the initial evolution scales given momenta
    */
-  virtual pair<Energy,Energy> calculateFinalFinalScales(const ShowerPPair &)=0;
+  pair<Energy,Energy> calculateFinalFinalScales(const Lorentz5Momentum & p1, 
+                                                const Lorentz5Momentum & p2,
+						bool colouredfirst);
 
   /**
-   *  Calculate the initial evolution scales for two initial-state particles
+   *  Calculate the initial evolution scales given momenta
    */
-  virtual pair<Energy,Energy> calculateInitialInitialScales(const ShowerPPair &)=0;
+  pair<Energy,Energy> calculateInitialInitialScales(const Lorentz5Momentum& p1, 
+						    const Lorentz5Momentum& p2);
 
   /**
-   *  Calculate the initial evolution scales for one initial 
-   *  and one final-state particles
+   *  Calculate the initial evolution scales given momenta
    */
-  virtual pair<Energy,Energy> calculateInitialFinalScales(const ShowerPPair &,
-							  const bool isDecayCase)=0;
+  pair<Energy,Energy> calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
+						  const bool isDecayCase);
+
+
   //@}
 
 private:
 
   /**
    * The assignment operator is private and must never be called.
    * In fact, it should not even be implemented.
    */
   PartnerFinder & operator=(const PartnerFinder &);
 
 private:
 
   /**
    *  Method for choosing colour partner
    */
    int partnerMethod_;
 
   /**
    *  Choice for the QED radiation partner
    */
   int QEDPartner_;
 
   /**
    *  Choice of the scale
    */
   int scaleChoice_;
+
+  /**
+   *  Flags controlling the initial conditions for the shower
+   */
+  //@{
+  /**
+   * Initial conditions for the shower with a final-final colour
+   * connection
+   */
+  unsigned int _finalFinalConditions; 
+
+  /**
+   * Initial conditions for the shower with an initial-final decay colour
+   * connection. This is done according to the top decay colour 
+   *  connection calculation in JHEP12(2003)_045. The options act as follows:
+   *  0: This is the default 'symmetric' choice which more or less divides
+   *     the phase space evenly between the parent and its charged child.
+   *  1: This 'maximal' choice maximises the phase space available for 
+   *     gluons emitted from the charged child.
+   *  2: This (experimental) 'smooth' choice does not suffer from
+   *     a discontinuity at the boundary between the region populated by
+   *     emissions from the charged child and the region populated by emissions
+   *     from the parent. This does, however, mean that the phase space 
+   *     available for emissions from the charged child is fairly minimal.
+   */
+  unsigned int _initialFinalDecayConditions;
+
+  /**
+   * Initial conditions for the shower with an initial-initial colour
+   * connection. This is done according to the colour connection 
+   * calculation in JHEP12(2003)_045. The options act as follows:
+   *  0: This is the default 'symmetric' choice which more or less divides
+   *     the phase space evenly between the two incoming partons.
+   *  1: This increases the phase space for emission from "parton b".
+   *  2: This increases the phase space for emission from "parton c".
+   */
+  unsigned int _initialInitialConditions;
+  //@}
 };
 
 }
 
 #endif /* HERWIG_PartnerFinder_H */
diff --git a/Shower/QTilde/Default/QTildeFinder.cc b/Shower/QTilde/Default/QTildeFinder.cc
deleted file mode 100644
--- a/Shower/QTilde/Default/QTildeFinder.cc
+++ /dev/null
@@ -1,273 +0,0 @@
-// -*- C++ -*-
-//
-// QTildeFinder.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 QTildeFinder class.
-//
-
-#include "QTildeFinder.h"
-#include "ThePEG/Interface/ClassDocumentation.h"
-#include "ThePEG/Persistency/PersistentOStream.h"
-#include "ThePEG/Persistency/PersistentIStream.h"
-#include "ThePEG/Interface/Switch.h"
-#include "ThePEG/Repository/EventGenerator.h"
-#include "ThePEG/EventRecord/Event.h"
-#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
-#include "ThePEG/Repository/UseRandom.h" 
-#include "ThePEG/Utilities/DescribeClass.h"
-
-using namespace Herwig;
-
-DescribeClass<QTildeFinder,Herwig::PartnerFinder>
-describeQTildeFinder ("Herwig::QTildeFinder","HwShower.so");
-
-void QTildeFinder::persistentOutput(PersistentOStream & os) const {
-  os << _finalFinalConditions << _initialFinalDecayConditions
-     << _initialInitialConditions;
-}
-
-void QTildeFinder::persistentInput(PersistentIStream & is, int) {
-  is >> _finalFinalConditions >> _initialFinalDecayConditions
-     >>_initialInitialConditions;
-}
-
-void QTildeFinder::Init() {
-
-  static ClassDocumentation<QTildeFinder> 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<QTildeFinder,unsigned int> interfaceFinalFinalConditions
-    ("FinalFinalConditions",
-     "The initial conditions for the shower of a final-final colour connection",
-     &QTildeFinder::_finalFinalConditions, 0, false, false);
-  static SwitchOption interfaceFinalFinalConditionsSymmetric
-    (interfaceFinalFinalConditions,
-     "Symmetric",
-     "The symmetric choice",
-     0);
-  static SwitchOption interfaceFinalFinalConditionsColoured
-    (interfaceFinalFinalConditions,
-     "Coloured",
-     "Maximal radiation from the coloured particle",
-     1);
-  static SwitchOption interfaceFinalFinalConditionsAntiColoured
-    (interfaceFinalFinalConditions,
-     "AntiColoured",
-     "Maximal emission from the anticoloured particle",
-     2);
-  static SwitchOption interfaceFinalFinalConditionsRandom
-    (interfaceFinalFinalConditions,
-     "Random",
-     "Randomly selected maximal emission from one of the particles",
-     3);
-
-  static Switch<QTildeFinder,unsigned int> interfaceInitialFinalDecayConditions
-    ("InitialFinalDecayConditions",
-     "The initial conditions for the shower of an initial-final"
-     " decay colour connection.",
-     &QTildeFinder::_initialFinalDecayConditions, 0, false, false);
-  static SwitchOption interfaceInitialFinalDecayConditionsSymmetric
-    (interfaceInitialFinalDecayConditions,
-     "Symmetric",
-     "The symmetric choice",
-     0);
-  static SwitchOption interfaceInitialFinalDecayConditionsMaximal
-    (interfaceInitialFinalDecayConditions,
-     "Maximal",
-     "Maximal radiation from the decay product",
-     1);
-  static SwitchOption interfaceInitialFinalDecayConditionsSmooth
-    (interfaceInitialFinalDecayConditions,
-     "Smooth",
-     "Smooth matching in the soft limit",
-     2);
-
-  static Switch<QTildeFinder,unsigned int> interfaceInitialInitialConditions
-    ("InitialInitialConditions",
-     "The initial conditions for the shower of an initial-initial"
-     " colour connection.",
-     &QTildeFinder::_initialInitialConditions, 0, false, false);
-  static SwitchOption interfaceInitialInitialConditionsSymmetric
-    (interfaceInitialInitialConditions,
-     "Symmetric",
-     "The symmetric choice",
-     0);
-  static SwitchOption interfaceInitialInitialConditionsMaximiseB
-    (interfaceInitialInitialConditions,
-     "MaximiseB",
-     "Maximal radiation from parton b",
-     1);
-  static SwitchOption interfaceInitialInitialConditionsMaximiseC
-    (interfaceInitialInitialConditions,
-     "MaximiseC",
-     "Maximal radiation from parton c",
-     2);
-}
-
-pair<Energy,Energy> QTildeFinder::
-calculateInitialFinalScales(const ShowerPPair &ppair, const bool isDecayCase) {
-  return
-    calculateInitialFinalScales(ppair.first->momentum(),ppair.second->momentum(),isDecayCase);
-}
-
-pair<Energy,Energy> QTildeFinder::
-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 
-    Energy2  mc2 = sqr(pc.mass());
-    Energy2  Q2  = -(pb-pc).m2();
-    return pair<Energy,Energy>(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.));
-    double PROD     = 0.25*sqr(1. - a + c + lambda);
-    double ktilde_b, ktilde_c,cosi(0.);
-    switch(initialFinalDecayConditions()) {
-    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 is a problem if very small here use 1GeV as minimum
-      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;
-    default:
-      throw Exception() << "Invalid option for decay shower's phase space"
-			<< " QTildeFinder::calculateInitialFinalScales"
-			<< Exception::abortnow;
-    }
-    return pair<Energy,Energy>(sqrt(mb2*ktilde_b),sqrt(mb2*ktilde_c));
-  }
-}
-
-pair<Energy,Energy> QTildeFinder::
-calculateInitialInitialScales(const ShowerPPair &ppair) {
-  return
-    calculateInitialInitialScales(ppair.first->momentum(),
-				  ppair.second->momentum());
-}
-
-pair<Energy,Energy> QTildeFinder::
-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
-  Energy Q = sqrt((p1+p2).m2());
-  if(_initialInitialConditions==1) {
-    return pair<Energy,Energy>(sqrt(2.0)*Q,sqrt(0.5)*Q);
-  } else if(_initialInitialConditions==2) {
-    return pair<Energy,Energy>(sqrt(0.5)*Q,sqrt(2.0)*Q);
-  } else {
-    return pair<Energy,Energy>(Q,Q);
-  }
-}
-
-pair<Energy,Energy> QTildeFinder::
-calculateFinalFinalScales(const ShowerPPair & pp) {
-  bool colouredFirst =
-    pp.first->colourLine()&&
-    pp.first->colourLine()==pp.second->antiColourLine();
-  return calculateFinalFinalScales(pp.first->momentum(),pp.second->momentum(),
-				   colouredFirst);
-}
-
-pair<Energy,Energy> QTildeFinder::
-calculateFinalFinalScales(Lorentz5Momentum p1, Lorentz5Momentum p2,
-			  bool colouredFirst) {
-  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 QTildeFinder::calculateFinalFinalScales()"
-			<< Exception::eventerror;
-    }
-    b = 0.;
-  }
-  if(c<0.) {
-    if(c<-eps) {
-      throw Exception() << "Negative Mass squared c = " << c
-			<< "in QTildeFinder::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. 
-  double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c))
-                 *(sqrt(b)-1.-sqrt(c))*(sqrt(c)-1.-sqrt(b)));
-  // symmetric case
-  unsigned int iopt=finalFinalConditions();
-  Energy firstQ,secondQ;
-  if(iopt==0) {
-    firstQ  = sqrt(0.5*Q2*(1.+b-c+lam));
-    secondQ = sqrt(0.5*Q2*(1.-b+c+lam));
-  }
-  // assymetric choice
-  else {
-    double kappab,kappac;
-    // calculate kappa with coloured line getting maximum
-    if((iopt==1&&colouredFirst)|| // first particle coloured+maximal for coloured
-       (iopt==2&&!colouredFirst)|| // first particle anticoloured+maximal for acoloured
-       (iopt==3&&UseRandom::rndbool(0.5))) { // random choice
-      kappab=4.*(1.-2.*sqrt(c)-b+c);
-      kappac=c+0.25*sqr(1.-b-c+lam)/(kappab-b);
-    }
-    else {
-      kappac=4.*(1.-2.*sqrt(b)-c+b);
-      kappab=b+0.25*sqr(1.-b-c+lam)/(kappac-c);
-    }
-    // calculate the scales
-    firstQ  = sqrt(Q2*kappab);
-    secondQ = sqrt(Q2*kappac);
-  }
-  return pair<Energy,Energy>(firstQ, secondQ);
-}
diff --git a/Shower/QTilde/Default/QTildeFinder.h b/Shower/QTilde/Default/QTildeFinder.h
deleted file mode 100644
--- a/Shower/QTilde/Default/QTildeFinder.h
+++ /dev/null
@@ -1,206 +0,0 @@
-// -*- C++ -*-
-//
-// QTildeFinder.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_QTildeFinder_H
-#define HERWIG_QTildeFinder_H
-//
-// This is the declaration of the QTildeFinder class.
-//
-
-#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
-#include "Herwig/Shower/QTilde/ShowerConfig.h"
-#include "ThePEG/Interface/Interfaced.h"
-
-namespace Herwig {
-using namespace ThePEG;
-
-/** \ingroup Shower
- *
- *  The QTildeFinder class is responsible for finding the partners and
- *  setting the initial evolution scales for the shower evolution described
- *  in JHEP 0312:045,2003.
- *
- * @see \ref QTildeFinderInterfaces "The interfaces"
- * defined for QTildeFinder.
- */
-class QTildeFinder: public PartnerFinder {
-
-public:
-
-  /**
-   * The default constructor.
-   */
-  QTildeFinder() :  _finalFinalConditions(0),
-		    _initialFinalDecayConditions(0),
-		    _initialInitialConditions(0) {}
-
-  /** @name Functions used by the persistent I/O system. */
-  //@{
-  /**
-   * Function used to write out object persistently.
-   * @param os the persistent output stream written to.
-   */
-  void persistentOutput(PersistentOStream & os) const;
-
-  /**
-   * Function used to read in object persistently.
-   * @param is the persistent input stream read from.
-   * @param version the version number of the object when written.
-   */
-  void persistentInput(PersistentIStream & is, int version);
-  //@}
-
-  /**
-   * The standard Init function used to initialize the interfaces.
-   * Called exactly once for each class by the class description system
-   * before the main function starts or
-   * when this class is dynamically loaded.
-   */
-  static void Init();
-
-public:
-
-  /**
-   *  Calculate the initial evolution scales given momenta
-   */
-  pair<Energy,Energy> calculateFinalFinalScales(Lorentz5Momentum p1, Lorentz5Momentum p2,
-						bool colouredfirst);
-
-  /**
-   *  Calculate the initial evolution scales given momenta
-   */
-  pair<Energy,Energy> calculateInitialInitialScales(const Lorentz5Momentum& p1, 
-						    const Lorentz5Momentum& p2);
-
-  /**
-   *  Calculate the initial evolution scales given momenta
-   */
-  pair<Energy,Energy> calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc,
-						  const bool isDecayCase);
-
-protected:
-
-  /**
-   * Given a pair of particles, supposedly partners w.r.t. an interaction,
-   * this method returns their initial evolution scales as a pair.
-   * If something wrong happens, it returns the null (ZERO,ZERO) pair. 
-   * This method is used by the above setXXXInitialEvolutionScales 
-   * methods.
-   */
-  //@{
-  /**
-   *  Calculate the initial evolution scales for two final-state particles
-   */
-  virtual pair<Energy,Energy> calculateFinalFinalScales(const ShowerPPair &);
-
-  /**
-   *  Calculate the initial evolution scales for two initial-state particles
-   */
-  virtual pair<Energy,Energy> calculateInitialInitialScales(const ShowerPPair &);
-
-  /**
-   *  Calculate the initial evolution scales for one initial 
-   *  and one final-state particles
-   */
-  virtual pair<Energy,Energy> calculateInitialFinalScales(const ShowerPPair &,
-							  const bool isDecayCase);
-  //@}
-
-  /**
-   * Access function for the initial conditions for the shower
-   */
-  //@{
-  /**
-   * Initial conditions for the shower of a final-final colour connection
-   * - 0 is the symmetric choice
-   * - 1 is maximal emmision from the coloured particle
-   * - 2 is maximal emmision from the anticoloured particle
-   * - 3 is randomly selected maximal emmision
-   */
-  unsigned int finalFinalConditions() const 
-  {return _finalFinalConditions;}
-
-  /**
-   * Initial conditions for the shower of an initial-final decay colour connection
-   * - 0 is the symmetric choice
-   * - 1 is maximal emission from the decay product
-   * - 2 is the smooth choice
-   */
-  unsigned int initialFinalDecayConditions() const
-  {return _initialFinalDecayConditions;}
-  //@}
-
-protected:
-
-  /** @name Clone Methods. */
-  //@{
-  /**
-   * Make a simple clone of this object.
-   * @return a pointer to the new object.
-   */
-  virtual IBPtr clone() const {return new_ptr(*this);}
-
-  /** Make a clone of this object, possibly modifying the cloned object
-   * to make it sane.
-   * @return a pointer to the new object.
-   */
-  virtual IBPtr fullclone() const {return new_ptr(*this);}
-  //@}
-
-private:
-
-  /**
-   * The assignment operator is private and must never be called.
-   * In fact, it should not even be implemented.
-   */
-  QTildeFinder & operator=(const QTildeFinder &);
-
-private:
-
-  /**
-   *  Flags controlling the initial conditions for the shower
-   */
-  //@{
-  /**
-   * Initial conditions for the shower with a final-final colour
-   * connection
-   */
-  unsigned int _finalFinalConditions; 
-
-  /**
-   * Initial conditions for the shower with an initial-final decay colour
-   * connection. This is done according to the top decay colour 
-   *  connection calculation in JHEP12(2003)_045. The options act as follows:
-   *  0: This is the default 'symmetric' choice which more or less divides
-   *     the phase space evenly between the parent and its charged child.
-   *  1: This 'maximal' choice maximises the phase space available for 
-   *     gluons emitted from the charged child.
-   *  2: This (experimental) 'smooth' choice does not suffer from
-   *     a discontinuity at the boundary between the region populated by
-   *     emissions from the charged child and the region populated by emissions
-   *     from the parent. This does, however, mean that the phase space 
-   *     available for emissions from the charged child is fairly minimal.
-   */
-  unsigned int _initialFinalDecayConditions;
-
-  /**
-   * Initial conditions for the shower with an initial-initial colour
-   * connection. This is done according to the colour connection 
-   * calculation in JHEP12(2003)_045. The options act as follows:
-   *  0: This is the default 'symmetric' choice which more or less divides
-   *     the phase space evenly between the two incoming partons.
-   *  1: This increases the phase space for emission from "parton b".
-   *  2: This increases the phase space for emission from "parton c".
-   */
-  unsigned int _initialInitialConditions;
-  //@}
-};
-
-}
-
-#endif /* HERWIG_QTildeFinder_H */
diff --git a/Shower/QTilde/Makefile.am b/Shower/QTilde/Makefile.am
--- a/Shower/QTilde/Makefile.am
+++ b/Shower/QTilde/Makefile.am
@@ -1,43 +1,42 @@
 SUBDIRS = Matching
 pkglib_LTLIBRARIES = HwShower.la
 HwShower_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 25:0:0
 HwShower_la_SOURCES =  \
 Couplings/ShowerAlphaQCD.h Couplings/ShowerAlphaQCD.cc \
 Couplings/ShowerAlphaQED.h Couplings/ShowerAlphaQED.cc\
 QTildeShowerHandler.h QTildeShowerHandler.fh QTildeShowerHandler.cc \
 SplittingFunctions/HalfHalfOneSplitFn.h SplittingFunctions/HalfHalfOneSplitFn.cc\
 SplittingFunctions/OneOneOneSplitFn.h SplittingFunctions/OneOneOneSplitFn.cc\
 SplittingFunctions/OneOneOneMassiveSplitFn.h SplittingFunctions/OneOneOneMassiveSplitFn.cc\
 SplittingFunctions/ZeroZeroOneSplitFn.h SplittingFunctions/ZeroZeroOneSplitFn.cc\
 SplittingFunctions/OneHalfHalfSplitFn.h SplittingFunctions/OneHalfHalfSplitFn.cc\
 SplittingFunctions/HalfOneHalfSplitFn.h SplittingFunctions/HalfOneHalfSplitFn.cc\
 SplittingFunctions/CMWOneOneOneSplitFn.h SplittingFunctions/CMWOneOneOneSplitFn.cc\
 SplittingFunctions/CMWHalfHalfOneSplitFn.h SplittingFunctions/CMWHalfHalfOneSplitFn.cc\
 Default/QTildeSudakov.cc Default/QTildeSudakov.h\
 Default/Decay_QTildeShowerKinematics1to2.cc \
 Default/Decay_QTildeShowerKinematics1to2.h  \
 Default/IS_QTildeShowerKinematics1to2.cc    Default/IS_QTildeShowerKinematics1to2.h  \
 Default/FS_QTildeShowerKinematics1to2.cc    Default/FS_QTildeShowerKinematics1to2.h  \
-Default/QTildeFinder.cc Default/QTildeFinder.h\
 Default/QTildeReconstructor.cc Default/QTildeReconstructor.h Default/QTildeReconstructor.tcc \
 Base/KinematicsReconstructor.cc \
 Base/KinematicsReconstructor.h \
 Base/KinematicsReconstructor.fh \
 Base/HardTree.cc Base/HardTree.h Base/HardTree.fh \
 Base/HardBranching.h Base/HardBranching.fh Base/HardBranching.cc\
 Base/PartnerFinder.h Base/PartnerFinder.fh Base/PartnerFinder.cc \
 Base/ShowerVeto.h Base/ShowerVeto.fh Base/ShowerVeto.cc \
 Base/FullShowerVeto.h Base/FullShowerVeto.fh Base/FullShowerVeto.cc \
 SplittingFunctions/SplittingGenerator.cc SplittingFunctions/SplittingGenerator.h\
 SplittingFunctions/SplittingGenerator.fh \
 Base/ShowerTree.h Base/ShowerTree.fh Base/ShowerTree.cc \
 ShowerConfig.h ShowerConfig.cc \
 Base/Branching.h \
 Base/ShowerParticle.cc  Base/ShowerParticle.fh  Base/ShowerParticle.h \
 Base/ShowerKinematics.fh  Base/ShowerKinematics.h Base/ShowerKinematics.cc \
 Base/ShowerBasis.fh  Base/ShowerBasis.h Base/ShowerBasis.cc \
 Base/ShowerProgenitor.fh Base/ShowerProgenitor.h \
 Base/SudakovFormFactor.cc Base/SudakovFormFactor.h Base/SudakovFormFactor.fh \
 SplittingFunctions/SplittingFunction.h SplittingFunctions/SplittingFunction.fh \
 SplittingFunctions/SplittingFunction.cc  \
 Base/ShowerVertex.cc Base/ShowerVertex.fh Base/ShowerVertex.h
diff --git a/src/defaults/Shower.in b/src/defaults/Shower.in
--- a/src/defaults/Shower.in
+++ b/src/defaults/Shower.in
@@ -1,310 +1,310 @@
 # -*- ThePEG-repository -*-
 
 ############################################################
 # Setup of default parton shower
 #
 # Useful switches for users are marked near the top of 
 # this file.
 #
 # Don't edit this file directly, but reset the switches
 # in your own input files!
 ############################################################
 library HwMPI.so
 library HwShower.so
 library HwMatching.so
 
 mkdir /Herwig/Shower
 cd /Herwig/Shower
 
 create Herwig::QTildeShowerHandler ShowerHandler
 newdef ShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
 newdef ShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
 # use LO PDFs for Shower, can be changed later
 newdef ShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
 newdef ShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
 newdef ShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
 newdef ShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
 
 #####################################
 # initial setup, don't change these!
 #####################################
 create Herwig::SplittingGenerator SplittingGenerator
 create Herwig::ShowerAlphaQCD AlphaQCD
 create Herwig::ShowerAlphaQED AlphaQED
-create Herwig::QTildeFinder PartnerFinder
+create Herwig::PartnerFinder PartnerFinder
 newdef PartnerFinder:PartnerMethod 1
 newdef PartnerFinder:ScaleChoice 1
 create Herwig::QTildeReconstructor KinematicsReconstructor
 newdef KinematicsReconstructor:ReconstructionOption Colour3
 newdef KinematicsReconstructor:InitialStateReconOption SofterFraction
 newdef KinematicsReconstructor:InitialInitialBoostOption LongTransBoost
 
 newdef /Herwig/Partons/RemnantDecayer:AlphaS  AlphaQCD
 newdef /Herwig/Partons/RemnantDecayer:AlphaEM AlphaQED
 
 newdef ShowerHandler:PartnerFinder PartnerFinder
 newdef ShowerHandler:KinematicsReconstructor KinematicsReconstructor
 newdef ShowerHandler:SplittingGenerator SplittingGenerator
 newdef ShowerHandler:SpinCorrelations Yes
 newdef ShowerHandler:SoftCorrelations Singular
 
 ##################################################################
 # Intrinsic pT
 #
 # Recommended: 
 # 1.9 GeV for Tevatron W/Z production.
 # 2.1 GeV for LHC W/Z production at 10 TeV
 # 2.2 GeV for LHC W/Z production at 14 TeV
 #
 # Set all parameters to 0 to disable
 ##################################################################
 newdef ShowerHandler:IntrinsicPtGaussian 1.3*GeV
 newdef ShowerHandler:IntrinsicPtBeta  	0
 newdef ShowerHandler:IntrinsicPtGamma 	0*GeV
 newdef ShowerHandler:IntrinsicPtIptmax   0*GeV
 #############################################################
 # Set up truncated shower handler.
 #############################################################
 
 create Herwig::PowhegShowerHandler PowhegShowerHandler
 set PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
 set PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
 newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
 newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
 newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
 newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
 newdef PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
 newdef PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
 newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
 newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
 newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
 newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
 newdef PowhegShowerHandler:PartnerFinder PartnerFinder
 newdef PowhegShowerHandler:KinematicsReconstructor KinematicsReconstructor
 newdef PowhegShowerHandler:SplittingGenerator SplittingGenerator
 newdef PowhegShowerHandler:Interactions QCDandQED
 newdef PowhegShowerHandler:SpinCorrelations Yes
 newdef PowhegShowerHandler:SoftCorrelations Singular
 newdef PowhegShowerHandler:IntrinsicPtGaussian 1.3*GeV
 newdef PowhegShowerHandler:IntrinsicPtBeta  	0
 newdef PowhegShowerHandler:IntrinsicPtGamma 	0*GeV
 newdef PowhegShowerHandler:IntrinsicPtIptmax   0*GeV
 newdef PowhegShowerHandler:ReconstructionOption OffShell5
 
 #############################################################
 # End of interesting user servicable section.
 #
 # Anything that follows below should only be touched if you 
 # know what you're doing. 
 #
 # Really.
 #############################################################
 #
 # a few default values
 newdef ShowerHandler:MECorrMode 1
 newdef ShowerHandler:ReconstructionOption OffShell5
 newdef AlphaQCD:ScaleFactor 1.0
 newdef AlphaQCD:NPAlphaS 2
 newdef AlphaQCD:Qmin 0.935
 newdef AlphaQCD:NumberOfLoops 2
 newdef AlphaQCD:InputOption 1
 newdef AlphaQCD:AlphaMZ 0.126234
 #
 #
 # Lets set up all the splittings
 create Herwig::HalfHalfOneSplitFn QtoQGammaSplitFn
 set QtoQGammaSplitFn:InteractionType QED
 set QtoQGammaSplitFn:ColourStructure ChargedChargedNeutral
 set QtoQGammaSplitFn:AngularOrdered Yes
 
 create Herwig::HalfHalfOneSplitFn QtoQGSplitFn
 newdef QtoQGSplitFn:InteractionType QCD
 newdef QtoQGSplitFn:ColourStructure TripletTripletOctet
 set QtoQGSplitFn:AngularOrdered Yes
 
 create Herwig::OneOneOneSplitFn GtoGGSplitFn
 newdef GtoGGSplitFn:InteractionType QCD
 newdef GtoGGSplitFn:ColourStructure OctetOctetOctet
 set GtoGGSplitFn:AngularOrdered Yes
 
 create Herwig::OneOneOneMassiveSplitFn WtoWGammaSplitFn
 newdef WtoWGammaSplitFn:InteractionType QED
 newdef WtoWGammaSplitFn:ColourStructure ChargedChargedNeutral
 set WtoWGammaSplitFn:AngularOrdered Yes
 
 create Herwig::OneHalfHalfSplitFn GtoQQbarSplitFn
 newdef GtoQQbarSplitFn:InteractionType QCD
 newdef GtoQQbarSplitFn:ColourStructure OctetTripletTriplet
 set GtoQQbarSplitFn:AngularOrdered Yes
 
 create Herwig::OneHalfHalfSplitFn GammatoQQbarSplitFn
 newdef GammatoQQbarSplitFn:InteractionType QED
 newdef GammatoQQbarSplitFn:ColourStructure NeutralChargedCharged
 set GammatoQQbarSplitFn:AngularOrdered Yes
 
 create Herwig::HalfOneHalfSplitFn QtoGQSplitFn
 newdef QtoGQSplitFn:InteractionType QCD
 newdef QtoGQSplitFn:ColourStructure TripletOctetTriplet
 set QtoGQSplitFn:AngularOrdered Yes
 
 create Herwig::HalfOneHalfSplitFn QtoGammaQSplitFn
 newdef QtoGammaQSplitFn:InteractionType QED
 newdef QtoGammaQSplitFn:ColourStructure ChargedNeutralCharged
 set QtoGammaQSplitFn:AngularOrdered Yes
 #
 # Now the Sudakovs
 create Herwig::QTildeSudakov SudakovCommon
 newdef SudakovCommon:Alpha AlphaQCD
 newdef SudakovCommon:cutoffKinScale 0.0*GeV
 newdef SudakovCommon:PDFmax 1.0
 newdef SudakovCommon:CutOffOption pT
 newdef SudakovCommon:pTmin 1.222798*GeV
 
 cp SudakovCommon QtoQGSudakov
 newdef QtoQGSudakov:SplittingFunction QtoQGSplitFn
 newdef QtoQGSudakov:PDFmax 1.9
 
 cp SudakovCommon QtoQGammaSudakov
 set QtoQGammaSudakov:SplittingFunction QtoQGammaSplitFn
 set QtoQGammaSudakov:Alpha AlphaQED
 set QtoQGammaSudakov:PDFmax 1.9
 
 cp QtoQGammaSudakov LtoLGammaSudakov
 # Technical parameter to stop evolution.
 set LtoLGammaSudakov:pTmin 0.000001
 
 cp SudakovCommon GtoGGSudakov
 newdef GtoGGSudakov:SplittingFunction GtoGGSplitFn
 newdef GtoGGSudakov:PDFmax 2.0
 
 cp SudakovCommon WtoWGammaSudakov
 newdef WtoWGammaSudakov:SplittingFunction WtoWGammaSplitFn
 set WtoWGammaSudakov:Alpha AlphaQED
 
 cp SudakovCommon GtoQQbarSudakov
 newdef GtoQQbarSudakov:SplittingFunction GtoQQbarSplitFn
 newdef GtoQQbarSudakov:PDFmax 120.0
 
 cp SudakovCommon GammatoQQbarSudakov
 newdef GammatoQQbarSudakov:SplittingFunction GammatoQQbarSplitFn
 set GammatoQQbarSudakov:Alpha AlphaQED
 newdef GammatoQQbarSudakov:PDFmax 120.0
 
 cp SudakovCommon GtobbbarSudakov
 newdef GtobbbarSudakov:SplittingFunction GtoQQbarSplitFn
 newdef GtobbbarSudakov:PDFmax 40000.0
 
 cp SudakovCommon GtoccbarSudakov
 newdef GtoccbarSudakov:SplittingFunction GtoQQbarSplitFn
 newdef GtoccbarSudakov:PDFmax 2000.0
 
 cp SudakovCommon QtoGQSudakov
 newdef QtoGQSudakov:SplittingFunction QtoGQSplitFn
 
 cp SudakovCommon QtoGammaQSudakov
 newdef QtoGammaQSudakov:SplittingFunction QtoGammaQSplitFn
 set QtoGammaQSudakov:Alpha AlphaQED
 
 cp SudakovCommon utoGuSudakov
 newdef utoGuSudakov:SplittingFunction QtoGQSplitFn
 newdef utoGuSudakov:PDFFactor OverOneMinusZ
 newdef utoGuSudakov:PDFmax 5.0
 
 cp SudakovCommon dtoGdSudakov
 newdef dtoGdSudakov:SplittingFunction QtoGQSplitFn
 newdef dtoGdSudakov:PDFFactor OverOneMinusZ
 
 
 #
 # Now add the final splittings
 #
 do SplittingGenerator:AddFinalSplitting u->u,g; QtoQGSudakov
 do SplittingGenerator:AddFinalSplitting d->d,g; QtoQGSudakov
 do SplittingGenerator:AddFinalSplitting s->s,g; QtoQGSudakov
 do SplittingGenerator:AddFinalSplitting c->c,g; QtoQGSudakov
 do SplittingGenerator:AddFinalSplitting b->b,g; QtoQGSudakov
 do SplittingGenerator:AddFinalSplitting t->t,g; QtoQGSudakov
 #
 do SplittingGenerator:AddFinalSplitting g->g,g; GtoGGSudakov
 #
 do SplittingGenerator:AddFinalSplitting g->u,ubar; GtoQQbarSudakov
 do SplittingGenerator:AddFinalSplitting g->d,dbar; GtoQQbarSudakov
 do SplittingGenerator:AddFinalSplitting g->s,sbar; GtoQQbarSudakov
 do SplittingGenerator:AddFinalSplitting g->c,cbar; GtoccbarSudakov
 do SplittingGenerator:AddFinalSplitting g->b,bbar; GtobbbarSudakov
 do SplittingGenerator:AddFinalSplitting g->t,tbar; GtoQQbarSudakov
 #
 do SplittingGenerator:AddFinalSplitting gamma->u,ubar; GammatoQQbarSudakov
 do SplittingGenerator:AddFinalSplitting gamma->d,dbar; GammatoQQbarSudakov
 do SplittingGenerator:AddFinalSplitting gamma->s,sbar; GammatoQQbarSudakov
 do SplittingGenerator:AddFinalSplitting gamma->c,cbar; GammatoQQbarSudakov
 do SplittingGenerator:AddFinalSplitting gamma->b,bbar; GammatoQQbarSudakov
 do SplittingGenerator:AddFinalSplitting gamma->t,tbar; GammatoQQbarSudakov
 do SplittingGenerator:AddFinalSplitting gamma->e-,e+; GammatoQQbarSudakov
 do SplittingGenerator:AddFinalSplitting gamma->mu-,mu+; GammatoQQbarSudakov
 do SplittingGenerator:AddFinalSplitting gamma->tau-,tau+; GammatoQQbarSudakov
 #
 do SplittingGenerator:AddFinalSplitting u->u,gamma; QtoQGammaSudakov
 do SplittingGenerator:AddFinalSplitting d->d,gamma; QtoQGammaSudakov
 do SplittingGenerator:AddFinalSplitting s->s,gamma; QtoQGammaSudakov
 do SplittingGenerator:AddFinalSplitting c->c,gamma; QtoQGammaSudakov
 do SplittingGenerator:AddFinalSplitting b->b,gamma; QtoQGammaSudakov
 do SplittingGenerator:AddFinalSplitting t->t,gamma; QtoQGammaSudakov
 
 do SplittingGenerator:AddFinalSplitting e-->e-,gamma; LtoLGammaSudakov
 do SplittingGenerator:AddFinalSplitting mu-->mu-,gamma; LtoLGammaSudakov
 do SplittingGenerator:AddFinalSplitting tau-->tau-,gamma; LtoLGammaSudakov
 
 do SplittingGenerator:AddFinalSplitting W+->W+,gamma; WtoWGammaSudakov
 #
 # Now lets add the initial splittings. Remember the form a->b,c; means
 # that the current particle b is given and we backward branch to new
 # particle a which is initial state and new particle c which is final state
 #
 do SplittingGenerator:AddInitialSplitting u->u,g; QtoQGSudakov
 do SplittingGenerator:AddInitialSplitting d->d,g; QtoQGSudakov
 do SplittingGenerator:AddInitialSplitting s->s,g; QtoQGSudakov
 do SplittingGenerator:AddInitialSplitting c->c,g; QtoQGSudakov
 do SplittingGenerator:AddInitialSplitting b->b,g; QtoQGSudakov
 do SplittingGenerator:AddInitialSplitting u->u,gamma; QtoQGammaSudakov
 do SplittingGenerator:AddInitialSplitting d->d,gamma; QtoQGammaSudakov
 do SplittingGenerator:AddInitialSplitting s->s,gamma; QtoQGammaSudakov
 do SplittingGenerator:AddInitialSplitting c->c,gamma; QtoQGammaSudakov
 do SplittingGenerator:AddInitialSplitting b->b,gamma; QtoQGammaSudakov
 do SplittingGenerator:AddInitialSplitting t->t,gamma; QtoQGammaSudakov
 
 do SplittingGenerator:AddInitialSplitting g->g,g; GtoGGSudakov
 #
 do SplittingGenerator:AddInitialSplitting g->d,dbar; GtoQQbarSudakov
 do SplittingGenerator:AddInitialSplitting g->u,ubar; GtoQQbarSudakov
 do SplittingGenerator:AddInitialSplitting g->s,sbar; GtoQQbarSudakov
 do SplittingGenerator:AddInitialSplitting g->c,cbar; GtoccbarSudakov
 do SplittingGenerator:AddInitialSplitting g->b,bbar; GtobbbarSudakov
 #
 do SplittingGenerator:AddInitialSplitting gamma->d,dbar; GammatoQQbarSudakov
 do SplittingGenerator:AddInitialSplitting gamma->u,ubar; GammatoQQbarSudakov
 do SplittingGenerator:AddInitialSplitting gamma->s,sbar; GammatoQQbarSudakov
 do SplittingGenerator:AddInitialSplitting gamma->c,cbar; GammatoQQbarSudakov
 do SplittingGenerator:AddInitialSplitting gamma->b,bbar; GammatoQQbarSudakov
 #
 do SplittingGenerator:AddInitialSplitting d->g,d; dtoGdSudakov
 do SplittingGenerator:AddInitialSplitting u->g,u; utoGuSudakov
 do SplittingGenerator:AddInitialSplitting s->g,s; QtoGQSudakov
 do SplittingGenerator:AddInitialSplitting c->g,c; QtoGQSudakov
 do SplittingGenerator:AddInitialSplitting b->g,b; QtoGQSudakov
 do SplittingGenerator:AddInitialSplitting dbar->g,dbar; dtoGdSudakov
 do SplittingGenerator:AddInitialSplitting ubar->g,ubar; utoGuSudakov
 do SplittingGenerator:AddInitialSplitting sbar->g,sbar; QtoGQSudakov
 do SplittingGenerator:AddInitialSplitting cbar->g,cbar; QtoGQSudakov
 do SplittingGenerator:AddInitialSplitting bbar->g,bbar; QtoGQSudakov
 #
 do SplittingGenerator:AddInitialSplitting d->gamma,d; QtoGammaQSudakov
 do SplittingGenerator:AddInitialSplitting u->gamma,u; QtoGammaQSudakov
 do SplittingGenerator:AddInitialSplitting s->gamma,s; QtoGammaQSudakov
 do SplittingGenerator:AddInitialSplitting c->gamma,c; QtoGammaQSudakov
 do SplittingGenerator:AddInitialSplitting b->gamma,b; QtoGammaQSudakov
 do SplittingGenerator:AddInitialSplitting dbar->gamma,dbar; QtoGammaQSudakov
 do SplittingGenerator:AddInitialSplitting ubar->gamma,ubar; QtoGammaQSudakov
 do SplittingGenerator:AddInitialSplitting sbar->gamma,sbar; QtoGammaQSudakov
 do SplittingGenerator:AddInitialSplitting cbar->gamma,cbar; QtoGammaQSudakov
 do SplittingGenerator:AddInitialSplitting bbar->gamma,bbar; QtoGammaQSudakov