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/PDF/HwRemDecayer.cc b/PDF/HwRemDecayer.cc
--- a/PDF/HwRemDecayer.cc
+++ b/PDF/HwRemDecayer.cc
@@ -1,1537 +1,2478 @@
 // -*- 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;
 
+//  cout << "the protons" <<endl;
+//  cout << P1/GeV <<endl;
+//  cout << P2/GeV <<endl;
+  
   ig1 = x_g1*P1;
   ig2 = x_g2*P2;
-  
+
+cout << "g1 = " << ig1/GeV <<endl;
+cout << "g2 = " <<ig2/GeV <<endl;
+
+
   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()){
+	//TODO:fix
     		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; 
 
+
+//cout << "Make " << N << "soft interactions" <<endl;
   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;
     }
-  
+ 
+//cout <<"g1: "<< g1/GeV<<endl;
+//cout <<"g2: "<< g2/GeV<<endl;
+ 
     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());
+
+
+//    int flow = UseRandom::rnd2(1,1);
+//    cout << "flow = " << flow <<endl;
+/*
+    flow=2;
+    switch (flow) {
+      case 0: 
+      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);
+      break;
+
+      case 1:
+
+      //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);
+      cl1->addColoured(gluons.first, !anti.first);
+      cl2->addColoured(softRems_.second, anti.second);
+      cl2->addColoured(gluons.second, !anti.second);
+*/
+    //  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_.first,anti.first);
+//      cl1->addColoured(gluons.first,!anti.first);
+//      cl2->addColoured(softRems_.second, anti.second);
+//      cl2->addColoured(gluons.second, !anti.second);
+
+      cl1->addColoured(softRems_.second, anti.first);
+      cl1->addColoured(gluons.first, !anti.first);
+      cl2->addColoured(softRems_.first, anti.second);
+      cl2->addColoured(gluons.second, !anti.second);
+
+  //  }
+/*
     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);
       cl1->addColoured(gluons.first, !anti.first);
       cl2->addColoured(softRems_.second, 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.;
+//    vector<Lorentz5Momentum>::const_iterator p = momenta.begin();
+//    vector<Energy>::const_iterator  m = masses.begin();
+    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;
+    }
+
+  }
+//  cout << "scaling factor = " << xi <<endl;
+  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_new(unsigned int N) {	
+  if(N == 0) return;
+  int Nmpi = N;
+
+  for(int j=0;j<Nmpi;j++){
+  ///////////////////////
+  // 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_;
+  // Parametrization of the ladder multiplicity
+  // 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
+  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){
+    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; 
+
+  // 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<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();
+ 
+  // 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);
+         cout << partons[0]->id() <<endl;
+         cl1->addColoured(partons[0]);     
+         // Connect old rems to the first gluon
+         cout << partons[1]->id() <<endl;
+         oldRems_.first->colourLine(anti.first)
+              ->addAntiColoured(partons[1]);
+
+         // Connect gluons with each other
+         for (unsigned int i=1; i<partons.size()-2; i++){
+            cout << i << endl;
+            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];
+    }
+
+  // generate the momenta
+  double eps = 1e-10/double(ncl);
+  vector<double> xi(ncl);
+  vector<Energy> tempEnergy(ncl);
+//  unsigned int its(0);
+  Energy sum1(ZERO);
+  double yy(0.);
+  while(its < _maxtries) {
+  //  cout << "try nr: " << its << endl;
+    ++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;
+  //  cout << "angle generation : " << iterations <<endl;
+    bool success=false;
+    Energy pTmax=ZERO;
+    Energy firstPt=ZERO;
+    for(unsigned int i = 0; i<ncl; ++i) {
+//      cout << i <<endl;
+      //To do: Different options for soft pt sampling
+      //1) pT1>pT2...pTN
+      //2) pT1>pT2>..>pTN
+      //3) flat
+      //4) y dependent
+      int triesPt=0;
+      Energy pt;
+      Energy ptTest;
+//      int pTgeneration=2;
+      switch(PtDistribution_) {
+        case 0: //default softPt()
+           pt=softPt();
+           break;
+        case 1: //pTordered
+           if(i==0){
+             pt=softPt();
+             pTmax=pt;
+  //           cout << "pTmax = " << pt/GeV <<endl;
+           }else{
+             do{
+              pt=softPt();
+    //          cout << "soft pT = " << pt/GeV << endl;
+             }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);
+            }
+            
+ 
+      }
+//      ofstream myfile3("softPt_in_Jadach.txt", ios::app );
+//      myfile3 << pt/GeV << endl;
+//      myfile3.close();
+ 
+      Energy2 ptp = pt*pt; //TODO: cannot be < ZERO (assert!)
+      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);
+//       cout << "px = " <<px/GeV << "py = "<<py/GeV <<endl;
+       sumxIt += px;
+       sumyIt += py;
+      }else{
+       //calculate azimuth angle s.t 
+      // double factor;
+      Energy pTdummy;
+      pTdummy = sqrt(sumxIt*sumxIt+sumyIt*sumyIt);
+//      cout << "pt = " << pTdummy/GeV <<endl;
+      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;
+      }
+      
+      }
+    }
+   if(success){
+//     cout << "calculation of angle was successful" <<endl;
+     break;
+   }
+   }
+    sumx /= ncl;
+    sumy /= ncl;
+  //  cout << sumy/GeV << " " <<sumx/GeV <<endl;
+    // find the sum of the transverse mass
+    Energy sumtm=ZERO;
+    for(unsigned int ix = 0; ix<ncl; ++ix) {
+    //  cout <<"sum x = " << sumx/GeV << " sum y = "<< sumy/GeV <<endl;
+      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
+      //TODO: vector mit Ncl komponenten statt temporaory 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) {
+//        cout <<"tries: "<< its <<endl;
+        yy += dyy;
+        suceeded=true;
+        break;
+      }
+      yy += dyy;
+    }
+    if(suceeded){
+ //      cout << "tries: " << its <<endl;
+       break;
+    }
+    if(its > 100) eps *= 10.;
+  }
+  if(its==_maxtries){
+//    cout << "maxtries reached" <<endl;
+    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];
+    //ofstream myfile3("softPt_Jadach_2.txt", ios::app );
+    //myfile3 << mom[i].vect().perp()/GeV << endl;
+    //myfile3.close();
+
+  }
+  return true;
+}
+
+
 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
-  
   // 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_;
 
   // Parametrization of the ladder multiplicity
   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
   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(!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;
  
   //momenta of the gluon pair generated in 
-  //each step
+  //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;
     }
 
     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){
+       } 
+       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));
-      	
-      }else{
+      }
+         if(i>1){
       	//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));  	
       }
-      
     }catch(ExtraSoftScatterVeto){
       tries++;
       i--;
       continue;
     }
     //reset counter
     tries = 1;
   }
   //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;
   }
   
   ////////
   }
   ////////
 }
 
 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,762 @@
 // -*- 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); //old multiperipheral model
+                doSoftInteractions_multiPeriph_new(N); //new multiperipheral model
   	}
   }
   
   /**
    * Create N soft gluon interactions (old version)
    */
   void doSoftInteractions_old(unsigned int N);
   
   /**
    * Create N soft gluon interactions - multiperhpheral kinematics
    */
   void doSoftInteractions_multiPeriph(unsigned int N);
 
+
+  void doSoftInteractions_multiPeriph_new(unsigned int N);
+
+  void generateMomentumGluons(PPtr remnant1, PPtr remnant2,
+                                    vector<Lorentz5Momentum> &softGluons, Energy &energy,
+                                    const Lorentz5Momentum &cm)const;
+
+  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/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