diff --git a/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.cc b/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.cc
--- a/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.cc
+++ b/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.cc
@@ -1,563 +1,571 @@
 // -*- C++ -*-
 //
 // PhasespaceHelpers.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.
 //
 
 #include "PhasespaceHelpers.h"
 #include "Herwig/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
 
 using namespace Herwig;
 using namespace Herwig::PhasespaceHelpers;
 using namespace Herwig::RandomHelpers;
 
 Energy PhasespaceInfo::generateMass(tcPDPtr data, 
                                     const pair<Energy,Energy>& range) {
 
   double xlow = sqr(range.first)/sHat;
   if ( range.first < ZERO )
     xlow = -xlow;
 
   double xup = sqr(range.second)/sHat;
   if ( range.second < ZERO )
     xup = -xup;
 
   double mu = sqr(data->hardProcessMass())/sHat;
   double gamma = sqr(data->hardProcessWidth())/sHat;
 
   if ( gamma < 1e-14 )
     gamma = 0.0;
 
   if ( M0 != ZERO )
     x0 = M0/sqrtSHat;
   if ( Mc != ZERO )
     xc = Mc/sqrtSHat;
 
   double r = rnd();
 
   pair<double,double> event;
   
   if ( gamma == 0. ) {
 
     if ( mu < xlow || mu > xup ) {
       if ( abs(xlow-mu) < xc )
         xlow = mu;
       if ( abs(xup-mu) < xc )
         xup = mu;
     }
 
     if ( mu < xlow || mu > xup ) {
       event =
       generate(inverse(mu,xlow,xup),r);
     } else {
       pair<double,double> pLeft(xlow,xlow < mu-x0 ? mu-x0 : xlow);
       pair<double,double> pRight(xup > mu+x0 ? mu+x0 : xup,xup);
       pair<double,double> fLeft(pLeft.second < mu-x0 ? mu-x0 : pLeft.second,mu-xc);
       if ( fLeft.first >= fLeft.second ) fLeft.first = fLeft.second;
       pair<double,double> fRight(mu+xc,pRight.first > mu+x0 ? mu+x0 : pRight.first);
       if ( fRight.first >= fRight.second ) fRight.second = fRight.first;
 
       
       const bool pL= abs( pLeft.first - pLeft.second ) < 1e-14;
       const bool fL= abs( fLeft.first - fLeft.second ) < 1e-14;
       const bool fR= abs( fRight.first - fRight.second ) < 1e-14;
       const bool pR= abs( pRight.first - pRight.second ) < 1e-14;
 
 
       if ( !pL && !fL && !fR && !pR ) {
         event =
         generate((piecewise(),
                   inverse(mu,pLeft.first,pLeft.second),
                   match(flat(fLeft.first,fLeft.second))) +
                  match((piecewise(),
                         flat(fRight.first,fRight.second),
                         match(inverse(mu,pRight.first,pRight.second)))),
                  r);
       } else if ( pL && !fL && !fR && !pR ) {
         event =
         generate(flat(fLeft.first,fLeft.second) +
                  match((piecewise(),
                         flat(fRight.first,fRight.second),
                         match(inverse(mu,pRight.first,pRight.second)))), r);
       } else if ( !pL && !fL && !fR && pR ) {
         event =
         generate((piecewise(),
                   inverse(mu,pLeft.first,pLeft.second),
                   match(flat(fLeft.first,fLeft.second))) +
                  match(flat(fRight.first,fRight.second)),
                  r);
       } else if ( pL && fL && !fR && !pR ) {
         event =
         generate((piecewise(),flat(fRight.first,fRight.second),
                   match(inverse(mu,pRight.first,pRight.second))), r);
       } else if ( !pL && !fL && fR && pR ) {
         event =
         generate((piecewise(),
                   inverse(mu,pLeft.first,pLeft.second),
                   match(flat(fLeft.first,fLeft.second))),r);
       } else if ( pL && !fL && !fR && pR ) {
         event = generate(flat(fLeft.first,fLeft.second) +
                  match(flat(fRight.first,fRight.second)),r);
       } else if ( pL && fL && !fR && pR ) {
         event = generate(flat(fRight.first,fRight.second),r);
       } else if ( pL && !fL && fR && pR ) {
         event = generate(flat(fLeft.first,fLeft.second),r);
       } else if ( pL && fL && fR && pR ) {
         throw Veto();
       } else assert(false);
     }
   } else {
     event = generate(breitWigner(mu,gamma,xlow,xup),r);
   }
 
   if ( abs(event.first) < xc )
     throw Veto();
 
   weight *= event.second;
 
   Energy res = sqrt(abs(event.first)*sHat);
   if ( event.first < 0. )
     res = -res;
 
   return res;
 }
 
 Lorentz5Momentum PhasespaceInfo::generateKt(const Lorentz5Momentum& p1,
                                             const Lorentz5Momentum& p2,
                                             Energy pt) {
 
   double phi = 2.*Constants::pi*rnd();
   weight *= 2.*Constants::pi;
 
   Lorentz5Momentum P = p1 + p2;
 
   Energy2 Q2 = abs(P.m2());
 
   Lorentz5Momentum Q =  Lorentz5Momentum(ZERO,ZERO,ZERO,sqrt(Q2),sqrt(Q2));
 
   bool boost =
     abs((P-Q).vect().mag2()/GeV2) > 1e-8 ||
     abs((P-Q).t()/GeV) > 1e-4;
   boost &= (P*Q-Q.mass2())/GeV2 > 1e-8;
 
   Lorentz5Momentum inFrame1;
   if ( boost )
     inFrame1 = p1 + ((P*p1-Q*p1)/(P*Q-Q.mass2()))*(P-Q);
   else
     inFrame1 = p1;
 
   Energy ptx = inFrame1.x();
   Energy pty = inFrame1.y();
   Energy q = 2.*inFrame1.z();
 
   Energy Qp = sqrt(4.*(sqr(ptx)+sqr(pty))+sqr(q));
   Energy Qy = sqrt(4.*sqr(pty)+sqr(q));
 
   double cPhi = cos(phi);
   double sPhi = sqrt(1.-sqr(cPhi));
   if ( phi > Constants::pi )
     sPhi = -sPhi;
 
   Lorentz5Momentum kt;
 
   kt.setT(ZERO);
   kt.setX(pt*Qy*cPhi/Qp);
   kt.setY(-pt*(4*ptx*pty*cPhi/Qp+q*sPhi)/Qy);
   kt.setZ(2.*pt*(-ptx*q*cPhi/Qp + pty*sPhi)/Qy);
 
   if ( boost )
     kt = kt + ((P*kt-Q*kt)/(P*Q-Q.mass2()))*(P-Q);
   kt.setMass(-pt);
   kt.rescaleRho();
 
   return kt;
 
 }
 
 void PhasespaceTree::setup(const Tree2toNDiagram& diag, 
                            int pos) {
   doMirror = false;
 
   pair<int,int> dchildren =  diag.children(pos);
 
   data = diag.allPartons()[pos];
 
   spacelike = pos < diag.nSpace();
 
   if ( pos == 0 )
     externalId = 0;
 
   if ( dchildren.first == -1 ) {
     externalId = diag.externalId(pos);
     leafs.insert(externalId);
     return;
   }
 
   children.push_back(PhasespaceTree());
   children.back().setup(diag,dchildren.first);
   children.push_back(PhasespaceTree());
   children.back().setup(diag,dchildren.second);
 
   if ( !children[0].children.empty() &&
       children[1].children.empty() &&
       !spacelike )
     swap(children[0],children[1]);
   if ( spacelike &&
       !children[0].spacelike )
     swap(children[0],children[1]);
 
   copy(children[0].leafs.begin(),children[0].leafs.end(),
        inserter(leafs,leafs.begin()));
   copy(children[1].leafs.begin(),children[1].leafs.end(),
        inserter(leafs,leafs.begin()));
 
 }
 
 void PhasespaceTree::setupMirrored(const Tree2toNDiagram& diag, 
                                    int pos) {
 
   doMirror = true;
 
   spacelike = pos < diag.nSpace();
 
   pair<int,int> dchildren;
   if (pos != 0 && spacelike)
     dchildren = {diag.parent(pos), diag.children(diag.parent(pos)).second};
   else if ( !spacelike ) dchildren = diag.children(pos);
   else                   dchildren = {-1,-1};
 
   data = diag.allPartons()[pos];
 
   if ( pos == diag.nSpace() - 1 )
     externalId = 1;
 
   if ( dchildren.first == -1 ) {
     externalId = diag.externalId(pos);
     leafs.insert(externalId);
     return;
   }
   
 
 
   children.push_back(PhasespaceTree());
   children.back().setupMirrored(diag,dchildren.first);
   children.push_back(PhasespaceTree());
   children.back().setupMirrored(diag,dchildren.second);
 
   if ( !children[0].children.empty() &&
       children[1].children.empty() &&
       !spacelike )
     swap(children[0],children[1]);
   if ( spacelike &&
       !children[0].spacelike ) {
     assert (false);
   }
 
   copy(children[0].leafs.begin(),children[0].leafs.end(),
        inserter(leafs,leafs.begin()));
   copy(children[1].leafs.begin(),children[1].leafs.end(),
        inserter(leafs,leafs.begin()));
 
 }
 
 void PhasespaceTree::print(int in) {
   for (int i = 0; i != in; i++)
     cerr << "   ";
   cerr << " |- "  << data->PDGName() << " " << externalId << "\n" << flush;
   if ( !children.empty() ) {
     children[1].print(in+1);
     children[0].print(in+int(!spacelike));
   }
   else {
     cerr << "\n";
   }
 }
 
 void PhasespaceTree::init(const vector<Lorentz5Momentum>& meMomenta) {
 
   if ( children.empty() ) {
     massRange.first = meMomenta[externalId].mass();
     massRange.second = massRange.first;
     if ( !doMirror && externalId == 1 )
       momentum = meMomenta[1];
     if ( doMirror && externalId == 0 )
       momentum = meMomenta[0];
     momentum.setMass(meMomenta[externalId].mass());
     return;
   }
 
   children[0].init(meMomenta);
   children[1].init(meMomenta);
 
   if ( !children[0].spacelike &&
       !children[1].spacelike ) {
     massRange.first =
     children[0].massRange.first +
     children[1].massRange.first;
   }
 
 }
 
 void PhasespaceTree::generateKinematics(PhasespaceInfo& info,
                                         vector<Lorentz5Momentum>& meMomenta) {
   
   if ( !doMirror && externalId == 0 ) {
     init(meMomenta);
-    momentum = meMomenta[0];
+    Energy2 s = (meMomenta[0]+meMomenta[1]).m2();
+    double sign = meMomenta[0].z() >= ZERO ? 1. : -1;
+    momentum = Lorentz5Momentum(ZERO,ZERO,sign*sqrt(s)/2.,sqrt(s)/2.,ZERO);
+    backwardMomentum = Lorentz5Momentum(ZERO,ZERO,-sign*sqrt(s)/2.,sqrt(s)/2.,ZERO);
   }
   else if ( doMirror && externalId == 1) {
     init(meMomenta);
-    momentum = meMomenta[1];
+    Energy2 s = (meMomenta[0]+meMomenta[1]).m2();
+    double sign = meMomenta[0].z() >= ZERO ? 1. : -1;
+    momentum = Lorentz5Momentum(ZERO,ZERO,-sign*sqrt(s)/2.,sqrt(s)/2.,ZERO);
+    backwardMomentum = Lorentz5Momentum(ZERO,ZERO,sign*sqrt(s)/2.,sqrt(s)/2.,ZERO);
   }
 
   if ( children.empty() ) {
     if ( ( !doMirror && externalId != 1 )
         || ( doMirror && externalId !=0 ) )
       meMomenta[externalId] = momentum;
     return;
   }
 
     // s-channel
   if ( ( !doMirror && externalId == 0 &&
           children[0].externalId == 1 )
      || ( doMirror && externalId == 1 &&
           children[0].externalId == 0 ) ) {
         children[1].momentum = meMomenta[0] + meMomenta[1];
         children[1].momentum.setMass(info.sqrtSHat);
         children[1].momentum.rescaleEnergy();
         children[1].generateKinematics(info,meMomenta);
         return;
       }
 
   if ( !spacelike ) {
 
     Energy mij = momentum.mass();
     Energy mi,mj;
     
       // work out the mass for the first child
     if ( !children[0].children.empty() ) {
       Energy sumOthers = ZERO;
       for ( size_t k = 2; k < meMomenta.size(); ++k )
         if ( children[1].leafs.find(k) != children[1].leafs.end() )
           sumOthers += meMomenta[k].mass();
       children[0].massRange.second = momentum.mass() - sumOthers;
       if ( children[0].massRange.second < children[0].massRange.first )
         throw Veto();
       if ( children[0].massRange.second > momentum.mass() )
         throw Veto();
       mi = info.generateMass(children[0].data,children[0].massRange);
       children[0].momentum.setMass(mi);
     } else {
       mi = children[0].momentum.mass();
     }
 
       // work out the mass for the second child
     if ( !children[1].children.empty() ) {
       children[1].massRange.second = momentum.mass()-children[0].momentum.mass();
       if ( children[1].massRange.second < children[1].massRange.first )
         throw Veto();
       mj = info.generateMass(children[1].data,children[1].massRange);
       children[1].momentum.setMass(mj);
     } else {
       mj = children[1].momentum.mass();
     }
 
     Energy2 mij2 = sqr(mij);
     Energy2 mi2 = sqr(mi);
     Energy2 mj2 = sqr(mj);
 
       // perform the decay
     Energy4 lambda2 = sqr(mij2-mi2-mj2)-4.*mi2*mj2;
     if ( lambda2 <= ZERO )
       throw Veto();
     Energy2 lambda = sqrt(lambda2);
     double phi = 2.*Constants::pi*info.rnd();
     double cosPhi = cos(phi);
     double sinPhi = sqrt(1.-sqr(cosPhi));
     if ( phi > Constants::pi )
       sinPhi = -sinPhi;
     info.weight *= Constants::pi*lambda/(2.*mij2);
     double cosTheta = 2.*info.rnd() - 1.;
     double sinTheta = sqrt(1.-sqr(cosTheta));
     Energy p = lambda/(2.*mij);
     children[0].momentum.setX(p*cosPhi*sinTheta);
     children[0].momentum.setY(p*sinPhi*sinTheta);
     children[0].momentum.setZ(p*cosTheta);
     children[0].momentum.rescaleEnergy();
     if ( momentum.m2() <= ZERO )
       throw Veto();
     Boost out = momentum.boostVector();
     if ( out.mag2() > Constants::epsilon ) {
       children[0].momentum.boost(out);
     }
     children[1].momentum = momentum - children[0].momentum;
     children[1].momentum.setMass(mj);
     children[1].momentum.rescaleEnergy();
 
       // go on with next branchings
     children[0].generateKinematics(info,meMomenta);
     children[1].generateKinematics(info,meMomenta);
 
     return;
 
   }
 
     // get the minimum mass of the `W' system
   Energy Wmin = ZERO;
   PhasespaceTree* current = &children[0];
   while ( !(current->children.empty()) ) {
     Wmin += current->children[1].massRange.first;
     current = &(current->children[0]);
   }
 
     // get the CM energy avaialble
-  Energy2 s = (momentum+meMomenta[int(!doMirror)]).m2();
+  Energy2 s = (momentum+backwardMomentum).m2();
   if ( s <= ZERO )
     throw Veto();
 
     // generate a mass for the timelike child
   Energy mi;
   if ( !children[1].children.empty() ) {
     children[1].massRange.second = sqrt(s)-Wmin;
     if ( children[1].massRange.second < children[1].massRange.first )
       throw Veto();
     mi = info.generateMass(children[1].data,children[1].massRange);
     children[1].momentum.setMass(mi);
   } else {
     mi = children[1].momentum.mass();
   }
   Energy2 mi2 = sqr(mi);
  
     // wether or not this is the last 2->2 scatter
   bool lastScatter = children[0].children[0].children.empty();
 
     // `W' mass relevant for the other boundaries
   Energy MW = Wmin;
 
     // generate a mass for second outgoing leg, if needed
   if ( lastScatter )
     if ( !children[0].children[1].children.empty() ) {
         // get the maximum `W' mass
       Energy Wmax = sqrt(s)-children[1].momentum.mass();
       children[0].children[1].massRange.second = Wmax;
       if ( children[0].children[1].massRange.second <
           children[0].children[1].massRange.first )
         throw Veto();
       MW = info.generateMass(children[0].children[1].data,
                              children[0].children[1].massRange);
       children[0].children[1].momentum.setMass(MW);
     }
   Energy2 MW2 = sqr(MW);
 
   Energy ma = momentum.mass();
   Energy2 ma2 = sqr(ma);
   if ( ma < ZERO )
     ma2 = -ma2;
-  Energy mb = meMomenta[int(!doMirror)].mass();
+  Energy mb = backwardMomentum.mass();
   Energy2 mb2 = sqr(mb);
+  if ( mb < ZERO )
+    mb2 = -mb2;
 
     // pick the ys variable
   Energy2 ys = ZERO;
   if ( !lastScatter ) {
     ys = info.rnd()*(sqr(sqrt(s)-mi)-MW2);
     info.weight *= (sqr(sqrt(s)-mi)-MW2)/info.sHat;
   }
 
   Energy4 lambda2 = sqr(s-ma2-mb2)-4.*ma2*mb2;
   if ( lambda2 <= ZERO ) {
     throw Veto();
   }
   Energy2 lambda = sqrt(lambda2);
   info.weight *= info.sHat/(4.*lambda);
 
     // get the boundaries on the momentum transfer
   Energy4 rho2 = sqr(s-ys-MW2-mi2)-4.*mi2*(ys+MW2);
   if ( rho2 < ZERO )
     throw Veto();
   Energy2 rho = sqrt(rho2);
   Energy4 tau2 =
   ys*(ma2-mb2+s)
   - sqr(s)+s*(ma2+mb2+mi2+MW2)-(mi2-MW2)*(ma2-mb2);
   pair<Energy2,Energy2> tBounds
   ((tau2-rho*lambda)/(2.*s),(tau2+rho*lambda)/(2.*s));
   children[0].massRange.first = sqrt(abs(tBounds.first));
   if ( tBounds.first < ZERO )
     children[0].massRange.first = -children[0].massRange.first;
   children[0].massRange.second = sqrt(abs(tBounds.second));
   if ( tBounds.second < ZERO )
     children[0].massRange.second = -children[0].massRange.second;
 
     // generate a momentum transfer
   Energy mai = info.generateMass(children[0].data,children[0].massRange);
   children[0].momentum.setMass(mai);
   Energy2 t = sqr(mai);
   if ( mai < ZERO )
     t = -t;
 
   Energy2 u = -s -t + ys + ma2 + mb2 + mi2 + MW2;
   Energy2 st = s - ma2 - mb2;
   Energy2 tt = t - mi2 - ma2;
   Energy2 ut = u - mi2 - mb2;
 
     // get the timelike momentum
   double xa = (-st*ut+2.*mb2*tt)/lambda2;
   double xb = (-st*tt+2.*ma2*ut)/lambda2;
   Energy2 pt2 = (st*tt*ut-ma2*sqr(ut)-mb2*sqr(tt)-mi2*sqr(st)+4.*ma2*mb2*mi2)/lambda2;
   if ( pt2 < ZERO )
     throw Veto();
   Energy pt = sqrt(pt2);
 
   children[1].momentum =
-  xa*momentum + xb*meMomenta[int(!doMirror)] 
-  + info.generateKt(momentum,meMomenta[int(!doMirror)],pt);
+  xa*momentum + xb*backwardMomentum 
+  + info.generateKt(momentum,backwardMomentum,pt);
   children[1].momentum.setMass(mi);
   children[1].momentum.rescaleEnergy();
 
   children[0].momentum =
   momentum - children[1].momentum;
   children[0].momentum.setMass(mai);
   bool changeSign = false;
   if ( children[0].momentum.t() < ZERO && mai > ZERO ) changeSign = true;
   if ( mai < ZERO )
     children[0].momentum.rescaleRho();
   else
     children[0].momentum.rescaleEnergy();
   if ( changeSign ) children[0].momentum.setT(-children[0].momentum.t());
  
   children[1].generateKinematics(info,meMomenta);
 
   if ( !lastScatter ) {
     children[0].generateKinematics(info,meMomenta);
   } else {
     children[0].children[1].momentum =
-    meMomenta[int(!doMirror)] + children[0].momentum;
+    backwardMomentum + children[0].momentum;
     children[0].children[1].momentum.setMass(MW);
     children[0].children[1].momentum.rescaleEnergy();
     children[0].children[1].generateKinematics(info,meMomenta);
   }
 
 }
 
 
 void PhasespaceTree::put(PersistentOStream& os) const {
   os << children.size();
   if ( !children.empty() ) {
     children[0].put(os);
     children[1].put(os);
   }
   os << data << externalId << leafs << spacelike << doMirror;
 }
 
 void PhasespaceTree::get(PersistentIStream& is) {
   size_t nc; is >> nc;
   assert(nc == 0 || nc == 2);
   if ( nc == 2 ) {
     children.resize(2,PhasespaceTree());
     children[0].get(is);
     children[1].get(is);
   }
   is >> data >> externalId >> leafs >> spacelike >> doMirror;
 }
diff --git a/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.h b/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.h
--- a/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.h
+++ b/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.h
@@ -1,189 +1,194 @@
 // -*- C++ -*-
 //
 // PhasespaceHelpers.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_PhasespaceHelpers_H
 #define HERWIG_PhasespaceHelpers_H
 
 #include "ThePEG/Config/ThePEG.h"
 #include "ThePEG/PDT/ParticleData.h"
 
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 #include "Herwig/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
 
 namespace Herwig {
 
 using namespace ThePEG;
 
   namespace PhasespaceHelpers {
 
     /**
      * \ingroup Matchbox
      * \author Simon Platzer
      * \brief General information for phase space generation
      */
     struct PhasespaceInfo {
 
       /**
        * The center of mass energy squared.
        */
       Energy2 sHat;
 
       /**
        * The center of mass energy.
        */
       Energy sqrtSHat;
 
       /**
        * The phase space weight.
        */
       double weight;
 
       /**
        * The random number generator.
        */
       StreamingRnd rnd;
 
       /**
        * Parameter steering from which on propagator virtualities are
        * sampled flat.
        */
       double x0;
 
       /**
        * Parameter steering at which virtuality singularities of
        * propagators are actually cut off.
        */
       double xc;
 
       /**
        * Parameter steering from which on propagator virtualities are
        * sampled flat.
        */
       Energy M0;
 
       /**
        * Parameter steering at which virtuality singularities of
        * propagators are actually cut off.
        */
       Energy Mc;
 
       /**
        * Generate a mass for the given particle type and mass range.
        */
       Energy generateMass(tcPDPtr, const pair<Energy,Energy>&);
 
       /**
        * Calculate a transverse momentum for the given momenta,
        * invariant pt and azimuth.
        */
       Lorentz5Momentum generateKt(const Lorentz5Momentum& p1,
 				  const Lorentz5Momentum& p2,
 				  Energy pt);
 
     };
 
     /**
      * \ingroup Matchbox
      * \author Simon Platzer, Ken Arnold
      * \brief A phase space tree.
      */
     struct PhasespaceTree {
 
       /**
        * Default constructor
        */
       PhasespaceTree()
 	: massRange(ZERO,ZERO), externalId(-1), 
 	  spacelike(false) {}
 
       /**
        * The particle running along this line.
        */
       tcPDPtr data;
 
       /**
        * The allowed mass range for this line.
        */
       pair<Energy,Energy> massRange;
 
       /**
        * The momentum running along this line.
        */
       Lorentz5Momentum momentum;
 
       /**
+       * A backward momentum, if needed
+       */
+      Lorentz5Momentum backwardMomentum;
+
+      /**
        * The external leg id of this line, if external.
        */
       int externalId;
 
       /**
        * The children lines; if empty this is an external line.
        */
       vector<PhasespaceTree> children;
 
       /**
        * External lines originating from this line.
        */
       set<int> leafs;
 
       /**
        * Wether or not this is a spacelike line.
        */
       bool spacelike;
 
       /**
        * Wether or not this is a mirrored channel.
        */
       bool doMirror;
       
       /**
        * Setup from diagram at given position.
        */
       void setup(const Tree2toNDiagram&, int pos = 0);
 
       /**
        * Setup mirror from diagram at given position.
        */
       void setupMirrored(const Tree2toNDiagram& diag, int pos);
 
       /**
        * Initialize using masses as given by mass() members of the
        * final state momenta
        */
       void init(const vector<Lorentz5Momentum>&);
 
       /**
        * Generate kinematics for the children
        */
       void generateKinematics(PhasespaceInfo&,
 			      vector<Lorentz5Momentum>&);
 
       /**
        * Write phase space tree to ostream
        */
       void put(PersistentOStream&) const;
 
       /**
        * Read phase space tree from istream
        */
       void get(PersistentIStream&);
 
       /**
        * Print tree, only for debugging purposes
        */
       void print(int in = 0);
 
     };
 
   }
 
 }
 
 #endif // HERWIG_PhasespaceHelpers_H