diff --git a/MatrixElement/Matchbox/Phasespace/FFLightInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FFLightInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FFLightInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FFLightInvertedTildeKinematics.cc
@@ -1,136 +1,137 @@
 // -*- C++ -*-
 //
 // FFLightInvertedTildeKinematics.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2012 The Herwig Collaboration
 //
 // Herwig is licenced under version 2 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 FFLightInvertedTildeKinematics class.
 //
 
 #include "FFLightInvertedTildeKinematics.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Repository/EventGenerator.h"
 
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 FFLightInvertedTildeKinematics::FFLightInvertedTildeKinematics() {}
 
 FFLightInvertedTildeKinematics::~FFLightInvertedTildeKinematics() {}
 
 IBPtr FFLightInvertedTildeKinematics::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr FFLightInvertedTildeKinematics::fullclone() const {
   return new_ptr(*this);
 }
 
 bool FFLightInvertedTildeKinematics::doMap(const double * r) {
 
   if ( ptMax() < ptCut() ) {
     jacobian(0.0);
     return false;
   }
 
   Lorentz5Momentum emitter = bornEmitterMomentum();
   Lorentz5Momentum spectator = bornSpectatorMomentum();
 
   double mapping = 1.0;
   pair<Energy,double> ptz = generatePtZ(mapping,r);
   if ( mapping == 0.0 ) {
     jacobian(0.0);
     return false;
   }
 
   Energy pt = ptz.first;
 
   double z = ptz.second;
   double y = sqr(pt/lastScale())/(z*(1.-z));
   if ( y < 0. || y > 1. ||
        z < 0. || z > 1. ) {
     jacobian(0.0);
     return false;
   }
 
   mapping *= (1.-y)/(z*(1.-z));
   jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
 
   double phi = 2.*Constants::pi*r[2];
   Lorentz5Momentum kt
     = getKt(emitter,spectator,pt,phi);
 
   subtractionParameters().resize(2);
   subtractionParameters()[0] = y;
   subtractionParameters()[1] = z;
 
   realEmitterMomentum() = z*emitter + y*(1.-z)*spectator + kt;
   realEmissionMomentum() = (1.-z)*emitter + y*z*spectator - kt;
   realSpectatorMomentum() = (1.-y)*spectator;
 
   realEmitterMomentum().setMass(ZERO);
   realEmitterMomentum().rescaleEnergy();
   realEmissionMomentum().setMass(ZERO);
   realEmissionMomentum().rescaleEnergy();
   realSpectatorMomentum().setMass(ZERO);
   realSpectatorMomentum().rescaleEnergy();
 
   return true;
 
 }
 
 Energy FFLightInvertedTildeKinematics::lastPt() const {
 
   Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
   double y = subtractionParameters()[0];
   double z = subtractionParameters()[1];
   return scale * sqrt(y*z*(1.-z));
 
 }
 
 double FFLightInvertedTildeKinematics::lastZ() const {
   return subtractionParameters()[1];
 }
 
 Energy FFLightInvertedTildeKinematics::ptMax() const {
   return lastScale()/2.;
 }
 
 pair<double,double> FFLightInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
   hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
+  if(pt>hardPt) return make_pair(0.5,0.5);
   double s = sqrt(1.-sqr(pt/hardPt));
   return make_pair(0.5*(1.-s),0.5*(1.+s));
 }
 
 // If needed, insert default implementations of virtual function defined
 // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
 
 
 void FFLightInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
 }
 
 void FFLightInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
 }
 
 void FFLightInvertedTildeKinematics::Init() {
 
   static ClassDocumentation<FFLightInvertedTildeKinematics> documentation
     ("FFLightInvertedTildeKinematics inverts the final-final tilde "
      "kinematics.");
 
 }
 
 // *** 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<FFLightInvertedTildeKinematics,InvertedTildeKinematics>
 describeHerwigFFLightInvertedTildeKinematics("Herwig::FFLightInvertedTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.cc
@@ -1,303 +1,304 @@
 // -*- C++ -*-
 //
 // FFMassiveInvertedTildeKinematics.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2012 The Herwig Collaboration
 //
 // Herwig is licenced under version 2 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 FFMassiveInvertedTildeKinematics class.
 //
 
 #include "FFMassiveInvertedTildeKinematics.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Repository/EventGenerator.h"
 
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 #include "Herwig/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
 
 using namespace Herwig;
 
 FFMassiveInvertedTildeKinematics::FFMassiveInvertedTildeKinematics() {}
 
 FFMassiveInvertedTildeKinematics::~FFMassiveInvertedTildeKinematics() {}
 
 IBPtr FFMassiveInvertedTildeKinematics::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr FFMassiveInvertedTildeKinematics::fullclone() const {
   return new_ptr(*this);
 }
 
 bool FFMassiveInvertedTildeKinematics::doMap(const double * r) {
 
   if ( ptMax() < ptCut() ) {
     jacobian(0.0);
     return false;
   }
 
   Lorentz5Momentum emitter = bornEmitterMomentum();
   Lorentz5Momentum spectator = bornSpectatorMomentum();
 
   double mapping = 1.0;
   pair<Energy,double> ptz = generatePtZ(mapping,r);
   if ( mapping == 0.0 ) {
     jacobian(0.0);
     return false;
   }
 
   Energy pt = ptz.first;
   double z = ptz.second;
   
   Energy scale = (emitter+spectator).m();
   // masses
   double mui2 = sqr( realEmitterData()->hardProcessMass() / scale );
   double mu2  = sqr( realEmissionData()->hardProcessMass() / scale );
   double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale );
   double Mui2 = sqr( bornEmitterData()->hardProcessMass() / scale );
   double Muj2 = sqr( bornSpectatorData()->hardProcessMass() / scale );
   
   double y = ( sqr( pt / scale ) + sqr(1.-z)*mui2 + z*z*mu2 ) /
       (z*(1.-z)*(1.-mui2-mu2-muj2));
   
   // now this is done in ptzAllowed!!
   // check (y,z) phasespace boundary
   // 2011-11-09: never occurred
   /**
   double bar = 1.-mui2-mu2-muj2;
   double ym = 2.*sqrt(mui2)*sqrt(mu2)/bar;
   double yp = 1. - 2.*sqrt(muj2)*(1.-sqrt(muj2))/bar;
   double zm = ( (2.*mui2+bar*y)*(1.-y) - sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
     ( 2.*(1.-y)*(mui2+mu2+bar*y) );
   double zp = ( (2.*mui2+bar*y)*(1.-y) + sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
     ( 2.*(1.-y)*(mui2+mu2+bar*y) );
   if ( y<ym || y>yp || z<zm || z>zp ) {
     cout << "A problem occurred in FFMassiveInvertedTildeKinematics::doMap: ";
     if(y<ym) cout << "y<ym " << endl;
     if(y>yp) cout << "y>yp " << endl;
     if(z<zm) cout << "z<zm " << endl;
     if(z>zp) cout << "z>zp " << endl;
     jacobian(0.0);
     return false;
   }
   **/
       
   mapping /= z*(1.-z);
   jacobian( mapping*(1.-y)*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)) *
     sqr(1.-mui2-mu2-muj2) / rootOfKallen(1.,Mui2,Muj2) );
 
   double phi = 2.*Constants::pi*r[2];
   /* // not used ???
   Lorentz5Momentum kt
     = getKt(emitter,spectator,pt,phi);
   */
 
   subtractionParameters().resize(2);
   subtractionParameters()[0] = y;
   subtractionParameters()[1] = z;
   
   // kinematics from FFMassiveKinematics.cc
   // updated 2011-08-23
   // updated 2011-11-08
   Energy2 sbar = sqr(scale) * (1.-mui2-mu2-muj2);
   // CMF: particle energies
   Energy Ei = ( sbar*(1.-(1.-z)*(1.-y)) + 2.*sqr(scale)*mui2 ) / (2.*scale);
   Energy E  = ( sbar*(1.-    z *(1.-y)) + 2.*sqr(scale)*mu2  ) / (2.*scale);
   Energy Ej = ( sbar*(1.-           y ) + 2.*sqr(scale)*muj2 ) / (2.*scale);
   // CMF: momenta in z-direction (axis of pEmitter &pEmitter pSpectator)  
   Energy qi3 = (2.*Ei*Ej-z*(1.-y)*sbar     ) / 2./sqrt(Ej*Ej-sqr(scale)*muj2);
   Energy q3  = (2.*E *Ej-(1.-z)*(1.-y)*sbar) / 2./sqrt(Ej*Ej-sqr(scale)*muj2);
   Energy qj3 = sqrt(   sqr(Ej) - sqr(scale)*muj2 );
   // get z axis in the dipole's CMF which is parallel to pSpectator
   Boost toCMF = (emitter+spectator).findBoostToCM();
   Lorentz5Momentum pjAux = spectator; pjAux.boost(toCMF);
   ThreeVector<double> pjAxis = pjAux.vect().unit();
   // set the momenta in this special reference frame
   // note that pt might in some cases differ from the physical pt!
   Energy ptResc = sqrt( sqr(Ei)-sqr(scale)*mui2-sqr(qi3) );
   Lorentz5Momentum em  ( ptResc*cos(phi), -ptResc*sin(phi), qi3, Ei );
   Lorentz5Momentum emm ( -ptResc*cos(phi), ptResc*sin(phi), q3, E );
   Lorentz5Momentum spe ( 0.*GeV, 0.*GeV, qj3, Ej );
   // rotate back
   em.rotateUz (pjAxis);
   emm.rotateUz(pjAxis);
   spe.rotateUz(pjAxis);
   // boost back
   em.boost (-toCMF);
   emm.boost(-toCMF);
   spe.boost(-toCMF);
   // mass shells, rescale energy
   em.setMass(scale*sqrt(mui2));
   em.rescaleEnergy();
   emm.setMass(scale*sqrt(mu2));
   emm.rescaleEnergy();
   spe.setMass(scale*sqrt(muj2));
   spe.rescaleEnergy();
 
   // book
   realEmitterMomentum() = em;
   realEmissionMomentum() = emm;
   realSpectatorMomentum() = spe;
   
   return true;
 
 }
 
 Energy FFMassiveInvertedTildeKinematics::lastPt() const {
   
   Energy scale = (bornEmitterMomentum()+bornSpectatorMomentum()).m();
   // masses
   double mui2 = sqr( realEmitterData()->hardProcessMass() / scale );
   double mu2  = sqr( realEmissionData()->hardProcessMass() / scale );
   double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale );
   
   double y = subtractionParameters()[0];
   double z = subtractionParameters()[1];
   
   Energy ret = scale * sqrt( y * (1.-mui2-mu2-muj2) * z*(1.-z) - sqr(1.-z)*mui2 - sqr(z)*mu2 );
 
   return ret;
 }
 
 double FFMassiveInvertedTildeKinematics::lastZ() const {
   return subtractionParameters()[1];
 }
     
 Energy FFMassiveInvertedTildeKinematics::ptMax() const {
   
   Energy scale = (bornEmitterMomentum()+bornSpectatorMomentum()).m();
   // masses
   double mui2 = sqr( realEmitterData()->hardProcessMass() / scale );
   double mu2  = sqr( realEmissionData()->hardProcessMass() / scale );
   double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale );
   
   Energy ptmax = rootOfKallen( mui2, mu2, sqr(1.-sqrt(muj2)) ) /
     ( 2.-2.*sqrt(muj2) ) * scale;
   
   return ptmax > 0.*GeV ? ptmax : 0.*GeV;
 }
 
 // NOTE: bounds calculated at this step may be too loose
 pair<double,double> FFMassiveInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
 
   hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
+  if(pt>hardPt) return make_pair(0.5,0.5);
   
   Energy scale = (bornEmitterMomentum()+bornSpectatorMomentum()).m();
   // masses
   double mui2 = sqr( realEmitterData()->hardProcessMass() / scale );
   double mu2  = sqr( realEmissionData()->hardProcessMass() / scale );
   double muj2 = sqr( realSpectatorData()->hardProcessMass() / scale );
   
   double zp = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) +
     rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
     sqrt( 1.-sqr(pt/hardPt) ) ) /
     ( 2.*sqr(1.-sqrt(muj2)) );
   double zm = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) -
     rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
     sqrt( 1.-sqr(pt/hardPt) ) ) /
     ( 2.*sqr(1.-sqrt(muj2)) );
     
   return make_pair(zm,zp);
 }
 
 bool FFMassiveInvertedTildeKinematics::ptzAllowed(pair<Energy,double> ptz) const {
   // masses
   double mui2 = sqr( realEmitterData()->hardProcessMass() / lastScale() );
   double mu2  = sqr( realEmissionData()->hardProcessMass() / lastScale() );
   double muj2 = sqr( realSpectatorData()->hardProcessMass() / lastScale() );
   
   Energy pt = ptz.first;
   double z = ptz.second;
   
   double y = ( sqr( pt / lastScale() ) + sqr(1.-z)*mui2 + z*z*mu2 ) /
       (z*(1.-z)*(1.-mui2-mu2-muj2));
   
   // check (y,z) phasespace boundary
   // TODO: is y boundary necessary?
   double bar = 1.-mui2-mu2-muj2;
   double ym = 2.*sqrt(mui2)*sqrt(mu2)/bar;
   double yp = 1. - 2.*sqrt(muj2)*(1.-sqrt(muj2))/bar;
   double zm = ( (2.*mui2+bar*y)*(1.-y) - sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
     ( 2.*(1.-y)*(mui2+mu2+bar*y) );
   double zp = ( (2.*mui2+bar*y)*(1.-y) + sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
     ( 2.*(1.-y)*(mui2+mu2+bar*y) );
   
   if ( y<ym || y>yp || z<zm || z>zp ) return false;
   return true;
 }
 
 pair<Energy,double> FFMassiveInvertedTildeKinematics::generatePtZ(double& jac, const double * r) const {
 
   double kappaMin = 
     ptCut() != ZERO ?
     sqr(ptCut()/ptMax()) :
     sqr(0.1*GeV/GeV);
 
   double kappa;
 
   using namespace RandomHelpers;
 
   if ( ptCut() > ZERO ) {
     pair<double,double> kw =
       generate(inverse(0.,kappaMin,1.),r[0]);
     kappa = kw.first;
     jac *= kw.second;
   } else {
     pair<double,double> kw =
       generate((piecewise(),
 		flat(1e-4,kappaMin),
 		match(inverse(0.,kappaMin,1.))),r[0]);
     kappa = kw.first;
     jac *= kw.second;
   }
 
   Energy pt = sqrt(kappa)*ptMax();
 
   pair<double,double> zLims = zBounds(pt);
 
   pair<double,double> zw =
     generate(inverse(0.,zLims.first,zLims.second)+
 	     inverse(1.,zLims.first,zLims.second),r[1]);
 
   double z = zw.first;
   jac *= zw.second;
 
   jac *= sqr(ptMax()/lastScale());
   
   if( !ptzAllowed(make_pair(pt,z)) ) jac = 0.;
 
   return make_pair(pt,z);
 
 }
 
 // If needed, insert default implementations of virtual function defined
 // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
 
 
 void FFMassiveInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
 }
 
 void FFMassiveInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
 }
 
 void FFMassiveInvertedTildeKinematics::Init() {
 
   static ClassDocumentation<FFMassiveInvertedTildeKinematics> documentation
     ("FFMassiveInvertedTildeKinematics inverts the final-final tilde "
      "kinematics.");
 
 }
 
 // *** 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<FFMassiveInvertedTildeKinematics,InvertedTildeKinematics>
 describeHerwigFFMassiveInvertedTildeKinematics("Herwig::FFMassiveInvertedTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc
@@ -1,138 +1,134 @@
 // -*- C++ -*-
 //
 // FILightInvertedTildeKinematics.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2012 The Herwig Collaboration
 //
 // Herwig is licenced under version 2 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 FILightInvertedTildeKinematics class.
 //
 
 #include "FILightInvertedTildeKinematics.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 FILightInvertedTildeKinematics::FILightInvertedTildeKinematics() {}
 
 FILightInvertedTildeKinematics::~FILightInvertedTildeKinematics() {}
 
 IBPtr FILightInvertedTildeKinematics::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr FILightInvertedTildeKinematics::fullclone() const {
   return new_ptr(*this);
 }
 
 bool FILightInvertedTildeKinematics::doMap(const double * r) {
 
   if ( ptMax() < ptCut() ) {
     jacobian(0.0);
     return false;
   }
 
   Lorentz5Momentum emitter = bornEmitterMomentum();
   Lorentz5Momentum spectator = bornSpectatorMomentum();
 
   double mapping = 1.0;
   pair<Energy,double> ptz = generatePtZ(mapping,r);
   if ( mapping == 0.0 ) {
     jacobian(0.0);
     return false;
   }
 
   Energy pt = ptz.first;
 
   double z = ptz.second;
   double y = sqr(pt/lastScale())/(z*(1.-z));
   double x = 1./(1.+y);
 
   if ( x < spectatorX() || x > 1. ||
        z < 0. || z > 1. ) {
     jacobian(0.0);
     return false;
   }
 
   mapping /= z*(1.-z);
   jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
 
   double phi = 2.*Constants::pi*r[2];
   Lorentz5Momentum kt
     = getKt(spectator,emitter,pt,phi,true);
 
   subtractionParameters().resize(2);
   subtractionParameters()[0] = x;
   subtractionParameters()[1] = z;
 
   realEmitterMomentum() = z*emitter + (1.-z)*((1.-x)/x)*spectator + kt;
   realEmissionMomentum() = (1.-z)*emitter + z*((1.-x)/x)*spectator - kt;
   realSpectatorMomentum() = (1./x)*spectator;
 
   realEmitterMomentum().setMass(ZERO);
   realEmitterMomentum().rescaleEnergy();
   realEmissionMomentum().setMass(ZERO);
   realEmissionMomentum().rescaleEnergy();
   realSpectatorMomentum().setMass(ZERO);
   realSpectatorMomentum().rescaleEnergy();
 
   return true;
 
 }
 
 Energy FILightInvertedTildeKinematics::lastPt() const {
 
   Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
   double x = subtractionParameters()[0];
   double z = subtractionParameters()[1];
   return scale * sqrt(z*(1.-z)*(1.-x)/x);
 
 }
 
 double FILightInvertedTildeKinematics::lastZ() const {
   return subtractionParameters()[1];
 }
 
 Energy FILightInvertedTildeKinematics::ptMax() const {
   double x = spectatorX();
   return sqrt((1.-x)/x)*lastScale()/2.;
 }
 
 pair<double,double> FILightInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
   hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
+  if(pt>hardPt) return make_pair(0.5,0.5);
   double s = sqrt(1.-sqr(pt/hardPt));
   return make_pair(0.5*(1.-s),0.5*(1.+s));
 }
 
-
-// If needed, insert default implementations of virtual function defined
-// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
-
-
 void FILightInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
 }
 
 void FILightInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
 }
 
 void FILightInvertedTildeKinematics::Init() {
 
   static ClassDocumentation<FILightInvertedTildeKinematics> documentation
     ("FILightInvertedTildeKinematics inverts the final-initial tilde "
      "kinematics.");
 
 }
 
 // *** 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<FILightInvertedTildeKinematics,InvertedTildeKinematics>
 describeHerwigFILightInvertedTildeKinematics("Herwig::FILightInvertedTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/FIMassiveInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FIMassiveInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FIMassiveInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FIMassiveInvertedTildeKinematics.cc
@@ -1,178 +1,179 @@
 // -*- C++ -*-
 //
 // FIMassiveInvertedTildeKinematics.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2012 The Herwig Collaboration
 //
 // Herwig is licenced under version 2 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 FIMassiveInvertedTildeKinematics class.
 //
 
 #include "FIMassiveInvertedTildeKinematics.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 FIMassiveInvertedTildeKinematics::FIMassiveInvertedTildeKinematics() {}
 
 FIMassiveInvertedTildeKinematics::~FIMassiveInvertedTildeKinematics() {}
 
 IBPtr FIMassiveInvertedTildeKinematics::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr FIMassiveInvertedTildeKinematics::fullclone() const {
   return new_ptr(*this);
 }
 
 bool FIMassiveInvertedTildeKinematics::doMap(const double * r) { 
   if ( ptMax() < ptCut() ) {
     jacobian(0.0);
     return false;
   }
 
   Lorentz5Momentum emitter = bornEmitterMomentum();
   Lorentz5Momentum spectator = bornSpectatorMomentum();
 
   double mapping = 1.0;
   pair<Energy,double> ptz = generatePtZ(mapping,r);
   if ( mapping == 0.0 ) {
     jacobian(0.0);
     return false;
   }
 
   Energy pt = ptz.first;
   double z = ptz.second;
 
   Energy2 mi2 = sqr(realEmitterData()->hardProcessMass());
   Energy2 m2  = sqr(realEmissionData()->hardProcessMass());
   Energy2 Mi2 = sqr(bornEmitterData()->hardProcessMass());
 
   Energy2 scale=2.*emitter*spectator;
   double y = (pt*pt+(1.-z)*mi2+z*m2-z*(1.-z)*Mi2) / (z*(1.-z)*scale);
   double x = 1./(1.+y);
 
   if ( x < spectatorX() ) {
     jacobian(0.0);
     return false;
   }
 
   // no additional massive factors
   mapping /= z*(1.-z);
   jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
 
   double phi = 2.*Constants::pi*r[2];
   Lorentz5Momentum kt = getKt(spectator,emitter,pt,phi,true);
 
   subtractionParameters().resize(2);
   subtractionParameters()[0] = x;
   subtractionParameters()[1] = z;
 
   realEmitterMomentum() = z*emitter +
     (sqr(pt)+mi2-z*z*Mi2)/(z*scale)*spectator + kt;
   realEmissionMomentum() = (1.-z)*emitter +
     (pt*pt+m2-sqr(1.-z)*Mi2)/((1.-z)*scale)*spectator - kt;
   realSpectatorMomentum() = (1.+y)*spectator;
 
   double mui2 = x*mi2/scale;
   double mu2  = x*m2/scale;
   double Mui2 = x*Mi2/scale;
   double xp = 1. + Mui2 - sqr(sqrt(mui2)+sqrt(mu2));
   double zm = .5*( 1.-x+Mui2+mui2-mui2 -
   		   sqrt( sqr(1.-x+Mui2-mui2-mu2)-4.*mui2*mu2 ) ) / (1.-x+Mui2);
   double zp = .5*( 1.-x+Mui2+mui2-mui2 +
   		   sqrt( sqr(1.-x+Mui2-mui2-mu2)-4.*mui2*mu2 ) ) / (1.-x+Mui2);
 
   if ( x > xp || z < zm || z > zp ) {
     jacobian(0.0);
     return false;
   }
   realEmitterMomentum().setMass(sqrt(mi2));
   realEmitterMomentum().rescaleEnergy();
   realEmissionMomentum().setMass(sqrt(m2));
   realEmissionMomentum().rescaleEnergy();
   realSpectatorMomentum().setMass(ZERO);
   realSpectatorMomentum().rescaleEnergy();
   return true;
 
 }
 
 Energy FIMassiveInvertedTildeKinematics::lastPt() const {
   Energy2 mi2 = sqr(realEmitterData()->hardProcessMass());
   Energy2 m2  = sqr(realEmissionData()->hardProcessMass());
   Energy2 Mi2 = sqr(bornEmitterData()->hardProcessMass());
 
   Energy2 scale = Mi2 - (realEmitterMomentum()+realEmissionMomentum()-realSpectatorMomentum()).m2();
   double x = subtractionParameters()[0];
   double z = subtractionParameters()[1];
 
   return sqrt( z*(1.-z)*(1.-x)/x*scale -
 	       ((1.-z)*mi2+z*m2-z*(1.-z)*Mi2) );
 
 }
 
 double FIMassiveInvertedTildeKinematics::lastZ() const {
   return subtractionParameters()[1];
 }
 
 Energy FIMassiveInvertedTildeKinematics::ptMax() const {
   Energy2 mi2 = sqr(realEmitterData()->hardProcessMass());
   Energy2 m2  = sqr(realEmissionData()->hardProcessMass());
   Energy2 Mi2 = sqr(bornEmitterData()->hardProcessMass());
   double x = spectatorX();
   // s^star/x
   Energy2 scale=2.*bornEmitterMomentum()*bornSpectatorMomentum();
   Energy2 s = scale * (1.-x)/x + Mi2;
   Energy ptmax = .5 * sqrt(s) * rootOfKallen( s/s, mi2/s, m2/s );
   return ptmax > 0.*GeV ? ptmax : 0.*GeV;
 }
 
 pair<double,double> FIMassiveInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
   hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
+  if(pt>hardPt) return make_pair(0.5,0.5);
   Energy2 mi2 = sqr(realEmitterData()->hardProcessMass());
   Energy2 m2  = sqr(realEmissionData()->hardProcessMass());
   Energy2 Mi2 = sqr(bornEmitterData()->hardProcessMass());
   // s^star/x
   Energy2 scale=2.*bornEmitterMomentum()*bornSpectatorMomentum();
   Energy2 s = scale * (1.-spectatorX())/spectatorX() +  Mi2;
 
   double zm = .5*( 1.+(mi2-m2)/s - rootOfKallen(s/s,mi2/s,m2/s) *
 		    sqrt( 1.-sqr(pt/hardPt) ) );
   double zp = .5*( 1.+(mi2-m2)/s + rootOfKallen(s/s,mi2/s,m2/s) *
 		    sqrt( 1.-sqr(pt/hardPt) ) );
   return make_pair(zm, zp);
 
 }
 
 
 // If needed, insert default implementations of virtual function defined
 // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
 
 
 void FIMassiveInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
 }
 
 void FIMassiveInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
 }
 
 void FIMassiveInvertedTildeKinematics::Init() {
 
   static ClassDocumentation<FIMassiveInvertedTildeKinematics> documentation
     ("FIMassiveInvertedTildeKinematics inverts the final-initial tilde "
      "kinematics.");
 
 }
 
 // *** 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<FIMassiveInvertedTildeKinematics,InvertedTildeKinematics>
 describeHerwigFIMassiveInvertedTildeKinematics("Herwig::FIMassiveInvertedTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc
@@ -1,148 +1,149 @@
 // -*- C++ -*-
 //
 // IFLightInvertedTildeKinematics.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2012 The Herwig Collaboration
 //
 // Herwig is licenced under version 2 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 IFLightInvertedTildeKinematics class.
 //
 
 #include "IFLightInvertedTildeKinematics.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Repository/EventGenerator.h"
 
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 IFLightInvertedTildeKinematics::IFLightInvertedTildeKinematics() {}
 
 IFLightInvertedTildeKinematics::~IFLightInvertedTildeKinematics() {}
 
 IBPtr IFLightInvertedTildeKinematics::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr IFLightInvertedTildeKinematics::fullclone() const {
   return new_ptr(*this);
 }
 
 bool IFLightInvertedTildeKinematics::doMap(const double * r) {
 
   if ( ptMax() < ptCut() ) {
     jacobian(0.0);
     return false;
   }
 
   Lorentz5Momentum emitter = bornEmitterMomentum();
   Lorentz5Momentum spectator = bornSpectatorMomentum();
 
   double mapping = 1.0;
   pair<Energy,double> ptz = generatePtZ(mapping,r);
   if ( mapping == 0.0 ) {
     jacobian(0.0);
     return false;
   }
 
   Energy pt = ptz.first;
   double z = ptz.second;
 
   double ratio = sqr(pt/lastScale());
   double rho = 1. - 4.*ratio*z*(1.-z) / sqr(1. - z + ratio);
   if ( rho < 0. ) {
     jacobian(0.0);
     return false;
   }
 
   double x = 0.5*(1./ratio)*(1.-z+ratio)*(1.-sqrt(rho));
   double u = 0.5*(1./(1.-z))*(1.-z+ratio)*(1.-sqrt(rho));
 
   if ( x < emitterX() || x > 1. || 
        u < 0. || u > 1. ) {
     jacobian(0.0);
     return false;
   }
 
   mapping *= (1.-x)/((1.-z)*(z*(1.-z)+sqr(x-z)));
   jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
 
   double phi = 2.*Constants::pi*r[2];
   Lorentz5Momentum kt = getKt(emitter,spectator,pt,phi,true);
 
   subtractionParameters().resize(2);
   subtractionParameters()[0] = x;
   subtractionParameters()[1] = u;
 
   realEmitterMomentum() = (1./x)*emitter;
   realEmissionMomentum() = ((1.-x)*(1.-u)/x)*emitter + u*spectator + kt;
   realSpectatorMomentum() = ((1.-x)*u/x)*emitter + (1.-u)*spectator - kt;
 
   realEmitterMomentum().setMass(ZERO);
   realEmitterMomentum().rescaleEnergy();
   realEmissionMomentum().setMass(ZERO);
   realEmissionMomentum().rescaleEnergy();
   realSpectatorMomentum().setMass(ZERO);
   realSpectatorMomentum().rescaleEnergy();
 
   return true;
 
 }
 
 Energy IFLightInvertedTildeKinematics::lastPt() const {
 
   Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
   double x = subtractionParameters()[0];
   double u = subtractionParameters()[1];
   return scale * sqrt(u*(1.-u)*(1.-x)/x);
 
 }
 
 Energy IFLightInvertedTildeKinematics::ptMax() const {
   double x = emitterX();
   return sqrt((1.-x)/x)*lastScale()/2.;
 }
 
 pair<double,double> IFLightInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
   hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
+  if(pt>hardPt) return make_pair(0.5,0.5);
   double s = sqrt(1.-sqr(pt/hardPt));
   double x = emitterX();
   return make_pair(0.5*(1.+x-(1.-x)*s),0.5*(1.+x+(1.-x)*s));
 }
 
 double IFLightInvertedTildeKinematics::lastZ() const {
   double x = subtractionParameters()[0];
   double u = subtractionParameters()[1];
   return 1. - (1.-x)*(1.-u);
 }
 
 // If needed, insert default implementations of virtual function defined
 // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
 
 
 void IFLightInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
 }
 
 void IFLightInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
 }
 
 void IFLightInvertedTildeKinematics::Init() {
 
   static ClassDocumentation<IFLightInvertedTildeKinematics> documentation
     ("IFLightInvertedTildeKinematics inverts the initial-final tilde "
      "kinematics.");
 
 }
 
 // *** 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<IFLightInvertedTildeKinematics,InvertedTildeKinematics>
 describeHerwigIFLightInvertedTildeKinematics("Herwig::IFLightInvertedTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/IFMassiveInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/IFMassiveInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/IFMassiveInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/IFMassiveInvertedTildeKinematics.cc
@@ -1,152 +1,153 @@
 // -*- C++ -*-
 //
 // IFMassiveInvertedTildeKinematics.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2012 The Herwig Collaboration
 //
 // Herwig is licenced under version 2 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 IFMassiveInvertedTildeKinematics class.
 //
 
 #include "IFMassiveInvertedTildeKinematics.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Repository/EventGenerator.h"
 
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 IFMassiveInvertedTildeKinematics::IFMassiveInvertedTildeKinematics() {}
 
 IFMassiveInvertedTildeKinematics::~IFMassiveInvertedTildeKinematics() {}
 
 IBPtr IFMassiveInvertedTildeKinematics::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr IFMassiveInvertedTildeKinematics::fullclone() const {
   return new_ptr(*this);
 }
 
 bool IFMassiveInvertedTildeKinematics::doMap(const double * r) { 
   if ( ptMax() < ptCut() ) {
     jacobian(0.0);
     return false;
   }
 
   Lorentz5Momentum emitter = bornEmitterMomentum();
   Lorentz5Momentum spectator = bornSpectatorMomentum();
   Energy2 scale = 2.*(spectator*emitter);
 
   double mapping = 1.0;
   pair<Energy,double> ptz = generatePtZ(mapping,r);
   if ( mapping == 0.0 ){
     jacobian(0.0);
     return false;
   }
 
   Energy pt = ptz.first;
   double z = ptz.second;
   double ratio = sqr(pt)/scale;
 
   // x
   double u = ratio/(1.-z);
   double x = (z*(1.-z)-ratio)/(1.-z-ratio);
   double up = (1.-x) / (1.-x+(x*sqr(bornSpectatorData()->hardProcessMass())/scale));
 
   if ( x < emitterX() || x > 1. || u > up) {
     jacobian(0.0);
     return false;
   }
 
   pt = sqrt(scale*u*(1.-u)*(1.-x));  
   Energy magKt = sqrt(scale*u*(1.-u)*(1.-x)/x - sqr(u*bornSpectatorData()->hardProcessMass()));
 
   // TODO: why not mapping /= sqr(z*(1.-z)-ratio)/(1.-z) ? (see phd thesis, (5.74))
   mapping /= sqr(z*(1.-z)-ratio)/(1.-z-ratio);
   // mapping *= (1.-u)/(1.-2.*u+u*u*alpha)/x;
   jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
 
   double phi = 2.*Constants::pi*r[2];
   Lorentz5Momentum kt = getKt(emitter,spectator,magKt,phi,true);
   subtractionParameters().resize(2);
   subtractionParameters()[0] = x;
   subtractionParameters()[1] = u;
   
   realEmitterMomentum() = (1./x)*emitter;
   realEmissionMomentum() = (-kt*kt-u*u*sqr(bornSpectatorData()->hardProcessMass()))/(u*scale)*emitter +  
     u*spectator - kt;
   realSpectatorMomentum() = (-kt*kt+sqr(bornSpectatorData()->hardProcessMass())*u*(2.-u))/((1.-u)*scale)*emitter + 
     (1.-u)*spectator + kt;
 
   realEmitterMomentum().setMass(ZERO);
   realEmitterMomentum().rescaleEnergy();
   realEmissionMomentum().setMass(ZERO);
   realEmissionMomentum().rescaleEnergy();
   realSpectatorMomentum().setMass(bornSpectatorData()->hardProcessMass());
   realSpectatorMomentum().rescaleEnergy();
   return true;
 
 }
 
 Energy IFMassiveInvertedTildeKinematics::lastPt() const {
   Energy2 scale = 2.*(bornEmitterMomentum()*bornSpectatorMomentum());
   double x = subtractionParameters()[0];
   double u = subtractionParameters()[1];
   // TODO: can't make sense of the following comment
   // there was no factor 1/x in massless case >> check
   //  return scale * sqrt(u*(1.-u)*(1.-x)/x);
   return sqrt(scale*u*(1.-u)*(1.-x));
 }
 
 double IFMassiveInvertedTildeKinematics::lastZ() const {
   double x = subtractionParameters()[0];
   double u = subtractionParameters()[1];
   return 1. - (1.-x)*(1.-u);
 }
 
 Energy IFMassiveInvertedTildeKinematics::ptMax() const {
   Energy2 scale = 2.*(bornEmitterMomentum()*bornSpectatorMomentum());
   double x = emitterX();
   return sqrt(scale*(1.-x))/2.;
 }
 
 pair<double,double> IFMassiveInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
   hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
+  if(pt>hardPt) return make_pair(0.5,0.5);
   double s = sqrt(1.-sqr(pt/hardPt));
   double x = emitterX();
   return make_pair(0.5*(1.+x-(1.-x)*s),0.5*(1.+x+(1.-x)*s));
 }
 
 
 
 // If needed, insert default implementations of virtual function defined
 // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
 
 
 void IFMassiveInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
 }
 
 void IFMassiveInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
 }
 
 void IFMassiveInvertedTildeKinematics::Init() {
 
   static ClassDocumentation<IFMassiveInvertedTildeKinematics> documentation
     ("IFMassiveInvertedTildeKinematics inverts the initial-final tilde "
      "kinematics.");
 
 }
 
 // *** 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<IFMassiveInvertedTildeKinematics,InvertedTildeKinematics>
 describeHerwigIFMassiveInvertedTildeKinematics("Herwig::IFMassiveInvertedTildeKinematics", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.cc
@@ -1,144 +1,145 @@
 // -*- C++ -*-
 //
 // IILightInvertedTildeKinematics.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2012 The Herwig Collaboration
 //
 // Herwig is licenced under version 2 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 IILightInvertedTildeKinematics class.
 //
 
 #include "IILightInvertedTildeKinematics.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Repository/EventGenerator.h"
 
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 IBPtr IILightInvertedTildeKinematics::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr IILightInvertedTildeKinematics::fullclone() const {
   return new_ptr(*this);
 }
 
 bool IILightInvertedTildeKinematics::doMap(const double * r) {
 
   if ( ptMax() < ptCut() ) {
     jacobian(0.0);
     return false;
   }
 
   Lorentz5Momentum emitter = bornEmitterMomentum();
   Lorentz5Momentum spectator = bornSpectatorMomentum();
 
   double mapping = 1.0;
   pair<Energy,double> ptz = generatePtZ(mapping,r,2.0);
   if ( mapping == 0.0 ) {
     jacobian(0.0);
     return false;
   }
 
   Energy pt = ptz.first;
   double z = ptz.second;
 
   double ratio = sqr(pt/lastScale());
   double x = z*(1.-z) / ( 1. - z + ratio );
   double v = ratio * z / ( 1. - z + ratio );
 
   if ( x < emitterX() || x > 1. ||
        v < 0. || v > 1.-x ) {
     jacobian(0.0);
     return false;
   }
 
   mapping /= z*(1.-z);
   jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
 
   subtractionParameters().resize(2);
   subtractionParameters()[0] = x;
   subtractionParameters()[1] = v;
 
   double phi = 2.*Constants::pi*r[2];
   Lorentz5Momentum kt = getKt(emitter,spectator,pt,phi);
 
   realEmitterMomentum() = (1./x)*emitter;
   realEmissionMomentum() = ((1.-x-v)/x)*emitter+v*spectator+kt;
   realSpectatorMomentum() = spectator;
 
   realEmitterMomentum().setMass(ZERO);
   realEmitterMomentum().rescaleEnergy();
   realEmissionMomentum().setMass(ZERO);
   realEmissionMomentum().rescaleEnergy();
   realSpectatorMomentum().setMass(ZERO);
   realSpectatorMomentum().rescaleEnergy();
 
   K = realEmitterMomentum() + realSpectatorMomentum() - realEmissionMomentum();
   K2 = K.m2();
 
   Ktilde = emitter + spectator;
   KplusKtilde = K + Ktilde;
 
   KplusKtilde2 = KplusKtilde.m2();
 
   return true;
 
 }
 
 Energy IILightInvertedTildeKinematics::lastPt() const {
 
   Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
   double x = subtractionParameters()[0];
   double v = subtractionParameters()[1];
   return scale * sqrt(v*(1.-x-v)/x);
 
 }
 
 double IILightInvertedTildeKinematics::lastZ() const {
   double x = subtractionParameters()[0];
   double v = subtractionParameters()[1];
   return x + v;
 }
 
 Energy IILightInvertedTildeKinematics::ptMax() const {
   return 0.5*(1.-emitterX())/sqrt(emitterX())*lastScale();
 }
 
-pair<double,double> IILightInvertedTildeKinematics::zBounds(Energy pt, Energy) const {
-  double root = sqr(1.-emitterX())-4*emitterX()*sqr(pt/lastScale());
-  root=sqrt(max(root,0.));
+pair<double,double> IILightInvertedTildeKinematics::zBounds(Energy pt, Energy hardPt) const {
+  hardPt = hardPt == ZERO ? ptMax() : min(hardPt,ptMax());
+  if(pt>hardPt) return make_pair(0.5,0.5);
+  double root = (1.-emitterX())*sqrt(1.-sqr(pt/hardPt));
   return make_pair(0.5*( 1.+emitterX() - root),0.5*( 1.+emitterX() + root));
 }
 
 void IILightInvertedTildeKinematics::persistentOutput(PersistentOStream & os) const {
   os << ounit(K,GeV) << ounit(K2,GeV2) << ounit(Ktilde,GeV)
      << ounit(KplusKtilde,GeV) << ounit(KplusKtilde2,GeV2);
 }
 
 void IILightInvertedTildeKinematics::persistentInput(PersistentIStream & is, int) {
   is >> iunit(K,GeV) >> iunit(K2,GeV2) >> iunit(Ktilde,GeV)
      >> iunit(KplusKtilde,GeV) >> iunit(KplusKtilde2,GeV2);
 }
 
 void IILightInvertedTildeKinematics::Init() {
 
   static ClassDocumentation<IILightInvertedTildeKinematics> documentation
     ("IILightInvertedTildeKinematics inverts the initial-initial tilde "
      "kinematics.");
 
 }
 
 // *** 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<IILightInvertedTildeKinematics,InvertedTildeKinematics>
 describeHerwigIILightInvertedTildeKinematics("Herwig::IILightInvertedTildeKinematics", "Herwig.so");