diff --git a/Hadronization/ColourReconnector.cc b/Hadronization/ColourReconnector.cc
--- a/Hadronization/ColourReconnector.cc
+++ b/Hadronization/ColourReconnector.cc
@@ -1,822 +1,822 @@
 // -*- C++ -*-
 //
 // ColourReconnector.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 ColourReconnector class.
 //
 
 #include "ColourReconnector.h"
 #include "Cluster.h"
 
 #include <ThePEG/Utilities/DescribeClass.h>
 #include <ThePEG/Repository/UseRandom.h>
 #include <ThePEG/PDT/StandardMatchers.h>
 #include <ThePEG/Persistency/PersistentOStream.h>
 #include <ThePEG/Persistency/PersistentIStream.h>
 
 #include <ThePEG/Interface/Switch.h>
 #include <ThePEG/Interface/Parameter.h>
 
 #include "Herwig/Utilities/Maths.h"
 
 using namespace Herwig;
 
 using CluVecIt = ColourReconnector::CluVecIt;
 using Constants::pi;
 using Constants::twopi;
 
 
 DescribeClass<ColourReconnector,Interfaced>
 describeColourReconnector("Herwig::ColourReconnector","");
 
 IBPtr ColourReconnector::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr ColourReconnector::fullclone() const {
   return new_ptr(*this);
 }
 
 void ColourReconnector::rearrange(ClusterVector & clusters) {
   if (_clreco == 0) return;
 
   // need at least two clusters
   if (clusters.size() < 2) return;
 
   // do the colour reconnection
   switch (_algorithm) {
     case 0: _doRecoPlain(clusters); break;
     case 1: _doRecoStatistical(clusters); break;
     case 2: _doRecoBaryonic(clusters); break;
   }
 }
 
 
 Energy2 ColourReconnector::_clusterMassSum(const PVector & q, 
                                            const PVector & aq) const {
   const size_t nclusters = q.size();
   assert (aq.size() == nclusters);
   Energy2 sum = ZERO;
   for (size_t i = 0; i < nclusters; i++)
     sum += ( q[i]->momentum() + aq[i]->momentum() ).m2();
   return sum;
 }
 
 
 bool ColourReconnector::_containsColour8(const ClusterVector & cv,
                                          const vector<size_t> & P) const {
   assert (P.size() == cv.size());
   for (size_t i = 0; i < cv.size(); i++) {
     tcPPtr p = cv[i]->colParticle();
     tcPPtr q = cv[P[i]]->antiColParticle();
     if (_isColour8(p, q)) return true;
   }
   return false;
 }
 
 
 void ColourReconnector::_doRecoStatistical(ClusterVector & cv) const {
 
   const size_t nclusters = cv.size();
 
   // initially, enumerate (anti)quarks as given in the cluster vector
   ParticleVector q, aq;
   for (size_t i = 0; i < nclusters; i++) {
     q.push_back( cv[i]->colParticle() );
     aq.push_back( cv[i]->antiColParticle() );
   }
 
   // annealing scheme
   Energy2 t, delta;
   Energy2 lambda = _clusterMassSum(q,aq);
   const unsigned _ntries = _triesPerStepFactor * nclusters;
   
   // find appropriate starting temperature by measuring the largest lambda
   // difference in some dry-run random rearrangements
   {
     vector<Energy2> typical;
     for (int i = 0; i < 10; i++) {
       const pair <int,int> toswap = _shuffle(q,aq,5);
       ParticleVector newaq = aq;
       swap (newaq[toswap.first], newaq[toswap.second]);
       Energy2 newlambda = _clusterMassSum(q,newaq);
       typical.push_back( abs(newlambda - lambda) );
     }
     t = _initTemp * Math::median(typical);
   }
 
   // anneal in up to _annealingSteps temperature steps
   for (unsigned step = 0; step < _annealingSteps; step++) {
 
     // For this temperature step, try to reconnect _ntries times. Stop the
     // algorithm if no successful reconnection happens.
     unsigned nSuccess = 0;
     for (unsigned it = 0; it < _ntries; it++) {
       
       // make a random rearrangement
       const unsigned maxtries = 10;
       const pair <int,int> toswap = _shuffle(q,aq,maxtries);
       const int i = toswap.first;
       const int j = toswap.second;
 
       // stop here if we cannot find any allowed reconfiguration
       if (i == -1) break;
 
       // create a new antiquark vector with the two partons swapped
       ParticleVector newaq = aq;
       swap (newaq[i], newaq[j]);
 
       // Check if lambda would decrease. If yes, accept the reconnection. If no,
       // accept it only with a probability given by the current Boltzmann
       // factor. In the latter case we set p = 0 if the temperature is close to
       // 0, to avoid division by 0.
       Energy2 newlambda = _clusterMassSum(q,newaq);
       delta = newlambda - lambda;
       double prob = 1.0;
       if (delta > ZERO) prob = ( abs(t) < 1e-8*MeV2 ) ? 0.0 : exp(-delta/t);
       if (UseRandom::rnd() < prob) {
         lambda = newlambda;
         swap (newaq, aq);
         nSuccess++;
       }
     }
     if (nSuccess == 0) break;
 
     // reduce temperature
     t *= _annealingFactor;
   }
 
   // construct the new cluster vector
   ClusterVector newclusters;
   for (size_t i = 0; i < nclusters; i++) {
     ClusterPtr cl = new_ptr( Cluster( q[i], aq[i] ) );
     newclusters.push_back(cl);
   }
   swap(newclusters,cv);
   return;
 }
 
 
 void ColourReconnector::_doRecoPlain(ClusterVector & cv) const {
 
   ClusterVector newcv = cv;
 
   // try to avoid systematic errors by randomising the reconnection order
   long (*p_irnd)(long) = UseRandom::irnd;
   random_shuffle( newcv.begin(), newcv.end(), p_irnd ); 
 
   // iterate over all clusters
   for (CluVecIt cit = newcv.begin(); cit != newcv.end(); cit++) {
     // find the cluster which, if reconnected with *cit, would result in the
     // smallest sum of cluster masses
     // NB this method returns *cit if no reconnection partner can be found
     CluVecIt candidate = _findRecoPartner(cit, newcv);
 
     // skip this cluster if no possible reshuffling partner can be found
     if (candidate == cit) continue;
 
     // accept the reconnection with probability _preco.
     if (UseRandom::rnd() < _preco) {
 
       pair <ClusterPtr,ClusterPtr> reconnected = _reconnect(*cit, *candidate);
 
       // Replace the clusters in the ClusterVector. The order of the
       // colour-triplet partons in the cluster vector is retained here.
 
       // replace *cit by reconnected.first
       *cit = reconnected.first;
 
       // replace candidate by reconnected.second
       *candidate = reconnected.second;
     }
   }
 
   swap(cv,newcv);
   return;
 }
 
 namespace {
   inline bool hasDiquark(CluVecIt cit) {
     for(int i = 0; i<(*cit)->numComponents(); i++) {
         if (DiquarkMatcher::Check(*((*cit)->particle(i)->dataPtr()))) 
             return true;
     }
     return false;
   }
 }
 
 
 // Implementation of the baryonic reconnection algorithm
 void ColourReconnector::_doRecoBaryonic(ClusterVector & cv) const {
 
     ClusterVector newcv = cv;
 
     ClusterVector deleted; deleted.reserve(cv.size());
 
     // try to avoid systematic errors by randomising the reconnection order
     long (*p_irnd)(long) = UseRandom::irnd;
     random_shuffle( newcv.begin(), newcv.end(), p_irnd ); 
 
     // iterate over all clusters
     for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) {
       //avoid clusters already containing diuarks
       if (hasDiquark(cit)) continue;
     
       //skip the cluster to be deleted later 3->2 cluster
       if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
         continue;
   
       // Skip all found baryonic clusters, this biases the algorithm but implementing
       // something like re-reconnection is ongoing work
       if ((*cit)->numComponents()==3) continue;
 
       // Find a candidate suitable for reconnection
       CluVecIt baryonic1, baryonic2;
       bool isBaryonicCandidate = false;
       CluVecIt candidate = _findPartnerBaryonic(cit, newcv, 
                                                 isBaryonicCandidate, 
                                                 deleted,
                                                 baryonic1, baryonic2);
 
       // skip this cluster if no possible reconnection partner can be found
       if ( !isBaryonicCandidate && candidate==cit ) 
         continue;
 
       if ( isBaryonicCandidate 
            && UseRandom::rnd() < _precoBaryonic ) {
         deleted.push_back(*baryonic2);
      
         // Function that does the reconnection from 3 -> 2 clusters
         ClusterPtr b1, b2;
         _makeBaryonicClusters(*cit,*baryonic1,*baryonic2, b1, b2);
 
         *cit = b1;
         *baryonic1 = b2;
 
         // Baryonic2 is easily skipped in the next loop 
       }
   
       // Normal 2->2 Colour reconnection
       if ( !isBaryonicCandidate 
            && UseRandom::rnd() < _preco ) {
         auto reconnected = _reconnectBaryonic(*cit, *candidate);
         *cit = reconnected.first;
         *candidate = reconnected.second;
       }
     }
 
     // create a new vector of clusters except for the ones which are "deleted" during
     // baryonic reconnection
     ClusterVector clustervector;
     for ( const auto & cluster : newcv )
       if ( find(deleted.begin(),
                 deleted.end(), cluster) == deleted.end() )
         clustervector.push_back(cluster);
 
     swap(cv,clustervector);
 }
 
 
 
 namespace {
 
 double calculateRapidityRF(const Lorentz5Momentum & q1, 
                            const Lorentz5Momentum & p2) {
   //calculate rapidity wrt the direction of q1
   //angle between the particles in the RF of cluster of q1
 
   // calculate the z component of p2 w.r.t the direction of q1
   const Energy pz = p2.vect() * q1.vect().unit();
   if ( pz == ZERO ) return 0.;
         
   // Transverse momentum of p2 w.r.t the direction of q1
   const Energy pt = sqrt(p2.vect().mag2() - sqr(pz));
   
   // Transverse mass pf p2 w.r.t to the direction of q1
   const Energy mtrans = sqrt(p2.mass()*p2.mass() + (pt*pt));
 
   // Correct formula
   const double y2 = log((p2.t() + abs(pz))/mtrans);
 
   return ( pz < ZERO ) ? -y2 : y2;
 }
 
 }
 
 
 CluVecIt ColourReconnector::_findPartnerBaryonic(
              CluVecIt cl, ClusterVector & cv,
              bool & baryonicCand,
              const ClusterVector& deleted,
              CluVecIt &baryonic1, 
              CluVecIt &baryonic2 ) const {
 
     using Constants::pi;
     using Constants::twopi;
 
     // Returns a candidate for possible reconnection
     CluVecIt candidate = cl;
 
     bool bcand = false;
 
     double maxrap = 0.0;
     double minrap = 0.0;
     double maxrapNormal = 0.0;  
     double minrapNormal = 0.0;
     double maxsumnormal = 0.0;
 
     double maxsum = 0.0;
     double secondsum = 0.0;
 
 
     // boost into RF of cl 
     Lorentz5Momentum cl1 = (*cl)->momentum();
     const Boost boostv(-cl1.boostVector());
     cl1.boost(boostv);
     // boost constituents of cl into RF of cl
     Lorentz5Momentum p1col = (*cl)->colParticle()->momentum();
     Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum(); 
     p1col.boost(boostv);
     p1anticol.boost(boostv);
 
  
     for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
       //avoid looping over clusters containing diquarks
       if ( hasDiquark(cit) ) continue;
       if ( (*cit)->numComponents()==3 ) continue;
       if ( cit==cl ) continue;
 
       //skip the cluster to be deleted later 3->2 cluster
       if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() )
         continue;
 
       if ( (*cl)->isBeamCluster() && (*cit)->isBeamCluster() )
         continue;
 
       // stop it putting far apart clusters together
       if ( ( (**cl).vertex()-(**cit).vertex() ).m() >_maxDistance )
         continue;
   
       const bool Colour8 = 
         _isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() )  
         ||
         _isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) ;
       if ( Colour8 ) continue;
 
 
       // boost constituents of cit into RF of cl
       Lorentz5Momentum p2col = (*cit)->colParticle()->momentum();
       Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum();
 
       p2col.boost(boostv);
       p2anticol.boost(boostv);
 
       // calculate the rapidity of the other constituents of the clusters
       // w.r.t axis of p1anticol.vect.unit
       const double rapq = calculateRapidityRF(p1anticol,p2col);
       const double rapqbar = calculateRapidityRF(p1anticol,p2anticol);
 
       // configuration for normal CR
       if ( rapq > 0.0 && rapqbar < 0.0 
            && rapq > maxrap 
            && rapqbar < minrap ) {
         maxrap = rapq;
         minrap = rapqbar;
         //sum of rapidities of quarks
         const double normalsum = abs(rapq) + abs(rapqbar);
         if ( normalsum > maxsumnormal ) {
           maxsumnormal = normalsum;
           maxrapNormal = rapq;
           minrapNormal = rapqbar;   
           bcand = false;
           candidate = cit;
         }
       } 
 
       if ( rapq < 0.0 && rapqbar >0.0 
           && rapqbar > maxrapNormal 
           && rapq < minrapNormal ) {
         maxrap = rapqbar;
         minrap = rapq;
         const double sumrap = abs(rapqbar) + abs(rapq);
         // first candidate gets here. If second baryonic candidate has higher Ysum than the first
         // one, the second candidate becomes the first one and the first the second.
         if (sumrap > maxsum) {
           if(maxsum != 0){
             baryonic2 = baryonic1;
             baryonic1 = cit;
             bcand = true;
           } else {
             baryonic1 = cit;
           }      
           maxsum = sumrap;
         } else {
           if (sumrap > secondsum && sumrap != maxsum) {
             secondsum = sumrap;
             bcand = true;
             baryonic2 = cit;
           }
         }  
       }
 
     }
 
     if(bcand == true){
       baryonicCand = true;
     }
 
     return candidate;
 }
 
 
 CluVecIt ColourReconnector::_findRecoPartner(CluVecIt cl,
                                              ClusterVector & cv) const {
 
   CluVecIt candidate = cl;
   Energy minMass = 1*TeV;
   for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
 
     // don't even look at original cluster
     if(cit==cl) continue;
 
     // don't allow colour octet clusters
     if ( _isColour8( (*cl)->colParticle(),
               (*cit)->antiColParticle() )  ||
          _isColour8( (*cit)->colParticle(),
               (*cl)->antiColParticle() ) ) {
       continue;
     }
 
     // stop it putting beam remnants together
     if((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue;
 
     // stop it putting far apart clusters together
     if(((**cl).vertex()-(**cit).vertex()).m()>_maxDistance) continue;
 
     // momenta of the old clusters
     Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() + 
                           (*cl)->antiColParticle()->momentum();
     Lorentz5Momentum p2 = (*cit)->colParticle()->momentum() + 
                           (*cit)->antiColParticle()->momentum();
 
     // momenta of the new clusters
     Lorentz5Momentum p3 = (*cl)->colParticle()->momentum() + 
                           (*cit)->antiColParticle()->momentum();
     Lorentz5Momentum p4 = (*cit)->colParticle()->momentum() + 
                           (*cl)->antiColParticle()->momentum();
 
     Energy oldMass = abs( p1.m() ) + abs( p2.m() );
     Energy newMass = abs( p3.m() ) + abs( p4.m() );
 
 
     if ( newMass < oldMass && newMass < minMass ) {
       minMass = newMass;
       candidate = cit;
     }
   }
 
   return candidate;
 }
 
 // forms two baryonic clusters from three clusters
 void ColourReconnector::_makeBaryonicClusters(
                 ClusterPtr &c1, ClusterPtr &c2, 
                 ClusterPtr &c3, 
                 ClusterPtr &newcluster1,
                 ClusterPtr &newcluster2) const{
 
     //make sure they all have 2 components
     assert(c1->numComponents()==2);
     assert(c2->numComponents()==2);
     assert(c3->numComponents()==2);
     //abandon children
     c1->colParticle()->abandonChild(c1);
     c1->antiColParticle()->abandonChild(c1);
     c2->colParticle()->abandonChild(c2);
     c2->antiColParticle()->abandonChild(c2);
     c3->colParticle()->abandonChild(c3);
     c3->antiColParticle()->abandonChild(c3);
 
     newcluster1 = new_ptr(Cluster(c1->colParticle(),c2->colParticle(), c3->colParticle()));
     c1->colParticle()->addChild(newcluster1);
     c2->colParticle()->addChild(newcluster1);
     c3->colParticle()->addChild(newcluster1);
     newcluster1->setVertex(LorentzPoint());
 
     newcluster2 = new_ptr(Cluster(c1->antiColParticle(), c2->antiColParticle(),
     c3->antiColParticle()));
     c1->antiColParticle()->addChild(newcluster2);
     c2->antiColParticle()->addChild(newcluster2);
     c3->antiColParticle()->addChild(newcluster2);
     newcluster2->setVertex(LorentzPoint());
 }
 
 pair <ClusterPtr,ClusterPtr>
 ColourReconnector::_reconnect(ClusterPtr &c1, ClusterPtr &c2) const {
 
   // choose the other possibility to form two clusters from the given
   // constituents
 
   assert(c1->numComponents()==2);
   assert(c2->numComponents()==2);
   int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1);
   for(unsigned int ix=0;ix<2;++ix) {
     if     (c1->particle(ix)->hasColour(false)) c1_col  = ix;
     else if(c1->particle(ix)->hasColour(true )) c1_anti = ix;
     if     (c2->particle(ix)->hasColour(false)) c2_col  = ix;
     else if(c2->particle(ix)->hasColour(true )) c2_anti = ix;
   }
   assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0);
 
   ClusterPtr newCluster1
     = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) );
 
   newCluster1->setVertex(0.5*( c1->colParticle()->vertex() + 
                                c2->antiColParticle()->vertex() ));
 
   if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true);
   if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true);
 
   ClusterPtr newCluster2
     = new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) );
 
   newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + 
                                c1->antiColParticle()->vertex() ));
 
   if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true);
   if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true);
 
   return pair <ClusterPtr,ClusterPtr> (newCluster1, newCluster2);
 }
 
 
 
 
 
 pair <ClusterPtr,ClusterPtr>
 ColourReconnector::_reconnectBaryonic(ClusterPtr &c1, ClusterPtr &c2) const {
 
   // choose the other possibility to form two clusters from the given
   // constituents
 
   assert(c1->numComponents()==2);
   assert(c2->numComponents()==2);
   int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1);
   for(unsigned int ix=0;ix<2;++ix) {
     if     (c1->particle(ix)->hasColour(false)) c1_col  = ix;
     else if(c1->particle(ix)->hasColour(true )) c1_anti = ix;
     if     (c2->particle(ix)->hasColour(false)) c2_col  = ix;
     else if(c2->particle(ix)->hasColour(true )) c2_anti = ix;
   }
   assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0);
 
 c1->colParticle()->abandonChild(c1);
 c2->antiColParticle()->abandonChild(c2);
 
   ClusterPtr newCluster1
     = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) );
 
   c1->colParticle()->addChild(newCluster1);
   c2->antiColParticle()->addChild(newCluster1);
 
   newCluster1->setVertex(0.5*( c1->colParticle()->vertex() + 
              c2->antiColParticle()->vertex() ));
 
   if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true);
   if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true);
 
   c1->antiColParticle()->abandonChild(c1);
   c2->colParticle()->abandonChild(c2);
 
   ClusterPtr newCluster2
     = new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) );
 
     c1->antiColParticle()->addChild(newCluster2);
     c2->colParticle()->addChild(newCluster2);
 
   newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + 
              c1->antiColParticle()->vertex() ));
 
   if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true);
   if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true);
 
   return pair <ClusterPtr,ClusterPtr> (newCluster1, newCluster2);
 }
 
 
 pair <int,int> ColourReconnector::_shuffle
   (const PVector & q, const PVector & aq, unsigned maxtries) const {
 
   const size_t nclusters = q.size();
   assert (nclusters > 1);
   assert (aq.size() == nclusters);
 
   int i, j;
   unsigned tries = 0;
   bool octet;
 
   do {
     // find two different random integers in the range [0, nclusters)
     i = UseRandom::irnd( nclusters );
     do { j = UseRandom::irnd( nclusters ); } while (i == j);
 
     // check if one of the two potential clusters would be a colour octet state
     octet = _isColour8( q[i], aq[j] ) || _isColour8( q[j], aq[i] ) ;
     tries++;
   } while (octet && tries < maxtries);
 
   if (octet) i = j = -1;
   return make_pair(i,j);
 }
 
 
 
 bool ColourReconnector::_isColour8(tcPPtr p, tcPPtr q) const {
   bool octet = false;
 
   // make sure we have a triplet and an anti-triplet
   if ( ( p->hasColour() && q->hasAntiColour() ) ||
        ( p->hasAntiColour() && q->hasColour() ) ) {
 
     // true if p and q are originated from a colour octet
     if ( !p->parents().empty() && !q->parents().empty() ) {
         octet = ( p->parents()[0] == q->parents()[0] ) &&
           ( p->parents()[0]->data().iColour() == PDT::Colour8 );
     }
 
     // (Final) option: check if same colour8 parent
     // or already found an octet.
     if(_octetOption==0||octet) return octet;
 
     // (All) option handling more octets
     // by browsing particle history/colour lines.
     tColinePtr cline,aline;
 
     // Get colourlines form final states.
     if(p->hasColour() && q->hasAntiColour()) {
       cline  = p->    colourLine();
       aline  = q->antiColourLine();
     }
     else {
       cline  = q->    colourLine();
       aline  = p->antiColourLine();
     }
     
     // Follow the colourline of p.
     if ( !p->parents().empty() ) {
       tPPtr parent = p->parents()[0];
       while (parent) {
         if(parent->data().iColour() == PDT::Colour8) {
             // Coulour8 particles should have a colour 
             // and an anticolour line. Currently the 
             // remnant has none of those. Since the children 
             // of the remnant are not allowed to emit currently, 
             // the colour octet remnant is handled by the return 
             // statement above. The assert also catches other 
             // colour octets without clines. If the children of 
             // a remnant should be allowed to emit, the remnant 
             // should get appropriate colour lines and 
             // colour states.
           //  See Ticket: #407 
           //  assert(parent->colourLine()&&parent->antiColourLine());
           octet = (parent->    colourLine()==cline &&
                    parent->antiColourLine()==aline);
         }
         if(octet||parent->parents().empty()) break;
         parent = parent->parents()[0];
       }
     }
   }
 
   return octet;
 }
 
 
 void ColourReconnector::persistentOutput(PersistentOStream & os) const {
   os << _clreco << _preco << _precoBaryonic << _algorithm << _initTemp << _annealingFactor
      << _annealingSteps << _triesPerStepFactor << ounit(_maxDistance,femtometer)
      << _octetOption;
 }
 
 void ColourReconnector::persistentInput(PersistentIStream & is, int) {
   is >> _clreco >> _preco >> _precoBaryonic >> _algorithm >> _initTemp >> _annealingFactor
      >> _annealingSteps >> _triesPerStepFactor >> iunit(_maxDistance,femtometer)
      >> _octetOption;
 }
 
 
 void ColourReconnector::Init() {
 
   static ClassDocumentation<ColourReconnector> documentation
     ("This class is responsible of the colour reconnection.");
 
 
   static Switch<ColourReconnector,int> interfaceColourReconnection
     ("ColourReconnection",
      "Colour reconnections",
      &ColourReconnector::_clreco, 0, true, false);
   static SwitchOption interfaceColourReconnectionNo
     (interfaceColourReconnection,
      "No",
      "Colour reconnections off",
      0);
   static SwitchOption interfaceColourReconnectionYes
     (interfaceColourReconnection,
      "Yes",
      "Colour reconnections on",
      1);
   
 
   static Parameter<ColourReconnector,double> interfaceMtrpAnnealingFactor
     ("AnnealingFactor",
      "The annealing factor is the ratio of the temperatures in two successive "
      "temperature steps.",
      &ColourReconnector::_annealingFactor, 0.9, 0.0, 1.0,
      false, false, Interface::limited);
 
   static Parameter<ColourReconnector,unsigned> interfaceMtrpAnnealingSteps
     ("AnnealingSteps",
      "Number of temperature steps in the statistical annealing algorithm",
      &ColourReconnector::_annealingSteps, 50, 1, 10000,
      false, false, Interface::limited);
 
   static Parameter<ColourReconnector,double> interfaceMtrpTriesPerStepFactor
     ("TriesPerStepFactor",
      "The number of reconnection tries per temperature steps is the number of "
      "clusters times this factor.",
      &ColourReconnector::_triesPerStepFactor, 5.0, 0.0, 100.0,
      false, false, Interface::limited);
 
   
   static Parameter<ColourReconnector,double> interfaceMtrpInitialTemp
     ("InitialTemperature",
      "Factor used to determine the initial temperature from the median of the "
      "energy change in a few random rearrangements.",
      &ColourReconnector::_initTemp, 0.1, 0.00001, 100.0,
      false, false, Interface::limited);
 
 
   static Parameter<ColourReconnector,double> interfaceRecoProb 
     ("ReconnectionProbability",
      "Probability that a found reconnection possibility is actually accepted",
      &ColourReconnector::_preco, 0.5, 0.0, 1.0,
      false, false, Interface::limited);
 
   static Parameter<ColourReconnector,double> interfaceRecoProbBaryonic
     ("ReconnectionProbabilityBaryonic",
      "Probability that a found reconnection possibility is actually accepted",
      &ColourReconnector::_precoBaryonic, 0.5, 0.0, 1.0,
      false, false, Interface::limited);
 
 
   static Switch<ColourReconnector,int> interfaceAlgorithm
     ("Algorithm",
      "Specifies the colour reconnection algorithm",
      &ColourReconnector::_algorithm, 0, true, false);
   static SwitchOption interfaceAlgorithmPlain
     (interfaceAlgorithm,
      "Plain",
      "Plain colour reconnection as in Herwig 2.5.0",
      0);
   static SwitchOption interfaceAlgorithmStatistical
     (interfaceAlgorithm,
      "Statistical",
      "Statistical colour reconnection using simulated annealing",
      1);
   static SwitchOption interfaceAlgorithmBaryonic
     (interfaceAlgorithm,
-     "BaryonicReco",
+     "Baryonic",
      "Baryonic cluster reconnection",
      2);
 
   static Parameter<ColourReconnector,Length> interfaceMaxDistance
     ("MaxDistance",
      "Maximum distance between the clusters at which to consider rearrangement"
      " to avoid colour reconneections of displaced vertices",
      &ColourReconnector::_maxDistance, femtometer, 1000.*femtometer, 0.0*femtometer, 1e100*femtometer,
      false, false, Interface::limited);
 
 
   static Switch<ColourReconnector,unsigned int> interfaceOctetTreatment
     ("OctetTreatment",
      "Which octets are not allowed to be reconnected",
      &ColourReconnector::_octetOption, 0, false, false);
   static SwitchOption interfaceOctetTreatmentFinal
     (interfaceOctetTreatment,
      "Final",
      "Only prevent for the final (usuaslly non-perturbative) g -> q qbar splitting",
      0);
   static SwitchOption interfaceOctetTreatmentAll
     (interfaceOctetTreatment,
      "All",
      "Prevent for all octets",
      1);
 
 }
 
diff --git a/MatrixElement/Hadron/MEMinBias.cc b/MatrixElement/Hadron/MEMinBias.cc
--- a/MatrixElement/Hadron/MEMinBias.cc
+++ b/MatrixElement/Hadron/MEMinBias.cc
@@ -1,167 +1,174 @@
 // -*- C++ -*-
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the MEMinBias class.
 //
 
 #include "MEMinBias.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/SimplePhaseSpace.h"
 //#include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/Handlers/StandardXComb.h"
 #include "ThePEG/Interface/Parameter.h"
 
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 #include "ThePEG/PDT/EnumParticles.h"
 #include "ThePEG/MatrixElement/Tree2toNDiagram.h"
 
 void MEMinBias::getDiagrams() const {
   int maxflav(2);
   // Pomeron data
   tcPDPtr pom = getParticleData(990);
 
   for ( int i = 1; i <= maxflav; ++i ) {
     for( int j=1; j <= i; ++j){
       tcPDPtr q1 = getParticleData(i);
       tcPDPtr q1b = q1->CC();
       tcPDPtr q2 = getParticleData(j);
       tcPDPtr q2b = q2->CC();
 
       // For each flavour we add:
       //qq -> qq
       add(new_ptr((Tree2toNDiagram(3), q1, pom, q2, 1, q1, 2, q2, -1)));
       //qqb -> qqb
       add(new_ptr((Tree2toNDiagram(3), q1, pom, q2b, 1, q1, 2, q2b, -2)));
       //qbqb -> qbqb
       add(new_ptr((Tree2toNDiagram(3), q1b, pom, q2b, 1, q1b, 2, q2b, -3)));
     }
   }
 }
 
 Energy2 MEMinBias::scale() const {
-  return sqr(10*GeV);
+  return sqr(Scale_);
 }
 
 int MEMinBias::nDim() const {
   return 0;
 }
 
 void MEMinBias::setKinematics() {
   HwMEBase::setKinematics(); // Always call the base class method first.
 }
 
 bool MEMinBias::generateKinematics(const double *) {
   // generate the masses of the particles
   for ( int i = 2, N = meMomenta().size(); i < N; ++i ) {
     meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->generateMass());
   }
 
   Energy q = ZERO;
   try {
     q = SimplePhaseSpace::
       getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass());
   } catch ( ImpossibleKinematics & e ) {
     return false;
   }
 
   Energy pt = ZERO;
   meMomenta()[2].setVect(Momentum3( pt,  pt, q));
   meMomenta()[3].setVect(Momentum3(-pt, -pt, -q));
 
   meMomenta()[2].rescaleEnergy();
   meMomenta()[3].rescaleEnergy();
 
   jacobian(1.0);
   return true;
 }
 
 double MEMinBias::me2() const {
   //tuned so it gives the correct normalization for xmin = 0.11
   return csNorm_*(sqr(generator()->maximumCMEnergy())/GeV2);
 }
 
 CrossSection MEMinBias::dSigHatDR() const {
   return me2()*jacobian()/sHat()*sqr(hbarc);
 }
 
 unsigned int MEMinBias::orderInAlphaS() const {
   return 2;
 }
 
 unsigned int MEMinBias::orderInAlphaEW() const {
   return 0;
 }
 
 Selector<MEBase::DiagramIndex>
 MEMinBias::diagrams(const DiagramVector & diags) const {
   Selector<DiagramIndex> sel;
   for ( DiagramIndex i = 0; i < diags.size(); ++i ) 
     sel.insert(1.0, i);
 
   return sel;
 }
 
 Selector<const ColourLines *>
 MEMinBias::colourGeometries(tcDiagPtr diag) const {
 
   static ColourLines qq("1 4, 3 5");
   static ColourLines qqb("1 4, -3 -5");
   static ColourLines qbqb("-1 -4, -3 -5");
 
   Selector<const ColourLines *> sel;
   
   switch(diag->id()){
   case -1:
     sel.insert(1.0, &qq);
     break;
   case -2:
     sel.insert(1.0, &qqb);
     break;
   case -3:
     sel.insert(1.0, &qbqb);
     break;
   }
   return sel;
 }
 
 
 IBPtr MEMinBias::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr MEMinBias::fullclone() const {
   return new_ptr(*this);
 }
 
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<MEMinBias,HwMEBase>
 describeHerwigMEMinBias("Herwig::MEMinBias", "HwMEHadron.so");
 
 void MEMinBias::persistentOutput(PersistentOStream & os) const {
   os << csNorm_;
 }
 
 void MEMinBias::persistentInput(PersistentIStream & is, int) {
   is >> csNorm_;
 }
 
 void MEMinBias::Init() {
 
   static ClassDocumentation<MEMinBias> documentation
     ("There is no documentation for the MEMinBias class");
      
   static Parameter<MEMinBias,double> interfacecsNorm
     ("csNorm",
      "Normalization of the min-bias cross section.",
      &MEMinBias::csNorm_, 
      1.0, 0.0, 100.0, 
      false, false, Interface::limited);
+  static Parameter<MEMinBias,Energy> interfaceScale
+    ("Scale",
+     "Scale for the Min Bias matrix element.",
+     &MEMinBias::Scale_,GeV,
+     1.0*GeV, 0.0*GeV, 100.0*GeV,
+     false, false, Interface::limited);
+
 }
 
diff --git a/MatrixElement/Hadron/MEMinBias.h b/MatrixElement/Hadron/MEMinBias.h
--- a/MatrixElement/Hadron/MEMinBias.h
+++ b/MatrixElement/Hadron/MEMinBias.h
@@ -1,190 +1,195 @@
 // -*- C++ -*-
 #ifndef HERWIG_MEMinBias_H
 #define HERWIG_MEMinBias_H
 //
 // This is the declaration of the MEMinBias class.
 //
 
 #include "Herwig/MatrixElement/HwMEBase.h"
 
 namespace Herwig {
 
 using namespace ThePEG;
 
 /**
  * The MEMinBias class provides a simple colour singlet exchange matrix element
  * to be used in the soft component of the multiple scattering model of the 
  * underlying event
  *
  * @see \ref MEMinBiasInterfaces "The interfaces"
  * defined for MEMinBias.
  */
 class MEMinBias: public HwMEBase {
 
 public:
 
   /**
    * The default constructor.
    */
-  MEMinBias() : csNorm_(1.) {}
+  MEMinBias() : csNorm_(1.), Scale_(2.*GeV) {}
 
 public:
 
   /** @name Virtual functions required by the MEBase class. */
   //@{
   /**
    * Return the order in \f$\alpha_S\f$ in which this matrix
    * element is given.
    */
   virtual unsigned int orderInAlphaS() const;
 
   /**
    * Return the order in \f$\alpha_{EW}\f$ in which this matrix
    * element is given.
    */
   virtual unsigned int orderInAlphaEW() const;
 
   /**
    * The matrix element for the kinematical configuration
    * previously provided by the last call to setKinematics(), suitably
    * scaled by sHat() to give a dimension-less number.
    * @return the matrix element scaled with sHat() to give a
    * dimensionless number.
    */
   virtual double me2() const;
 
   /**
    * Return the scale associated with the last set phase space point.
    */
   virtual Energy2 scale() const;
 
   /**
    * Set the typed and momenta of the incoming and outgoing partons to
    * be used in subsequent calls to me() and colourGeometries()
    * according to the associated XComb object. If the function is
    * overridden in a sub class the new function must call the base
    * class one first.
    */
   virtual void setKinematics();
 
   /**
    * The number of internal degrees of freedom used in the matrix
    * element.
    */
   virtual int nDim() const;
 
   /**
    * Generate internal degrees of freedom given nDim() uniform
    * random numbers in the interval \f$ ]0,1[ \f$. To help the phase space
    * generator, the dSigHatDR should be a smooth function of these
    * numbers, although this is not strictly necessary.
    * @param r a pointer to the first of nDim() consecutive random numbers.
    * @return true if the generation succeeded, otherwise false.
    */
   virtual bool generateKinematics(const double * r);
 
   /**
    * Return the matrix element squared differential in the variables
    * given by the last call to generateKinematics().
    */
   virtual CrossSection dSigHatDR() const;
 
   /**
    * Add all possible diagrams with the add() function.
    */
   virtual void getDiagrams() const;
 
   /**
    * Get diagram selector. With the information previously supplied with the
    * setKinematics method, a derived class may optionally
    * override this method to weight the given diagrams with their
    * (although certainly not physical) relative probabilities.
    * @param dv the diagrams to be weighted.
    * @return a Selector relating the given diagrams to their weights.
    */
   virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
 
   /**
    * Return a Selector with possible colour geometries for the selected
    * diagram weighted by their relative probabilities.
    * @param diag the diagram chosen.
    * @return the possible colour geometries weighted by their
    * relative probabilities.
    */
   virtual Selector<const ColourLines *>
   colourGeometries(tcDiagPtr diag) 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.
    */
 
 
   /**
    * 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).
 
 
 private:
   /**
   *  Normalization of the min-bias cross section
   */
   double csNorm_;	
 
   /**
+   * Scale for the Min Bias matrix element
+   */
+  Energy Scale_;
+
+  /**
    * The assignment operator is private and must never be called.
    * In fact, it should not even be implemented.
    */
   MEMinBias & operator=(const MEMinBias &) = delete;
 
 };
 
 }
 
 #endif /* HERWIG_MEMinBias_H */
diff --git a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
--- a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
+++ b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
@@ -1,1658 +1,1675 @@
 // -*- C++ -*-
 //
 // MatchboxMEBase.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 MatchboxMEBase class.
 //
 
 #include "MatchboxMEBase.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/PDF/PDF.h"
 #include "ThePEG/PDT/PDT.h"
 #include "ThePEG/StandardModel/StandardModelBase.h"
 #include "ThePEG/Cuts/Cuts.h"
 #include "ThePEG/Handlers/StdXCombGroup.h"
 #include "ThePEG/EventRecord/SubProcess.h"
 #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
 #include "Herwig/MatrixElement/Matchbox/Utility/DiagramDrawer.h"
 #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
 #include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h"
 #include "Herwig/API/RunDirectories.h"
 #include "Herwig/MatrixElement/ProductionMatrixElement.h"
 #include "Herwig/MatrixElement/HardVertex.h"
 
 #include <cctype>
 
 #include <iterator>
 using std::ostream_iterator;
 
 using namespace Herwig;
 
 MatchboxMEBase::MatchboxMEBase() 
   : MEBase(), 
     theOneLoop(false),
     theOneLoopNoBorn(false),
     theOneLoopNoLoops(false),
     theNoCorrelations(false),
     theHavePDFs(false,false), checkedPDFs(false) {}
 
 MatchboxMEBase::~MatchboxMEBase() {}
 
 Ptr<MatchboxFactory>::tptr MatchboxMEBase::factory() const {
   return MatchboxFactory::currentFactory();
 }
 
 Ptr<Tree2toNGenerator>::tptr MatchboxMEBase::diagramGenerator() const { return factory()->diagramGenerator(); }
 
 Ptr<ProcessData>::tptr MatchboxMEBase::processData() const { return factory()->processData(); }
 
 unsigned int MatchboxMEBase::getNLight() const { return factory()->nLight(); }
 
 vector<long> MatchboxMEBase::getNLightJetVec() const { return factory()->nLightJetVec(); }
 
 vector<long> MatchboxMEBase::getNHeavyJetVec() const { return factory()->nHeavyJetVec(); }
 
 vector<long> MatchboxMEBase::getNLightProtonVec() const { return factory()->nLightProtonVec(); }
 
 double MatchboxMEBase::factorizationScaleFactor() const { return factory()->factorizationScaleFactor(); }
 
 double MatchboxMEBase::renormalizationScaleFactor() const { return factory()->renormalizationScaleFactor(); }
 
 bool MatchboxMEBase::fixedCouplings() const { return factory()->fixedCouplings(); }
 
 bool MatchboxMEBase::fixedQEDCouplings() const { return factory()->fixedQEDCouplings(); }
 
 bool MatchboxMEBase::checkPoles() const { return factory()->checkPoles(); }
 
 bool MatchboxMEBase::verbose() const { return factory()->verbose(); }
 
 bool MatchboxMEBase::initVerbose() const { return factory()->initVerbose(); }
 
 void MatchboxMEBase::getDiagrams() const {
 
   if ( diagramGenerator() && processData() ) {
 
     vector<Ptr<Tree2toNDiagram>::ptr> diags;
 
     vector<Ptr<Tree2toNDiagram>::ptr>& res =
       processData()->diagramMap()[subProcess().legs];
     if ( res.empty() ) {
       res = diagramGenerator()->generate(subProcess().legs,orderInAlphaS(),orderInAlphaEW());
     }
     copy(res.begin(),res.end(),back_inserter(diags));
     processData()->fillMassGenerators(subProcess().legs);
 
     if ( diags.empty() )
       return;
 
     for (auto const & d :  diags )
       add(d);
     
     return;
 
   }
 
   throw Exception()
     << "MatchboxMEBase::getDiagrams() expects a Tree2toNGenerator and ProcessData object.\n"
     << "Please check your setup." << Exception::runerror;
 
 }
 
 Selector<MEBase::DiagramIndex> 
 MatchboxMEBase::diagrams(const DiagramVector & diags) const {
 
   if ( phasespace() ) {
     return phasespace()->selectDiagrams(diags);
   }
 
   throw Exception()
     << "MatchboxMEBase::diagrams() expects a MatchboxPhasespace object.\n"
     << "Please check your setup." << Exception::runerror;
   return Selector<MEBase::DiagramIndex>();
 
 }
 
 Selector<const ColourLines *>
 MatchboxMEBase::colourGeometries(tcDiagPtr diag) const {
 
   if ( matchboxAmplitude() ) {
     if ( matchboxAmplitude()->haveColourFlows() ) {
       if ( matchboxAmplitude()->treeAmplitudes() )
 	matchboxAmplitude()->prepareAmplitudes(this);
       return matchboxAmplitude()->colourGeometries(diag);
     }
   }
 
   Ptr<Tree2toNDiagram>::tcptr tdiag =
     dynamic_ptr_cast<Ptr<Tree2toNDiagram>::tcptr>(diag);
   assert(diag && processData());
 
   vector<ColourLines*>& flows = processData()->colourFlowMap()[tdiag];
 
   if ( flows.empty() ) {
 
     list<list<list<pair<int,bool> > > > cflows =
       ColourBasis::colourFlows(tdiag);
 
     for ( auto const & fit : cflows)
       flows.push_back(new ColourLines(ColourBasis::cfstring(fit)));
 
   }
 
   Selector<const ColourLines *> res;
   for ( auto const & f : flows ) res.insert(1.0,f);
 
   return res;
 
 }
 
 void MatchboxMEBase::constructVertex(tSubProPtr sub, const ColourLines* cl) {
 
   if ( !canFillRhoMatrix() || !factory()->spinCorrelations() )
     return;
 
   assert(matchboxAmplitude());
   assert(matchboxAmplitude()->colourBasis());
 
   // get the colour structure for the selected colour flow
   size_t cStructure = 
     matchboxAmplitude()->colourBasis()->tensorIdFromFlow(lastXComb().lastDiagram(),cl);
 
   // hard process for processing the spin info
   tPVector hard;
   hard.push_back(sub->incoming().first);
   hard.push_back(sub->incoming().second);
   
   vector<PDT::Spin> out;
   for ( auto const & p  : sub->outgoing() ) {
     out.push_back(p->data().iSpin());
     hard.push_back(p);
   }
 
   // calculate dummy wave functions to fill the spin info
   static vector<VectorWaveFunction> dummyPolarizations;
   static vector<SpinorWaveFunction> dummySpinors;
   static vector<SpinorBarWaveFunction> dummyBarSpinors;
   for ( size_t k = 0; k < hard.size(); ++k ) {
     if ( hard[k]->data().iSpin() == PDT::Spin1Half ) {
       if ( hard[k]->id() > 0 && k > 1 ) {
 	SpinorBarWaveFunction(dummyBarSpinors,hard[k],
 			      outgoing, true);
       } else if ( hard[k]->id() < 0 && k > 1 ) {
 	SpinorWaveFunction(dummySpinors,hard[k],
 			   outgoing, true);
       } else if ( hard[k]->id() > 0 && k < 2 ) {
 	SpinorWaveFunction(dummySpinors,hard[k],
 			   incoming, false);
       } else if ( hard[k]->id() < 0 && k < 2 ) {
 	SpinorBarWaveFunction(dummyBarSpinors,hard[k],
 			      incoming, false);
       }
     } 
     else if ( hard[k]->data().iSpin() == PDT::Spin1 ) {
       VectorWaveFunction(dummyPolarizations,hard[k],
 			 k > 1 ? outgoing : incoming,
 			 k > 1 ? true : false,
 			 hard[k]->data().hardProcessMass() == ZERO);
     }
     else if (hard[k]->data().iSpin() == PDT::Spin0 ) {
       ScalarWaveFunction(hard[k],k > 1 ? outgoing : incoming,
 			 k > 1 ? true : false);
     }
     else
       assert(false);
   }
 
   // fill the production matrix element
   ProductionMatrixElement pMe(mePartonData()[0]->iSpin(),
 			      mePartonData()[1]->iSpin(),
 			      out);
   for ( map<vector<int>,CVector>::const_iterator lamp = lastLargeNAmplitudes().begin();
 	lamp != lastLargeNAmplitudes().end(); ++lamp ) {
     vector<unsigned int> pMeHelicities
       = matchboxAmplitude()->physicalHelicities(lamp->first);
     pMe(pMeHelicities) = lamp->second[cStructure];
   }
 
   // set the spin information
   HardVertexPtr hardvertex = new_ptr(HardVertex());
   hardvertex->ME(pMe);
   if ( sub->incoming().first->spinInfo() )
     sub->incoming().first->spinInfo()->productionVertex(hardvertex);
   if ( sub->incoming().second->spinInfo() )
     sub->incoming().second->spinInfo()->productionVertex(hardvertex);
   for ( auto const & p : sub->outgoing() )
     if ( p->spinInfo() )
       p->spinInfo()->productionVertex(hardvertex);
   
 
 }
 
 unsigned int MatchboxMEBase::orderInAlphaS() const {
   return subProcess().orderInAlphaS;
 }
 
 unsigned int MatchboxMEBase::orderInAlphaEW() const {
   return subProcess().orderInAlphaEW;
 }
 
 void MatchboxMEBase::setXComb(tStdXCombPtr xc) {
 
   MEBase::setXComb(xc);
   lastMatchboxXComb(xc);
   if ( phasespace() )
     phasespace()->setXComb(xc);
   if ( scaleChoice() )
     scaleChoice()->setXComb(xc);
   if ( matchboxAmplitude() )
     matchboxAmplitude()->setXComb(xc);
   if (theMerger){
     theMerger->setME(this);
     theMerger->setXComb( xc );
   }
 
 }
 
 double MatchboxMEBase::generateIncomingPartons(const double* r1, const double* r2) {
   
   // shamelessly stolen from PartonExtractor.cc
 
   Energy2 shmax = lastCuts().sHatMax();
   Energy2 shmin = lastCuts().sHatMin();
   Energy2 sh = shmin*pow(shmax/shmin, *r1);
   double ymax = lastCuts().yHatMax();
   double ymin = lastCuts().yHatMin();
   double km = log(shmax/shmin);
   ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/sh)));
   ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/sh)));
 
   double y = ymin + (*r2)*(ymax - ymin);
   double x1 = exp(-0.5*log(lastS()/sh) + y);
   double x2 = exp(-0.5*log(lastS()/sh) - y);
 
   Lorentz5Momentum P1 = lastParticles().first->momentum();
   LorentzMomentum p1 = lightCone((P1.rho() + P1.e())*x1, Energy());
   p1.rotateY(P1.theta());
   p1.rotateZ(P1.phi());
   meMomenta()[0] = p1;
 
   Lorentz5Momentum P2 = lastParticles().second->momentum();
   LorentzMomentum p2 = lightCone((P2.rho() + P2.e())*x2, Energy());
   p2.rotateY(P2.theta());
   p2.rotateZ(P2.phi());
   meMomenta()[1] = p2;
 
   lastXCombPtr()->lastX1X2(make_pair(x1,x2));
   lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
 
   return km*(ymax - ymin);
 
 }
 
 bool MatchboxMEBase::generateKinematics(const double * r) {
 
   if ( phasespace() ) {
 
     jacobian(phasespace()->generateKinematics(r,meMomenta()));
     if ( jacobian() == 0.0 )
       return false;
 
     setScale();
     if (theMerger&&!theMerger->generateKinematics(r)){
       return false;
     }
 
     logGenerateKinematics(r);
 
     assert(lastMatchboxXComb());
 
     if ( nDimAmplitude() > 0 ) {
       amplitudeRandomNumbers().resize(nDimAmplitude());
       copy(r + nDimPhasespace(), 
 	   r + nDimPhasespace() + nDimAmplitude(),
 	   amplitudeRandomNumbers().begin());
     }
 
     if ( nDimInsertions() > 0 ) {
       insertionRandomNumbers().resize(nDimInsertions());
       copy(r + nDimPhasespace() + nDimAmplitude(), 
 	   r + nDimPhasespace() + nDimAmplitude() + nDimInsertions(),
 	   insertionRandomNumbers().begin());
     }
 
     return true;
 
   }
 
   throw Exception()
     << "MatchboxMEBase::generateKinematics() expects a MatchboxPhasespace object.\n"
     << "Please check your setup." << Exception::runerror;
 
   return false;
 
 }
 
 int MatchboxMEBase::nDim() const { 
 
   if ( lastMatchboxXComb() )
     return nDimPhasespace() + nDimAmplitude() + nDimInsertions();
 
   int ampAdd = 0;
   if ( matchboxAmplitude() ) {
     ampAdd = matchboxAmplitude()->nDimAdditional();
   }
 
   int insertionAdd = 0;
   for ( auto const & v : virtuals() ) {
     insertionAdd = max(insertionAdd,v->nDimAdditional());
   }
 
   return nDimBorn() + ampAdd + insertionAdd;
 
 }
 
 int MatchboxMEBase::nDimBorn() const { 
 
   if ( lastMatchboxXComb() )
     return nDimPhasespace();
 
   if ( phasespace() )
     return phasespace()->nDim(diagrams().front()->partons());
 
   throw Exception()
     << "MatchboxMEBase::nDim() expects a MatchboxPhasespace object.\n"
     << "Please check your setup." << Exception::runerror;
 
   return 0;
 
 }
 
 void MatchboxMEBase::setScale(Energy2 ren, Energy2 fac) const {
   if ( haveX1X2() ) {
     lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
   }
   Energy2 fcscale = (fac == ZERO) ? factorizationScale() : fac;
   Energy2 fscale = fcscale*sqr(factorizationScaleFactor());
   Energy2 rscale = (ren == ZERO ? renormalizationScale() : ren)*sqr(renormalizationScaleFactor());
   Energy2 ewrscale = renormalizationScaleQED();
   lastXCombPtr()->lastScale(fscale);
   lastXCombPtr()->lastCentralScale(fcscale);
   lastXCombPtr()->lastShowerScale(showerScale());
   lastMatchboxXComb()->lastRenormalizationScale(rscale);
   if ( !fixedCouplings() ) {
     if ( rscale > lastCuts().scaleMin() )
       lastXCombPtr()->lastAlphaS(SM().alphaS(rscale));
     else
       lastXCombPtr()->lastAlphaS(SM().alphaS(lastCuts().scaleMin()));
   } else {
     lastXCombPtr()->lastAlphaS(SM().alphaS());
   }
   if ( !fixedQEDCouplings() ) {
     lastXCombPtr()->lastAlphaEM(SM().alphaEMME(ewrscale));
   } else {
     lastXCombPtr()->lastAlphaEM(SM().alphaEMMZ());
   }
   logSetScale();
 }
 
 Energy2 MatchboxMEBase::factorizationScale() const {
   if ( scaleChoice() ) {
     return scaleChoice()->factorizationScale();
   }
 
   throw Exception()
     << "MatchboxMEBase::factorizationScale() expects a MatchboxScaleChoice object.\n"
     << "Please check your setup." << Exception::runerror;
 
   return ZERO;
 
 }
 
 Energy2 MatchboxMEBase::renormalizationScale() const {
   if ( scaleChoice() ) {
     return scaleChoice()->renormalizationScale();
   }
 
   throw Exception()
     << "MatchboxMEBase::renormalizationScale() expects a MatchboxScaleChoice object.\n"
     << "Please check your setup." << Exception::runerror;
 
   return ZERO;
 
 }
 
 Energy2 MatchboxMEBase::renormalizationScaleQED() const {
   if ( scaleChoice() ) {
     return scaleChoice()->renormalizationScaleQED();
   }
   return renormalizationScale();
 }
 
 Energy2 MatchboxMEBase::showerScale() const {
   if ( scaleChoice() ) {
     return scaleChoice()->showerScale();
   }
 
   throw Exception()
     << "MatchboxMEBase::showerScale() expects a MatchboxScaleChoice object.\n"
     << "Please check your setup." << Exception::runerror;
 
   return ZERO;
 
 }
 
 void MatchboxMEBase::setVetoScales(tSubProPtr) const {}
 
 bool MatchboxMEBase::havePDFWeight1() const {
   if ( checkedPDFs )
     return theHavePDFs.first;
   theHavePDFs.first = 
     factory()->isIncoming(mePartonData()[0]) && 
     lastXCombPtr()->partonBins().first->pdf();
   theHavePDFs.second = 
     factory()->isIncoming(mePartonData()[1]) &&
     lastXCombPtr()->partonBins().second->pdf();
   checkedPDFs = true;
   return theHavePDFs.first;
 }
 
 bool MatchboxMEBase::havePDFWeight2() const {
   if ( checkedPDFs )
     return theHavePDFs.second;
   theHavePDFs.first = 
     factory()->isIncoming(mePartonData()[0]) && 
     lastXCombPtr()->partonBins().first->pdf();
   theHavePDFs.second = 
     factory()->isIncoming(mePartonData()[1]) &&
     lastXCombPtr()->partonBins().second->pdf();
   checkedPDFs = true;
   return theHavePDFs.second;
 }
 
 void MatchboxMEBase::getPDFWeight(Energy2 factorizationScale) const {
 
   if ( !havePDFWeight1() && !havePDFWeight2() ) {
     lastMEPDFWeight(1.0);
     logPDFWeight();
     return;
   }
 
   double w = 1.;
 
   if ( havePDFWeight1() )
     w *= pdf1(factorizationScale);
 
   if ( havePDFWeight2() )
     w *= pdf2(factorizationScale);
 
   lastMEPDFWeight(w);
 
   logPDFWeight();
 
 }
 
 double MatchboxMEBase::pdf1(Energy2 fscale, double xEx, double xFactor) const {
 
   assert(lastXCombPtr()->partonBins().first->pdf());
 
   if ( xEx < 1. && lastX1()*xFactor >= xEx ) {
     return
       ( ( 1. - lastX1()*xFactor ) / ( 1. - xEx ) ) *
       lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
 						     lastPartons().first->dataPtr(),
 						     fscale == ZERO ? lastScale() : fscale,
 						     xEx)/xEx;
   }
 
   return lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
 							lastPartons().first->dataPtr(),
 							fscale == ZERO ? lastScale() : fscale,
 							lastX1()*xFactor)/lastX1()/xFactor;
 }
 
 double MatchboxMEBase::pdf2(Energy2 fscale, double xEx, double xFactor) const {
 
   assert(lastXCombPtr()->partonBins().second->pdf());
 
   if ( xEx < 1. && lastX2()*xFactor >= xEx ) {
     return
       ( ( 1. - lastX2()*xFactor ) / ( 1. - xEx ) ) *
       lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
 						      lastPartons().second->dataPtr(),
 						      fscale == ZERO ? lastScale() : fscale,
 						      xEx)/xEx;
   }
 
   return lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
 							 lastPartons().second->dataPtr(),
 							 fscale == ZERO ? lastScale() : fscale,
 							 lastX2()*xFactor)/lastX2()/xFactor;
 
 }
 
 double MatchboxMEBase::me2() const {
 
   if ( matchboxAmplitude() ) {
 
     if ( matchboxAmplitude()->treeAmplitudes() )
       matchboxAmplitude()->prepareAmplitudes(this);
 
     double res = 
       matchboxAmplitude()->me2()*
       me2Norm();
 
     return res;
 
   }
 
   throw Exception()
     << "MatchboxMEBase::me2() expects a MatchboxAmplitude object.\n"
     << "Please check your setup." << Exception::runerror;
   return 0.;
 
 }
 
 double MatchboxMEBase::largeNME2(Ptr<ColourBasis>::tptr largeNBasis) const {
 
   if ( matchboxAmplitude() ) {
 
     if ( matchboxAmplitude()->treeAmplitudes() ) {
       largeNBasis->prepare(mePartonData(),false);
       matchboxAmplitude()->prepareAmplitudes(this);
     }
 
     double res = 
       matchboxAmplitude()->largeNME2(largeNBasis)*
       me2Norm();
 
     return res;
 
   }
 
   throw Exception()
     << "MatchboxMEBase::largeNME2() expects a MatchboxAmplitude object.\n"
     << "Please check your setup." << Exception::runerror;
   return 0.;
 
 }
 
 double MatchboxMEBase::finalStateSymmetry() const {
 
   if ( symmetryFactor() > 0.0 )
     return symmetryFactor();
 
   double sFactor = 1.;
 
   map<long,int> counts;
 
   cPDVector checkData;
   copy(mePartonData().begin()+2,mePartonData().end(),back_inserter(checkData));
 
   cPDVector::iterator p = checkData.begin();
   while ( !checkData.empty() ) {
     if ( counts.find((**p).id()) != counts.end() ) {
       counts[(**p).id()] += 1;
     } else {
       counts[(**p).id()] = 1;
     }
     checkData.erase(p);
     p = checkData.begin();
     continue;
   }
 
   for ( auto const & c : counts) {
     if ( c.second == 1 )
       continue;
     if ( c.second == 2 )
       sFactor /= 2.;
     else if ( c.second == 3 )
       sFactor /= 6.;
     else if ( c.second == 4 )
       sFactor /= 24.;
   }
 
   symmetryFactor(sFactor);
 
   return symmetryFactor();
 
 }
 
 double MatchboxMEBase::me2Norm(unsigned int addAlphaS) const {
 
-  // assume that we always have incoming
-  // spin-1/2 or massless spin-1 particles
-  double fac = 1./4.;
+  double fac = 1.;
 
-  if ( hasInitialAverage() )
-    fac = 1.;
+  if ( !hasInitialAverage() ) {
+    for ( size_t k = 0; k < 2; ++k ) {
+      // Spin 0 does not involve any additional factors
+      if ( mePartonData()[k]->iSpin() == PDT::Spin1Half )
+	fac *= 1./2.;
+      else if ( mePartonData()[k]->iSpin() == PDT::Spin1 &&
+		mePartonData()[k]->hardProcessMass() == ZERO )
+	fac *= 1./2.;
+      else if ( mePartonData()[k]->iSpin() == PDT::Spin1 &&
+		mePartonData()[k]->hardProcessMass() > ZERO )
+	fac *= 1./3.;
+      else if ( mePartonData()[k]->iSpin() == PDT::Spin3Half )
+	fac *= 1./4.;
+      else if ( mePartonData()[k]->iSpin() == PDT::Spin2 &&
+		mePartonData()[k]->hardProcessMass() == ZERO )
+	fac *= 1./2.;
+      else if ( mePartonData()[k]->iSpin() == PDT::Spin2 &&
+		mePartonData()[k]->hardProcessMass() > ZERO )
+	fac *= 1./5.;
+    }
+  }  
 
   double couplings = 1.0;
   if ( (orderInAlphaS() > 0 || addAlphaS != 0) && !hasRunningAlphaS() ) {
     fac *= pow(lastAlphaS()/SM().alphaS(),double(orderInAlphaS()+addAlphaS));
     couplings *= pow(lastAlphaS(),double(orderInAlphaS()+addAlphaS));
   }
   if ( orderInAlphaEW() > 0 && !hasRunningAlphaEW() ) {
     fac *= pow(lastAlphaEM()/SM().alphaEMMZ(),double(orderInAlphaEW()));
     couplings *= pow(lastAlphaEM(),double(orderInAlphaEW()));
   }
   lastMECouplings(couplings);
 
   if ( !hasInitialAverage() ) {
     if ( mePartonData()[0]->iColour() == PDT::Colour3 || 
 	 mePartonData()[0]->iColour() == PDT::Colour3bar )
       fac /= SM().Nc();
     else if ( mePartonData()[0]->iColour() == PDT::Colour8 )
       fac /= (SM().Nc()*SM().Nc()-1.);
 
     if ( mePartonData()[1]->iColour() == PDT::Colour3 || 
 	 mePartonData()[1]->iColour() == PDT::Colour3bar )
       fac /= SM().Nc();
     else if ( mePartonData()[1]->iColour() == PDT::Colour8 )
       fac /= (SM().Nc()*SM().Nc()-1.);
   }
 
   return !hasFinalStateSymmetry() ? finalStateSymmetry()*fac : fac;
 
 }
 
 
 CrossSection MatchboxMEBase::prefactor()const{
   return (sqr(hbarc)/(2.*lastSHat())) *jacobian()* lastMEPDFWeight();
 }
 
 CrossSection MatchboxMEBase::dSigHatDRB() const {
   getPDFWeight();
   lastME2(me2());
   return oneLoopNoBorn()?ZERO:prefactor() * lastME2();
 }
 
 CrossSection MatchboxMEBase::dSigHatDRV() const {
   getPDFWeight();
   lastME2(me2());
   return ( oneLoop() && !oneLoopNoLoops() )?(prefactor() * oneLoopInterference()):ZERO;
 }
 
 CrossSection MatchboxMEBase::dSigHatDRI() const {
   getPDFWeight();
   lastME2(me2());
   CrossSection res=ZERO;
   if  (oneLoop() &&!onlyOneLoop())  {
     for ( auto const & v : virtuals()) {
       v->setXComb(lastXCombPtr());
       res += v->dSigHatDR();
     }
     if ( checkPoles() && oneLoop() )
       logPoles();
   }
   return res;
 }
 
 CrossSection MatchboxMEBase::dSigHatDRAlphaDiff(double alpha) const {
   getPDFWeight();
   lastME2(me2());
   CrossSection res=ZERO;
   for ( auto const & v: virtuals() ) {
     v->setXComb(lastXCombPtr());
     res+=v->dSigHatDRAlphaDiff( alpha);
   }
   return res;
 }
 
 CrossSection MatchboxMEBase::dSigHatDR() const {
   getPDFWeight();
   
   if (theMerger){
     lastMECrossSection(theMerger->MergingDSigDR());
     return lastMECrossSection();
   }
   else if (lastXCombPtr()->willPassCuts() ) {
 	lastME2(me2());
 	CrossSection _dSigHatDRB, _dSigHatDRV, _dSigHatDRI, res = ZERO;
 	// ----- dSigHatDRB -----
 	_dSigHatDRB = oneLoopNoBorn()?ZERO:prefactor() * lastME2();
 	// ----- dSigHatDRV -----
 	_dSigHatDRV = ( oneLoop() && !oneLoopNoLoops() )?(prefactor() * oneLoopInterference()):ZERO;
 	// ----- dSigHatDRI -----
     if  (oneLoop() &&!onlyOneLoop())  {
     for ( auto const & v : virtuals()) {
       v->setXComb(lastXCombPtr());
       res += v->dSigHatDR();
     }
     if ( checkPoles() && oneLoop() )
       logPoles();
   }
 	_dSigHatDRI = res;
 	// ----- finalizing -----
     lastMECrossSection(_dSigHatDRB + _dSigHatDRV + _dSigHatDRI);
     return lastMECrossSection();
   }
   else
   {
   lastME2(ZERO);
   lastMECrossSection(ZERO);
   return lastMECrossSection();
   }
 }
 
 double MatchboxMEBase::oneLoopInterference() const {
 
   if ( matchboxAmplitude() ) {
 
     if ( matchboxAmplitude()->oneLoopAmplitudes() )
       matchboxAmplitude()->prepareOneLoopAmplitudes(this);
 
     double res = 
       matchboxAmplitude()->oneLoopInterference()*
       me2Norm(1);
 
     return res;
 
   }
 
   throw Exception()
     << "MatchboxMEBase::oneLoopInterference() expects a MatchboxAmplitude object.\n"
     << "Please check your setup." << Exception::runerror;
   return 0.;
 
 }
 
 MatchboxMEBase::AccuracyHistogram::AccuracyHistogram(double low,
 						     double up,
 						     unsigned int nbins) 
   : lower(low), upper(up), 
     sameSign(0), oppositeSign(0), nans(0),
     overflow(0), underflow(0) {
 
   double step = (up-low)/nbins;
 
   for ( unsigned int k = 1; k <= nbins; ++k )
     bins[lower + k*step] = 0.0;
 
 }
 
 void MatchboxMEBase::AccuracyHistogram::book(double a, double b) {
   if ( ! (isfinite(a) && isfinite(b)) ) {
     ++nans;
     return;
   }
   if ( a*b >= 0. )
     ++sameSign;
   if ( a*b < 0. )
     ++oppositeSign;
   double r = 1.;
   if ( abs(a) != 0.0 )
     r = abs(1.-abs(b/a));
   else if ( abs(b) != 0.0 )
     r = abs(b);
   if ( log10(r) < lower || r == 0.0 ) {
     ++underflow;
     return;
   }
   if ( log10(r) > upper ) {
     ++overflow;
     return;
   }
   map<double,double>::iterator bin =
     bins.upper_bound(log10(r));
   if ( bin == bins.end() )
     return;
   bin->second += 1.;
 }
 
 void MatchboxMEBase::AccuracyHistogram::dump(const std::string& folder, const std::string& prefix,
 					     const cPDVector& proc) const {
   ostringstream fname("");
   for ( cPDVector::const_iterator p = proc.begin();
 	p != proc.end(); ++p )
     fname << (**p).PDGName();
   ofstream out((folder+"/"+prefix+fname.str()+".dat").c_str());
   out << "# same sign : " << sameSign << " opposite sign : "
       << oppositeSign << " nans : " << nans 
       << " overflow : " << overflow
       << " underflow : " << underflow << "\n";
   for ( map<double,double>::const_iterator b = bins.begin();
 	b != bins.end(); ++b ) {
     map<double,double>::const_iterator bp = b; --bp;
     if ( b->second != 0. ) {
       if ( b != bins.begin() )
 	out << bp->first;
       else
 	out << lower;
       out << " " << b->first
 	  << " " << b->second
 	  << "\n" << flush;
     }
   }
   ofstream gpout((folder+"/"+prefix+fname.str()+".gp").c_str());
   gpout << "set terminal png\n"
       << "set xlabel 'accuracy of pole cancellation [decimal places]'\n"
       << "set ylabel 'counts\n"
       << "set xrange [-20:0]\n"
       << "set output '" << prefix << fname.str() << ".png'\n"
       << "plot '" << prefix << fname.str() << ".dat' using (0.5*($1+$2)):3 with linespoints pt 7 ps 1 not";
 }
 
 void MatchboxMEBase::AccuracyHistogram::persistentOutput(PersistentOStream& os) const {
   os << lower << upper << bins
      << sameSign << oppositeSign << nans
      << overflow << underflow;
 }
 
 void MatchboxMEBase::AccuracyHistogram::persistentInput(PersistentIStream& is) {
   is >> lower >> upper >> bins
      >> sameSign >> oppositeSign >> nans
      >> overflow >> underflow;
 }
 
 void MatchboxMEBase::logPoles() const {
   double res2me = oneLoopDoublePole();
   double res1me = oneLoopSinglePole();
   double res2i = 0.;
   double res1i = 0.;
   for ( auto const & v : virtuals()) {
     res2i += v->oneLoopDoublePole();
     res1i += v->oneLoopSinglePole();
   }
   if (res2me != 0.0 || res2i != 0.0) epsilonSquarePoleHistograms[mePartonData()].book(res2me,res2i);
   if (res1me != 0.0 || res1i != 0.0) epsilonPoleHistograms[mePartonData()].book(res1me,res1i);
 }
 
 bool MatchboxMEBase::haveOneLoop() const {
   if ( matchboxAmplitude() )
     return matchboxAmplitude()->haveOneLoop();
   return false;
 }
 
 bool MatchboxMEBase::onlyOneLoop() const {
   if ( matchboxAmplitude() )
     return matchboxAmplitude()->onlyOneLoop();
   return false;
 }
 
 bool MatchboxMEBase::isDRbar() const {
   if ( matchboxAmplitude() )
     return matchboxAmplitude()->isDRbar();
   return false;
 }
 
 bool MatchboxMEBase::isDR() const {
   if ( matchboxAmplitude() )
     return matchboxAmplitude()->isDR();
   return false;
 }
 
 bool MatchboxMEBase::isCS() const {
   if ( matchboxAmplitude() )
     return matchboxAmplitude()->isCS();
   return false;
 }
 
 bool MatchboxMEBase::isBDK() const {
   if ( matchboxAmplitude() )
     return matchboxAmplitude()->isBDK();
   return false;
 }
 
 bool MatchboxMEBase::isExpanded() const {
   if ( matchboxAmplitude() )
     return matchboxAmplitude()->isExpanded();
   return false;
 }
 
 Energy2 MatchboxMEBase::mu2() const {
   if ( matchboxAmplitude() )
     return matchboxAmplitude()->mu2();
   return 0*GeV2;
 }
 
 double MatchboxMEBase::oneLoopDoublePole() const {
 
   if ( matchboxAmplitude() ) {
 
     return
       matchboxAmplitude()->oneLoopDoublePole()*
       me2Norm(1);
 
   }
 
   return 0.;
 
 }
 
 double MatchboxMEBase::oneLoopSinglePole() const {
 
   if ( matchboxAmplitude() ) {
 
     return 
       matchboxAmplitude()->oneLoopSinglePole()*
       me2Norm(1);
 
   }
 
   return 0.;
 
 }
 
 vector<SubtractionDipolePtr> 
 MatchboxMEBase::getDipoles(const vector<SubtractionDipolePtr>& dipoles,
 			   const vector<MatchboxMEBasePtr> & borns,bool slim) const {
 
   vector<SubtractionDipolePtr> res;
 
   // keep track of the dipoles we already did set up
   set<pair<pair<pair<int,int>,int>,pair<Ptr<MatchboxMEBase>::tptr,Ptr<SubtractionDipole>::tptr> > > done;
 
   cPDVector rep = diagrams().front()->partons();
   int nreal = rep.size();
 
   // now loop over configs
   for ( int emitter = 0; emitter < nreal; ++emitter ) {
 
     list<SubtractionDipolePtr> matchDipoles;
     for ( auto const & d : dipoles ) {
       if ( ! d->canHandleEmitter(rep,emitter) )
 	continue;
       matchDipoles.push_back(d);
     }
     if ( matchDipoles.empty() )
       continue;
 
     for ( int emission = 2; emission < nreal; ++emission ) {
       if ( emission == emitter )
 	continue;
 
       list<SubtractionDipolePtr> matchDipoles2;
       for ( auto const & d : matchDipoles ) {
 	if ( !d->canHandleSplitting(rep,emitter,emission) )
 	  continue;
 	matchDipoles2.push_back(d);
       }
       if ( matchDipoles2.empty() )
 	continue;
 
       map<Ptr<DiagramBase>::ptr,SubtractionDipole::MergeInfo> mergeInfo;
 
       for ( auto const & d : diagrams() ) {
 
 	Ptr<Tree2toNDiagram>::ptr check =
         new_ptr(Tree2toNDiagram(*dynamic_ptr_cast<Ptr<Tree2toNDiagram>::ptr>(d)));
 
 	map<int,int> theMergeLegs;
 
 	for ( unsigned int i = 0; i < check->external().size(); ++i )
 	  theMergeLegs[i] = -1;
 	int theEmitter = check->mergeEmission(emitter,emission,theMergeLegs);
 
 	// no underlying Born
 	if ( theEmitter == -1 )
 	  continue;
 
 	SubtractionDipole::MergeInfo info;
 	info.diagram = check;
 	info.emitter = theEmitter;
 	info.mergeLegs = theMergeLegs;
 	mergeInfo[d] = info;
 
       }
 
       if ( mergeInfo.empty() )
 	continue;
 
       for ( int spectator = 0; spectator < nreal; ++spectator ) {
 	if ( spectator == emitter || spectator == emission )
 	  continue;
 
 	list<SubtractionDipolePtr> matchDipoles3;
     for ( auto const & d : matchDipoles2 ) {
 	  if ( ! d->canHandleSpectator(rep,spectator) )
 	    continue;
 	  matchDipoles3.push_back(d);
 	}
 	if ( matchDipoles3.empty() )
 	  continue;
 
 	if ( noDipole(emitter,emission,spectator) )
 	  continue;
 
         for ( auto const & d : matchDipoles3 ) {
 
 	  if ( !d->canHandle(rep,emitter,emission,spectator) )
 	    continue;
 
           for ( auto const & b : borns ) {
 	    if ( b->onlyOneLoop() )
 	      continue;
 	    if ( done.find(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(b,d)))
 		 != done.end() )
 	      continue;
 	    // now get to work
 	    d->clearBookkeeping();
 	    d->realEmitter(emitter);
 	    d->realEmission(emission);
 	    d->realSpectator(spectator);
 	    d->realEmissionME(const_cast<MatchboxMEBase*>(this));
 	    d->underlyingBornME(b);
 	    d->setupBookkeeping(mergeInfo,slim);
 	    if ( ! d->empty() )  {
 	      res.push_back( d->cloneMe() );
 	      Ptr<SubtractionDipole>::tptr nDipole = res.back();
 	      done.insert(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(b,d)));
 	      if ( nDipole->isSymmetric() )
 		done.insert(make_pair(make_pair(make_pair(emission,emitter),spectator),make_pair(b,d)));
 	      ostringstream dname;
               if ( theMerger) {
                 dname << fullName();
                 if (theOneLoopNoBorn)  dname <<  ".virtual" << "." ;
                 dname   << b->name() << "."
                         << d->name() << ".[("
                         << emitter << "," << emission << ")," << spectator << "]";
               } else {
                 dname << fullName() << "." << b->name() << "."
 		    << d->name() << ".[("
 		    << emitter << "," << emission << ")," << spectator << "]";
               }
 	      if ( ! (generator()->preinitRegister(nDipole,dname.str()) ) )
 		throw Exception() << "MatchboxMEBase::getDipoles(): Dipole " << dname.str() << " already existing." << Exception::runerror;
 	      if ( !factory()->reweighters().empty() ) {
             for ( auto const & rw : factory()->reweighters())
 		  nDipole->addReweighter(rw);
 	      }
 	      if ( !factory()->preweighters().empty() ) {
             for ( auto const & rw : factory()->preweighters() )
 		  nDipole->addPreweighter(rw);
 	      }
 	      nDipole->cloneDependencies(dname.str(),slim);
 	    }
 	  }
 	}
       }
     }
   }
 
   vector<Ptr<SubtractionDipole>::tptr> partners;
   copy(res.begin(),res.end(),back_inserter(partners));
   for ( auto const & d : res )
     d->partnerDipoles(partners);
 
   return res;
 
 }
 
 double MatchboxMEBase::colourCorrelatedME2(pair<int,int> ij) const {
 
   if ( matchboxAmplitude() ) {
 
     if ( matchboxAmplitude()->treeAmplitudes() )
       matchboxAmplitude()->prepareAmplitudes(this);
 
     double res = 
       matchboxAmplitude()->colourCorrelatedME2(ij)*
       me2Norm();
 
     return res;
 
   }
 
   throw Exception()
     << "MatchboxMEBase::colourCorrelatedME2() expects a MatchboxAmplitude object.\n"
     << "Please check your setup." << Exception::runerror;
   return 0.;
 
 }
 
 double MatchboxMEBase::largeNColourCorrelatedME2(pair<int,int> ij,
 						 Ptr<ColourBasis>::tptr largeNBasis) const {
 
   if ( matchboxAmplitude() ) {
 
     if ( matchboxAmplitude()->treeAmplitudes() ) {
       largeNBasis->prepare(mePartonData(),false);
       matchboxAmplitude()->prepareAmplitudes(this);
     }
 
     double res = 
       matchboxAmplitude()->largeNColourCorrelatedME2(ij,largeNBasis)*
       me2Norm();
 
     return res;
 
   }
 
   throw Exception()
     << "MatchboxMEBase::largeNColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
     << "Please check your setup." << Exception::runerror;
   return 0.;
 
 }
 
 double MatchboxMEBase::spinColourCorrelatedME2(pair<int,int> ij,
 					       const SpinCorrelationTensor& c) const {
 
   if ( matchboxAmplitude() ) {
 
     if ( matchboxAmplitude()->treeAmplitudes() )
       matchboxAmplitude()->prepareAmplitudes(this);
 
     double res = 
       matchboxAmplitude()->spinColourCorrelatedME2(ij,c)*
       me2Norm();
 
     return res;
 
   }
 
   throw Exception()
     << "MatchboxMEBase::spinColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
     << "Please check your setup." << Exception::runerror;
   return 0.;
 
 }
 
 double MatchboxMEBase::spinCorrelatedME2(pair<int,int> ij,
 					 const SpinCorrelationTensor& c) const {
 
   if ( matchboxAmplitude() ) {
 
     if ( matchboxAmplitude()->treeAmplitudes() )
       matchboxAmplitude()->prepareAmplitudes(this);
 
     double res = 
       matchboxAmplitude()->spinCorrelatedME2(ij,c)*
       me2Norm();
 
     return res;
 
   }
 
   throw Exception()
     << "MatchboxMEBase::spinCorrelatedME2() expects a MatchboxAmplitude object.\n"
     << "Please check your setup." << Exception::runerror;
   return 0.;
 
 }
 
 
 void MatchboxMEBase::flushCaches() {
   if ( theMerger )theMerger->flushCaches(); 
   MEBase::flushCaches();
   if ( matchboxAmplitude() )
     matchboxAmplitude()->flushCaches();
   for ( auto const & r : reweights() ) 
     r->flushCaches();
   for ( auto const &  v : virtuals()) 
     v->flushCaches();
 }
 
 
 void MatchboxMEBase::setKinematics() {
   MEBase::setKinematics();
   if ( theMerger ) 
     theMerger->setKinematics();
 }
 
 void MatchboxMEBase::clearKinematics() {
   MEBase::clearKinematics();
   if ( theMerger )
     theMerger->clearKinematics();
 }
 
 const MergerBasePtr MatchboxMEBase::merger() const {
   return theMerger;
 }
 
 MergerBasePtr MatchboxMEBase::merger() {
   return theMerger;
 }
 
 void MatchboxMEBase::merger(MergerBasePtr v) {
   theMerger = v;
 }
 
 
 
 
 void MatchboxMEBase::print(ostream& os) const {
 
   os << "--- MatchboxMEBase setup -------------------------------------------------------\n";
 
   os << " '" << name() << "' for subprocess:\n";
 
   os << "  ";
   for ( PDVector::const_iterator pp = subProcess().legs.begin();
 	pp != subProcess().legs.end(); ++pp ) {
     os << (**pp).PDGName() << " ";
     if ( pp == subProcess().legs.begin() + 1 )
       os << "-> ";
   }
   os << "\n";
 
   os << " including " << (oneLoop() ? "" : "no ") << "virtual corrections";
   if ( oneLoopNoBorn() )
     os << " without Born contributions";
   if ( oneLoopNoLoops() )
     os << " without loop contributions";
   os << "\n";
 
   if ( oneLoop() && !onlyOneLoop() ) {
     os << " using insertion operators\n";
     for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
 	    virtuals().begin(); v != virtuals().end(); ++v ) {
       os << " '" << (**v).name() << "' with " 
 	 << ((**v).isDR() ? "" : "C") << "DR/";
       if ( (**v).isCS() )
 	os << "CS";
       if ( (**v).isBDK() )
 	os << "BDK";
       if ( (**v).isExpanded() )
 	os << "expanded";
       os << " conventions\n";
     }
   }
 
   os << "--------------------------------------------------------------------------------\n";
 
   os << flush;
 
 }
 
 void MatchboxMEBase::printLastEvent(ostream& os) const {
 
   os << "--- MatchboxMEBase last event information --------------------------------------\n";
 
   os << " for matrix element '" << name() << "'\n";
 
   os << " process considered:\n ";
 
   int in = 0;
   for ( cPDVector::const_iterator p = mePartonData().begin();
 	p != mePartonData().end(); ++p ) {
     os << (**p).PDGName() << " ";
     if ( ++in == 2 )
       os << " -> ";
   }
 
   os << " kinematic environment as set by the XComb " << lastXCombPtr() << ":\n"
      << " sqrt(shat)/GeV = " << sqrt(lastSHat()/GeV2)
      << " x1 = " << lastX1() << " x2 = " << lastX2() 
      << " alphaS = " << lastAlphaS() << "\n";
 
   os << " momenta/GeV generated from random numbers\n ";
   copy(lastXComb().lastRandomNumbers().begin(),
        lastXComb().lastRandomNumbers().end(),ostream_iterator<double>(os," "));
   os << ":\n ";
 
   for ( vector<Lorentz5Momentum>::const_iterator p = meMomenta().begin();
 	p != meMomenta().end(); ++p ) {
     os << (*p/GeV) << "\n ";
   }
 
   os << "last cross section/nb calculated was:\n "
      << (lastMECrossSection()/nanobarn) << " (pdf weight " << lastMEPDFWeight() << ")\n";
 
   os << "--------------------------------------------------------------------------------\n";
 
   os << flush;
 
 }
 
 void MatchboxMEBase::logGenerateKinematics(const double * r) const {
 
   if ( !verbose() )
     return;
 
   generator()->log() << "'" << name() << "' generated kinematics\nfrom "
 		     << nDim() << " random numbers:\n";
   copy(r,r+nDim(),ostream_iterator<double>(generator()->log()," "));
   generator()->log() << "\n";
 
   generator()->log() << "storing phase space information in XComb "
 		     << lastXCombPtr() << "\n";
 
   generator()->log() << "generated phase space point (in GeV):\n";
 
   vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
   cPDVector::const_iterator dit = mePartonData().begin();
 
   for ( ; pit != meMomenta().end() ; ++pit, ++dit )
     generator()->log() << (**dit).PDGName() << " : "
 		       << (*pit/GeV) << "\n";
 
   generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
 		     << "and Jacobian = " << jacobian() << " sHat/GeV2 = "
 		     << (lastSHat()/GeV2) << "\n" << flush;
 
 }
 
 void MatchboxMEBase::logSetScale() const {
 
   if ( !verbose() )
     return;
 
   generator()->log() << "'" << name() << "' set scales using XComb " << lastXCombPtr() << ":\n"
 		     << "scale/GeV2 = " << (scale()/GeV2) << " xi_R = "
 		     << renormalizationScaleFactor() << " xi_F = "
 		     << factorizationScaleFactor() << "\n"
 		     << "alpha_s = " << lastAlphaS() << "\n" << flush;
 
 }
 
 void MatchboxMEBase::logPDFWeight() const {
 
   if ( !verbose() )
     return;
 
   generator()->log() << "'" << name() << "' calculated pdf weight = "
 		     << lastMEPDFWeight() << " from XComb "
 		     << lastXCombPtr() << "\n"
 		     << "x1 = " << lastX1() << " (" << (mePartonData()[0]->coloured() ? "" : "not ") << "used) "
 		     << "x2 = " << lastX2() << " (" << (mePartonData()[1]->coloured() ? "" : "not ") << "used)\n"
 		     << flush;
 
 }
 
 void MatchboxMEBase::logME2() const {
 
   if ( !verbose() )
     return;
 
   generator()->log() << "'" << name() << "' evaluated me2 using XComb "
 		     << lastXCombPtr() << "\n"
 		     << "and phase space point (in GeV):\n";
 
   vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
   cPDVector::const_iterator dit = mePartonData().begin();
 
   for ( ; pit != meMomenta().end() ; ++pit, ++dit )
     generator()->log() << (**dit).PDGName() << " : "
 		       << (*pit/GeV) << "\n";
 
   generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
 		     << "sHat/GeV2 = " << (lastSHat()/GeV2) << "\n" << flush;
 
 }
 
 void MatchboxMEBase::logDSigHatDR() const {
 
   if ( !verbose() )
     return;
 
   generator()->log() << "'" << name() << "' evaluated cross section using XComb "
 		     << lastXCombPtr() << "\n"
 		     << "Jacobian = " << jacobian() << " sHat/GeV2 = "
 		     << (lastSHat()/GeV2) << " dsig/nb = "
 		     << (lastMECrossSection()/nanobarn) << "\n" << flush;
 
 }
 
 void MatchboxMEBase::cloneDependencies(const std::string& prefix,bool slim) {
 
   if ( phasespace() && !slim ) {
     Ptr<MatchboxPhasespace>::ptr myPhasespace = phasespace()->cloneMe();
     ostringstream pname;
     pname << (prefix == "" ? fullName() : prefix) << "/" << myPhasespace->name();
     if ( ! (generator()->preinitRegister(myPhasespace,pname.str()) ) )
       throw Exception() << "MatchboxMEBase::cloneDependencies(): Phasespace generator " << pname.str() << " already existing." << Exception::runerror;
     myPhasespace->cloneDependencies(pname.str());
     phasespace(myPhasespace);
   }
 
   theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());
 
   if ( matchboxAmplitude() ) {
     Ptr<MatchboxAmplitude>::ptr myAmplitude = matchboxAmplitude()->cloneMe();
     ostringstream pname;
     pname << (prefix == "" ? fullName() : prefix) << "/" << myAmplitude->name();
     if ( ! (generator()->preinitRegister(myAmplitude,pname.str()) ) ){
       throw Exception() << "MatchboxMEBase::cloneDependencies(): Amplitude " << pname.str() << " already existing." << Exception::runerror;
     }
     myAmplitude->cloneDependencies(pname.str(),slim);
     matchboxAmplitude(myAmplitude);
     amplitude(myAmplitude);
     matchboxAmplitude()->orderInGs(orderInAlphaS());
     matchboxAmplitude()->orderInGem(orderInAlphaEW());
   }
 
   if ( scaleChoice() &&!slim ) {
     Ptr<MatchboxScaleChoice>::ptr myScaleChoice = scaleChoice()->cloneMe();
     ostringstream pname;
     pname << (prefix == "" ? fullName() : prefix) << "/" << myScaleChoice->name();
     if ( ! (generator()->preinitRegister(myScaleChoice,pname.str()) ) )
       throw Exception() << "MatchboxMEBase::cloneDependencies(): Scale choice " << pname.str() << " already existing." << Exception::runerror;
     scaleChoice(myScaleChoice);
   }
 
   for ( auto &  rw : theReweights ) {
     Ptr<MatchboxReweightBase>::ptr myReweight = rw->cloneMe();
     ostringstream pname;
     pname << (prefix == "" ? fullName() : prefix) << "/" << rw->name();
     if ( ! (generator()->preinitRegister(myReweight,pname.str()) ) )
       throw Exception() << "MatchboxMEBase::cloneDependencies(): Reweight " << pname.str() << " already existing." << Exception::runerror;
     myReweight->cloneDependencies(pname.str());
     rw = myReweight;
   }
 
   for ( auto & v : virtuals()) {
     Ptr<MatchboxInsertionOperator>::ptr myIOP = v->cloneMe();
     ostringstream pname;
     pname << (prefix == "" ? fullName() : prefix) << "/" << v->name();
     if ( ! (generator()->preinitRegister(myIOP,pname.str()) ) )
       throw Exception() << "MatchboxMEBase::cloneDependencies(): Insertion operator " << pname.str() << " already existing." << Exception::runerror;
     v = myIOP;
   }
 
 }
 
 void MatchboxMEBase::prepareXComb(MatchboxXCombData& xc) const {
 
   // fixme We need to pass on the partons from the xcmob here, not
   // assuming one subprocess per matrix element
 
   if ( phasespace() )
     xc.nDimPhasespace(phasespace()->nDim(diagrams().front()->partons()));
 
   if ( matchboxAmplitude() ) {
     xc.nDimAmplitude(matchboxAmplitude()->nDimAdditional());
     if ( matchboxAmplitude()->colourBasis() ) {
       size_t cdim = 
  	matchboxAmplitude()->colourBasis()->prepare(diagrams(),noCorrelations());
       xc.colourBasisDim(cdim);
     }
     if ( matchboxAmplitude()->isExternal() ) {
       xc.externalId(matchboxAmplitude()->externalId(diagrams().front()->partons()));
     }
   }
 
   int insertionAdd = 0;
   for ( auto const &  v : virtuals() )
     insertionAdd = max(insertionAdd,v->nDimAdditional());
   
   xc.nDimInsertions(insertionAdd);
 
   xc.nLight(getNLight());
   if(xc.nLightJetVec().empty())
   for (auto const & id : getNLightJetVec())
     xc.nLightJetVec( id );
   
   if(xc.nHeavyJetVec().empty())
   for (auto const & id :getNHeavyJetVec())
     xc.nHeavyJetVec(id);
 
   if(xc.nLightProtonVec().empty())
   for (auto const & id : getNLightProtonVec())
     xc.nLightProtonVec(id);
 
   xc.olpId(olpProcess());
 
   if ( initVerbose() ) {
     ostringstream fname_strm;
     // only allow alphanumeric, / and _ in filename
     for (const char c : name()) {
         switch (c) {
           case '+' : fname_strm << "+"; break;
           case '-' : fname_strm << "-"; break;
           case '~' : fname_strm << "_tilde"; break;
           case ']' : break;
           case ',' : fname_strm << "__"; break;
           default  : fname_strm << (isalnum(c) ? c : '_'); break;
         }
     }
     fname_strm << ".diagrams";
     const string fname = fname_strm.str();
     ifstream test(fname.c_str());
     if ( !test ) {
       test.close();
       ofstream out(fname.c_str());
       for ( vector<Ptr<DiagramBase>::ptr>::const_iterator d = diagrams().begin();
 	    d != diagrams().end(); ++d ) {
 	DiagramDrawer::drawDiag(out,dynamic_cast<const Tree2toNDiagram&>(**d));
 	out << "\n";
       }
     }
   }
 
 }
 
 StdXCombPtr MatchboxMEBase::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
 				      tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
 				      tPExtrPtr newExtractor,	tCascHdlPtr newCKKW,
 				      const PBPair & newPartonBins, tCutsPtr newCuts,
 				      const DiagramVector & newDiagrams, bool mir,
 				      const PartonPairVec&,
 				      tStdXCombPtr newHead,
 				      tMEPtr newME) {
 
   if ( !newME )
     newME = this;
 
   Ptr<MatchboxXComb>::ptr xc =
     new_ptr(MatchboxXComb(newMaxEnergy, inc,
 			  newEventHandler, newSubProcessHandler,
 			  newExtractor, newCKKW,
 			  newPartonBins, newCuts, newME,
 			  newDiagrams, mir,
 			  newHead));
 
   prepareXComb(*xc);
 
   return xc;
 
 }
 
 StdXCombPtr MatchboxMEBase::makeXComb(tStdXCombPtr newHead,
 				      const PBPair & newPartonBins,
 				      const DiagramVector & newDiagrams,
 				      tMEPtr newME) {
   if ( !newME )
     newME = this;
 
   Ptr<MatchboxXComb>::ptr xc = 
     new_ptr(MatchboxXComb(newHead, newPartonBins, newME, newDiagrams));
 
   prepareXComb(*xc);
 
   return xc;
 }
 
 void MatchboxMEBase::persistentOutput(PersistentOStream & os) const {
   os << theLastXComb << thePhasespace 
      << theAmplitude << theScaleChoice << theVirtuals 
      << theReweights << theSubprocess << theOneLoop 
      << theOneLoopNoBorn << theOneLoopNoLoops
      << epsilonSquarePoleHistograms << epsilonPoleHistograms
      << theMerger
      << theOLPProcess << theNoCorrelations
      << theHavePDFs << checkedPDFs;
 }
 
 void MatchboxMEBase::persistentInput(PersistentIStream & is, int) {
   is >> theLastXComb >> thePhasespace 
      >> theAmplitude >> theScaleChoice >> theVirtuals 
      >> theReweights >> theSubprocess >> theOneLoop 
      >> theOneLoopNoBorn >> theOneLoopNoLoops
      >> epsilonSquarePoleHistograms >> epsilonPoleHistograms
      >> theMerger
      >> theOLPProcess >> theNoCorrelations
      >> theHavePDFs >> checkedPDFs;
   lastMatchboxXComb(theLastXComb);
 }
 
 void MatchboxMEBase::Init() {
 
   static ClassDocumentation<MatchboxMEBase> documentation
     ("MatchboxMEBase is the base class for matrix elements "
      "in the context of the matchbox NLO interface.");
 
 }
 
 IBPtr MatchboxMEBase::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr MatchboxMEBase::fullclone() const {
   return new_ptr(*this);
 }
 
 void MatchboxMEBase::doinit() {
   MEBase::doinit();
   if ( !theAmplitude )
     theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());
   if ( matchboxAmplitude() )
     matchboxAmplitude()->init();
   if ( phasespace() ) {
     phasespace()->init();
     matchboxAmplitude()->checkReshuffling(phasespace());
   }
   if ( scaleChoice() ) {
     scaleChoice()->init();
   }
   for (auto const & rw : theReweights)
     rw->init();
   for (auto const & v :  virtuals() )
     v->init();
 }
 
 
 void MatchboxMEBase::doinitrun() {
   MEBase::doinitrun();
   if ( matchboxAmplitude() )
     matchboxAmplitude()->initrun();
   
   if ( phasespace() )
     phasespace()->initrun();
   
   if ( scaleChoice() )
     scaleChoice()->initrun();
   
   for (auto const & rw : theReweights)
     rw->initrun();
   for (auto const & v :  virtuals() )
     v->initrun();
 }
   
 void MatchboxMEBase::dofinish() {
   MEBase::dofinish();
   for (auto const & b : epsilonSquarePoleHistograms ) {
     b.second.dump(factory()->poleData(),"epsilonSquarePoles-",b.first);
   }
   for (auto const & b :  epsilonPoleHistograms ) {
     b.second.dump(factory()->poleData(),"epsilonPoles-",b.first);
   }
 }
 
 // *** 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<MatchboxMEBase,MEBase>
 describeHerwigMatchboxMEBase("Herwig::MatchboxMEBase", "Herwig.so");
 
diff --git a/PDF/HwRemDecayer.cc b/PDF/HwRemDecayer.cc
--- a/PDF/HwRemDecayer.cc
+++ b/PDF/HwRemDecayer.cc
@@ -1,1537 +1,2183 @@
 // -*- C++ -*-
 //
 // HwRemDecayer.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 HwRemDecayer class.
 //
 
 #include "HwRemDecayer.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Interface/Reference.h"  
 #include "ThePEG/Interface/Parameter.h"  
 #include "ThePEG/Interface/Switch.h"  
 #include "ThePEG/Utilities/UtilityBase.h"
 #include "ThePEG/Utilities/SimplePhaseSpace.h"
 #include "ThePEG/Utilities/Throw.h"
 
 #include "Herwig/Shower/ShowerHandler.h"
 
 using namespace Herwig;
 
 namespace{
 
 const bool dbg = false;
 
 
   void reShuffle(Lorentz5Momentum &p1, Lorentz5Momentum &p2, Energy m1, Energy m2){
 
     Lorentz5Momentum ptotal(p1+p2);
     ptotal.rescaleMass();
 
     if( ptotal.m() < m1+m2 ) {
       if(dbg)
 	cerr << "Not enough energy to perform reshuffling \n";
       throw HwRemDecayer::ExtraSoftScatterVeto();
     }
 
     Boost boostv = -ptotal.boostVector();
     ptotal.boost(boostv);
 
     p1.boost(boostv);
     // set the masses and energies,
     p1.setMass(m1);
     p1.setE(0.5/ptotal.m()*(ptotal.m2()+sqr(m1)-sqr(m2)));
     p1.rescaleRho();
     // boost back to the lab
     p1.boost(-boostv);
 
     p2.boost(boostv);
     // set the masses and energies,
     p2.setMass(m2);
     p2.setE(0.5/ptotal.m()*(ptotal.m2()+sqr(m2)-sqr(m1)));
     p2.rescaleRho();
     // boost back to the lab
     p2.boost(-boostv);
   }
 }
 
 void HwRemDecayer::initialize(pair<tRemPPtr, tRemPPtr> rems, tPPair beam, Step & step,
 			      Energy forcedSplitScale) {
   // the step
   thestep = &step;
   // valence content of the hadrons
   theContent.first  = getHadronContent(beam.first);
   theContent.second = getHadronContent(beam.second);
   // momentum extracted from the hadrons
   theUsed.first  = Lorentz5Momentum();
   theUsed.second = Lorentz5Momentum();
   theMaps.first.clear();
   theMaps.second.clear();
   theX.first = 0.0;
   theX.second = 0.0;
   theRems = rems;
   _forcedSplitScale = forcedSplitScale;
   // check remnants attached to the right hadrons
   if( (theRems.first  && parent(theRems.first ) != beam.first ) ||
       (theRems.second && parent(theRems.second) != beam.second) )
     throw Exception() << "Remnant order wrong in "
 		      << "HwRemDecayer::initialize(...)"
 		      << Exception::runerror; 
   return;
 }
 
 void HwRemDecayer::split(tPPtr parton, HadronContent & content, 
 			 tRemPPtr rem, Lorentz5Momentum & used, 
 			 PartnerMap &partners, tcPDFPtr pdf, bool first) {
   theBeam = parent(rem);
   theBeamData = dynamic_ptr_cast<Ptr<BeamParticleData>::const_pointer>
     (theBeam->dataPtr());
   double currentx = parton->momentum().rho()/theBeam->momentum().rho();
   double check = rem==theRems.first ? theX.first : theX.second;
   check += currentx;
   if(1.0-check < 1e-3) throw ShowerHandler::ExtraScatterVeto();
   bool anti;
   Lorentz5Momentum lastp(parton->momentum());
   int lastID(parton->id());
   Energy oldQ(_forcedSplitScale);
   _pdf = pdf;
   //do nothing if already valence quark
   if(first && content.isValenceQuark(parton)) { 
     //set the extracted value, because otherwise no RemID could be generated.
     content.extract(lastID);
     // add the particle to the colour partners
     partners.push_back(make_pair(parton, tPPtr()));
     //set the sign
     anti = parton->hasAntiColour() && parton->id()!=ParticleID::g;
     if(rem==theRems.first) theanti.first  = anti;
     else                   theanti.second = anti;
     // add the x and return
     if(rem==theRems.first) theX.first  += currentx;
     else                   theX.second += currentx;
     return;
   }
   //or gluon for secondaries
   else if(!first && lastID == ParticleID::g) {
     partners.push_back(make_pair(parton, tPPtr()));
     // add the x and return
     if(rem==theRems.first) theX.first  += currentx;
     else                   theX.second += currentx;
     return; 
   }
   // if a sea quark.antiquark forced splitting to a gluon
   // Create the new parton with its momentum and parent/child relationship set
   PPtr newSea;
   if( !(lastID == ParticleID::g || 
 	lastID == ParticleID::gamma) ) {
     newSea = forceSplit(rem, -lastID, oldQ, currentx, lastp, used,content);
     ColinePtr cl = new_ptr(ColourLine());
     if(newSea->id() > 0) cl->    addColoured(newSea);
     else                 cl->addAntiColoured(newSea);
     // if a secondard scatter finished so return
     if(!first || content.isValenceQuark(ParticleID::g) ){
       partners.push_back(make_pair(parton, newSea));
       // add the x and return
       if(rem==theRems.first) theX.first  += currentx;
       else                   theX.second += currentx;
       if(first) content.extract(ParticleID::g);
       return;
     }
   }
   // otherwise evolve back to valence
   // final valence splitting
   PPtr newValence = forceSplit(rem, 
 			       lastID!=ParticleID::gamma ? 
 			       ParticleID::g : ParticleID::gamma,
 			       oldQ, currentx , lastp, used, content);
   // extract from the hadron to allow remnant to be determined
   content.extract(newValence->id());
   // case of a gluon going into the hard subprocess
   if( lastID == ParticleID::g ) {
     partners.push_back(make_pair(parton, tPPtr()));
     anti = newValence->hasAntiColour();
     if(rem==theRems.first) theanti.first  = anti;
     else                   theanti.second = anti;
     parton->colourLine(!anti)->addColoured(newValence, anti);
     return;
   }
   else if( lastID == ParticleID::gamma) {
     partners.push_back(make_pair(parton, newValence));
     anti = newValence->hasAntiColour();
     ColinePtr newLine(new_ptr(ColourLine()));
     newLine->addColoured(newValence, anti);
     if(rem==theRems.first) theanti.first  = anti;
     else                   theanti.second = anti;
     // add the x and return
     if(rem==theRems.first) theX.first  += currentx;
     else                   theX.second += currentx;
     return;
   }
   //The valence quark will always be connected to the sea quark with opposite sign
   tcPPtr particle;
   if(lastID*newValence->id() < 0){
     particle = parton;
     partners.push_back(make_pair(newSea, tPPtr()));
   } 
   else {
     particle = newSea;
     partners.push_back(make_pair(parton, tPPtr()));
   }
   anti = newValence->hasAntiColour();
   if(rem==theRems.first) theanti.first  = anti;
   else                   theanti.second = anti;
   
   if(particle->colourLine()) 
     particle->colourLine()->addAntiColoured(newValence);
   if(particle->antiColourLine()) 
     particle->antiColourLine()->addColoured(newValence);
   // add the x and return
   if(rem==theRems.first) theX.first  += currentx;
   else                   theX.second += currentx;
   return;
 }
 
 void HwRemDecayer::doSplit(pair<tPPtr, tPPtr> partons,
 			   pair<tcPDFPtr, tcPDFPtr> pdfs,
 			   bool first) {
   if(theRems.first) {
     ParticleVector children=theRems.first->children();
     for(unsigned int ix=0;ix<children.size();++ix) {
       if(children[ix]->dataPtr()==theRems.first->dataPtr()) 
 	theRems.first = dynamic_ptr_cast<RemPPtr>(children[ix]);
     }
   }
   if(theRems.second) {
     ParticleVector children=theRems.second->children();
     for(unsigned int ix=0;ix<children.size();++ix) {
       if(children[ix]->dataPtr()==theRems.second->dataPtr()) 
 	theRems.second = dynamic_ptr_cast<RemPPtr>(children[ix]);
     }
   }
   // forced splitting for first parton
   if(isPartonic(partons.first )) { 
     try {
       split(partons.first, theContent.first, theRems.first, 
 	    theUsed.first, theMaps.first, pdfs.first, first);
     }
     catch(ShowerHandler::ExtraScatterVeto) {
       throw ShowerHandler::ExtraScatterVeto();
     }
   }
   // forced splitting for second parton
   if(isPartonic(partons.second)) { 
     try {
       split(partons.second, theContent.second, theRems.second, 
 	    theUsed.second, theMaps.second, pdfs.second, first);
       // additional check for the remnants
       // if can't do the rescale veto the emission
       if(!first&&partons.first->data().coloured()&&
 	 partons.second->data().coloured()) {
 	Lorentz5Momentum pnew[2]=
 	  {theRems.first->momentum()  - theUsed.first  - partons.first->momentum(),
 	   theRems.second->momentum() - theUsed.second - partons.second->momentum()};
 
 	pnew[0].setMass(getParticleData(theContent.first.RemID())->constituentMass());
 	pnew[0].rescaleEnergy();
 	pnew[1].setMass(getParticleData(theContent.second.RemID())->constituentMass());
 	pnew[1].rescaleEnergy();
 	
 	for(unsigned int iy=0; iy<theRems.first->children().size(); ++iy)
 	  pnew[0] += theRems.first->children()[iy]->momentum();
 	
 	for(unsigned int iy=0; iy<theRems.second->children().size(); ++iy)
 	  pnew[1] += theRems.second->children()[iy]->momentum();
 	
 	Lorentz5Momentum ptotal=
 	  theRems.first ->momentum()-partons.first ->momentum()+
 	  theRems.second->momentum()-partons.second->momentum();
 	// add x limits
 	if(ptotal.m() < (pnew[0].m() + pnew[1].m()) ) {
 	  if(partons.second->id() != ParticleID::g){
 	    if(partons.second==theMaps.second.back().first) 
 	      theUsed.second -= theMaps.second.back().second->momentum();
 	    else
 	      theUsed.second -= theMaps.second.back().first->momentum();
 	    
 	    thestep->removeParticle(theMaps.second.back().first);
 	    thestep->removeParticle(theMaps.second.back().second);
 	  }
 	  theMaps.second.pop_back();
 	  theX.second -= partons.second->momentum().rho()/
 	    parent(theRems.second)->momentum().rho();
 	  throw ShowerHandler::ExtraScatterVeto();
 	}
       }
     }
     catch(ShowerHandler::ExtraScatterVeto){
       if(!partons.first||!partons.second||
 	 !theRems.first||!theRems.second) 
 	throw ShowerHandler::ExtraScatterVeto();
       //case of the first forcedSplitting worked fine
       theX.first -= partons.first->momentum().rho()/
 	parent(theRems.first)->momentum().rho();
       //case of the first interaction
       //throw veto immediately, because event get rejected anyway.
       if(first) throw ShowerHandler::ExtraScatterVeto();
       //secondary interactions have to end on a gluon, if parton 
       //was NOT a gluon, the forced splitting particles must be removed
       if(partons.first->id() != ParticleID::g) {
 	if(partons.first==theMaps.first.back().first) 
 	  theUsed.first -= theMaps.first.back().second->momentum();
 	else
 	  theUsed.first -= theMaps.first.back().first->momentum();
 	
 	thestep->removeParticle(theMaps.first.back().first);
 	thestep->removeParticle(theMaps.first.back().second);
       }
       theMaps.first.pop_back();
       throw ShowerHandler::ExtraScatterVeto();
     }
   }
   // veto if not enough energy for extraction
   if( !first &&(theRems.first ->momentum().e() - 
 		partons.first ->momentum().e() < 1.0e-3*MeV ||
 		theRems.second->momentum().e() - 
 		partons.second->momentum().e() < 1.0e-3*MeV )) {
     if(partons.first->id() != ParticleID::g) {
       if(partons.first==theMaps.first.back().first) 
 	theUsed.first -= theMaps.first.back().second->momentum();
       else
 	theUsed.first -= theMaps.first.back().first->momentum();
       thestep->removeParticle(theMaps.first.back().first);
       thestep->removeParticle(theMaps.first.back().second);
     }
     theMaps.first.pop_back();
     if(partons.second->id() != ParticleID::g) {
       if(partons.second==theMaps.second.back().first) 
 	theUsed.second -= theMaps.second.back().second->momentum();
       else
 	theUsed.second -= theMaps.second.back().first->momentum();
       thestep->removeParticle(theMaps.second.back().first);
       thestep->removeParticle(theMaps.second.back().second);
     }
     theMaps.second.pop_back();
     throw ShowerHandler::ExtraScatterVeto();
   }
 }
 
 void HwRemDecayer::mergeColour(tPPtr pold, tPPtr pnew, bool anti) const {
   ColinePtr clnew, clold;
   
   //save the corresponding colour lines
   clold = pold->colourLine(anti);
   clnew = pnew->colourLine(!anti);
 
   assert(clold);
 
   // There is already a colour line (not the final diquark)
   if(clnew){
 
     if( (clnew->coloured().size() + clnew->antiColoured().size()) > 1 ){
       if( (clold->coloured().size() + clold->antiColoured().size()) > 1 ){
 	//join the colour lines
 	//I don't use the join method, because potentially only (anti)coloured
 	//particles belong to one colour line
 	if(clold!=clnew){//procs are not already connected
 	  while ( !clnew->coloured().empty() ) {
 	    tPPtr p = clnew->coloured()[0];
 	    clnew->removeColoured(p);
 	    clold->addColoured(p);
 	  }
 	  while ( !clnew->antiColoured().empty() ) {
 	    tPPtr p = clnew->antiColoured()[0];
 	    clnew->removeAntiColoured(p);
 	    clold->addAntiColoured(p);
 	  }
 	}
 
       }else{
 	//if pold is the only member on it's 
 	//colour line, remove it.
 	clold->removeColoured(pold, anti);
 	//and add it to clnew
 	clnew->addColoured(pold, anti);
       }    
     } else{//pnnew is the only member on it's colour line.
       clnew->removeColoured(pnew, !anti);
       clold->addColoured(pnew, !anti);
     }
   } else {//there is no coline at all for pnew
     clold->addColoured(pnew, !anti);
   }
 }
 
 void HwRemDecayer::fixColours(PartnerMap partners, bool anti, 
 			      double colourDisrupt) const {
   PartnerMap::iterator prev;
   tPPtr pnew, pold;
   assert(partners.size()>=2);
 
   PartnerMap::iterator it=partners.begin();
   while(it != partners.end()) {
     //skip the first one to have a partner
     if(it==partners.begin()){
       it++;
       continue;
     }
     
     prev = it - 1;
     //determine the particles to work with
     pold = prev->first;
     if(prev->second) {
       if(!pold->coloured())
 	pold = prev->second;
       else if(pold->hasAntiColour() != anti)
 	pold = prev->second;
     }
     assert(pold);
 
     pnew = it->first;
     if(it->second) {
       if(it->second->colourLine(!anti)) //look for the opposite colour
 	pnew = it->second;
     }
     assert(pnew);
 
     // Implement the disruption of colour connections
     if( it != partners.end()-1 ) {//last one is diquark-has to be connected
       
       //has to be inside the if statement, so that the probability is
       //correctly counted:
       if( UseRandom::rnd() < colourDisrupt ){
 
 	if(!it->second){//check, whether we have a gluon
 	  mergeColour(pnew, pnew, anti);
 	}else{
 	  if(pnew==it->first)//be careful about the order
 	    mergeColour(it->second, it->first, anti);
 	  else
 	    mergeColour(it->first, it->second, anti);
 	}
 
 	it = partners.erase(it);
 	continue;
       }
 
     }
     // regular merging
     mergeColour(pold, pnew, anti);
     //end of loop
     it++;
   }
   return;
 }
 
 PPtr HwRemDecayer::forceSplit(const tRemPPtr rem, long child, Energy &lastQ, 
 			      double &lastx, Lorentz5Momentum &pf, 
 			      Lorentz5Momentum &p, 
 			      HadronContent & content) const {
   static const double eps=1e-6;
   // beam momentum
   Lorentz5Momentum beam = theBeam->momentum();
   // the last scale is minimum of last value and upper limit
   Energy minQ=_range*_kinCutoff*sqrt(lastx)/(1-lastx);
   if(minQ>lastQ) lastQ=minQ;
   // generate the new value of qtilde
   // weighted towards the lower value: dP/dQ = 1/Q -> Q(R) =
   // Q0 (Qmax/Q0)^R
   Energy q;
   unsigned int ntry=0,maxtry=100;
   double xExtracted = rem==theRems.first ? theX.first : theX.second;
   double zmin=  lastx/(1.-xExtracted) ,zmax,yy;
 
   if(1-lastx<eps) throw ShowerHandler::ExtraScatterVeto();
   do {
     q = minQ*pow(lastQ/minQ,UseRandom::rnd());
     yy   = 1.+0.5*sqr(_kinCutoff/q);
     zmax = yy - sqrt(sqr(yy)-1.); 
     ++ntry;
   }
   while(zmax<zmin&&ntry<maxtry);
   if(ntry==maxtry) throw ShowerHandler::ExtraScatterVeto();
   if(zmax-zmin<eps) throw ShowerHandler::ExtraScatterVeto();
   // now generate z as in FORTRAN HERWIG
   // use y = ln(z/(1-z)) as integration variable
   double ymin=log(zmin/(1.-zmin));
   double ymax=log(zmax/(1.-zmax));
   double dely=ymax-ymin;
   unsigned int nz=_nbinmax;
   dely/=nz;
   yy=ymin+0.5*dely;
   vector<int> ids;
   if(child==21||child==22) {
     ids=content.flav;
     for(unsigned int ix=0;ix<ids.size();++ix) ids[ix] *= content.sign;
   }
   else {
     ids.push_back(ParticleID::g);
   }
   // probabilities of the different types of possible splitting
   map<long,pair<double,vector<double> > > partonprob;
   double ptotal(0.);
   for(unsigned int iflav=0;iflav<ids.size();++iflav) {
     // only do each parton once
     if(partonprob.find(ids[iflav])!=partonprob.end()) continue;
     // particle data object
     tcPDPtr in = getParticleData(ids[iflav]);
     double psum(0.);
     vector<double> prob;
     for(unsigned int iz=0;iz<nz;++iz) {
       double ez=exp(yy);
       double wr=1.+ez;
       double zr=wr/ez;
       double wz=1./wr;
       double zz=wz*ez;
       double coup = child!=22 ? 
 	_alphaS ->value(sqr(max(wz*q,_kinCutoff))) :
 	_alphaEM->value(sqr(max(wz*q,_kinCutoff)));
       double az=wz*zz*coup;
       // g -> q qbar
       if(ids[iflav]==ParticleID::g) {
 	// calculate splitting function   
 	// SP as q is always less than forcedSplitScale, the pdf scale is fixed 
 	// pdfval = _pdf->xfx(theBeamData,in,sqr(q),lastx*zr); 
 	double pdfval=_pdf->xfx(theBeamData,in,sqr(_forcedSplitScale),lastx*zr);
 	if(pdfval>0.) psum += pdfval*az*0.5*(sqr(zz)+sqr(wz));
       }
       // q -> q g
       else {
 	// calculate splitting function
 	// SP as q is always less than forcedSplitScale, the pdf scale is fixed 
  	// pdfval = _pdf->xfx(theBeamData,in,sqr(q),lastx*zr);
  	double pdfval=_pdf->xfx(theBeamData,in,sqr(_forcedSplitScale),lastx*zr);
 	if(pdfval>0.) psum += pdfval*az*4./3.*(1.+sqr(wz))*zr;
       }
       if(psum>0.) prob.push_back(psum);
       yy+=dely;
     }
     if(psum>0.) partonprob[ids[iflav]] = make_pair(psum,prob);
     ptotal+=psum;
   }
   // select the flavour
   if(ptotal==0.) throw ShowerHandler::ExtraScatterVeto();
   ptotal *= UseRandom::rnd();
   map<long,pair<double,vector<double> > >::const_iterator pit;
   for(pit=partonprob.begin();pit!=partonprob.end();++pit) {
     if(pit->second.first>=ptotal) break;
     else                          ptotal -= pit->second.first;
   }
   if(pit==partonprob.end()) 
     throw Exception() << "Can't select parton for forced backward evolution in "
 		      << "HwRemDecayer::forceSplit" << Exception::eventerror;
   // select z 
   unsigned int iz=0;
   for(;iz<pit->second.second.size();++iz) {
     if(pit->second.second[iz]>ptotal) break;
   }
   if(iz==pit->second.second.size()) --iz;
   double ey=exp(ymin+dely*(float(iz+1)-UseRandom::rnd()));
   double z=ey/(1.+ey);
   Energy2 pt2=sqr((1.-z)*q)- z*sqr(_kinCutoff);
   // create the particle
   if(pit->first!=ParticleID::g) child=pit->first;
   PPtr parton = getParticleData(child)->produceParticle();
   Energy2 emittedm2 = sqr(parton->dataPtr()->constituentMass());
   // Now boost pcm and pf to z only frame
   Lorentz5Momentum p_ref = Lorentz5Momentum(ZERO,  beam.vect());
   Lorentz5Momentum n_ref = Lorentz5Momentum(ZERO, -beam.vect());
   // generate phi and compute pt of branching
   double phi = Constants::twopi*UseRandom::rnd();
   Energy pt=sqrt(pt2);
   Lorentz5Momentum qt   = LorentzMomentum(pt*cos(phi), pt*sin(phi), ZERO, ZERO);
   Axis axis(p_ref.vect().unit());
   if(axis.perp2()>0.) {
     LorentzRotation rot;
     double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
     rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
     qt.transform(rot);
   }
   // compute alpha for previous particle
   Energy2 p_dot_n  = p_ref*n_ref;
   double lastalpha =    pf*n_ref/p_dot_n;
   Lorentz5Momentum qtout=qt;
   Energy2 qtout2=-qt*qt;
   double alphaout=(1.-z)/z*lastalpha;
   double betaout=0.5*(emittedm2+qtout2)/alphaout/p_dot_n;
   Lorentz5Momentum k=alphaout*p_ref+betaout*n_ref+qtout;
   k.rescaleMass();
   parton->set5Momentum(k);
   pf+=k;
   lastQ=q;
   lastx/=z;
   p += parton->momentum();
   thestep->addDecayProduct(rem,parton,false);
   return parton;
 }
 
 void HwRemDecayer::setRemMasses() const {
   // get the masses of the remnants
   Energy mrem[2];
   Lorentz5Momentum ptotal,pnew[2];
   vector<tRemPPtr> theprocessed;
   theprocessed.push_back(theRems.first);
   theprocessed.push_back(theRems.second);
   // one remnant in e.g. DIS
   if(!theprocessed[0]||!theprocessed[1]) {
     tRemPPtr rem = theprocessed[0] ? theprocessed[0] : theprocessed[1];
-    int remid = theprocessed[0] ? theContent.first.RemID() : theContent.second.RemID();
     Lorentz5Momentum deltap(rem->momentum());
     // find the diquark and momentum we still need in the energy
     tPPtr diquark;
     vector<PPtr> progenitors;
     for(unsigned int ix=0;ix<rem->children().size();++ix) {
-      if(rem->children()[ix]->id()!=remid) {
+      if(!DiquarkMatcher::Check(rem->children()[ix]->data())) {
 	progenitors.push_back(rem->children()[ix]);
 	deltap -= rem->children()[ix]->momentum();
       }
       else 
 	diquark = rem->children()[ix];
     }
     // now find the total momentum of the hadronic final-state to
     // reshuffle against
     // find the hadron for this remnant
     tPPtr hadron=rem;
     do hadron=hadron->parents()[0];
     while(!hadron->parents().empty());
     // find incoming parton to hard process from this hadron
     tPPtr hardin = 
       generator()->currentEvent()->primaryCollision()->incoming().first==hadron ?
       generator()->currentEvent()->primarySubProcess()->incoming().first :
       generator()->currentEvent()->primarySubProcess()->incoming().second;
     tPPtr parent=hardin;
     vector<PPtr> tempprog;
     // find the outgoing particles emitted from the backward shower
     do {
       assert(!parent->parents().empty());
       tPPtr newparent=parent->parents()[0];
       if(newparent==hadron) break;
       for(unsigned int ix=0;ix<newparent->children().size();++ix) {
 	if(newparent->children()[ix]!=parent)
 	  findChildren(newparent->children()[ix],tempprog);
       }
       parent=newparent;
     }
     while(parent!=hadron);
     // add to list of potential particles to reshuffle against in right order
     for(unsigned int ix=tempprog.size();ix>0;--ix) progenitors.push_back(tempprog[ix-1]);
     // final-state particles which are colour connected
     tColinePair lines = make_pair(hardin->colourLine(),hardin->antiColourLine());
     vector<PPtr> others;
     for(ParticleVector::const_iterator 
 	  cit = generator()->currentEvent()->primarySubProcess()->outgoing().begin();
 	cit!= generator()->currentEvent()->primarySubProcess()->outgoing().end();++cit) {
       // colour connected
       if(lines.first&&lines.first==(**cit).colourLine()) {
 	findChildren(*cit,progenitors);
 	continue;
       }
       // anticolour connected
       if(lines.second&&lines.second==(**cit).antiColourLine()) {
 	findChildren(*cit,progenitors);
 	continue;
       }
       // not connected
       for(unsigned int ix=0;ix<(**cit).children().size();++ix)
 	others.push_back((**cit).children()[ix]);
     }
     // work out how much of the system needed for rescaling
     unsigned int iloc=0;
     Lorentz5Momentum psystem,ptotal;
     do {
       psystem+=progenitors[iloc]->momentum();
       ptotal = psystem + deltap;
       ptotal.rescaleMass();
       psystem.rescaleMass();
       ++iloc;
       if(ptotal.mass() > psystem.mass() + diquark->mass() &&
 	 psystem.mass()>1*MeV && DISRemnantOpt_<2 && ptotal.e() > 0.*GeV ) break;
     }
     while(iloc<progenitors.size());
     if(ptotal.mass() > psystem.mass() + diquark->mass()) --iloc;
     if(iloc==progenitors.size()) {
       // try touching the lepton in dis as a last restort
       for(unsigned int ix=0;ix<others.size();++ix) {
 	progenitors.push_back(others[ix]);
 	psystem+=progenitors[iloc]->momentum();
 	ptotal = psystem + deltap;
 	ptotal.rescaleMass();
 	psystem.rescaleMass();
 	++iloc;
       }
       --iloc;
       if(ptotal.mass() > psystem.mass() + diquark->mass()) {
 	if(DISRemnantOpt_==0||DISRemnantOpt_==2)
 	  Throw<Exception>() << "Warning had to adjust the momentum of the"
 			     << " non-colour connected"
 			     << " final-state, e.g. the scattered lepton in DIS"
 			     << Exception::warning;
 	else
 	  throw Exception() << "Can't set remnant momentum without adjusting "
 			    << "the momentum of the"
 			    << " non-colour connected"
 			    << " final-state, e.g. the scattered lepton in DIS"
 			    << " vetoing event"
 			    << Exception::eventerror;
       }
       else {
 	throw Exception() << "Can't put the remnant on-shell in HwRemDecayer::setRemMasses()"
 			  << Exception::eventerror;
       }
     }
     psystem.rescaleMass();
     LorentzRotation R = Utilities::getBoostToCM(make_pair(psystem, deltap));
     Energy pz = SimplePhaseSpace::getMagnitude(sqr(ptotal.mass()), 
 					       psystem.mass(), diquark->mass());
     LorentzRotation Rs(-(R*psystem).boostVector());
     Rs.boost(0.0, 0.0, pz/sqrt(sqr(pz) + sqr(psystem.mass())));
     Rs = Rs*R;
     // put remnant on shell
     deltap.transform(R);
     deltap.setMass(diquark->mass());
     deltap.setE(sqrt(sqr(diquark->mass())+sqr(pz)));
     deltap.rescaleRho();
     R.invert();
     deltap.transform(R);
     Rs = R*Rs;
     // apply transformation to required particles to absorb recoil
     for(unsigned int ix=0;ix<=iloc;++ix) {
       progenitors[ix]->deepTransform(Rs);
     }
     diquark->set5Momentum(deltap);
   }
   // two remnants
   else {
     for(unsigned int ix=0;ix<2;++ix) {
       if(!theprocessed[ix]) continue;
       pnew[ix]=Lorentz5Momentum();
       for(unsigned int iy=0;iy<theprocessed[ix]->children().size();++iy) {
 	pnew[ix]+=theprocessed[ix]->children()[iy]->momentum();
       }
       mrem[ix]=sqrt(pnew[ix].m2());
     }
     // now find the remnant remnant cmf frame
     Lorentz5Momentum prem[2]={theprocessed[0]->momentum(),
 			      theprocessed[1]->momentum()};
     ptotal=prem[0]+prem[1];
     ptotal.rescaleMass();
     // boost momenta to this frame
     if(ptotal.m()< (pnew[0].m()+pnew[1].m())) 
       throw Exception() << "Not enough energy in both remnants in " 
 			<< "HwRemDecayer::setRemMasses() " 
 			<< Exception::eventerror;
     
     Boost boostv(-ptotal.boostVector());
     ptotal.boost(boostv);
     for(unsigned int ix=0;ix<2;++ix) {
       prem[ix].boost(boostv);
       // set the masses and energies,
       prem[ix].setMass(mrem[ix]);
       prem[ix].setE(0.5/ptotal.m()*(sqr(ptotal.m())+sqr(mrem[ix])-sqr(mrem[1-ix])));
       prem[ix].rescaleRho();
       // boost back to the lab
       prem[ix].boost(-boostv);
       // set the momenta of the remnants
       theprocessed[ix]->set5Momentum(prem[ix]);
     }
     // boost the decay products
     Lorentz5Momentum ptemp;
     for(unsigned int ix=0;ix<2;++ix) { 
       Boost btorest(-pnew[ix].boostVector()); 
       Boost bfmrest( prem[ix].boostVector()); 
       for(unsigned int iy=0;iy<theprocessed[ix]->children().size();++iy) { 
 	ptemp=theprocessed[ix]->children()[iy]->momentum(); 
 	ptemp.boost(btorest); 
 	ptemp.boost(bfmrest); 
 	theprocessed[ix]->children()[iy]->set5Momentum(ptemp); 
       }
     }
   }
 }
 
 void HwRemDecayer::initSoftInteractions(Energy ptmin, InvEnergy2 beta){
   ptmin_ = ptmin;
   beta_ = beta;
 }
 
+
 Energy HwRemDecayer::softPt() const {
   Energy2 pt2(ZERO);
   double xmin(0.0), xmax(1.0), x(0);
 
   if(beta_ == ZERO){
     return UseRandom::rnd(0.0,(double)(ptmin_/GeV))*GeV;
   }
-
   if(beta_ < ZERO){
     xmin = 1.0;
     xmax = exp( -beta_*sqr(ptmin_) );    
   }else{
     xmin = exp( -beta_*sqr(ptmin_) );    
     xmax = 1.0;
   }
   x = UseRandom::rnd(xmin, xmax);
   pt2 = 1.0/beta_ * log(1/x);
   
   if( pt2 < ZERO || pt2 > sqr(ptmin_) )
     throw Exception() << "HwRemDecayer::softPt generation of pt "
                       << "outside allowed range [0," << ptmin_/GeV << "]."
-                      << Exception::runerror; 
+                      << Exception::runerror;
+ 
+    //ofstream myfile2("softPt.txt", ios::app );
+    //myfile2 << pt2/GeV2 <<" "<<sqrt(pt2)/GeV<< endl;
+    //myfile2.close();
+
 
   return sqrt(pt2);
 }
 
 void HwRemDecayer::softKinematics(Lorentz5Momentum &r1, Lorentz5Momentum &r2, 
 				  Lorentz5Momentum &g1, Lorentz5Momentum &g2) const {
 
   g1 = Lorentz5Momentum();
   g2 = Lorentz5Momentum();
   //All necessary variables for the two soft gluons
   Energy pt(softPt()), pz(ZERO);
   Energy2 pz2(ZERO);
   double phi(UseRandom::rnd(2.*Constants::pi));
   double x_g1(0.0), x_g2(0.0);
   
   //Get the external momenta
   tcPPair beam(generator()->currentEventHandler()->currentCollision()->incoming());
   Lorentz5Momentum P1(beam.first->momentum()), P2(beam.second->momentum());
 
   if(dbg){
     cerr << "new event --------------------\n"
 	 << *(beam.first) << *(softRems_.first) 
 	 << "-------------------\n"
 	 << *(beam.second) << *(softRems_.second) << endl;
   }
   
   //parton mass
   Energy mp;
   if(quarkPair_){
   	mp = getParticleData(ParticleID::u)->constituentMass();
   }else{
   	mp = mg_;
   }
   
   //Get x_g1 and x_g2
   //first limits
   double xmin  = sqr(ptmin_)/4.0/(P1+P2).m2();
   double x1max = (r1.e()+abs(r1.z()))/(P1.e() + abs(P1.z()));
   double x2max = (r2.e()+abs(r2.z()))/(P2.e() + abs(P2.z()));
   double x1;
   
   if(!multiPeriph_){
   	//now generate according to 1/x
   	x_g1 = xmin * exp(UseRandom::rnd(log(x1max/xmin)));
   	x_g2 = xmin * exp(UseRandom::rnd(log(x2max/xmin)));
   }else{
   	if(valOfN_==0) return;
   	
   	double param = (1/(2*valOfN_+1))*initTotRap_;
   	do{
   		// need 1-x instead of x to get the proper final momenta
     		x1 = UseRandom::rndGauss(gaussWidth_, 1 - (exp(param)-1)/exp(param));
     	}while(x1 < 0 || x1>=1.0);
   	x_g1 = x1max*x1;
   	x_g2 = x2max*x1;
   }
   
 
   if(dbg)
     cerr << x1max << " " << x_g1 << endl << x2max << " " << x_g2 << endl;
 
   Lorentz5Momentum ig1, ig2, cmf;
 
   ig1 = x_g1*P1;
   ig2 = x_g2*P2;
-  
+
   ig1.setMass(mp);
   ig2.setMass(mp);
   ig1.rescaleEnergy();
   ig2.rescaleEnergy();
 
   cmf = ig1 + ig2;
   //boost vector from cmf to lab
   Boost boostv(cmf.boostVector());
-
   //outgoing gluons in cmf
   g1.setMass(mp);
   g2.setMass(mp);
 
   g1.setX(pt*cos(phi));
   g2.setX(-pt*cos(phi));
   g1.setY(pt*sin(phi));
   g2.setY(-pt*sin(phi));
-  
-  pz2 = cmf.m2()/4 - sqr(mp) - sqr(pt);
+  pz2 = cmf.m2()/4 - sqr(mp) - (pt*pt);
 
-  if(pz2/GeV2 < 0.0){
+  if( pz2/GeV2 < 0.0 ){
     if(dbg)
       cerr << "EXCEPTION not enough energy...." << endl;
     throw ExtraSoftScatterVeto();
   }
   
   if(!multiPeriph_){
-  	if(UseRandom::rndbool())
+  	if(UseRandom::rndbool()){
     		pz = sqrt(pz2);
-  	else
-    		pz = -sqrt(pz2);
+
+  	}else
+                pz = -sqrt(pz2);
   }else{
-  	pz = pz2 > ZERO ? sqrt(pz2) : ZERO;
+        pz = pz2 > ZERO ? sqrt(pz2) : ZERO;
+
   }
   
 
   if(dbg)
-    cerr << "pz has been calculated to: " << pz/GeV << endl;
+    cerr << "pz1 has been calculated to: " << pz/GeV << endl;
 
 
   g1.setZ(pz);
   g2.setZ(-pz);
   g1.rescaleEnergy();
   g2.rescaleEnergy();
 
   if(dbg){
     cerr << "check inv mass in cmf frame: " << (g1+g2).m()/GeV 
 	 << " vs. lab frame: " << (ig1+ig2).m()/GeV << endl;
   }
 
   g1.boost(boostv);
   g2.boost(boostv);
 
   //recalc the remnant momenta
   Lorentz5Momentum r1old(r1), r2old(r2);
-  
-  r1 -= ig1;
-  r2 -= ig2;
-  
+  r1 -= g1;
+  r2 -= g2;
+
   try{
     reShuffle(r1, r2, r1old.m(), r2old.m());
   }catch(ExtraSoftScatterVeto){
     r1 = r1old;
     r2 = r2old;
     throw ExtraSoftScatterVeto();
   }
 
   if(dbg){
     cerr << "remnant 1,2 momenta: " << r1/GeV << "--" << r2/GeV << endl;
     cerr << "remnant 1,2 masses: " << r1.m()/GeV << " " << r2.m()/GeV << endl;
     cerr << "check momenta in the lab..." << (-r1old-r2old+r1+r2+g1+g2)/GeV << endl;
   }
 }
 
 void HwRemDecayer::doSoftInteractions_old(unsigned int N) {
   if(N == 0) return;
   if(!softRems_.first || !softRems_.second)
     throw Exception() << "HwRemDecayer::doSoftInteractions: no "
                       << "Remnants available."
                       << Exception::runerror; 
 
   if( ptmin_ == -1.*GeV )
     throw Exception() << "HwRemDecayer::doSoftInteractions: init "
                       << "code has not been called! call initSoftInteractions."
                       << Exception::runerror; 
 
+
   Lorentz5Momentum g1, g2;
   Lorentz5Momentum r1(softRems_.first->momentum()), r2(softRems_.second->momentum());
   unsigned int tries(1), i(0);
 
   for(i=0; i<N; i++){
     //check how often this scattering has been regenerated
     if(tries > maxtrySoft_) break;
 
     if(dbg){
       cerr << "new try \n" << *softRems_.first << *softRems_.second << endl;
     }
 
     try{
       softKinematics(r1, r2, g1, g2);
 
     }catch(ExtraSoftScatterVeto){
       tries++;
       i--;
       continue;
     }
-  
+ 
+ 
     PPair oldrems = softRems_;
     PPair gluons = make_pair(addParticle(softRems_.first, ParticleID::g, g1), 
 			     addParticle(softRems_.second, ParticleID::g, g2));
     //now reset the remnants with the new ones
     softRems_.first = addParticle(softRems_.first, softRems_.first->id(), r1);
     softRems_.second = addParticle(softRems_.second, softRems_.second->id(), r2);
 
     //do the colour connections
     pair<bool, bool> anti = make_pair(oldrems.first->hasAntiColour(),
 				      oldrems.second->hasAntiColour());
 
     ColinePtr cl1 = new_ptr(ColourLine());
     ColinePtr cl2 = new_ptr(ColourLine());
-    if( UseRandom::rnd() < colourDisrupt_ ){//this is the member variable, i.e. SOFT colour disruption
-      //connect the remnants independent of the gluons
-      oldrems.first->colourLine(anti.first)->addColoured(softRems_.first, anti.first);
-      oldrems.second->colourLine(anti.second)->addColoured(softRems_.second, anti.second);
-      //connect the gluons to each other
-      cl1->addColoured(gluons.first);
-      cl1->addAntiColoured(gluons.second);
-      cl2->addColoured(gluons.second);
-      cl2->addAntiColoured(gluons.first);	
-    }else{
-      //connect the remnants to the gluons
-      oldrems.first->colourLine(anti.first)->addColoured(gluons.first, anti.first);
-      oldrems.second->colourLine(anti.second)->addColoured(gluons.second, anti.second);
-      //and the remaining colour line to the final remnant
-      cl1->addColoured(softRems_.first, anti.first);
+
+   //  case 2:
+oldrems.first->colourLine(anti.first)
+              ->addColoured(gluons.second,anti.second);
+cl2->addColoured(softRems_.first, anti.second);
+      cl2->addColoured(gluons.second, !anti.second);
+
+
+      oldrems.first->colourLine(anti.first)
+              ->addColoured(gluons.second,anti.second);
+      oldrems.second->colourLine(anti.second)
+              ->addColoured(gluons.first,anti.first);
+
+      cl1->addColoured(softRems_.second, anti.first);
       cl1->addColoured(gluons.first, !anti.first);
-      cl2->addColoured(softRems_.second, anti.second);
+      cl2->addColoured(softRems_.first, anti.second);
       cl2->addColoured(gluons.second, !anti.second);
-    }
-
     //reset counter
     tries = 1;
   }
 
   if(dbg)
     cerr << "generated " << i << "th soft scatters\n";
 
 }
 
+// Solve the reshuffling equation to rescale the remnant momenta
+double bisectReshuffling(const vector<PPtr>& particles,
+                         Energy w,
+                         double target = -16., double maxLevel = 80.) {
+
+  double level = 0;
+  double left = 0;
+  double right = 1;
+
+  double check = -1.;
+  double xi = -1;
+
+  while ( level < maxLevel ) {
+
+    xi = (left+right)*pow(0.5,level+1.);
+    check = 0.;
+    for (vector<PPtr>::const_iterator p = particles.begin(); p != particles.end(); ++p){
+      check += sqrt(sqr(xi)*((*p)->momentum().vect().mag2())+sqr((*p)->mass()))/w;
+    }
+
+    if ( log10(abs(1.-check)) <= target )
+      break;
+
+    left *= 2.;
+    right *= 2.;
+
+    if ( check >= 1. ) {
+      right -= 1.;
+      ++level;
+    }
+
+    if ( check < 1. ) {
+      left += 1.;
+      ++level;
+    }
+
+  }
+  return xi;
+
+}
+
+LorentzRotation HwRemDecayer::rotate(const LorentzMomentum &p) const {
+  LorentzRotation R;
+  static const double ptcut = 1e-20;
+  Energy2 pt2 = sqr(p.x())+sqr(p.y());
+  Energy2 pp2 = sqr(p.z())+pt2;
+  double phi, theta;
+  if(pt2 <= pp2*ptcut) {
+     if(p.z() > ZERO) theta = 0.;
+     else theta = Constants::pi;
+     phi = 0.;
+  } else {
+     Energy pp = sqrt(pp2);
+     Energy pt = sqrt(pt2);
+     double ct = p.z()/pp;
+     double cf = p.x()/pt;
+     phi = -acos(cf);
+     theta = acos(ct);
+  }
+  // Rotate first around the z axis to put p in the x-z plane
+  // Then rotate around the Y axis to put p on the z axis
+  R.rotateZ(phi).rotateY(theta);
+  return R;
+}
+
+struct vectorSort{
+  bool operator() (Lorentz5Momentum i,Lorentz5Momentum j) {return(i.rapidity() < j.rapidity());}
+} ySort;
+
+
+
 void HwRemDecayer::doSoftInteractions_multiPeriph(unsigned int N) {	
   if(N == 0) return;
   int Nmpi = N;
+
   for(int j=0;j<Nmpi;j++){
   ///////////////////////
-  // Average number of partons in the ladder
-  
+  // TODO: parametrization of the ladder multiplicity (need to tune to 900GeV, 7Tev and 13Tev) 
   // Parameterize the ladder multiplicity to: ladderMult_ = A_0 * (s/1TeV^2)^alpha  
   // with the two tunable parameters A_0 =ladderNorm_ and alpha = ladderPower_
 
   // Get the collision energy
   Energy energy(generator()->maximumCMEnergy());
 
   double reference = sqr(energy/TeV);
-  double ladderMult_;
 
+  // double ladderMult_;
   // Parametrization of the ladder multiplicity
-  ladderMult_ = ladderNorm_ * pow( ( reference ) , ladderPower_ );
+  // ladderMult_ = ladderNorm_ * pow( ( reference ) , ladderPower_ );
 
- 
   int avgN = floor(ladderMult_*log((softRems_.first->momentum()
   		+softRems_.second->momentum()).m()/mg_) + ladderbFactor_);
   initTotRap_ = abs(softRems_.first->momentum().rapidity())
   		+abs(softRems_.second->momentum().rapidity());
     
-  //generate the poisson distribution with mean avgN
+  // Generate the poisson distribution with mean avgN
   double L = exp(-double(avgN));
   int k = 0;
   double p = 1;
   do {
       k++;
       p *= UseRandom::rnd();
   } while( p > L);
   N=k-1;
   
   valOfN_=N;
-  if(N == 0) return;
+  if(N == 0){
+    continue;
+  }
   
   if(!softRems_.first || !softRems_.second)
     throw Exception() << "HwRemDecayer::doSoftInteractions: no "
                       << "Remnants available."
                       << Exception::runerror; 
 
   if( ptmin_ == -1.*GeV )
     throw Exception() << "HwRemDecayer::doSoftInteractions: init "
                       << "code has not been called! call initSoftInteractions."
                       << Exception::runerror; 
 
-  Lorentz5Momentum g1, g2, q1, q2;
-  //intermediate gluons; erased in if there is 
-  //another step
-  Lorentz5Momentum gint1, gint2;
+  // The remnants
+  PPtr rem1 = softRems_.first;
+  PPtr rem2 = softRems_.second;
+  // Vector for the ladder particles
+  vector<Lorentz5Momentum> ladderMomenta;
+  // Remnant momenta
+  Lorentz5Momentum r1(softRems_.first->momentum()), r2(softRems_.second->momentum());
+  Lorentz5Momentum cm =r1+r2;
+
+  // Initialize partons in the ladder
+  // The toy masses are needed for the correct calculation of the available energy
+  Lorentz5Momentum sumMomenta;
+  for(unsigned int i = 0; i<2*N; i++) {
+      // choose constituents
+      Energy newMass = ZERO;
+      Energy toyMass;
+      if(i<2){
+        // u and d have the same mass so its enough to use u 
+        toyMass = getParticleData(ParticleID::u)->constituentMass();
+      }
+      else{
+        toyMass = getParticleData(ParticleID::g)->constituentMass();
+      }
+      Lorentz5Momentum cp(ZERO,ZERO,ZERO,newMass,newMass);
+      // dummy container for the momentum that is used for momentum conservation
+      Lorentz5Momentum dummy(ZERO,ZERO,ZERO,toyMass,toyMass);
+      ladderMomenta.push_back(cp);
+      sumMomenta+=dummy;
+  }
+
+  // Get the beam energy
+  tcPPair beam(generator()->currentEventHandler()->currentCollision()->incoming());
+  Lorentz5Momentum P1(beam.first->momentum()), P2(beam.second->momentum());
+
+  // Calculate available energy for the partons
+  double x_1(0.0), x_2(0.0);
+  double x1max = (r1.e()+abs(r1.z()))/(P1.e() + abs(P1.z()));
+  double x2max = (r2.e()+abs(r2.z()))/(P2.e() + abs(P2.z()));
+  double x1;
+  double param = (1/(2*valOfN_+1))*initTotRap_;
+  do{
+       // Need 1-x instead of x to get the proper final momenta
+       x1 = UseRandom::rndGauss(gaussWidth_, 1 - (exp(param)-1)/exp(param));
+  }while(x1 < 0 || x1>=1.0);
+       x_1 = x1max*x1;
+       x_2 = x2max*x1;
+
+  // Remnants 1 and 2 need to be rescaled later by this amount
+  Lorentz5Momentum ig1 = x1*r1;
+  Lorentz5Momentum ig2 = x1*r2;
+
+  // The available energy that is used to generate the ladder
+  // sumMomenta is the the sum of rest masses of the ladder partons
+  // the available energy goes all into the kinematics
+  Energy availableEnergy = (ig1+ig2).m() - sumMomenta.m();
  
-  //momenta of the gluon pair generated in 
-  //each step
-  vector< pair<Lorentz5Momentum,Lorentz5Momentum> > gluonMomPairs;
-  //momenta of remnants
-  Lorentz5Momentum r1(softRems_.first->momentum()), r2(softRems_.second->momentum());
-  
-  unsigned int tries(1), i(0);
-  //generate the momenta of particles in the ladder
-  for(i = 0; i < N; i++){
-    //check how often this scattering has been regenerated
-    //and break if exeeded maximal number
-    if(tries > maxtrySoft_) break;
-    	
-    if(dbg){
-      cerr << "new try \n" << *softRems_.first << *softRems_.second << endl;
+  // If not enough energy then continue
+  if ( availableEnergy < ZERO ) {
+    continue;
+  }
+
+ unsigned int its(0); 
+  // Generate the momenta of the partons in the ladder
+  if ( !(doPhaseSpaceGenerationGluons(ladderMomenta,availableEnergy,its)) ){
+    continue;
+  }
+ // Add gluon mass and rescale
+ Lorentz5Momentum totalMomPartons;
+ Lorentz5Momentum totalMassLessPartons;
+
+ // Sort the ladder partons according to their rapidity and then choose which ones will be the quarks
+ sort(ladderMomenta.begin(),ladderMomenta.end(),ySort);
+
+  int countPartons=0;
+  long quarkID=0;
+  // Choose between up and down quarks
+  int choice = UseRandom::rnd2(1,1);
+  switch (choice) {
+    case 0: quarkID = ParticleID::u; break;
+    case 1: quarkID = ParticleID::d; break;
+  }
+
+  for (auto &p:ladderMomenta){
+    totalMomPartons+=p;
+    // Set the mass of the gluons and the two quarks in the ladder
+    if(countPartons==0 || countPartons==(ladderMomenta.size()-1)){
+      p.setMass( getParticleData(quarkID)->constituentMass() );
+    }else{
+      p.setMass( getParticleData(ParticleID::g)->constituentMass() );
+    }
+    p.rescaleEnergy();
+    countPartons++;
+  }
+
+  // Continue if energy conservation is violated 
+  if ( abs(availableEnergy - totalMomPartons.m()) > 1e-8*GeV){
+    continue;
+  }
+
+  // Boost momenta into CM frame
+  const Boost boostv(-totalMomPartons.boostVector());
+  Lorentz5Momentum totalMomentumAfterBoost;
+  for ( unsigned int i=0; i<ladderMomenta.size(); i++){
+    ladderMomenta[i].boost(boostv);
+    totalMomentumAfterBoost +=ladderMomenta[i];
+  }
+
+  const Boost boostvR(-cm.boostVector());
+  r1.boost(boostvR);
+  r2.boost(boostvR);
+
+  // Remainig energy after generation of soft ladder
+  Energy remainingEnergy = cm.m() - totalMomentumAfterBoost.m();
+
+  // Continue if not enough energy
+  if(remainingEnergy<ZERO) {
+    continue;
+  }
+
+  vector<PPtr> remnants;
+  rem1->set5Momentum(r1);
+  rem2->set5Momentum(r2);
+  remnants.push_back(rem1);
+  remnants.push_back(rem2);
+
+  vector<PPtr> reshuffledRemnants;
+  Lorentz5Momentum totalMomentumAll;
+
+  // Bisect reshuffling for rescaling of remnants
+  double xi_remnants = bisectReshuffling(remnants,remainingEnergy);
+
+  // Rescale remnants
+  for ( auto &rems: remnants ) {
+    Lorentz5Momentum reshuffledMomentum;
+    reshuffledMomentum = xi_remnants*rems->momentum();
+
+    reshuffledMomentum.setMass(getParticleData(softRems_.first->id())->constituentMass());
+    reshuffledMomentum.rescaleEnergy();
+    reshuffledMomentum.boost(-boostvR);
+    rems->set5Momentum(reshuffledMomentum);
+    totalMomentumAll+=reshuffledMomentum;
+  }
+  // Then the other particles   
+  for ( auto &p:ladderMomenta ) {
+        p.boost(-boostvR);
+        totalMomentumAll+=p;
+  }
+
+  // Do the colour connections
+  // Original rems are the ones which are connected to other parts of the event
+  PPair oldRems_ = softRems_;
+
+  pair<bool, bool> anti = make_pair(oldRems_.first->hasAntiColour(),
+                                      oldRems_.second->hasAntiColour());
+
+  // Replace first remnant
+  softRems_.first = addParticle(softRems_.first, softRems_.first->id(),
+                               remnants[0]->momentum());
+
+  // Switch for random connections
+  if (randomRemnantConnection_){
+    disrupt_ = UseRandom::rnd3(1,1,1);
+  }
+
+  // Connect the old remnant to the new remnant
+  if(disrupt_==0){
+    oldRems_.first->colourLine(anti.first)->addColoured(softRems_.first, anti.first);
+  }
+  // Replace second remnant
+  softRems_.second = addParticle(softRems_.second, softRems_.second->id(),
+                                remnants[1]->momentum());
+  // This connects the old remnants to the new remnants
+  if(disrupt_==0){
+    oldRems_.second->colourLine(anti.second)->addColoured(softRems_.second, anti.second);
+  }
+  // Add all partons to the first remnant for the event record
+  vector<PPtr> partons;
+  int count=0;
+
+  // Choose the colour connections and position of quark antiquark
+  // Choose between R1-q-g..g-qbar-R2 or R1-qbar-g...g-q-R2 
+  // (place of quark antiquarks in the ladder)
+  int quarkPosition = UseRandom::rnd2(1,1);
+
+  for (auto &p:ladderMomenta){
+
+    if(p.mass()==getParticleData(ParticleID::u)->constituentMass()){ 
+      if(count==0){
+        if(quarkPosition==0){
+          partons.push_back(addParticle(softRems_.first, quarkID, p));
+          count++;
+        }else{
+          partons.push_back(addParticle(softRems_.first, -quarkID, p));
+          count++;          
+        }
+      }else{
+        if(quarkPosition==0){
+          partons.push_back(addParticle(softRems_.first, -quarkID, p));
+        }else{
+          partons.push_back(addParticle(softRems_.first, quarkID, p));
+        }
+      }
+  }else{
+        partons.push_back(addParticle(softRems_.first, ParticleID::g, p));
+  }
+  softRems_.first = addParticle(softRems_.first, softRems_.first->id(),
+                                        softRems_.first->momentum());
+
+
+  if (disrupt_==0) {
+   oldRems_.first->colourLine(anti.first)->addColoured(softRems_.first, anti.first);
+}
+  }
+
+
+  // Need to differenciate between the two quark positions, this defines the 
+  // colour connections to the new remnants and old remnants
+  if(quarkPosition==0){
+  switch(disrupt_){
+    case 0:
+        {
+        // ladder self contained
+        if(partons.size()==2){
+          ColinePtr clq =  new_ptr(ColourLine());
+          clq->addColoured(partons[0]);
+          clq->addAntiColoured(partons[1]);
+         }
+
+         ColinePtr clfirst =  new_ptr(ColourLine());
+         ColinePtr cllast =  new_ptr(ColourLine());
+
+         if(partons.size()>2){
+           clfirst->addColoured(partons[0]);
+    	   clfirst->addAntiColoured(partons[1]);
+    	   cllast->addAntiColoured(partons[partons.size()-1]);
+           cllast->addColoured(partons[partons.size()-2]);
+           //now the remaining gluons
+           for (unsigned int i=1; i<partons.size()-2; i++){
+             ColinePtr cl = new_ptr(ColourLine());
+             cl->addColoured(partons[i]);
+             cl->addAntiColoured(partons[i+1]);
+           }
+         }
+         break;
+         }
+    case 1:
+         {
+         // Connect remnants with ladder
+         // Case 1 (only 2 quarks)
+         if(partons.size()==2){
+           ColinePtr clq =  new_ptr(ColourLine());
+           clq->addColoured(partons[0]);
+           clq->addAntiColoured(partons[1]);
+           // and connect remnants to old rems
+              oldRems_.first->colourLine(anti.first)
+                          ->addColoured(softRems_.first, anti.first);
+              oldRems_.second->colourLine(anti.second)
+                           ->addColoured(softRems_.second, anti.second);
+         } else {
+         // Case 2
+         // New remnant 1 with first quark in ladder
+         ColinePtr cl1 = new_ptr(ColourLine());
+         cl1->addColoured(softRems_.first,anti.first);
+         cl1->addColoured(partons[0]);     
+         // Connect old rems to the first gluon
+         oldRems_.first->colourLine(anti.first)
+              ->addAntiColoured(partons[1]);
+
+         // Connect gluons with each other
+         for (unsigned int i=1; i<partons.size()-2; i++){
+            ColinePtr cl = new_ptr(ColourLine());
+            cl->addColoured(partons[i]);
+            cl->addAntiColoured(partons[i+1]);
+
+         }
+         // Connect remnant 2 with with last gluon
+         ColinePtr cl2 = new_ptr(ColourLine());
+
+         cl2->addColoured(softRems_.second,anti.second);
+         cl2->addColoured(partons[partons.size()-2]);
+
+         // Connect old remnant 2 with antiquark
+         oldRems_.second->colourLine(anti.first)
+              ->addAntiColoured(partons[partons.size()-1]);
+
+         }
+         break;
+         }
+      case 2:
+         {
+         if(partons.size()==2){
+           ColinePtr clq =  new_ptr(ColourLine());
+           clq->addColoured(partons[0]);
+           clq->addAntiColoured(partons[1]);
+              oldRems_.first->colourLine(anti.first)
+                          ->addColoured(softRems_.first, anti.first);
+              oldRems_.second->colourLine(anti.second)
+                           ->addColoured(softRems_.second, anti.second);
+         }else{
+         // NNew remnant 1 with first gluon in ladder
+         ColinePtr cl1 = new_ptr(ColourLine());
+         cl1->addColoured(softRems_.first,anti.first);
+         cl1->addColoured(partons[1],!anti.first);
+         // Connect old rems to the first gluon as well
+         oldRems_.first->colourLine(anti.first)
+              ->addColoured(partons[1],anti.first);
+         // Gluons with each other
+         // Connect quark with second gluon
+         ColinePtr clq1 = new_ptr(ColourLine());
+         clq1->addColoured(partons[0]);
+         clq1->addAntiColoured(partons[2]);
+         // First gluon already is handled
+         for (unsigned int i=2; i<partons.size()-2; i++){
+            ColinePtr cl = new_ptr(ColourLine());
+            cl->addColoured(partons[i]);
+            cl->addAntiColoured(partons[i+1]);
+         }
+         // Connect remnant 2 with with last gluon
+         ColinePtr cl2 = new_ptr(ColourLine());
+         cl2->addColoured(softRems_.second,anti.second);
+         cl2->addColoured(partons[partons.size()-2],!anti.second);
+
+         // Connect original remnant 2 with antiquark
+         oldRems_.second->colourLine(anti.first)
+              ->addColoured(partons[partons.size()-1],anti.second);
+         }
+         break;
+         }
+    }
+  } else {
+   switch(disrupt_){
+     case 0:
+        {
+        if(partons.size()==2){
+          ColinePtr clq =  new_ptr(ColourLine());
+          clq->addAntiColoured(partons[0]);
+          clq->addColoured(partons[1]);
+         }
+
+         ColinePtr clfirst =  new_ptr(ColourLine());
+         ColinePtr cllast =  new_ptr(ColourLine());
+
+         if(partons.size()>2){
+           clfirst->addAntiColoured(partons[0]);
+    	   clfirst->addColoured(partons[1]);
+    	   cllast->addColoured(partons[partons.size()-1]);
+           cllast->addAntiColoured(partons[partons.size()-2]);
+           //now the remaining gluons
+           for (unsigned int i=1; i<partons.size()-2; i++){
+             ColinePtr cl = new_ptr(ColourLine());
+             cl->addAntiColoured(partons[i]);
+             cl->addColoured(partons[i+1]);
+           }
+         }
+         break;
+         }
+    case 1:
+         {
+         //Case 1 (only 2 quarks)
+         if(partons.size()==2){
+           ColinePtr clq =  new_ptr(ColourLine());
+           clq->addAntiColoured(partons[0]);
+           clq->addColoured(partons[1]);
+           // and connect remnants to old rems
+              oldRems_.first->colourLine(anti.first)
+                          ->addColoured(softRems_.first, anti.first);
+              oldRems_.second->colourLine(anti.second)
+                           ->addColoured(softRems_.second, anti.second);
+         } else {
+         //Case 2
+         //new remnant 1 with first gluon in ladder
+         ColinePtr cl1 = new_ptr(ColourLine());
+         cl1->addColoured(softRems_.first,anti.first);
+         cl1->addColoured(partons[1],!anti.first);     
+         //Connect old rems to the first antiquark 
+         oldRems_.first->colourLine(anti.first)
+              ->addColoured(partons[0],anti.first);
+         //gluons with each other
+         for (unsigned int i=1; i<partons.size()-2; i++){
+            ColinePtr cl = new_ptr(ColourLine());
+            cl->addAntiColoured(partons[i]);
+            cl->addColoured(partons[i+1]);
+         }
+         //Connect remnant 2 with with last quark
+         ColinePtr cl2 = new_ptr(ColourLine());
+
+         cl2->addColoured(softRems_.second,anti.second);
+         cl2->addColoured(partons[partons.size()-1],!anti.second);
+        
+         //Connect original remnant with first gluon
+         oldRems_.second->colourLine(anti.first)
+              ->addColoured(partons[partons.size()-2],anti.second);
+         }
+         break;
+         }
+      case 2:
+         {
+	 // only two quarks case
+         if(partons.size()==2){
+           ColinePtr clq =  new_ptr(ColourLine());
+           clq->addAntiColoured(partons[0]);
+           clq->addColoured(partons[1]);
+           // and connect remnants to old rems
+              oldRems_.first->colourLine(anti.first)
+                          ->addColoured(softRems_.first, anti.first);
+              oldRems_.second->colourLine(anti.second)
+                           ->addColoured(softRems_.second, anti.second);
+         }else{
+         //new remnant 1 with first gluon in ladder
+         ColinePtr cl1 = new_ptr(ColourLine());
+         cl1->addColoured(softRems_.first,anti.first);
+         cl1->addColoured(partons[1],!anti.first);
+         //Connect old rems to the first antiquark
+         oldRems_.first->colourLine(anti.first)
+              ->addColoured(partons[0],anti.first);
+         //gluons with each other
+         // connect quark with second gluon
+         ColinePtr clq1 = new_ptr(ColourLine());
+         clq1->addColoured(partons[partons.size()-1]);
+         clq1->addAntiColoured(partons[partons.size()-2]);
+         // first gluon already is handled
+         for (unsigned int i=1; i<partons.size()-2; i++){
+            ColinePtr cl = new_ptr(ColourLine());
+            cl->addAntiColoured(partons[i]);
+            cl->addColoured(partons[i+1]);
+         }
+         //Connect remnant 2 with with last gluon
+         ColinePtr cl2 = new_ptr(ColourLine());
+         cl2->addColoured(softRems_.second,anti.second);
+         cl2->addColoured(partons[partons.size()-2],!anti.second);
+
+         //Connect original remnant 2 with antiquark
+         oldRems_.second->colourLine(anti.first)
+              ->addColoured(partons[partons.size()-3],anti.second);
+         }
+         break;
+         }
+  }
+
+  } 
+
+ }
+
+ 
+}
+// Do the phase space generation here is 1 to 1 the same from UA5 model
+bool HwRemDecayer::doPhaseSpaceGenerationGluons(vector<Lorentz5Momentum> &softGluons, Energy CME, unsigned int &its)
+    const{
+
+  // Define the parameters
+  unsigned int _maxtries = 300;
+
+  double alog = log(CME*CME/GeV2);
+  unsigned int ncl = softGluons.size();
+  // calculate the slope parameters for the different clusters
+  // outside loop to save time
+  vector<Lorentz5Momentum> mom(ncl);
+
+  // Sets the slopes depending on the constituent quarks of the cluster
+  for(unsigned int ix=0;ix<ncl;++ix)
+    { 
+      mom[ix]=softGluons[ix];
     }
 
-    try{
-      if(i==0){
-      	//generated partons in the ladder are quark-antiquark
-      	quarkPair_ = true;
-      	//first splitting: remnant -> remnant + quark
-      	softKinematics(r1, r2, q1, q2);
-      }else if(i==1){
-        
-        //generated pair is gluon-gluon
-      	quarkPair_ = false;
-      	//second splitting; quark -> quark + gluon
-      	softKinematics(q1, q2, g1, g2);
-      	//record generated gluon pair
-      	gluonMomPairs.push_back(make_pair(g1,g2));
-      	
+  // generate the momenta
+  double eps = 1e-10/double(ncl);
+  vector<double> xi(ncl);
+  vector<Energy> tempEnergy(ncl);
+  Energy sum1(ZERO);
+  double yy(0.);
+  while(its < _maxtries) {
+    ++its;
+    Energy sumx = ZERO;
+    Energy sumy = ZERO;
+   unsigned int iterations(0);
+   unsigned int _maxtriesNew = 100;
+
+   while(iterations < _maxtriesNew) {
+    iterations++;
+    Energy sumxIt = ZERO;
+    Energy sumyIt = ZERO;
+    bool success=false;
+    Energy pTmax=ZERO;
+    Energy firstPt=ZERO;
+    for(unsigned int i = 0; i<ncl; ++i) {
+      // Different options for soft pt sampling
+      //1) pT1>pT2...pTN
+      //2) pT1>pT2>..>pTN
+      //3) flat
+      //4) y dependent
+      //5) Frist then flat
+      int triesPt=0;
+      Energy pt;
+      Energy ptTest;
+      switch(PtDistribution_) {
+        case 0: //default softPt()
+           pt=softPt();
+           break;
+        case 1: //pTordered
+           if(i==0){
+             pt=softPt();
+             pTmax=pt;
+           }else{
+             do{
+              pt=softPt();
+             }while(pt>pTmax);
+           }
+           break;
+         case 2: //strongly pT ordered
+             if ( i==0 ) {
+               pt=softPt();
+               pTmax=pt;
+              } else {
+                do {
+                  if ( triesPt==20 ) { 
+                    pt=pTmax;
+                    break;
+                  }
+                  pt=softPt();
+                  triesPt++;
+                } while ( pt>pTmax );
+                pTmax=pt;
+           }
+           break;
+         case 3: //flat
+            pt = UseRandom::rnd(0.0,(double)(ptmin_/GeV))*GeV;
+            break;
+         case 4: //flat below first pT
+            if ( i==0 ) {
+              pt = softPt();
+              firstPt = pt;
+            } else {
+              pt = firstPt * UseRandom::rnd();
+            }
+            break;
+         case 5: //flat but rising below first pT
+            if ( i==0 ) {
+		pt=softPt();
+                firstPt = pt;
+            } else {
+                  pt = firstPt * pow(UseRandom::rnd(),1/2);
+            }
+            
+ 
+      }
+ 
+      Energy2 ptp = pt*pt;
+      if(ptp <= ZERO) pt = - sqrt(-ptp);
+      else pt = sqrt(ptp);
+      // randomize azimuth
+      Energy px,py;
+      //randomize the azimuth, but the last one should cancel all others
+      if(i<ncl-1){
+       randAzm(pt,px,py);
+       // set transverse momentum
+       mom[i].setX(px);
+       mom[i].setY(py);
+       sumxIt += px;
+       sumyIt += py;
       }else{
-      	//consequent splittings gluon -> gluon+gluon
-      	//but, the previous gluon gets deleted
-      	quarkPair_ = false;
-      	//save first the previous gluon momentum
-      	g1=gluonMomPairs.back().first;
-      	g2=gluonMomPairs.back().second;
-      	//split gluon momentum
-      	softKinematics(g1, g2, gint1, gint2);
-      	
-      	//erase the last entry
-      	gluonMomPairs.pop_back();
-      	
-      	//add new gluons
-      	gluonMomPairs.push_back(make_pair(g1,g2));
-      	gluonMomPairs.push_back(make_pair(gint1,gint2));  	
+       //calculate azimuth angle s.t 
+      // double factor;
+      Energy pTdummy;
+      pTdummy = sqrt(sumxIt*sumxIt+sumyIt*sumyIt);
+      if( pTdummy < ptmin_ ){
+        px=-sumxIt;
+        py=-sumyIt;
+        mom[i].setX(px);
+        mom[i].setY(py);
+      
+        sumxIt+=px;
+        sumyIt+=py;
+        sumx = sumxIt;
+        sumy = sumyIt;
+        success=true;
       }
       
-    }catch(ExtraSoftScatterVeto){
-      tries++;
-      i--;
-      continue;
+      }
     }
-    //reset counter
-    tries = 1;
+   if(success){
+     break;
+   }
+   }
+    sumx /= ncl;
+    sumy /= ncl;
+    // find the sum of the transverse mass
+    Energy sumtm=ZERO;
+    for(unsigned int ix = 0; ix<ncl; ++ix) {
+      mom[ix].setX(mom[ix].x()-sumx);
+      mom[ix].setY(mom[ix].y()-sumy);
+      Energy2 pt2 = mom[ix].perp2();
+      // Use the z component of the clusters momentum for temporary storage
+      mom[ix].setZ(sqrt(pt2+mom[ix].mass2()));
+      sumtm += mom[ix].z();
+    }
+    // if transverse mass greater the CMS try again
+    if(sumtm > CME) continue;
+    for(unsigned int i = 0; i<ncl; i++) xi[i] = randUng(0.6,1.0);
+    // sort into ascending order
+    sort(xi.begin(), xi.end());
+    double ximin = xi[0];
+    double ximax = xi.back()-ximin;
+    xi[0] = 0.;
+    for(unsigned int i = ncl-2; i>=1; i--) xi[i+1] = (xi[i]-ximin)/ximax;
+    xi[1] = 1.;
+    yy= log(CME*CME/(mom[0].z()*mom[1].z()));
+    bool suceeded=false;
+    Energy sum2,sum3,sum4;
+    for(unsigned int j = 0; j<10; j++) {
+      sum1 = sum2 = sum3 = sum4 = ZERO;
+      for(unsigned int i = 0; i<ncl; i++) {
+        Energy tm = mom[i].z();
+        double ex = exp(yy*xi[i]);
+        sum1 += tm*ex;
+        sum2 += tm/ex;
+        sum3 += (tm*ex)*xi[i];
+        sum4 += (tm/ex)*xi[i];
+      }
+      double fy = alog-log(sum1*sum2/GeV2);
+      double dd = (sum3*sum2 - sum1*sum4)/(sum1*sum2);
+      double dyy = fy/dd;
+      if(abs(dyy/yy) < eps) {
+        yy += dyy;
+        suceeded=true;
+        break;
+      }
+      yy += dyy;
+    }
+    if(suceeded){
+       break;
+    }
+    if(its > 100) eps *= 10.;
   }
-  //if no gluons discard the ladder
-  if(gluonMomPairs.size()==0) return;
-  
-  if(dbg)
-    cerr << "generated " << i << "th soft scatters\n";
-  
-  //gluons from previous step
-  PPair oldgluons;  				      
-  //quark-antiquark pair
-  PPair quarks;
-  //colour direction of quark-antiquark (true if anticolour)
-  pair<bool, bool> anti_q;
-  //generate particles (quark-antiquark and gluons) in the 
-  //ladder from momenta generated above
-  for(i = 0; i <= gluonMomPairs.size(); i++){
-    //remnants before splitting
-    PPair oldrems = softRems_;
-    //current gluon pair
-    PPair gluons;
-    
-    //add quarks
-    if(i==0){
-    	quarks = make_pair(addParticle(softRems_.first, ParticleID::u, q1), 
-    			     addParticle(softRems_.second, ParticleID::ubar, q2));
-    	anti_q = make_pair(quarks.first->hasAntiColour(),
-    				      quarks.second->hasAntiColour());
-    }
-    			     
-    
-    
-    if(i>0){
-    	
-    	gluons = make_pair(addParticle(softRems_.first, ParticleID::g, 
-    				gluonMomPairs[i-1].first), 
-			     addParticle(softRems_.second, ParticleID::g, 
-			     	gluonMomPairs[i-1].second));
-    }
-    
-    //now reset the remnants with the new ones
-    softRems_.first = addParticle(softRems_.first, softRems_.first->id(), r1);
-    softRems_.second = addParticle(softRems_.second, softRems_.second->id(), r2);
-    
-    
-    //colour direction of remnats
-    pair<bool, bool> anti = make_pair(oldrems.first->hasAntiColour(),
-    				      oldrems.second->hasAntiColour());
-    
-    
-    
-    ColinePtr cl1 = new_ptr(ColourLine());
-    ColinePtr cl2 = new_ptr(ColourLine());
-    //first colour connect remnants and if no 
-    //gluons the quark-antiquark pair 
-      if(i==0){
-      		oldrems.first->colourLine(anti.first)->addColoured(softRems_.first, 
-      								anti.first);
-      		oldrems.second->colourLine(anti.second)->addColoured(softRems_.second, 
-      								anti.second);
-      		if(gluonMomPairs.size()==0){
-      			ColinePtr clq = new_ptr(ColourLine());
-      			clq->addColoured(quarks.first, anti_q.first);
-      			clq->addColoured(quarks.second, anti_q.second);
-      		}	
-      }//colour connect remnants from previous step and gluons to quarks or gluons
-      if(i!=0){
-      		oldrems.first->colourLine(anti.first)->addColoured(softRems_.first, 
-      								anti.first);
-      		oldrems.second->colourLine(anti.second)->addColoured(softRems_.second, 
-      								anti.second);
-      		if(i==1){
-      			cl1->addColoured(quarks.first, anti_q.first);
-      			cl1->addColoured(gluons.first, !anti_q.first);
-      			cl2->addColoured(quarks.second, anti_q.second);
-      			cl2->addColoured(gluons.second, !anti_q.second);
-      		}else{
-      			cl1->addColoured(oldgluons.first, anti_q.first);
-      			cl1->addColoured(gluons.first, !anti_q.first);
-      			cl2->addColoured(oldgluons.second, anti_q.second);
-      			cl2->addColoured(gluons.second, !anti_q.second);
-      		}
-      		
-      } //last step; connect last gluons 
-      if(i > 0 && i == gluonMomPairs.size()){
-      		ColinePtr clg = new_ptr(ColourLine());
-      		clg->addColoured(gluons.first, anti_q.first);
-      		clg->addColoured(gluons.second, anti_q.second);
-      }
-      
-      
-      //save gluons for the next step	    
-      if(i < gluonMomPairs.size()) oldgluons = gluons;
+  if(its==_maxtries){
+    return false;
   }
-  
-  ////////
+//    throw Exception() << "Can't generate soft underlying event in "
+//                      << "UA5Handler::generateCylindricalPS"
+//                      << Exception::eventerror;
+  double zz = log(CME/sum1);
+  for(unsigned int i = 0; i<ncl; i++) {
+    Energy tm = mom[i].z();
+    double E1 = exp(zz + yy*xi[i]);
+    mom[i].setZ(0.5*tm*(1./E1-E1));
+    mom[i].setE( 0.5*tm*(1./E1+E1));
+    softGluons[i]=mom[i];
+
   }
-  ////////
+  return true;
 }
 
+
 void HwRemDecayer::finalize(double colourDisrupt, unsigned int softInt){
   PPair diquarks;
   //Do the final Rem->Diquark or Rem->quark "decay"
   if(theRems.first) {
     diquarks.first = finalSplit(theRems.first, theContent.first.RemID(), 
 				theUsed.first);
     theMaps.first.push_back(make_pair(diquarks.first, tPPtr()));
   }
   if(theRems.second) {
     diquarks.second = finalSplit(theRems.second, theContent.second.RemID(), 
 				 theUsed.second);
     theMaps.second.push_back(make_pair(diquarks.second, tPPtr()));
   }
   setRemMasses();
   if(theRems.first) {
     fixColours(theMaps.first, theanti.first, colourDisrupt);
     if(theContent.first.hadron->id()==ParticleID::pomeron&&
        pomeronStructure_==0) fixColours(theMaps.first, !theanti.first, colourDisrupt);
   }
   if(theRems.second) {
     fixColours(theMaps.second, theanti.second, colourDisrupt);
     if(theContent.second.hadron->id()==ParticleID::pomeron&&
        pomeronStructure_==0) fixColours(theMaps.second, !theanti.second, colourDisrupt);
   }
 
   if( !theRems.first || !theRems.second ) return;
   //stop here if we don't have two remnants
   softRems_ = diquarks;
   doSoftInteractions(softInt);
 }
 
 
 HwRemDecayer::HadronContent
 HwRemDecayer::getHadronContent(tcPPtr hadron) const {
   HadronContent hc;
   hc.hadron = hadron->dataPtr();
   long id(hadron->id());
   // baryon
   if(BaryonMatcher::Check(hadron->data())) {
     hc.sign = id < 0? -1: 1;
     hc.flav.push_back((id = abs(id)/10)%10);
     hc.flav.push_back((id /= 10)%10);
     hc.flav.push_back((id /= 10)%10);
     hc.extracted = -1;
   }
   else if(hadron->data().id()==ParticleID::gamma ||
 	  (hadron->data().id()==ParticleID::pomeron && pomeronStructure_==1)) {
     hc.sign = 1;
     for(int ix=1;ix<6;++ix) {
       hc.flav.push_back( ix);
       hc.flav.push_back(-ix);
     }
   }
   else if(hadron->data().id()==ParticleID::pomeron ) {
     hc.sign = 1;
     hc.flav.push_back(ParticleID::g);
     hc.flav.push_back(ParticleID::g);
   }
   else if(hadron->data().id()==ParticleID::reggeon ) {
     hc.sign = 1;
     for(int ix=1;ix<3;++ix) {
       hc.flav.push_back( ix);
       hc.flav.push_back(-ix);
     }
   }
   hc.pomeronStructure = pomeronStructure_;
   return hc;
 }
 
 long HwRemDecayer::HadronContent::RemID() const{
   if(extracted == -1)
     throw Exception() << "Try to build a Diquark id without "
 		      << "having extracted something in "
 		      << "HwRemDecayer::RemID(...)"
 		      << Exception::runerror;
   //the hadron was a meson or photon
   if(flav.size()==2) return sign*flav[(extracted+1)%2];
  
   long remId;
   int id1(sign*flav[(extracted+1)%3]), 
     id2(sign*flav[(extracted+2)%3]),
     sign(0), spin(0);
 
   if (abs(id1) > abs(id2)) swap(id1, id2);
   sign = (id1 < 0) ? -1 : 1; // Needed for the spin 0/1 part
   remId = id2*1000+id1*100;
   // Now decide if we have spin 0 diquark or spin 1 diquark 
   if(id1 == id2) spin = 3; // spin 1                  
   else spin = 1; // otherwise spin 0  
   remId += sign*spin;  
   return remId;
 }
 
 
 tPPtr HwRemDecayer::addParticle(tcPPtr parent, long id, Lorentz5Momentum p) const {
 
   PPtr newp = new_ptr(Particle(getParticleData(id)));
   newp->set5Momentum(p);
   // Add the new remnant to the step, but don't do colour connections
   thestep->addDecayProduct(parent,newp,false);
   return newp;
 }
 
 void HwRemDecayer::findChildren(tPPtr part,vector<PPtr> & particles) const {
   if(part->children().empty()) particles.push_back(part);
   else {
     for(unsigned int ix=0;ix<part->children().size();++ix)
       findChildren(part->children()[ix],particles);
   }
 }
 
 ParticleVector HwRemDecayer::decay(const DecayMode &, 
 				   const Particle &, Step &) const {
   throw Exception() << "HwRemDecayer::decay(...) "
 		    << "must not be called explicitely."
 		    << Exception::runerror;
 }
 
 void HwRemDecayer::persistentOutput(PersistentOStream & os) const {
   os << ounit(_kinCutoff, GeV) << _range << _zbin << _ybin 
      << _nbinmax << _alphaS << _alphaEM << DISRemnantOpt_
-     << maxtrySoft_ << colourDisrupt_ << ladderPower_<< ladderNorm_ << ladderbFactor_ << pomeronStructure_
+     << maxtrySoft_ << colourDisrupt_ << ladderPower_<< ladderNorm_ << ladderMult_ << ladderbFactor_ << pomeronStructure_
      << ounit(mg_,GeV) << ounit(ptmin_,GeV) << ounit(beta_,sqr(InvGeV))
-     << allowTop_ << multiPeriph_ << valOfN_ << initTotRap_;
+     << allowTop_ << multiPeriph_ << valOfN_ << initTotRap_ << PtDistribution_
+     << disrupt_ << randomRemnantConnection_;
 }
 
 void HwRemDecayer::persistentInput(PersistentIStream & is, int) {
   is >> iunit(_kinCutoff, GeV) >> _range >> _zbin >> _ybin 
      >> _nbinmax >> _alphaS >> _alphaEM >> DISRemnantOpt_
-     >> maxtrySoft_ >> colourDisrupt_ >> ladderPower_ >> ladderNorm_ >> ladderbFactor_ >> pomeronStructure_
+     >> maxtrySoft_ >> colourDisrupt_ >> ladderPower_ >> ladderNorm_ >> ladderMult_ >> ladderbFactor_ >> pomeronStructure_
      >> iunit(mg_,GeV) >> iunit(ptmin_,GeV) >> iunit(beta_,sqr(InvGeV))
-     >> allowTop_ >> multiPeriph_ >> valOfN_ >> initTotRap_;
+     >> allowTop_ >> multiPeriph_ >> valOfN_ >> initTotRap_ >> PtDistribution_
+     >> disrupt_ >> randomRemnantConnection_;
 }
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<HwRemDecayer,RemnantDecayer>
 describeHerwigHwRemDecayer("Herwig::HwRemDecayer", "HwShower.so");
 
 void HwRemDecayer::Init() {
 
   static ClassDocumentation<HwRemDecayer> documentation
     ("The HwRemDecayer class decays the remnant for Herwig");
 
   static Parameter<HwRemDecayer,double> interfaceZBinSize
     ("ZBinSize",
      "The size of the vbins in z for the interpolation of the splitting function.",
      &HwRemDecayer::_zbin, 0.05, 0.001, 0.1,
      false, false, Interface::limited);
 
   static Parameter<HwRemDecayer,int> interfaceMaxBin
     ("MaxBin",
      "Maximum number of z bins",
      &HwRemDecayer::_nbinmax, 100, 10, 1000,
      false, false, Interface::limited);
 
   static Reference<HwRemDecayer,ShowerAlpha> interfaceAlphaS
     ("AlphaS",
      "Pointer to object to calculate the strong coupling",
      &HwRemDecayer::_alphaS, false, false, true, false, false);
 
   static Reference<HwRemDecayer,ShowerAlpha> interfaceAlphaEM
     ("AlphaEM",
      "Pointer to object to calculate the electromagnetic coupling",
      &HwRemDecayer::_alphaEM, false, false, true, false, false);
 
   static Parameter<HwRemDecayer,Energy> interfaceKinCutoff
     ("KinCutoff",
      "Parameter kinCutoff used to constrain qtilde",
      &HwRemDecayer::_kinCutoff, GeV, 0.75*GeV, 0.5*GeV, 10.0*GeV,
      false, false, Interface::limited);
 
   static Parameter<HwRemDecayer,double> interfaceEmissionRange
     ("EmissionRange",
      "Factor above the minimum possible value in which the forced splitting is allowed.",
      &HwRemDecayer::_range, 1.1, 1.0, 10.0,
      false, false, Interface::limited);
 
 
   static Switch<HwRemDecayer,unsigned int> interfaceDISRemnantOption
     ("DISRemnantOption",
      "Options for the treatment of the remnant in DIS",
      &HwRemDecayer::DISRemnantOpt_, 0, false, false);
   static SwitchOption interfaceDISRemnantOptionDefault
     (interfaceDISRemnantOption,
      "Default",
      "Use the minimum number of particles needed to take the recoil"
      " and allow the lepton to be used if needed",
      0);
   static SwitchOption interfaceDISRemnantOptionNoLepton
     (interfaceDISRemnantOption,
      "NoLepton",
      "Use the minimum number of particles needed to take the recoil but"
      " veto events where the lepton kinematics would need to be altered",
      1);
   static SwitchOption interfaceDISRemnantOptionAllParticles
     (interfaceDISRemnantOption,
      "AllParticles",
      "Use all particles in the colour connected system to take the recoil"
      " and use the lepton if needed.",
      2);
   static SwitchOption interfaceDISRemnantOptionAllParticlesNoLepton
     (interfaceDISRemnantOption,
      "AllParticlesNoLepton",
      "Use all the particles in the colour connected system to take the"
      " recoil but don't use the lepton.",
      3);
 
   static Parameter<HwRemDecayer,unsigned int> interfaceMaxTrySoft
     ("MaxTrySoft",
      "The maximum number of regeneration attempts for an additional soft scattering",
      &HwRemDecayer::maxtrySoft_, 10, 0, 100,
      false, false, Interface::limited);
 
   static Parameter<HwRemDecayer,double> interfacecolourDisrupt
     ("colourDisrupt",
      "Fraction of connections to additional soft subprocesses, which are colour disrupted.",
      &HwRemDecayer::colourDisrupt_, 
      1.0, 0.0, 1.0, 
      false, false, Interface::limited);
-  
+ 
+
+ 
   static Parameter<HwRemDecayer,double> interaceladderPower
     ("ladderPower",
      "The power factor in the ladder parameterization.",
      &HwRemDecayer::ladderPower_, 
      1.0, -5.0, 10.0, 
      false, false, Interface::limited);
 
   static Parameter<HwRemDecayer,double> interfaceladderNorm
     ("ladderNorm",
      "The normalization factor in the ladder parameterization",
      &HwRemDecayer::ladderNorm_,
      1.0, 0.0, 10.0,
      false, false, Interface::limited);
 
-  
+    static Parameter<HwRemDecayer,double> interfaceladderMult
+    ("ladderMult",
+     "The ladder multiplicity factor ",
+     &HwRemDecayer::ladderMult_,
+     1.0, 0.0, 10.0,
+     false, false, Interface::limited);
+
   static Parameter<HwRemDecayer,double> interfaceladderbFactor
     ("ladderbFactor",
      "The additive factor in the multiperipheral ladder multiplicity.",
      &HwRemDecayer::ladderbFactor_, 
      1.0, 0.0, 10.0, 
      false, false, Interface::limited);   
 
   static Parameter<HwRemDecayer,double> interfacegaussWidth
     ("gaussWidth",
      "The gaussian width of the fluctuation of longitudinal momentum fraction.",
      &HwRemDecayer::gaussWidth_, 
      0.1, 0.0, 1.0, 
      false, false, Interface::limited);   
 
 
   static Switch<HwRemDecayer,unsigned int> interfacePomeronStructure
     ("PomeronStructure",
      "Option for the treatment of the valance structure of the pomeron",
      &HwRemDecayer::pomeronStructure_, 0, false, false);
   static SwitchOption interfacePomeronStructureGluon
     (interfacePomeronStructure,
      "Gluon",
      "Assume the pomeron is a two gluon state",
      0);
   static SwitchOption interfacePomeronStructureQQBar
     (interfacePomeronStructure,
      "QQBar",
      "Assumne the pomeron is q qbar as for the photon,"
      " this option is not recommended and is provide for compatiblity with POMWIG",
      1);
 
   static Switch<HwRemDecayer,bool> interfaceAllowTop
     ("AllowTop",
      "Allow top quarks in the hadron",
      &HwRemDecayer::allowTop_, false, false, false);
   static SwitchOption interfaceAllowTopNo
     (interfaceAllowTop,
      "No",
      "Don't allow them",
      false);
   static SwitchOption interfaceAllowTopYes
     (interfaceAllowTop,
      "Yes",
      "Allow them",
      true);
    
    static Switch<HwRemDecayer,bool> interfaceMultiPeriph
     ("MultiPeriph",
      "Use multiperipheral kinematics",
      &HwRemDecayer::multiPeriph_, false, false, false);
   static SwitchOption interfaceMultiPeriphNo
     (interfaceMultiPeriph,
      "No",
      "Don't use multiperipheral",
      false);
   static SwitchOption interfaceMultiPeriphYes
     (interfaceMultiPeriph,
      "Yes",
      "Use multiperipheral kinematics",
      true);    
+  static Switch<HwRemDecayer,unsigned int> interfacePtDistribution
+    ("PtDistribution",
+     "Options for different pT generation methods",
+     &HwRemDecayer::PtDistribution_, 0, false, false);
+  static SwitchOption interfacePtDistributionDefault
+    (interfacePtDistribution,
+     "Default",
+     "Default generation of pT",
+     0);
+  static SwitchOption interfacePtDistributionOrdered
+    (interfacePtDistribution,
+     "Ordered",
+     "Ordered generation of pT,where the first pT is the hardest",
+     1);
+  static SwitchOption interfacePtDistributionStronglyOrdered
+    (interfacePtDistribution,
+     "StronglyOrdered",
+     "Strongly ordered generation of pT",
+     2);
+  static SwitchOption interfacePtDistributionFlat
+    (interfacePtDistribution,
+     "Flat",
+     "Sample from a flat pT distribution",
+     3);
+  static SwitchOption interfacePtDistributionFlatOrdered
+    (interfacePtDistribution,
+     "FlatOrdered",
+     "First pT normal, then flat",
+     4);
+  static SwitchOption interfacePtDistributionFlatOrdered2
+    (interfacePtDistribution,
+     "FlatOrdered2",
+     "First pT normal, then flat but steep",
+     5);
+
+static Switch<HwRemDecayer,unsigned int> interfaceRemnantConnection
+    ("RemnantConnection",
+     "Options for different colour connections between the remnant and the ladder",
+     &HwRemDecayer::disrupt_, 0, false, false);
+  static SwitchOption interfaceRemnantConnectionDisrupt
+    (interfaceRemnantConnection,
+     "Default",
+     "Connect new remnants to old remnants and the ladder is connected itself",
+     0);
+  static SwitchOption interfaceRemnantConnectionConnected
+    (interfaceRemnantConnection,
+     "Connected",
+     "Connect the remnants to the gluon ladder",
+     1);
+  static SwitchOption interfaceRemnantConnectionConnected2
+    (interfaceRemnantConnection,
+     "Connected2",
+     "Connect the remnants to the gluon ladder, possibility 2",
+     2);
+
+
+   static Switch<HwRemDecayer,bool> interfaceRandomConnection
+    ("RandomConnection",
+     "Choose between differnet remnant connections",
+     &HwRemDecayer::randomRemnantConnection_, false, false, false);
+  static SwitchOption interfaceRandomConnectionNo
+    (interfaceRandomConnection,
+     "No",
+     "Don't use random connections",
+     false);
+  static SwitchOption interfaceRandomConnectionYes
+    (interfaceRandomConnection,
+     "Yes",
+     "Use random connections ",
+     true);
+
 
 }
 
 bool HwRemDecayer::canHandle(tcPDPtr particle, tcPDPtr parton) const {
   if(! (StandardQCDPartonMatcher::Check(*parton) || parton->id()==ParticleID::gamma) ) {
     if(abs(parton->id())==ParticleID::t) {
       if(!allowTop_)
 	throw Exception() << "Top is not allow as a parton in the remant handling, please "
 			  << "use a PDF which does not contain top for the remnant"
 			  << " handling (preferred) or allow top in the remnant using\n"
 			  << " set " << fullName() << ":AllowTop Yes\n"
 			  << Exception::runerror;
     }
     else
       return false;
   }
   return HadronMatcher::Check(*particle) || particle->id()==ParticleID::gamma 
     || particle->id()==ParticleID::pomeron || particle->id()==ParticleID::reggeon;
 }
 
 bool HwRemDecayer::isPartonic(tPPtr parton) const {
   if(parton->parents().empty()) return false;
   tPPtr parent = parton->parents()[0];
   bool partonic = false;
   for(unsigned int ix=0;ix<parent->children().size();++ix) {
     if(dynamic_ptr_cast<tRemPPtr>(parent->children()[ix])) {
       partonic = true;
       break;
     }
   }
   return partonic;
 }
diff --git a/PDF/HwRemDecayer.h b/PDF/HwRemDecayer.h
--- a/PDF/HwRemDecayer.h
+++ b/PDF/HwRemDecayer.h
@@ -1,665 +1,757 @@
 // -*- C++ -*-
 //
 // HwRemDecayer.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_HwRemDecayer_H
 #define HERWIG_HwRemDecayer_H
 //
 // This is the declaration of the HwRemDecayer class.
 //
 
 #include "ThePEG/PDT/RemnantDecayer.h"
 #include "ThePEG/Handlers/EventHandler.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/EventRecord/SubProcess.h"
 #include "ThePEG/PDF/BeamParticleData.h"
 #include "Herwig/Shower/ShowerAlpha.h"
 #include "Herwig/PDT/StandardMatchers.h"
 #include "ThePEG/PDT/StandardMatchers.h"
 #include "HwRemDecayer.fh"
 
 namespace Herwig {
 using namespace ThePEG;
 /**
  * The HwRemDecayer class is responsible for the decay of the remnants. Additional 
  * secondary scatters have to be evolved backwards to a gluon, the
  * first/hard interaction has to be evolved back to a valence quark.
  * This is all generated inside this class,
  * which main methods are then called by the ShowerHandler.
  *
  * A simple forced splitting algorithm is used.
  * This takes the Remnant object produced from the PDF and backward
  * evolution (hadron - parton) and produce partons with the remaining 
  * flavours and with the correct colour connections.
  *
  * The algorithim operates by starting with the parton which enters the hard process.
  * If this is from the sea there is a forced branching to produce the antiparticle
  * from a gluon branching. If the parton entering the hard process was a gluon, or
  * a gluon was produced from the first step of the algorithm, there is then a further
  * branching back to a valence parton. After these partons have been produced a quark or
  * diquark is produced to give the remaining valence content of the incoming hadron.
  *
  * The forced branching are generated using a scale between QSpac and EmissionRange times
  * the minimum scale. The energy fractions are then distributed using
  * \f[\frac{\alpha_S}{2\pi}\frac{P(z)}{z}f(x/z,\tilde{q})\f]
  * with the massless splitting functions.
  *
  * \author Manuel B\"ahr
  *
  * @see \ref HwRemDecayerInterfaces "The interfaces"
  * defined for HwRemDecayer.
  */
 class HwRemDecayer: public RemnantDecayer {
 
 public:
 
   /** Typedef to store information about colour partners */
   typedef vector<pair<tPPtr, tPPtr> > PartnerMap;
 
 public:
 
   /**
    * The default constructor.
    */
   HwRemDecayer() : allowTop_(false), multiPeriph_(false), quarkPair_(false),
                    ptmin_(-1.*GeV), beta_(ZERO),
 		   maxtrySoft_(10), 
 		   colourDisrupt_(1.0),
 		   ladderbFactor_(0.0),
 		   ladderPower_(-0.08),
 		   ladderNorm_(1.0),
+ 		   ladderMult_(1.0),
 		   gaussWidth_(0.1),
 		   valOfN_(0), 
 		   initTotRap_(0),
 		   _kinCutoff(0.75*GeV), 
 		   _forcedSplitScale(2.5*GeV),
 		   _range(1.1), _zbin(0.05),_ybin(0.),
 		   _nbinmax(100), DISRemnantOpt_(0),
+                   PtDistribution_(0),
+                   disrupt_(0),
+		   randomRemnantConnection_(0),
 		   pomeronStructure_(0), mg_(ZERO) {}
 
   /** @name Virtual functions required by the Decayer class. */
   //@{
   /**
    * Check if this decayer can perfom the decay specified by the
    * given decay mode.
    * @return true if this decayer can handle the given mode, otherwise false.
    */
   virtual bool accept(const DecayMode &) const {
     return true;
   }
 
   /**
    * Return true if this decayer can handle the extraction of the \a   
    * extracted parton from the given \a particle.   
    */  
   virtual bool canHandle(tcPDPtr particle, tcPDPtr parton) const;
   
   /**   
    * Return true if this decayed can extract more than one parton from   
    * a particle.   
    */  
   virtual bool multiCapable() const {  
     return true;
   }
   
   /**
    * Perform a decay for a given DecayMode and a given Particle instance.
    * @param dm the DecayMode describing the decay.
    * @param p the Particle instance to be decayed.
    * @param step the step we are working on.
    * @return a ParticleVector containing the decay products.
    */
   virtual ParticleVector decay(const DecayMode & dm, const Particle & p, Step & step) const;
   //@}
 
 public:
 
   /** 
    * struct that is used to catch exceptions which are thrown
    * due to energy conservation issues of additional soft scatters
    */
   struct ExtraSoftScatterVeto {};
 
   /** @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();
 
   /**
    * Do several checks and initialization, for remnantdecay inside ShowerHandler.
    */
   void initialize(pair<tRemPPtr, tRemPPtr> rems, tPPair beam, Step & step,
 		  Energy forcedSplitScale);
 
   /**
    * Initialize the soft scattering machinery.
    * @param ptmin = the pt cutoff used in the UE model
    * @param beta = slope of the soft pt-spectrum
    */
   void initSoftInteractions(Energy ptmin, InvEnergy2 beta);
 
   /**
    * Perform the acual forced splitting.
    * @param partons is a pair of ThePEG::Particle pointers which store the final 
    * partons on which the shower ends.
    * @param pdfs are pointers to the pdf objects for both beams
    * @param first is a flage wether or not this is the first or a secondary interation
    */
   void doSplit(pair<tPPtr, tPPtr> partons, pair<tcPDFPtr, tcPDFPtr> pdfs, bool first);
 
   /**
    * Perform the final creation of the diquarks. Set the remnant masses and do 
    * all colour connections.
    * @param colourDisrupt = variable to control how many "hard" scatters
    * are colour isolated
    * @param softInt = parameter for the number of soft scatters
    */
   void finalize(double colourDisrupt=0.0, unsigned int softInt=0);
 
   /**
    *  Find the children
    */
   void findChildren(tPPtr,vector<PPtr> &) const;
 
 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);}
   //@}
 
 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() {
     Interfaced::doinit();
     _ybin=0.25/_zbin;
     mg_ = getParticleData(ParticleID::g)->constituentMass();
   }
   //@}
 
 private:
 
   /**
    * The assignment operator is private and must never be called.
    * In fact, it should not even be implemented.
    */
-  HwRemDecayer & operator=(const HwRemDecayer &) = delete;
+  HwRemDecayer & operator=(const HwRemDecayer &);
 
 public:
   
   /**                                                                           
    * Simple struct to store info about baryon quark and di-quark                
    * constituents.                                                              
    */                                                                           
   struct HadronContent {
 
     /**
      * manually extract the valence flavour \a id.
      */
     inline void extract(int id) {
       for(unsigned int i=0; i<flav.size(); i++) {
 	if(id == sign*flav[i]){
 	  if(hadron->id() == ParticleID::gamma || 
 	     (hadron->id() == ParticleID::pomeron && pomeronStructure==1) ||
 	     hadron->id() == ParticleID::reggeon) {
 	    flav[0] =  id;
 	    flav[1] = -id;
 	    extracted = 0;
 	    flav.resize(2);
 	  }
 	  else if (hadron->id() == ParticleID::pomeron && pomeronStructure==0) {
 	    extracted = 0;
 	  }
 	  else {
 	    extracted = i;
 	  }
 	  break;
 	}
       }      
     }
     
     /**
      * Return a proper particle ID assuming that \a id has been removed
      * from the hadron.
      */
     long RemID() const;
 
     /**
      * Method to determine whether \a parton is a quark from the sea.
      * @return TRUE if \a parton is neither a valence quark nor a gluon.
      */
     bool isSeaQuark(tcPPtr parton) const {
       return ((parton->id() != ParticleID::g) && ( !isValenceQuark(parton) ) );
     }
 
     /**
      * Method to determine whether \a parton is a valence quark.
      */
     bool isValenceQuark(tcPPtr parton) const {
       return isValenceQuark(parton->id());
     }
 
     /**
      * Method to determine whether \a parton is a quark from the sea.
      * @return TRUE if \a parton is neither a valence quark nor a gluon.
      */
     bool isSeaQuarkData(tcPDPtr partonData) const {
       return ((partonData->id() != ParticleID::g) && ( !isValenceQuarkData(partonData) ) );
     }
 
     /**
      * Method to determine whether \a parton is a valence quark.
      */
     bool isValenceQuarkData(tcPDPtr partonData) const {
       int id(sign*partonData->id());
       return find(flav.begin(),flav.end(),id) != flav.end();
     }
 
     /**
      * Method to determine whether \a parton is a valence quark.
      */
     bool isValenceQuark(int id) const {
       return find(flav.begin(),flav.end(),sign*id) != flav.end();
     }
 
     /** The valence flavours of the corresponding baryon. */                    
     vector<int> flav;                                                           
 
     /** The array index of the extracted particle. */
     int extracted;
 
     /** -1 if the particle is an anti-particle. +1 otherwise. */                
     int sign;
 
     /** The ParticleData objects of the hadron */
     tcPDPtr hadron;
 
     /** Pomeron treatment */
     unsigned int pomeronStructure;
   }; 
 
   /**
    * Return the hadron content objects for the incoming particles.
    */
   const pair<HadronContent, HadronContent>& content() const {
     return theContent;
   }
 
   /**
    * Return a HadronContent struct from a PPtr to a hadron.
    */
   HadronContent getHadronContent(tcPPtr hadron) const;
 
   /**
    * Set the hadron contents.
    */
   void setHadronContent(tPPair beam) {
     theContent.first  = getHadronContent(beam.first);
     theContent.second = getHadronContent(beam.second);
   }
 
 private:
 
   /**
    * Do the forced Splitting of the Remnant with respect to the 
    * extracted parton \a parton.
    * @param parton = PPtr to the parton going into the subprocess.
    * @param content = HadronContent struct to keep track of flavours.
    * @param rem = Pointer to the ThePEG::RemnantParticle.
    * @param used = Momentum vector to keep track of remaining momenta.
    * @param partners = Vector of pairs filled with tPPtr to the particles 
    * which should be colour connected.
    * @param pdf pointer to the PDF Object which is used for this particle
    * @param first = Flag for the first interaction.
    */
   void split(tPPtr parton, HadronContent & content, tRemPPtr rem, 
 	     Lorentz5Momentum & used, PartnerMap & partners, tcPDFPtr pdf, bool first);
 
   /**
    * Merge the colour lines of two particles
    * @param p1 = Pointer to particle 1
    * @param p2 = Pointer to particle 2
    * @param anti = flag to indicate, if (anti)colour was extracted as first parton.
    */
   void mergeColour(tPPtr p1, tPPtr p2, bool anti) const;
 
   /**
    * Set the colour connections.
    * @param partners = Object that holds the information which particles to connect.
    * @param anti = flag to indicate, if (anti)colour was extracted as first parton.
    * @param disrupt parameter for disruption of the colour structure
    */
   void fixColours(PartnerMap partners, bool anti, double disrupt) const;
 
   /**
    * Set the momenta of the Remnants properly and boost the decay particles.
    */
   void setRemMasses() const;
 
   /**
    * This creates a parton from the remaining flavours of the hadron. The
    * last parton used was a valance parton, so only 2 (or 1, if meson) flavours
    * remain to be used.
    */
   PPtr finalSplit(const tRemPPtr rem, long remID,
 		  Lorentz5Momentum usedMomentum) const {
     // Create the remnant and set its momentum, also reset all of the decay 
     // products from the hadron
     PPtr remnant = new_ptr(Particle(getParticleData(remID)));
     Lorentz5Momentum prem(rem->momentum()-usedMomentum);
     prem.setMass(getParticleData(remID)->constituentMass());
     prem.rescaleEnergy();
     remnant->set5Momentum(prem);
     // Add the remnant to the step, but don't do colour connections
     thestep->addDecayProduct(rem,remnant,false);
     return remnant;
   }
   
 
   /**
    * This takes the particle and find a splitting for np -> p + child and 
    * creates the correct kinematics and connects for such a split. This
    * Splitting has an upper bound on qtilde given by the energy argument
    * @param rem The Remnant
    * @param child The PDG code for the outgoing particle
    * @param oldQ  The maximum scale for the evolution
    * @param oldx  The fraction of the hadron's momentum carried by the last parton
    * @param pf    The momentum of the last parton at input and after branching at output
    * @param p     The total emitted momentum
    * @param content The content of the hadron
    */
   PPtr forceSplit(const tRemPPtr rem, long child, Energy &oldQ, double &oldx, 
 		  Lorentz5Momentum &pf, Lorentz5Momentum &p,
 		  HadronContent & content) const;
 
   /**
    *  Check if a particle is a parton from a hadron or not
    * @param parton The parton to be tested
    */
   bool isPartonic(tPPtr parton) const;
 
   /** @name Soft interaction methods. */
   //@{
 
   /**
    * Produce pt values according to dN/dp_T = N p_T exp(-beta_*p_T^2)
    */
   Energy softPt() const;
 
   /**
    * Get the 2 pairs of 5Momenta for the scattering. Needs calling of
    * initSoftInteractions.
    */
   void softKinematics(Lorentz5Momentum &r1, Lorentz5Momentum &r2, 
 		      Lorentz5Momentum &g1, Lorentz5Momentum &g2) const;
 
   /**
    * Create N soft gluon interactions
    */
   void doSoftInteractions(unsigned int N){
   	if(!multiPeriph_){
-  		doSoftInteractions_old(N);}
+  		doSoftInteractions_old(N);} //outdated model for soft interactions
   	else{
-  		doSoftInteractions_multiPeriph(N);
+                doSoftInteractions_multiPeriph(N); // Multiperipheral model
   	}
   }
   
   /**
    * Create N soft gluon interactions (old version)
    */
   void doSoftInteractions_old(unsigned int N);
   
   /**
-   * Create N soft gluon interactions - multiperhpheral kinematics
+   * Create N soft gluon interactions with multiperhpheral kinematics
    */
   void doSoftInteractions_multiPeriph(unsigned int N);
 
   /**
+   * Phase space generation for the ladder partons
+   */
+  bool doPhaseSpaceGenerationGluons(vector<Lorentz5Momentum> &softGluons, Energy energy, unsigned int &its)
+    const;
+
+  /**
+   * This returns the rotation matrix needed to rotate p into the z axis
+   */
+  LorentzRotation rotate(const LorentzMomentum &p) const;
+
+   /**
+    * Methods to generate random distributions also all stolen form UA5Handler
+    **/
+
+  template <typename T>
+  inline T gaussDistribution(T mean, T stdev) const{
+    double x = rnd();
+    x = sqrt(-2.*log(x));
+    double y;
+    randAzm(x,x,y);
+    return mean + stdev*x;
+  }
+
+
+  /**
+   * This returns a random number with a flat distribution
+   * [-A,A] plus gaussian tail with stdev B 
+   * TODO: Should move this to Utilities
+   * @param A The width of the flat part
+   * @param B The standard deviation of the gaussian tail
+   * @return the randomly generated value
+   */
+  inline double randUng(double A, double B) const{
+    double prun;
+    if(A == 0.) prun = 0.;
+    else prun = 1./(1.+B*1.2533/A);
+    if(rnd() < prun) return 2.*(rnd()-0.5)*A;
+    else {
+      double temp = gaussDistribution(0.,B);
+      if(temp < 0) return temp - abs(A);
+      else return temp + abs(A);
+    }
+  }
+  template <typename T>
+  inline void randAzm(T pt, T &px, T &py) const{
+    double c,s,cs;
+    while(true) {
+      c = 2.*rnd()-1.;
+      s = 2.*rnd()-1.;
+      cs = c*c+s*s;
+      if(cs <= 1.&&cs!=0.) break;
+    }
+    T qt = pt/cs;
+    px = (c*c-s*s)*qt;
+    py = 2.*c*s*qt;
+  }
+
+  inline Energy randExt(Energy AM0,InvEnergy B) const{
+    double r = rnd();
+    // Starting value
+    Energy am = AM0-log(r)/B;
+    for(int i = 1; i<20; ++i) {
+      double a = exp(-B*(am-AM0))/(1.+B*AM0);
+      double f = (1.+B*am)*a-r;
+      InvEnergy df = -B*B*am*a;
+      Energy dam = -f/df;
+      am += dam;
+      if(am<AM0) am = AM0 + .001*MeV;
+      if(abs(dam) < .001*MeV) break;
+    }
+    return am;
+  }
+
+  /**
    * Method to add a particle to the step
    * @param parent = pointer to the parent particle
    * @param id = Particle ID of the newly created particle
    * @param p = Lorentz5Momentum of the new particle
    */
   tPPtr addParticle(tcPPtr parent, long id, Lorentz5Momentum p) const;
   //@}
 
   /**
    * A flag which indicates, whether the extracted valence quark was a 
    * anti particle.
    */
   pair<bool, bool> theanti;
 
   /**
    * variable to sum up the x values of the extracted particles
    */
   pair<double, double> theX;
 
   /**Pair of HadronContent structs to know about the quark content of the beams*/
   pair<HadronContent, HadronContent> theContent;
 
   /**Pair of Lorentz5Momentum to keep track of the forced splitting product momenta*/
   pair<Lorentz5Momentum, Lorentz5Momentum> theUsed;
 
   /**
    * Pair of PartnerMap's to store the particles, which will be colour
    * connected in the end.
    */
   pair<PartnerMap, PartnerMap> theMaps;
 
   /**
    * Variable to hold a pointer to the current step. The variable is used to 
    * determine, wether decay(const DecayMode & dm, const Particle & p, Step & step) 
    * has been called in this event or not.
    */
   StepPtr thestep;
 
   /**
    * Pair of Remnant pointers. This is needed to boost
    * in the Remnant-Remnant CMF after all have been decayed.
    */
   pair<RemPPtr, RemPPtr> theRems;
 
   /**
    *  The beam particle data for the current incoming hadron
    */
   mutable tcPPtr theBeam;
 
   /**
    *  the beam data
    */
   mutable Ptr<BeamParticleData>::const_pointer theBeamData;
 
   /** 
    *  The PDF for the current initial-state shower 
    */ 
   mutable tcPDFPtr _pdf; 
   
 private:
 
   /**
    *  Switch to control handling of top quarks in proton
    */
   bool allowTop_;
   
   /**
    *  Switch to control using multiperipheral kinemaics
    */
   bool multiPeriph_;
   
   /**
    *  True if kinematics is to be calculated for quarks
    */
   bool quarkPair_;
 
   /** @name Soft interaction variables. */
   //@{
 
   /**
    * Pair of soft Remnant pointers, i.e. Diquarks.
    */
   tPPair softRems_;
 
   /**
    * ptcut of the UE model
    */
   Energy ptmin_;
 
   /**
    * slope of the soft pt-spectrum: dN/dp_T = N p_T exp(-beta*p_T^2)
    */
   InvEnergy2 beta_;
 
   /**
    *  Maximum number of attempts for the regeneration of an additional
    *  soft scattering, before the number of scatters is reduced.
    */
   unsigned int maxtrySoft_;
 
   /**
    * Variable to store the relative number of colour disrupted
    * connections to additional soft subprocesses.
    */
   double colourDisrupt_;
   
   /**
    * Variable to store the additive factor of the 
    multiperipheral ladder multiplicity.
    */
   double ladderbFactor_;
   
   /**
    * Variable of the parameterization of the ladder multiplicity.
    */
   double ladderPower_;
 
   /**
    * Variable of the parameterization of the ladder multiplicity.
    */
   double ladderNorm_;
 
+  double ladderMult_;
   /**
    * Variable to store the gaussian width of the 
    * fluctuation of the longitudinal momentum
    * fraction.
    */
   double gaussWidth_;
   
   /**
    * Variable to store the current total multiplicity 
    of a ladder.
    */
   double valOfN_;
   
   /**
    * Variable to store the initial total rapidity between 
    of the remnants.
    */
   double initTotRap_;
 
   //@}
 
   /** @name Forced splitting variables. */
   //@{
 
   /**
    *  The kinematic cut-off
    */
   Energy _kinCutoff;
   
   /**
    * The PDF freezing scale as set in ShowerHandler
    */
   Energy _forcedSplitScale;
 
   /**
    *  Range for emission
    */
   double _range;
 
   /**
    *  Size of the bins in z for the interpolation
    */
   double _zbin;
 
   /**
    *  Size of the bins in y for the interpolation
    */
   double _ybin;
 
   /**
    *  Maximum number of bins for the z interpolation
    */
   int _nbinmax;
 
   /**
    *  Pointer to the object calculating the QCD coupling
    */
   ShowerAlphaPtr _alphaS;
 
   /**
    *  Pointer to the object calculating the QED coupling
    */
   ShowerAlphaPtr _alphaEM; 
 
   /**
    *  Option for the DIS remnant
    */
   unsigned int DISRemnantOpt_;
 
   /**
+   *  Option for the pT generation
+   */
+  unsigned int PtDistribution_;
+
+
+  /**
+   * Option for remnant connection
+   */
+  unsigned int disrupt_;
+
+
+  bool randomRemnantConnection_;
+  /**
    *  Option for the treatment of the pomeron structure
    */
   unsigned int pomeronStructure_;
   //@}
 
   /**
    * The gluon constituent mass.
    */
   Energy mg_;
 
 };
 
 
 }
 
 #endif /* HERWIG_HwRemDecayer_H */
diff --git a/UnderlyingEvent/MPIHandler.cc b/UnderlyingEvent/MPIHandler.cc
--- a/UnderlyingEvent/MPIHandler.cc
+++ b/UnderlyingEvent/MPIHandler.cc
@@ -1,857 +1,857 @@
 // -*- C++ -*-
 //
 // MPIHandler.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 MPIHandler class.
 //
 
 #include "MPIHandler.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Handlers/StandardXComb.h"
 #include "ThePEG/Handlers/SubProcessHandler.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/Interface/RefVector.h"
 #include "ThePEG/Interface/ParVector.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/MatrixElement/MEBase.h"
 #include "ThePEG/Handlers/CascadeHandler.h"
 #include "ThePEG/Cuts/Cuts.h"
 #include "ThePEG/Cuts/JetCuts.h"
 #include "ThePEG/Cuts/SimpleKTCut.h"
 #include "ThePEG/PDF/PartonExtractor.h"
 
 #include "gsl/gsl_sf_bessel.h"
 
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 MPIHandler * MPIHandler::currentHandler_ = 0;
 
 bool MPIHandler::beamOK() const {
   return (HadronMatcher::Check(*eventHandler()->incoming().first)  &&
 	  HadronMatcher::Check(*eventHandler()->incoming().second) );
 }
 
 
 tStdXCombPtr MPIHandler::generate(unsigned int sel) {
   //generate a certain process
   if(sel+1 > processHandlers().size())
     throw Exception() << "MPIHandler::generate called with argument out of range"
                       << Exception::runerror;
 
   return processHandlers()[sel]->generate();
 }
 
 IBPtr MPIHandler::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr MPIHandler::fullclone() const {
   return new_ptr(*this);
 }
 
 void MPIHandler::finalize() {
   if( beamOK() ){
     statistics();
   }
 }
 
 void MPIHandler::initialize() {
   currentHandler_ = this;
   useMe();
   theHandler = generator()->currentEventHandler(); 
   //stop if the EventHandler is not present:
   assert(theHandler);
 
   //check if MPI is wanted
   if( !beamOK() ){
     throw Exception()  << "You have requested multiple parton-parton scattering,\n"
 		       << "but the model is not forseen for the beam setup you chose.\n" 
 		       << "You should therefore disable that by setting XXXGenerator:EventHandler:"
 		       << "CascadeHandler:MPIHandler to NULL"
                        << Exception::runerror;
   }
 
   numSubProcs_ = subProcesses().size();
 
   if( numSubProcs_ != cuts().size() ) 
     throw Exception() << "MPIHandler::each SubProcess needs a Cuts Object"
 		      << "ReferenceVectors are not equal in size"
 		      << Exception::runerror;
 
   if( additionalMultiplicities_.size()+1 != numSubProcs_ )
     throw Exception() << "MPIHandler: each additional SubProcess needs "
 		      << "a multiplicity assigned. This can be done in with "
 		      << "insert MPIHandler:additionalMultiplicities 0 1"
 		      << Exception::runerror;
 
   //identicalToUE_ = 0 hard process is identical to ue, -1 no one
   if( identicalToUE_ > (int)numSubProcs_ || identicalToUE_ < -1 )
     throw Exception() << "MPIHandler:identicalToUE has disallowed value"
 		      << Exception::runerror;
 
   // override the cuts for the additional scatters if energyExtrapolation_ is
   // set
   if (energyExtrapolation_ != 0 ) {
     overrideUECuts();
   }
 
   tcPDPtr gluon=getParticleData(ParticleID::g);
   //determine ptmin
   Ptmin_ = cuts()[0]->minKT(gluon);
 
   if(identicalToUE_ == -1){
     algorithm_ = 2;
   }else{
     if(identicalToUE_ == 0){
       //Need to work a bit, in case of LesHouches events for QCD2to2
       Ptr<StandardEventHandler>::pointer eH =
 	dynamic_ptr_cast<Ptr<StandardEventHandler>::pointer>(eventHandler());
 
       if( eH ) {
 	PtOfQCDProc_ = eH->cuts()->minKT(gluon);
 	// find the jet cut in the new style cuts
 	for(unsigned int ix=0;ix<eH->cuts()->multiCuts().size();++ix) {
 	  Ptr<JetCuts>::pointer jetCuts =
 	    dynamic_ptr_cast<Ptr<JetCuts>::pointer>(eH->cuts()->multiCuts()[ix]);
 	  if(jetCuts) {
 	    Energy ptMin=1e30*GeV;
 	    for(unsigned int iy=0;iy<jetCuts->jetRegions().size();++iy) {
 	      ptMin = min(ptMin,jetCuts->jetRegions()[iy]->ptMin());
 	    }
 	    if(ptMin<1e29*GeV&&ptMin>PtOfQCDProc_) PtOfQCDProc_ = ptMin;
 	  }
 	}
       }
       else {
 	if(PtOfQCDProc_ == -1.0*GeV)
 	  throw Exception() << "MPIHandler: You need to specify the pt cutoff "
 			    << "used to in the LesHouches file for QCD2to2 events"
 			    << Exception::runerror;
       }
     }
     else {
       PtOfQCDProc_ = cuts()[identicalToUE_]->minKT(gluon);
     }
 
     if(PtOfQCDProc_ > 2*Ptmin_)
       algorithm_ = 1;
     else
       algorithm_ = 0;
 
     if(PtOfQCDProc_ == ZERO)//pure MinBias mode
       algorithm_ = -1;
   }
 
   //Init all subprocesses
   for(unsigned int i=0; i<numSubProcs_; i++){
     theProcessHandlers.push_back(new_ptr(ProcessHandler()));
     processHandlers().back()->initialize(subProcesses()[i], 
 					 cuts()[i], eventHandler());
     processHandlers().back()->initrun();
   }
 
   //now calculate the individual Probabilities
   XSVector UEXSecs;
   UEXSecs.push_back(processHandlers()[0]->integratedXSec());
   //save the hard cross section
   hardXSec_ = UEXSecs.front();
 
   //determine sigma_soft and beta
   if(softInt_){//check that soft ints are requested
     GSLBisection rootFinder;
 
     if(twoComp_){
       //two component model
       /*
       GSLMultiRoot eqSolver;
       slopeAndTotalXSec eq(this);
       pair<CrossSection, Energy2> res = eqSolver.value(eq, 10*millibarn, 0.6*GeV2);
       softXSec_ = res.first;
       softMu2_ = res.second;
       */
       slopeBisection fs(this);
       try{
 	softMu2_ = rootFinder.value(fs, 0.3*GeV2, 1.*GeV2);
 	softXSec_ = fs.softXSec();
       }catch(GSLBisection::IntervalError){
         try{
 	// Very low energies (e.g. 200 GeV) need this.
 	// Very high energies (e.g. 100 TeV) produce issues with
 	// this choice. This is a temp. fix.  
           softMu2_ = rootFinder.value(fs, 0.25*GeV2, 1.3*GeV2);
           softXSec_ = fs.softXSec();
         }catch(GSLBisection::IntervalError){
 	throw Exception() <<
         "\n**********************************************************\n"	  "* Inconsistent MPI parameter choice for this beam energy *\n"
 	  "**********************************************************\n"
 	  "MPIHandler parameter choice is unable to reproduce\n"
 	  "the total cross section. Please check arXiv:0806.2949\n"
 	  "for the allowed parameter space."
 			  << Exception::runerror;
         }
       }
 
     }else{
       //single component model
       TotalXSecBisection fn(this);
       try{
 	softXSec_ = rootFinder.value(fn, 0*millibarn, 5000*millibarn);
       }catch(GSLBisection::IntervalError){
 	throw Exception() <<
 	"\n**********************************************************\n"
 	  "* Inconsistent MPI parameter choice for this beam energy *\n"
 	  "**********************************************************\n"
 	  "MPIHandler parameter choice is unable to reproduce\n"
 	  "the total cross section. Please check arXiv:0806.2949\n"
 	  "for the allowed parameter space."
 			  << Exception::runerror;      
       }
     }
 
     //now get the differential cross section at ptmin
     ProHdlPtr qcd = new_ptr(ProcessHandler());
     Energy eps = 0.1*GeV;
     Energy ptminPlus = Ptmin_ + eps;
     
     Ptr<SimpleKTCut>::pointer ktCut = new_ptr(SimpleKTCut(ptminPlus));
     ktCut->init();
     ktCut->initrun();
 
     CutsPtr qcdCut = new_ptr(Cuts(2*ptminPlus));
     qcdCut->add(dynamic_ptr_cast<tOneCutPtr>(ktCut));
     qcdCut->init();
     qcdCut->initrun();
       
     qcd->initialize(subProcesses()[0], qcdCut, eventHandler());
     qcd->initrun();
     
     // ds/dp_T^2 = 1/2/p_T ds/dp_T
     DiffXSec hardPlus = (hardXSec_-qcd->integratedXSec())/(2*Ptmin_*eps);
     
     betaBisection fn2(softXSec_, hardPlus, Ptmin_);
     try{
       beta_ = rootFinder.value(fn2, -10/GeV2, 2/GeV2);
     }catch(GSLBisection::IntervalError){
       try{
 	// Very low energies (e.g. 200 GeV) need this.
 	// Very high energies (e.g. 100 TeV) produce issues with
 	// this choice. This is a temp. fix.
         beta_ = rootFinder.value(fn2, -5/GeV2, 8./GeV2);
       }catch(GSLBisection::IntervalError){
       throw Exception() << "MPIHandler: slope of soft pt spectrum couldn't be "
 			<< "determined."
 			<< Exception::runerror;    
       }
     }
   }
 
   Probs(UEXSecs);
   //MultDistribution("probs.test");
 
   UEXSecs.clear();
 }
 
 void MPIHandler::MultDistribution(string filename) const {
   ofstream file;
   double p(0.0), pold(0.0);
   file.open(filename.c_str());
   //theMultiplicities  
   Selector<MPair>::const_iterator it = theMultiplicities.begin();
 
   while(it != theMultiplicities.end()){
     p = it->first;
     file << it->second.first << " " << it->second.second
 	 << " " << p-pold << '\n';
     it++;
     pold = p;
   }
 
   file << "sum of all probabilities: " << theMultiplicities.sum() 
        << endl;
 
   file.close();
 }
 
 void MPIHandler::statistics() const {
   ostream & file = generator()->misc();
   
   string line = "======================================="
     "=======================================\n";
 
   for(unsigned int i=0; i<cuts().size(); i++){
     Stat tot;
     if(i == 0)
       file << "Statistics for the UE process: \n";
     else
       file << "Statistics for additional hard Process " << i << ": \n";
 
     processHandlers()[i]->statistics(file, tot);
     file << "\n";
   }
 
   if(softInt_){
     file << line
 	 << "Eikonalized and soft cross sections:\n\n"
 	 << "Model parameters:                    "
 	 << "ptmin:   " << Ptmin_/GeV << " GeV"
       	 << ", mu2: " << invRadius_/sqr(1.*GeV) << " GeV2\n"
 	 << "                                     "
 	 << "DL mode: " << DLmode_
 	 << ", CMenergy: " << generator()->maximumCMEnergy()/GeV
 	 << " GeV" << '\n'
 	 << "hard inclusive cross section (mb):   "
 	 << hardXSec_/millibarn << '\n'
 	 << "soft inclusive cross section (mb):   "
 	 << softXSec_/millibarn << '\n'
 	 << "total cross section (mb):            "
 	 << totalXSecExp()/millibarn << '\n'
 	 << "inelastic cross section (mb):        "
 	 << inelXSec_/millibarn << '\n'
 	 << "soft inv radius (GeV2):              "
 	 << softMu2_/GeV2 << '\n'
 	 << "slope of soft pt spectrum (1/GeV2):  "
 	 << beta_*sqr(1.*GeV) << '\n'
 	 << "Average hard multiplicity:           "
 	 << avgNhard_ << '\n'
 	 << "Average soft multiplicity:           "
 	 << avgNsoft_ << '\n' << line << endl;
   }else{
     file << line
 	 << "Eikonalized and soft cross sections:\n\n"
 	 << "Model parameters:                    "
 	 << "ptmin:   " << Ptmin_/GeV << " GeV"
       	 << ", mu2: " << invRadius_/sqr(1.*GeV) << " GeV2\n"
 	 << "                                     "
 	 << ", CMenergy: " << generator()->maximumCMEnergy()/GeV
 	 << " GeV" << '\n'
 	 << "hard inclusive cross section (mb):   "
 	 << hardXSec_/millibarn << '\n'
 	 << "Average hard multiplicity:           "
 	 << avgNhard_ << '\n' << line << endl;
   }
 }
 
 unsigned int MPIHandler::multiplicity(unsigned int sel){
   if(sel==0){//draw from the pretabulated distribution
     MPair m = theMultiplicities.select(UseRandom::rnd());
     softMult_ = m.second;
     return m.first;
   } else{ //fixed multiplicities for the additional hard scatters
     if(additionalMultiplicities_.size() < sel)
       throw Exception() << "MPIHandler::multiplicity: process index "
 			<< "is out of range"
 			<< Exception::runerror;
 
     return additionalMultiplicities_[sel-1];
   }
 }
 
 
 void MPIHandler::Probs(XSVector UEXSecs) {
   GSLIntegrator integrator;
   unsigned int iH(1), iS(0);
   double P(0.0);
   double P0(0.0);//the probability for i hard and zero soft scatters
   Length bmax(500.0*sqrt(millibarn));
   //only one UE process will be drawn from a probability distribution,
   //so check that.
   assert(UEXSecs.size() == 1);
 
   for ( XSVector::const_iterator it = UEXSecs.begin();
         it != UEXSecs.end(); ++it ) {
 
     iH = 0; 
 
     //get the inel xsec
     Eikonalization inelint(this, -1, *it, softXSec_, softMu2_);
     inelXSec_ = integrator.value(inelint, ZERO, bmax);
 
     avgNhard_ = 0.0;
     avgNsoft_ = 0.0;
     bmax = 10.0*sqrt(millibarn);
     do{//loop over hard ints
       if(Algorithm()>-1 && iH==0){
 	iH++;
 	continue;
       }
       iS = 0;
       do{//loop over soft ints
 	if( ( Algorithm() == -1 && iS==0 && iH==0 ) ){
 	  iS++;
 	  continue;
 	}
 
 	Eikonalization integrand(this, iH*100+iS, *it, softXSec_, softMu2_);
       
 	if(Algorithm() > 0){
 	  P = integrator.value(integrand, ZERO, bmax) / *it;
 	}else{
 	  P = integrator.value(integrand, ZERO, bmax) / inelXSec_;
 	}
 	//store the probability
 	if(Algorithm()>-1){
 	  theMultiplicities.insert(P, make_pair(iH-1, iS));
 	  avgNhard_ += P*(iH-1);
 	}else{
 	  theMultiplicities.insert(P, make_pair(iH, iS));
 	  avgNhard_ += P*(iH);
 	}
 	avgNsoft_ += P*iS;
 	if(iS==0)
 	  P0 = P;
 
 	iS++;
       } while ( (iS < maxScatters_) && (iS < 5 || P > 1.e-15 ) && softInt_ );
       iH++;
     } while ( (iH < maxScatters_) && (iH < 5 || P0 > 1.e-15) );
   }
 }
 
 
 // calculate the integrand
 Length Eikonalization::operator() (Length b) const {
   unsigned int Nhard(0), Nsoft(0);
   //fac is just: d^2b=fac*db despite that large number
   Length fac(Constants::twopi*b);
   
   double chiTot(( theHandler->OverlapFunction(b)*hardXSec_ + 
 		  theHandler->OverlapFunction(b, softMu2_)*softXSec_ ) / 2.0);
 
   //total cross section wanted
   if(theoption == -2) return 2 * fac * ( 1 - exp(-chiTot) );
   //inelastic cross section
   if(theoption == -1) return fac * ( 1 - exp(- 2.0 * chiTot) );
 
   if(theoption >= 0){
     //encode multiplicities as: N_hard*100 + N_soft   
     Nhard = theoption/100;
     Nsoft = theoption%100;
 
     if(theHandler->Algorithm() > 0){
       //P_n*sigma_hard: n-1 extra scatters + 1 hard scatterer != hardXSec_
       return fac * Nhard * theHandler->poisson(b, hardXSec_, Nhard) *
 	theHandler->poisson(b, softXSec_, Nsoft, softMu2_);
     }else{
       //P_n*sigma_inel: n extra scatters
       return fac * theHandler->poisson(b, hardXSec_, Nhard) *
 	theHandler->poisson(b, softXSec_, Nsoft, softMu2_);
     }
 
   }else{
     throw Exception() << "Parameter theoption in Struct Eikonalization in " 
 		      << "MPIHandler.cc has not allowed value"
                       << Exception::runerror;
     return 0.0*meter;
   }
 }
 
 InvEnergy2 slopeBisection::operator() (Energy2 softMu2) const {
   GSLBisection root;
   //single component model
   TotalXSecBisection fn(handler_, softMu2);
   try{
     softXSec_ = root.value(fn, 0*millibarn, 5000*millibarn);
   }catch(GSLBisection::IntervalError){
     throw Exception() << "MPIHandler 2-Component model didn't work out."
 		      << Exception::runerror;
   }
   return handler_->slopeDiff(softXSec_, softMu2);
 }
 
 LengthDiff slopeInt::operator() (Length b) const {
   //fac is just: d^2b=fac*db
   Length fac(Constants::twopi*b);
   
   double chiTot(( handler_->OverlapFunction(b)*hardXSec_ + 
 		  handler_->OverlapFunction(b, softMu2_)*softXSec_ ) / 2.0);
 
   InvEnergy2 b2 = sqr(b/hbarc);
   //B*sigma_tot
   return fac * b2 * ( 1 - exp(-chiTot) );
 }
 
 double MPIHandler::factorial (unsigned int n) const {
 
   static double f[] = {1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880.,3.6288e6,
 		3.99168e7,4.790016e8,6.2270208e9,8.71782912e10,1.307674368e12,
 		2.0922789888e13,3.55687428096e14,6.402373705728e15,1.21645100408832e17,
 		2.43290200817664e18,5.10909421717094e19,1.12400072777761e21,
 		2.5852016738885e22,6.20448401733239e23,1.5511210043331e25,
 		4.03291461126606e26,1.08888694504184e28,3.04888344611714e29,
 		8.8417619937397e30,2.65252859812191e32,8.22283865417792e33,
 		2.63130836933694e35,8.68331761881189e36,2.95232799039604e38,
 		1.03331479663861e40,3.71993326789901e41,1.37637530912263e43,
 		5.23022617466601e44,2.03978820811974e46,8.15915283247898e47,
 		3.34525266131638e49,1.40500611775288e51,6.04152630633738e52,
 		2.65827157478845e54,1.1962222086548e56,5.50262215981209e57,
 		2.58623241511168e59,1.24139155925361e61,6.08281864034268e62,
 		3.04140932017134e64,1.55111875328738e66,8.06581751709439e67,
 		4.27488328406003e69,2.30843697339241e71,1.26964033536583e73,
 		7.10998587804863e74,4.05269195048772e76,2.35056133128288e78,
 		1.3868311854569e80,8.32098711274139e81,5.07580213877225e83,
 		3.14699732603879e85,1.98260831540444e87,1.26886932185884e89,
 		8.24765059208247e90,5.44344939077443e92,3.64711109181887e94,
 		2.48003554243683e96,1.71122452428141e98,1.19785716699699e100,
 		8.50478588567862e101,6.12344583768861e103,4.47011546151268e105,
 		3.30788544151939e107,2.48091408113954e109,1.88549470166605e111,
 		1.45183092028286e113,1.13242811782063e115,8.94618213078298e116,
 		7.15694570462638e118,5.79712602074737e120,4.75364333701284e122,
 		3.94552396972066e124,3.31424013456535e126,2.81710411438055e128,
 		2.42270953836727e130,2.10775729837953e132,1.85482642257398e134,
 		1.65079551609085e136,1.48571596448176e138,1.3520015276784e140,
 		1.24384140546413e142,1.15677250708164e144,1.08736615665674e146,
 		1.03299784882391e148,9.9167793487095e149,9.61927596824821e151,
 		9.42689044888325e153,9.33262154439442e155,9.33262154439442e157};
 
   if(n > maxScatters_) 
         throw Exception() << "MPIHandler::factorial called with too large argument"
                       << Exception::runerror;
   else
     return f[n];
 }
 
 InvArea MPIHandler::OverlapFunction(Length b, Energy2 mu2) const {
   if(mu2 == ZERO)
     mu2 = invRadius_;
 
   InvLength mu = sqrt(mu2)/hbarc;
   return (sqr(mu)/96/Constants::pi)*pow(mu*b, 3)*(gsl_sf_bessel_Kn(3, mu*b));
 }
 
 double MPIHandler::poisson(Length b, CrossSection sigma, unsigned int N, Energy2 mu2) const {
   if(sigma > 0*millibarn){
     return pow(OverlapFunction(b, mu2)*sigma, (double)N)/factorial(N)
       *exp(-OverlapFunction(b, mu2)*sigma);
   }else{
     return (N==0) ? 1.0 : 0.0;
   }
 }
 
 CrossSection MPIHandler::totalXSecDiff(CrossSection softXSec, 
 				       Energy2 softMu2) const {
   GSLIntegrator integrator;
   Eikonalization integrand(this, -2, hardXSec_, softXSec, softMu2);
   Length bmax = 500.0*sqrt(millibarn);
 
   CrossSection tot = integrator.value(integrand, ZERO, bmax);
   return (tot-totalXSecExp());
 }
 
 InvEnergy2 MPIHandler::slopeDiff(CrossSection softXSec, 
 				 Energy2 softMu2) const {
   GSLIntegrator integrator;
   Eikonalization integrand(this, -2, hardXSec_, softXSec, softMu2);
   Length bmax = 500.0*sqrt(millibarn);
 
   CrossSection tot = integrator.value(integrand, ZERO, bmax);
   
   slopeInt integrand2(this, hardXSec_, softXSec, softMu2);
   
   return integrator.value(integrand2, ZERO, bmax)/tot - slopeExp();
 }
 
 CrossSection MPIHandler::totalXSecExp() const {
   if(totalXSecExp_ != 0*millibarn)
     return totalXSecExp_;
 
   double pom_old = 0.0808;
   CrossSection coef_old = 21.7*millibarn;
   double pom_new_hard = 0.452;
   CrossSection coef_new_hard = 0.0139*millibarn;
   double pom_new_soft = 0.0667;
   CrossSection coef_new_soft = 24.22*millibarn;
 
   Energy energy(generator()->maximumCMEnergy());
 
   switch(DLmode_){
   case 1://old DL extrapolation
     return coef_old * pow(energy/GeV, 2*pom_old);
     break;
 
   case 2://old DL extrapolation fixed to CDF
     return 81.8*millibarn * pow(energy/1800.0/GeV, 2*pom_old);
     break;
     
   case 3://new DL extrapolation
-    return coef_new_hard * pow(energy/GeV, 2*pom_new_hard) + 
+    return 0.8*coef_new_hard * pow(energy/GeV, 2*pom_new_hard) + 
       coef_new_soft * pow(energy/GeV, 2*pom_new_soft);
     break;
     
   default:
     throw Exception() << "MPIHandler::totalXSecExp non-existing mode selected"
                       << Exception::runerror;   
   }
 }
 
 InvEnergy2 MPIHandler::slopeExp() const{
   //Currently return the slope as calculated in the pomeron fit by
   //Donnachie & Landshoff
   Energy energy(generator()->maximumCMEnergy());
   //slope at
   Energy e_0 = 1800*GeV;
   InvEnergy2 b_0 = 17/GeV2;
   return b_0 + log(energy/e_0)/GeV2;
 }
 
 void MPIHandler::overrideUECuts() {
   if(energyExtrapolation_==1)
     Ptmin_ = EEparamA_ * log(generator()->maximumCMEnergy() / EEparamB_);
   else if(energyExtrapolation_==2)
     Ptmin_ = pT0_*pow(double(generator()->maximumCMEnergy()/refScale_),b_);
   else
     assert(false);
   // create a new SimpleKTCut object with the calculated ptmin value
   Ptr<SimpleKTCut>::pointer newUEktCut = new_ptr(SimpleKTCut(Ptmin_));
   newUEktCut->init();
   newUEktCut->initrun();
 
   // create a new Cuts object with MHatMin = 2 * Ptmin_
   CutsPtr newUEcuts = new_ptr(Cuts(2*Ptmin_));
   newUEcuts->add(dynamic_ptr_cast<tOneCutPtr>(newUEktCut));
   newUEcuts->init();
   newUEcuts->initrun();
 
   // replace the old Cuts object
   cuts()[0] = newUEcuts;
 }
 
 void MPIHandler::persistentOutput(PersistentOStream & os) const {
   os << theMultiplicities << theHandler
      << theSubProcesses << theCuts << theProcessHandlers
      << additionalMultiplicities_ << identicalToUE_ 
      << ounit(PtOfQCDProc_, GeV) << ounit(Ptmin_, GeV) 
      << ounit(hardXSec_, millibarn) << ounit(softXSec_, millibarn)
      << ounit(beta_, 1/GeV2)
      << algorithm_ << ounit(invRadius_, GeV2)
      << numSubProcs_ << colourDisrupt_ << softInt_ << twoComp_ 
      << DLmode_ << ounit(totalXSecExp_, millibarn)
      << energyExtrapolation_ << ounit(EEparamA_, GeV) << ounit(EEparamB_, GeV)
      << ounit(refScale_,GeV) << ounit(pT0_,GeV) << b_
      << avgNhard_ << avgNsoft_ << softMult_ 
      << ounit(inelXSec_, millibarn) 
      << ounit(softMu2_, GeV2);
 }
 
 void MPIHandler::persistentInput(PersistentIStream & is, int) {
   is >> theMultiplicities >> theHandler
      >> theSubProcesses >> theCuts >> theProcessHandlers
      >> additionalMultiplicities_ >> identicalToUE_ 
      >> iunit(PtOfQCDProc_, GeV) >> iunit(Ptmin_, GeV)
      >> iunit(hardXSec_, millibarn) >> iunit(softXSec_, millibarn)
      >> iunit(beta_, 1/GeV2)
      >> algorithm_ >> iunit(invRadius_, GeV2)
      >> numSubProcs_ >> colourDisrupt_ >> softInt_ >> twoComp_ 
      >> DLmode_ >> iunit(totalXSecExp_, millibarn)
      >> energyExtrapolation_ >> iunit(EEparamA_, GeV) >> iunit(EEparamB_, GeV)
      >> iunit(refScale_,GeV) >> iunit(pT0_,GeV) >> b_
      >> avgNhard_ >> avgNsoft_ >> softMult_ 
      >> iunit(inelXSec_, millibarn) 
      >> iunit(softMu2_, GeV2);
   currentHandler_ = this;
 }
 
 void MPIHandler::clean() {
   // ThePEG's event handler's usual event cleanup doesn't reach these
   // XCombs. Need to do it by hand here.
   for ( size_t i = 0; i < theSubProcesses.size(); ++i ) {
     theSubProcesses[i]->pExtractor()->lastXCombPtr()->clean();
   }
 }
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<MPIHandler,Interfaced>
 describeHerwigMPIHandler("Herwig::MPIHandler",
 			 "JetCuts.so SimpleKTCut.so HwMPI.so");
 
 void MPIHandler::Init() {
 
   static ClassDocumentation<MPIHandler> documentation
     ("The MPIHandler class is the main administrator of the multiple interaction model", 
      "The underlying event was simulated with an eikonal model for multiple partonic interactions."
      "Details can be found in Ref.~\\cite{Bahr:2008dy,Bahr:2009ek}.", 
      "%\\cite{Bahr:2008dy}\n"
      "\\bibitem{Bahr:2008dy}\n"
      "  M.~Bahr, S.~Gieseke and M.~H.~Seymour,\n"
      "  ``Simulation of multiple partonic interactions in Herwig,''\n"
      "  JHEP {\\bf 0807}, 076 (2008)\n"
      "  [arXiv:0803.3633 [hep-ph]].\n"
      "  %%CITATION = JHEPA,0807,076;%%\n"
     "\\bibitem{Bahr:2009ek}\n"
      "  M.~Bahr, J.~M.~Butterworth, S.~Gieseke and M.~H.~Seymour,\n"
      "  ``Soft interactions in Herwig,''\n"
      "  arXiv:0905.4671 [hep-ph].\n"
      "  %%CITATION = ARXIV:0905.4671;%%\n"
      );
   
   static RefVector<MPIHandler,SubProcessHandler> interfaceSubhandlers
     ("SubProcessHandlers",
      "The list of sub-process handlers used in this EventHandler. ",
      &MPIHandler::theSubProcesses, -1, false, false, true, false, false);
 
   static RefVector<MPIHandler,Cuts> interfaceCuts
     ("Cuts",
      "List of cuts used for the corresponding list of subprocesses. These cuts "
      "should not be overidden in individual sub-process handlers.",
      &MPIHandler::theCuts, -1, false, false, true, false, false);
 
   static Parameter<MPIHandler,Energy2> interfaceInvRadius
     ("InvRadius",
      "The inverse hadron radius squared used in the overlap function",
      &MPIHandler::invRadius_, GeV2, 2.0*GeV2, 0.2*GeV2, 4.0*GeV2,
      true, false, Interface::limited);
 
   static ParVector<MPIHandler,int> interfaceadditionalMultiplicities
     ("additionalMultiplicities",
      "specify the multiplicities of secondary hard processes (multiple parton scattering)",
      &MPIHandler::additionalMultiplicities_, 
      -1, 0, 0, 3,
      false, false, true);
 
   static Parameter<MPIHandler,int> interfaceIdenticalToUE
     ("IdenticalToUE",
      "Specify which of the hard processes is identical to the UE one (QCD dijets)",
      &MPIHandler::identicalToUE_, -1, 0, 0,
      false, false, Interface::nolimits);
 
   static Parameter<MPIHandler,Energy> interfacePtOfQCDProc
     ("PtOfQCDProc",
      "Specify the value of the pt cutoff for the process that is identical to the UE one",
      &MPIHandler::PtOfQCDProc_, GeV, -1.0*GeV, ZERO, ZERO,
      false, false, Interface::nolimits);
 
   static Parameter<MPIHandler,double> interfacecolourDisrupt
     ("colourDisrupt",
      "Fraction of connections to additional subprocesses, which are colour disrupted.",
      &MPIHandler::colourDisrupt_, 
      0.0, 0.0, 1.0, 
      false, false, Interface::limited);
 
   
   static Switch<MPIHandler,bool> interfacesoftInt
     ("softInt",
      "Switch to enable soft interactions",
      &MPIHandler::softInt_, true, false, false);
 
   static SwitchOption interfacesoftIntYes
     (interfacesoftInt,
      "Yes",
      "enable the two component model",
      true);
   static SwitchOption interfacesoftIntNo
     (interfacesoftInt,
      "No",
      "disable the model",
      false);
 
 
   static Switch<MPIHandler,unsigned int> interEnergyExtrapolation
     ("EnergyExtrapolation",
      "Switch to ignore the cuts object at MPIHandler:Cuts[0]. "
      "Instead, extrapolate the pt cut.",
      &MPIHandler::energyExtrapolation_, 2, false, false);
   static SwitchOption interEnergyExtrapolationLog
     (interEnergyExtrapolation,
      "Log",
      "Use logarithmic dependence, ptmin = A * log (sqrt(s) / B).",
      1);
   static SwitchOption interEnergyExtrapolationPower
     (interEnergyExtrapolation,
      "Power",
      "Use power law, ptmin = pt_0 * (sqrt(s) / E_0)^b.",
      2);
   static SwitchOption interEnergyExtrapolationNo
     (interEnergyExtrapolation,
      "No",
      "Use manually set value for the minimal pt, "
      "specified in MPIHandler:Cuts[0]:OneCuts[0]:MinKT.",
      0);
 
   static Parameter<MPIHandler,Energy> interfaceEEparamA
     ("EEparamA",
      "Parameter A in the empirical parametrization "
      "ptmin = A * log (sqrt(s) / B)",
      &MPIHandler::EEparamA_, GeV, 0.6*GeV, ZERO, Constants::MaxEnergy,
      false, false, Interface::limited);
 
   static Parameter<MPIHandler,Energy> interfaceEEparamB
     ("EEparamB",
      "Parameter B in the empirical parametrization "
      "ptmin = A * log (sqrt(s) / B)",
      &MPIHandler::EEparamB_, GeV, 39.0*GeV, ZERO, Constants::MaxEnergy,
      false, false, Interface::limited);
 
   static Switch<MPIHandler,bool> interfacetwoComp
     ("twoComp",
      "switch to enable the model with a different radius for soft interactions",
      &MPIHandler::twoComp_, true, false, false);
 
   static SwitchOption interfacetwoCompYes
     (interfacetwoComp,
      "Yes",
      "enable the two component model",
      true);
   static SwitchOption interfacetwoCompNo
     (interfacetwoComp,
      "No",
      "disable the model",
      false);
 
 
   static Parameter<MPIHandler,CrossSection> interfaceMeasuredTotalXSec
     ("MeasuredTotalXSec",
      "Value for the total cross section (assuming rho=0). If non-zero, this "
      "overwrites the Donnachie-Landshoff parametrizations.",
      &MPIHandler::totalXSecExp_, millibarn, 0.0*millibarn, 0.0*millibarn, 0*millibarn,
      false, false, Interface::lowerlim);
 
   
   static Switch<MPIHandler,unsigned int> interfaceDLmode
     ("DLmode",
      "Choice of Donnachie-Landshoff parametrization for the total cross section.",
      &MPIHandler::DLmode_, 2, false, false);
   static SwitchOption interfaceDLmodeStandard
     (interfaceDLmode,
      "Standard",
      "Standard parametrization with s**0.08",
      1);
   static SwitchOption interfaceDLmodeCDF
     (interfaceDLmode,
      "CDF",
      "Standard parametrization but normalization fixed to CDF's measured value",
      2);
   static SwitchOption interfaceDLmodeNew
     (interfaceDLmode,
      "New",
      "Parametrization taking hard and soft pomeron contributions into account",
      3);
 
   static Parameter<MPIHandler,Energy> interfaceReferenceScale
     ("ReferenceScale",
      "The reference energy for power law energy extrapolation of pTmin",
      &MPIHandler::refScale_, GeV, 7000.0*GeV, 0.0*GeV, 20000.*GeV,
      false, false, Interface::limited);
 
   static Parameter<MPIHandler,Energy> interfacepTmin0
     ("pTmin0",
      "The pTmin at the reference scale for power law extrapolation of pTmin.",
      &MPIHandler::pT0_, GeV, 3.11*GeV, 0.0*GeV, 10.0*GeV,
      false, false, Interface::limited);
 
   static Parameter<MPIHandler,double> interfacePower
     ("Power",
      "The power for power law extrapolation of the pTmin cut-off.",
      &MPIHandler::b_, 0.21, 0.0, 10.0,
      false, false, Interface::limited);
 }
diff --git a/src/LHC-MB.in b/src/LHC-MB.in
--- a/src/LHC-MB.in
+++ b/src/LHC-MB.in
@@ -1,69 +1,103 @@
 # -*- ThePEG-repository -*-
 
 ################################################################################
 # This file contains our best tune to UE data from ATLAS at 7 TeV. More recent
 # tunes and tunes for other centre-of-mass energies as well as more usage
 # instructions can be obtained from this Herwig wiki page:
 # http://projects.hepforge.org/herwig/trac/wiki/MB_UE_tunes
 # The model for soft interactions and diffractions is explained in
 # [S. Gieseke, P. Kirchgaesser, F. Loshaj, arXiv:1612.04701]
 ################################################################################
 
 read snippets/PPCollider.in
 
 ##################################################
 # Technical parameters for this run
 ##################################################
 cd /Herwig/Generators
 ##################################################
 # LHC physics parameters (override defaults here) 
 ##################################################
 set EventGenerator:EventHandler:LuminosityFunction:Energy 7000.0
 
+# Intrinsic pT tune extrapolated to LHC energy
+set /Herwig/Shower/ShowerHandler:IntrinsicPtGaussian 2.2*GeV
 
 # Minimum Bias
 read snippets/MB.in
 
 # Read in parameters of the soft model recommended for MB/UE simulations
 read snippets/SoftTune.in
 
-# Diffraction model
+
+#Diffraction model
 read snippets/Diffraction.in
 
-# Read in snippet in order to use baryonic reconnection model with modified gluon splitting (uds)
-# For more details see [S. Gieseke, P. Kirchgaeßer, S. Plätzer. arXiv:1710.10906]]
-##############################################################################################
+#Turn on Baryonic reconnection
+read snippets/BaryonicReconnection.in
 
-# read snippets/BaryonicReconnection.in
 
+# Normalization of the Min bias cross section for correct diffractive cross section
+set /Herwig/MatrixElements/MEMinBias:csNorm 0.01
 
-##################################################
+
+# Use LHC parametrization of the cross section
+set /Herwig/UnderlyingEvent/MPIHandler:DLmode 3
+
+
+#some preliminary parameters for the MPI model which need to be tuned
+#TODO
+set /Herwig/UnderlyingEvent/MPIHandler:pTmin0 3.02
+set /Herwig/UnderlyingEvent/MPIHandler:InvRadius 0.9
+set /Herwig/UnderlyingEvent/MPIHandler:Power 0.308
+set /Herwig/Partons/RemnantDecayer:ladderMult 0.45
+set /Herwig/Partons/RemnantDecayer:ladderbFactor 1.0
+
+#ordered pTs
+set /Herwig/Partons/RemnantDecayer:PtDistribution 5
+#use random remnant connections
+set /Herwig/Partons/RemnantDecayer:RandomConnection Yes
+
+# set the correct PDFs
+cd /Herwig/Partons
+set HardLOPDF:PDFName CT14lo
+set ShowerLOPDF:PDFName CT14lo
+set MPIPDF:PDFName CT14lo
+set RemnantPDF:PDFName CT14lo
+
+#set /Herwig/Shower/ShowerHandler:CascadeHandler NULL     # Switch off parton shower
+#set /Herwig/Hadronization/HadronizationHandler:HadronizationHandler  NULL
+# do Soft interactions
+
+#set /Herwig/Hadronization/ColourReconnector:ColourReconnection No
+#################################################
 # Analyses
 ##################################################
 
+#Comment these lines out in order to use rivet analyses
+#cd /Herwig/Analysis
+#create ThePEG::RivetAnalysis RivetAnalysis RivetAnalysis.so
 
-# cd /Herwig/Analysis
-# create ThePEG::RivetAnalysis RivetAnalysis RivetAnalysis.so
+#cd /Herwig/Generators
+#insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
 
-# cd /Herwig/Generators
-# insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
-
-# insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_20XX_XXXXXXX
-
+#Some example analyses
+#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2012_I1084540
+#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2010_S8918562
 
 #set /Herwig/Analysis/Plot:EventNumber 54
 #cd /Herwig/Generators
 #insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/Plot
 
 #insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
 #set /Herwig/Analysis/HepMCFile:PrintEvent 1000000
 #set /Herwig/Analysis/HepMCFile:Format GenEvent
 #set /Herwig/Analysis/HepMCFile:Units GeV_mm
 #set /Herwig/Analysis/HepMCFile:Filename events.fifo
 
 ##################################################
 # Save run for later usage with 'Herwig run'
 ##################################################
 cd /Herwig/Generators
 saverun LHC-MB EventGenerator
 
diff --git a/src/snippets/BaryonicReconnection.in b/src/snippets/BaryonicReconnection.in
--- a/src/snippets/BaryonicReconnection.in
+++ b/src/snippets/BaryonicReconnection.in
@@ -1,17 +1,17 @@
 # Set strange quark mass to 0.45 in order to allow alternative gluon splitting
 set /Herwig/Particles/s:ConstituentMass 0.45*GeV
 set /Herwig/Particles/sbar:ConstituentMass 0.45*GeV
 
 # Use Baryonic Colour Reconnection Model
-set /Herwig/Hadronization/ColourReconnector:Algorithm BaryonicReco
+set /Herwig/Hadronization/ColourReconnector:Algorithm Baryonic
 
 # Allow alternative gluon splitting
 set /Herwig/Hadronization/PartonSplitter:Split uds
 
 # Parameters for the Baryonic Reconnection Model
 set /Herwig/Hadronization/ColourReconnector:ReconnectionProbability 0.772606     
 set /Herwig/Hadronization/ColourReconnector:ReconnectionProbabilityBaryonic 0.477612
 set /Herwig/UnderlyingEvent/MPIHandler:pTmin0 3.053252     
 set /Herwig/UnderlyingEvent/MPIHandler:InvRadius 1.282032
 set /Herwig/Hadronization/HadronSelector:PwtSquark 0.291717
 set /Herwig/Hadronization/PartonSplitter:SplitPwtSquark 0.824135
diff --git a/src/snippets/MB.in b/src/snippets/MB.in
--- a/src/snippets/MB.in
+++ b/src/snippets/MB.in
@@ -1,43 +1,55 @@
 ##################################################
 # MEMinBias Matrix Element
 ##################################################
 
 # MPI model settings
 set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0
 
 ## Report the correct cross section
 cd /Herwig/Generators
 create Herwig::MPIXSecReweighter MPIXSecReweighter
 insert EventGenerator:EventHandler:PostSubProcessHandlers 0 MPIXSecReweighter
 set EventGenerator:EventHandler:CascadeHandler NULL
 
 clear EventGenerator:EventHandler:SubProcessHandlers[0]
 
 ##################################################
 # Create separate SubProcessHandler for MinBias
 ##################################################
 cd /Herwig/MatrixElements/
 cp SubProcess QCDMinBias
 
 set QCDMinBias:CascadeHandler /Herwig/Shower/ShowerHandler
 set QCDMinBias:CascadeHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
 set QCDMinBias:DecayHandler /Herwig/Decays/DecayHandler
 
 # Due to numerics the pomeron could be seen as timelike.
 set /Herwig/Shower/ShowerHandler:SplitHardProcess No
 set /Herwig/DipoleShower/DipoleShowerHandler:SplitHardProcess No
 
 insert QCDMinBias:MatrixElements[0] MEMinBias
 
 cd /Herwig/Generators
 # MinBias parameters used for the new kinematics of soft MPI
-set /Herwig/Cuts/MinBiasCuts:X1Min 0.11
-set /Herwig/Cuts/MinBiasCuts:X2Min 0.11
+set /Herwig/Cuts/MinBiasCuts:X1Min 0.011
+set /Herwig/Cuts/MinBiasCuts:X2Min 0.011
+
+
+#PDFs for MPI,Underlying Event, Min Bias
+# set the correct PDFs
+cd /Herwig/Partons
+set HardLOPDF:PDFName CT14lo
+set ShowerLOPDF:PDFName CT14lo
+set MPIPDF:PDFName CT14lo
+set RemnantPDF:PDFName CT14lo
+
 
 # Needed to get the correct fraction of diffractive events
-set /Herwig/MatrixElements/MEMinBias:csNorm 4.5584
+set /Herwig/MatrixElements/MEMinBias:csNorm 0.01
+
+set /Herwig/MatrixElements/MEMinBias:Scale 2.0
 
 set EventGenerator:EventHandler:Cuts /Herwig/Cuts/MinBiasCuts
 
 cd /Herwig/MatrixElements/
 insert /Herwig/Generators/EventGenerator:EventHandler:SubProcessHandlers[0] QCDMinBias