diff --git a/MatrixElement/Hadron/MEMinBias.cc b/MatrixElement/Hadron/MEMinBias.cc
--- a/MatrixElement/Hadron/MEMinBias.cc
+++ b/MatrixElement/Hadron/MEMinBias.cc
@@ -1,167 +1,174 @@
 // -*- C++ -*-
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the MEMinBias class.
 //
 
 #include "MEMinBias.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/SimplePhaseSpace.h"
 //#include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/Handlers/StandardXComb.h"
 #include "ThePEG/Interface/Parameter.h"
 
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 #include "ThePEG/PDT/EnumParticles.h"
 #include "ThePEG/MatrixElement/Tree2toNDiagram.h"
 
 void MEMinBias::getDiagrams() const {
   int maxflav(2);
   // Pomeron data
   tcPDPtr pom = getParticleData(990);
 
   for ( int i = 1; i <= maxflav; ++i ) {
     for( int j=1; j <= i; ++j){
       tcPDPtr q1 = getParticleData(i);
       tcPDPtr q1b = q1->CC();
       tcPDPtr q2 = getParticleData(j);
       tcPDPtr q2b = q2->CC();
 
       // For each flavour we add:
       //qq -> qq
       add(new_ptr((Tree2toNDiagram(3), q1, pom, q2, 1, q1, 2, q2, -1)));
       //qqb -> qqb
       add(new_ptr((Tree2toNDiagram(3), q1, pom, q2b, 1, q1, 2, q2b, -2)));
       //qbqb -> qbqb
       add(new_ptr((Tree2toNDiagram(3), q1b, pom, q2b, 1, q1b, 2, q2b, -3)));
     }
   }
 }
 
 Energy2 MEMinBias::scale() const {
-  return sqr(10*GeV);
+  return sqr(Scale_*GeV);
 }
 
 int MEMinBias::nDim() const {
   return 0;
 }
 
 void MEMinBias::setKinematics() {
   HwMEBase::setKinematics(); // Always call the base class method first.
 }
 
 bool MEMinBias::generateKinematics(const double *) {
   // generate the masses of the particles
   for ( int i = 2, N = meMomenta().size(); i < N; ++i ) {
     meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->generateMass());
   }
 
   Energy q = ZERO;
   try {
     q = SimplePhaseSpace::
       getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass());
   } catch ( ImpossibleKinematics & e ) {
     return false;
   }
 
   Energy pt = ZERO;
   meMomenta()[2].setVect(Momentum3( pt,  pt, q));
   meMomenta()[3].setVect(Momentum3(-pt, -pt, -q));
 
   meMomenta()[2].rescaleEnergy();
   meMomenta()[3].rescaleEnergy();
 
   jacobian(1.0);
   return true;
 }
 
 double MEMinBias::me2() const {
   //tuned so it gives the correct normalization for xmin = 0.11
   return csNorm_*(sqr(generator()->maximumCMEnergy())/GeV2);
 }
 
 CrossSection MEMinBias::dSigHatDR() const {
   return me2()*jacobian()/sHat()*sqr(hbarc);
 }
 
 unsigned int MEMinBias::orderInAlphaS() const {
   return 2;
 }
 
 unsigned int MEMinBias::orderInAlphaEW() const {
   return 0;
 }
 
 Selector<MEBase::DiagramIndex>
 MEMinBias::diagrams(const DiagramVector & diags) const {
   Selector<DiagramIndex> sel;
   for ( DiagramIndex i = 0; i < diags.size(); ++i ) 
     sel.insert(1.0, i);
 
   return sel;
 }
 
 Selector<const ColourLines *>
 MEMinBias::colourGeometries(tcDiagPtr diag) const {
 
   static ColourLines qq("1 4, 3 5");
   static ColourLines qqb("1 4, -3 -5");
   static ColourLines qbqb("-1 -4, -3 -5");
 
   Selector<const ColourLines *> sel;
   
   switch(diag->id()){
   case -1:
     sel.insert(1.0, &qq);
     break;
   case -2:
     sel.insert(1.0, &qqb);
     break;
   case -3:
     sel.insert(1.0, &qbqb);
     break;
   }
   return sel;
 }
 
 
 IBPtr MEMinBias::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr MEMinBias::fullclone() const {
   return new_ptr(*this);
 }
 
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<MEMinBias,HwMEBase>
 describeHerwigMEMinBias("Herwig::MEMinBias", "HwMEHadron.so");
 
 void MEMinBias::persistentOutput(PersistentOStream & os) const {
   os << csNorm_;
 }
 
 void MEMinBias::persistentInput(PersistentIStream & is, int) {
   is >> csNorm_;
 }
 
 void MEMinBias::Init() {
 
   static ClassDocumentation<MEMinBias> documentation
     ("There is no documentation for the MEMinBias class");
      
   static Parameter<MEMinBias,double> interfacecsNorm
     ("csNorm",
      "Normalization of the min-bias cross section.",
      &MEMinBias::csNorm_, 
      1.0, 0.0, 100.0, 
      false, false, Interface::limited);
+  static Parameter<MEMinBias,double> interfaceScale
+    ("Scale",
+     "Scale for the Min Bias matrix element.",
+     &MEMinBias::Scale_,
+     1.0, 0.0, 100.0,
+     false, false, Interface::limited);
+ 
 }
 
diff --git a/MatrixElement/Hadron/MEMinBias.h b/MatrixElement/Hadron/MEMinBias.h
--- a/MatrixElement/Hadron/MEMinBias.h
+++ b/MatrixElement/Hadron/MEMinBias.h
@@ -1,190 +1,195 @@
 // -*- C++ -*-
 #ifndef HERWIG_MEMinBias_H
 #define HERWIG_MEMinBias_H
 //
 // This is the declaration of the MEMinBias class.
 //
 
 #include "Herwig/MatrixElement/HwMEBase.h"
 
 namespace Herwig {
 
 using namespace ThePEG;
 
 /**
  * The MEMinBias class provides a simple colour singlet exchange matrix element
  * to be used in the soft component of the multiple scattering model of the 
  * underlying event
  *
  * @see \ref MEMinBiasInterfaces "The interfaces"
  * defined for MEMinBias.
  */
 class MEMinBias: public HwMEBase {
 
 public:
 
   /**
    * The default constructor.
    */
-  MEMinBias() : csNorm_(1.) {}
+  MEMinBias() : csNorm_(1.), Scale_(2.) {}
 
 public:
 
   /** @name Virtual functions required by the MEBase class. */
   //@{
   /**
    * Return the order in \f$\alpha_S\f$ in which this matrix
    * element is given.
    */
   virtual unsigned int orderInAlphaS() const;
 
   /**
    * Return the order in \f$\alpha_{EW}\f$ in which this matrix
    * element is given.
    */
   virtual unsigned int orderInAlphaEW() const;
 
   /**
    * The matrix element for the kinematical configuration
    * previously provided by the last call to setKinematics(), suitably
    * scaled by sHat() to give a dimension-less number.
    * @return the matrix element scaled with sHat() to give a
    * dimensionless number.
    */
   virtual double me2() const;
 
   /**
    * Return the scale associated with the last set phase space point.
    */
   virtual Energy2 scale() const;
 
   /**
    * Set the typed and momenta of the incoming and outgoing partons to
    * be used in subsequent calls to me() and colourGeometries()
    * according to the associated XComb object. If the function is
    * overridden in a sub class the new function must call the base
    * class one first.
    */
   virtual void setKinematics();
 
   /**
    * The number of internal degrees of freedom used in the matrix
    * element.
    */
   virtual int nDim() const;
 
   /**
    * Generate internal degrees of freedom given nDim() uniform
    * random numbers in the interval \f$ ]0,1[ \f$. To help the phase space
    * generator, the dSigHatDR should be a smooth function of these
    * numbers, although this is not strictly necessary.
    * @param r a pointer to the first of nDim() consecutive random numbers.
    * @return true if the generation succeeded, otherwise false.
    */
   virtual bool generateKinematics(const double * r);
 
   /**
    * Return the matrix element squared differential in the variables
    * given by the last call to generateKinematics().
    */
   virtual CrossSection dSigHatDR() const;
 
   /**
    * Add all possible diagrams with the add() function.
    */
   virtual void getDiagrams() const;
 
   /**
    * Get diagram selector. With the information previously supplied with the
    * setKinematics method, a derived class may optionally
    * override this method to weight the given diagrams with their
    * (although certainly not physical) relative probabilities.
    * @param dv the diagrams to be weighted.
    * @return a Selector relating the given diagrams to their weights.
    */
   virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
 
   /**
    * Return a Selector with possible colour geometries for the selected
    * diagram weighted by their relative probabilities.
    * @param diag the diagram chosen.
    * @return the possible colour geometries weighted by their
    * relative probabilities.
    */
   virtual Selector<const ColourLines *>
   colourGeometries(tcDiagPtr diag) const;
   //@}
 
 
 public:
   /** @name Functions used by the persistent I/O system. */
   //@{
   /**
    * Function used to write out object persistently.
    * @param os the persistent output stream written to.
    */
   void persistentOutput(PersistentOStream & os) const;
 
   /**
    * Function used to read in object persistently.
    * @param is the persistent input stream read from.
    * @param version the version number of the object when written.
    */
   void persistentInput(PersistentIStream & is, int version);
   //@}
 
   /**
    * The standard Init function used to initialize the interfaces.
    * Called exactly once for each class by the class description system
    * before the main function starts or
    * when this class is dynamically loaded.
    */
 
 
   /**
    * The standard Init function used to initialize the interfaces.
    * Called exactly once for each class by the class description system
    * before the main function starts or
    * when this class is dynamically loaded.
    */
   static void Init();
 
 protected:
 
   /** @name Clone Methods. */
   //@{
   /**
    * Make a simple clone of this object.
    * @return a pointer to the new object.
    */
   virtual IBPtr clone() const;
 
   /** Make a clone of this object, possibly modifying the cloned object
    * to make it sane.
    * @return a pointer to the new object.
    */
   virtual IBPtr fullclone() const;
   //@}
 
 
 // If needed, insert declarations of virtual function defined in the
 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
 
 
 private:
   /**
   *  Normalization of the min-bias cross section
   */
   double csNorm_;	
 
   /**
+   * Scale for the Min Bias matrix element
+   */
+  double Scale_;
+
+  /**
    * The assignment operator is private and must never be called.
    * In fact, it should not even be implemented.
    */
   MEMinBias & operator=(const MEMinBias &) = delete;
 
 };
 
 }
 
 #endif /* HERWIG_MEMinBias_H */
diff --git a/PDF/HwRemDecayer.cc b/PDF/HwRemDecayer.cc
--- a/PDF/HwRemDecayer.cc
+++ b/PDF/HwRemDecayer.cc
@@ -1,2276 +1,2183 @@
 // -*- C++ -*-
 //
 // HwRemDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2017 The Herwig Collaboration
 //
 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
 //
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the HwRemDecayer class.
 //
 
 #include "HwRemDecayer.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/Interface/Reference.h"  
 #include "ThePEG/Interface/Parameter.h"  
 #include "ThePEG/Interface/Switch.h"  
 #include "ThePEG/Utilities/UtilityBase.h"
 #include "ThePEG/Utilities/SimplePhaseSpace.h"
 #include "ThePEG/Utilities/Throw.h"
 
 #include "Herwig/Shower/ShowerHandler.h"
 
 using namespace Herwig;
 
 namespace{
 
 const bool dbg = false;
 
 
   void reShuffle(Lorentz5Momentum &p1, Lorentz5Momentum &p2, Energy m1, Energy m2){
 
     Lorentz5Momentum ptotal(p1+p2);
     ptotal.rescaleMass();
 
     if( ptotal.m() < m1+m2 ) {
       if(dbg)
 	cerr << "Not enough energy to perform reshuffling \n";
       throw HwRemDecayer::ExtraSoftScatterVeto();
     }
 
     Boost boostv = -ptotal.boostVector();
     ptotal.boost(boostv);
 
     p1.boost(boostv);
     // set the masses and energies,
     p1.setMass(m1);
     p1.setE(0.5/ptotal.m()*(ptotal.m2()+sqr(m1)-sqr(m2)));
     p1.rescaleRho();
     // boost back to the lab
     p1.boost(-boostv);
 
     p2.boost(boostv);
     // set the masses and energies,
     p2.setMass(m2);
     p2.setE(0.5/ptotal.m()*(ptotal.m2()+sqr(m2)-sqr(m1)));
     p2.rescaleRho();
     // boost back to the lab
     p2.boost(-boostv);
   }
 }
 
 void HwRemDecayer::initialize(pair<tRemPPtr, tRemPPtr> rems, tPPair beam, Step & step,
 			      Energy forcedSplitScale) {
   // the step
   thestep = &step;
   // valence content of the hadrons
   theContent.first  = getHadronContent(beam.first);
   theContent.second = getHadronContent(beam.second);
   // momentum extracted from the hadrons
   theUsed.first  = Lorentz5Momentum();
   theUsed.second = Lorentz5Momentum();
   theMaps.first.clear();
   theMaps.second.clear();
   theX.first = 0.0;
   theX.second = 0.0;
   theRems = rems;
   _forcedSplitScale = forcedSplitScale;
   // check remnants attached to the right hadrons
   if( (theRems.first  && parent(theRems.first ) != beam.first ) ||
       (theRems.second && parent(theRems.second) != beam.second) )
     throw Exception() << "Remnant order wrong in "
 		      << "HwRemDecayer::initialize(...)"
 		      << Exception::runerror; 
   return;
 }
 
 void HwRemDecayer::split(tPPtr parton, HadronContent & content, 
 			 tRemPPtr rem, Lorentz5Momentum & used, 
 			 PartnerMap &partners, tcPDFPtr pdf, bool first) {
   theBeam = parent(rem);
   theBeamData = dynamic_ptr_cast<Ptr<BeamParticleData>::const_pointer>
     (theBeam->dataPtr());
   double currentx = parton->momentum().rho()/theBeam->momentum().rho();
   double check = rem==theRems.first ? theX.first : theX.second;
   check += currentx;
   if(1.0-check < 1e-3) throw ShowerHandler::ExtraScatterVeto();
   bool anti;
   Lorentz5Momentum lastp(parton->momentum());
   int lastID(parton->id());
   Energy oldQ(_forcedSplitScale);
   _pdf = pdf;
   //do nothing if already valence quark
   if(first && content.isValenceQuark(parton)) { 
     //set the extracted value, because otherwise no RemID could be generated.
     content.extract(lastID);
     // add the particle to the colour partners
     partners.push_back(make_pair(parton, tPPtr()));
     //set the sign
     anti = parton->hasAntiColour() && parton->id()!=ParticleID::g;
     if(rem==theRems.first) theanti.first  = anti;
     else                   theanti.second = anti;
     // add the x and return
     if(rem==theRems.first) theX.first  += currentx;
     else                   theX.second += currentx;
     return;
   }
   //or gluon for secondaries
   else if(!first && lastID == ParticleID::g) {
     partners.push_back(make_pair(parton, tPPtr()));
     // add the x and return
     if(rem==theRems.first) theX.first  += currentx;
     else                   theX.second += currentx;
     return; 
   }
   // if a sea quark.antiquark forced splitting to a gluon
   // Create the new parton with its momentum and parent/child relationship set
   PPtr newSea;
   if( !(lastID == ParticleID::g || 
 	lastID == ParticleID::gamma) ) {
     newSea = forceSplit(rem, -lastID, oldQ, currentx, lastp, used,content);
     ColinePtr cl = new_ptr(ColourLine());
     if(newSea->id() > 0) cl->    addColoured(newSea);
     else                 cl->addAntiColoured(newSea);
     // if a secondard scatter finished so return
     if(!first || content.isValenceQuark(ParticleID::g) ){
       partners.push_back(make_pair(parton, newSea));
       // add the x and return
       if(rem==theRems.first) theX.first  += currentx;
       else                   theX.second += currentx;
       if(first) content.extract(ParticleID::g);
       return;
     }
   }
   // otherwise evolve back to valence
   // final valence splitting
   PPtr newValence = forceSplit(rem, 
 			       lastID!=ParticleID::gamma ? 
 			       ParticleID::g : ParticleID::gamma,
 			       oldQ, currentx , lastp, used, content);
   // extract from the hadron to allow remnant to be determined
   content.extract(newValence->id());
   // case of a gluon going into the hard subprocess
   if( lastID == ParticleID::g ) {
     partners.push_back(make_pair(parton, tPPtr()));
     anti = newValence->hasAntiColour();
     if(rem==theRems.first) theanti.first  = anti;
     else                   theanti.second = anti;
     parton->colourLine(!anti)->addColoured(newValence, anti);
     return;
   }
   else if( lastID == ParticleID::gamma) {
     partners.push_back(make_pair(parton, newValence));
     anti = newValence->hasAntiColour();
     ColinePtr newLine(new_ptr(ColourLine()));
     newLine->addColoured(newValence, anti);
     if(rem==theRems.first) theanti.first  = anti;
     else                   theanti.second = anti;
     // add the x and return
     if(rem==theRems.first) theX.first  += currentx;
     else                   theX.second += currentx;
     return;
   }
   //The valence quark will always be connected to the sea quark with opposite sign
   tcPPtr particle;
   if(lastID*newValence->id() < 0){
     particle = parton;
     partners.push_back(make_pair(newSea, tPPtr()));
   } 
   else {
     particle = newSea;
     partners.push_back(make_pair(parton, tPPtr()));
   }
   anti = newValence->hasAntiColour();
   if(rem==theRems.first) theanti.first  = anti;
   else                   theanti.second = anti;
   
   if(particle->colourLine()) 
     particle->colourLine()->addAntiColoured(newValence);
   if(particle->antiColourLine()) 
     particle->antiColourLine()->addColoured(newValence);
   // add the x and return
   if(rem==theRems.first) theX.first  += currentx;
   else                   theX.second += currentx;
   return;
 }
 
 void HwRemDecayer::doSplit(pair<tPPtr, tPPtr> partons,
 			   pair<tcPDFPtr, tcPDFPtr> pdfs,
 			   bool first) {
   if(theRems.first) {
     ParticleVector children=theRems.first->children();
     for(unsigned int ix=0;ix<children.size();++ix) {
       if(children[ix]->dataPtr()==theRems.first->dataPtr()) 
 	theRems.first = dynamic_ptr_cast<RemPPtr>(children[ix]);
     }
   }
   if(theRems.second) {
     ParticleVector children=theRems.second->children();
     for(unsigned int ix=0;ix<children.size();++ix) {
       if(children[ix]->dataPtr()==theRems.second->dataPtr()) 
 	theRems.second = dynamic_ptr_cast<RemPPtr>(children[ix]);
     }
   }
   // forced splitting for first parton
   if(isPartonic(partons.first )) { 
     try {
       split(partons.first, theContent.first, theRems.first, 
 	    theUsed.first, theMaps.first, pdfs.first, first);
     }
     catch(ShowerHandler::ExtraScatterVeto) {
       throw ShowerHandler::ExtraScatterVeto();
     }
   }
   // forced splitting for second parton
   if(isPartonic(partons.second)) { 
     try {
       split(partons.second, theContent.second, theRems.second, 
 	    theUsed.second, theMaps.second, pdfs.second, first);
       // additional check for the remnants
       // if can't do the rescale veto the emission
       if(!first&&partons.first->data().coloured()&&
 	 partons.second->data().coloured()) {
 	Lorentz5Momentum pnew[2]=
 	  {theRems.first->momentum()  - theUsed.first  - partons.first->momentum(),
 	   theRems.second->momentum() - theUsed.second - partons.second->momentum()};
 
 	pnew[0].setMass(getParticleData(theContent.first.RemID())->constituentMass());
 	pnew[0].rescaleEnergy();
 	pnew[1].setMass(getParticleData(theContent.second.RemID())->constituentMass());
 	pnew[1].rescaleEnergy();
 	
 	for(unsigned int iy=0; iy<theRems.first->children().size(); ++iy)
 	  pnew[0] += theRems.first->children()[iy]->momentum();
 	
 	for(unsigned int iy=0; iy<theRems.second->children().size(); ++iy)
 	  pnew[1] += theRems.second->children()[iy]->momentum();
 	
 	Lorentz5Momentum ptotal=
 	  theRems.first ->momentum()-partons.first ->momentum()+
 	  theRems.second->momentum()-partons.second->momentum();
 	// add x limits
 	if(ptotal.m() < (pnew[0].m() + pnew[1].m()) ) {
 	  if(partons.second->id() != ParticleID::g){
 	    if(partons.second==theMaps.second.back().first) 
 	      theUsed.second -= theMaps.second.back().second->momentum();
 	    else
 	      theUsed.second -= theMaps.second.back().first->momentum();
 	    
 	    thestep->removeParticle(theMaps.second.back().first);
 	    thestep->removeParticle(theMaps.second.back().second);
 	  }
 	  theMaps.second.pop_back();
 	  theX.second -= partons.second->momentum().rho()/
 	    parent(theRems.second)->momentum().rho();
 	  throw ShowerHandler::ExtraScatterVeto();
 	}
       }
     }
     catch(ShowerHandler::ExtraScatterVeto){
       if(!partons.first||!partons.second||
 	 !theRems.first||!theRems.second) 
 	throw ShowerHandler::ExtraScatterVeto();
       //case of the first forcedSplitting worked fine
       theX.first -= partons.first->momentum().rho()/
 	parent(theRems.first)->momentum().rho();
       //case of the first interaction
       //throw veto immediately, because event get rejected anyway.
       if(first) throw ShowerHandler::ExtraScatterVeto();
       //secondary interactions have to end on a gluon, if parton 
       //was NOT a gluon, the forced splitting particles must be removed
       if(partons.first->id() != ParticleID::g) {
 	if(partons.first==theMaps.first.back().first) 
 	  theUsed.first -= theMaps.first.back().second->momentum();
 	else
 	  theUsed.first -= theMaps.first.back().first->momentum();
 	
 	thestep->removeParticle(theMaps.first.back().first);
 	thestep->removeParticle(theMaps.first.back().second);
       }
       theMaps.first.pop_back();
       throw ShowerHandler::ExtraScatterVeto();
     }
   }
   // veto if not enough energy for extraction
   if( !first &&(theRems.first ->momentum().e() - 
 		partons.first ->momentum().e() < 1.0e-3*MeV ||
 		theRems.second->momentum().e() - 
 		partons.second->momentum().e() < 1.0e-3*MeV )) {
     if(partons.first->id() != ParticleID::g) {
       if(partons.first==theMaps.first.back().first) 
 	theUsed.first -= theMaps.first.back().second->momentum();
       else
 	theUsed.first -= theMaps.first.back().first->momentum();
       thestep->removeParticle(theMaps.first.back().first);
       thestep->removeParticle(theMaps.first.back().second);
     }
     theMaps.first.pop_back();
     if(partons.second->id() != ParticleID::g) {
       if(partons.second==theMaps.second.back().first) 
 	theUsed.second -= theMaps.second.back().second->momentum();
       else
 	theUsed.second -= theMaps.second.back().first->momentum();
       thestep->removeParticle(theMaps.second.back().first);
       thestep->removeParticle(theMaps.second.back().second);
     }
     theMaps.second.pop_back();
     throw ShowerHandler::ExtraScatterVeto();
   }
 }
 
 void HwRemDecayer::mergeColour(tPPtr pold, tPPtr pnew, bool anti) const {
   ColinePtr clnew, clold;
   
   //save the corresponding colour lines
   clold = pold->colourLine(anti);
   clnew = pnew->colourLine(!anti);
 
   assert(clold);
 
   // There is already a colour line (not the final diquark)
   if(clnew){
 
     if( (clnew->coloured().size() + clnew->antiColoured().size()) > 1 ){
       if( (clold->coloured().size() + clold->antiColoured().size()) > 1 ){
 	//join the colour lines
 	//I don't use the join method, because potentially only (anti)coloured
 	//particles belong to one colour line
 	if(clold!=clnew){//procs are not already connected
 	  while ( !clnew->coloured().empty() ) {
 	    tPPtr p = clnew->coloured()[0];
 	    clnew->removeColoured(p);
 	    clold->addColoured(p);
 	  }
 	  while ( !clnew->antiColoured().empty() ) {
 	    tPPtr p = clnew->antiColoured()[0];
 	    clnew->removeAntiColoured(p);
 	    clold->addAntiColoured(p);
 	  }
 	}
 
       }else{
 	//if pold is the only member on it's 
 	//colour line, remove it.
 	clold->removeColoured(pold, anti);
 	//and add it to clnew
 	clnew->addColoured(pold, anti);
       }    
     } else{//pnnew is the only member on it's colour line.
       clnew->removeColoured(pnew, !anti);
       clold->addColoured(pnew, !anti);
     }
   } else {//there is no coline at all for pnew
     clold->addColoured(pnew, !anti);
   }
 }
 
 void HwRemDecayer::fixColours(PartnerMap partners, bool anti, 
 			      double colourDisrupt) const {
   PartnerMap::iterator prev;
   tPPtr pnew, pold;
   assert(partners.size()>=2);
 
   PartnerMap::iterator it=partners.begin();
   while(it != partners.end()) {
     //skip the first one to have a partner
     if(it==partners.begin()){
       it++;
       continue;
     }
     
     prev = it - 1;
     //determine the particles to work with
     pold = prev->first;
     if(prev->second) {
       if(!pold->coloured())
 	pold = prev->second;
       else if(pold->hasAntiColour() != anti)
 	pold = prev->second;
     }
     assert(pold);
 
     pnew = it->first;
     if(it->second) {
       if(it->second->colourLine(!anti)) //look for the opposite colour
 	pnew = it->second;
     }
     assert(pnew);
 
     // Implement the disruption of colour connections
     if( it != partners.end()-1 ) {//last one is diquark-has to be connected
       
       //has to be inside the if statement, so that the probability is
       //correctly counted:
       if( UseRandom::rnd() < colourDisrupt ){
 
 	if(!it->second){//check, whether we have a gluon
 	  mergeColour(pnew, pnew, anti);
 	}else{
 	  if(pnew==it->first)//be careful about the order
 	    mergeColour(it->second, it->first, anti);
 	  else
 	    mergeColour(it->first, it->second, anti);
 	}
 
 	it = partners.erase(it);
 	continue;
       }
 
     }
     // regular merging
     mergeColour(pold, pnew, anti);
     //end of loop
     it++;
   }
   return;
 }
 
 PPtr HwRemDecayer::forceSplit(const tRemPPtr rem, long child, Energy &lastQ, 
 			      double &lastx, Lorentz5Momentum &pf, 
 			      Lorentz5Momentum &p, 
 			      HadronContent & content) const {
   static const double eps=1e-6;
   // beam momentum
   Lorentz5Momentum beam = theBeam->momentum();
   // the last scale is minimum of last value and upper limit
   Energy minQ=_range*_kinCutoff*sqrt(lastx)/(1-lastx);
   if(minQ>lastQ) lastQ=minQ;
   // generate the new value of qtilde
   // weighted towards the lower value: dP/dQ = 1/Q -> Q(R) =
   // Q0 (Qmax/Q0)^R
   Energy q;
   unsigned int ntry=0,maxtry=100;
   double xExtracted = rem==theRems.first ? theX.first : theX.second;
   double zmin=  lastx/(1.-xExtracted) ,zmax,yy;
 
   if(1-lastx<eps) throw ShowerHandler::ExtraScatterVeto();
   do {
     q = minQ*pow(lastQ/minQ,UseRandom::rnd());
     yy   = 1.+0.5*sqr(_kinCutoff/q);
     zmax = yy - sqrt(sqr(yy)-1.); 
     ++ntry;
   }
   while(zmax<zmin&&ntry<maxtry);
   if(ntry==maxtry) throw ShowerHandler::ExtraScatterVeto();
   if(zmax-zmin<eps) throw ShowerHandler::ExtraScatterVeto();
   // now generate z as in FORTRAN HERWIG
   // use y = ln(z/(1-z)) as integration variable
   double ymin=log(zmin/(1.-zmin));
   double ymax=log(zmax/(1.-zmax));
   double dely=ymax-ymin;
   unsigned int nz=_nbinmax;
   dely/=nz;
   yy=ymin+0.5*dely;
   vector<int> ids;
   if(child==21||child==22) {
     ids=content.flav;
     for(unsigned int ix=0;ix<ids.size();++ix) ids[ix] *= content.sign;
   }
   else {
     ids.push_back(ParticleID::g);
   }
   // probabilities of the different types of possible splitting
   map<long,pair<double,vector<double> > > partonprob;
   double ptotal(0.);
   for(unsigned int iflav=0;iflav<ids.size();++iflav) {
     // only do each parton once
     if(partonprob.find(ids[iflav])!=partonprob.end()) continue;
     // particle data object
     tcPDPtr in = getParticleData(ids[iflav]);
     double psum(0.);
     vector<double> prob;
     for(unsigned int iz=0;iz<nz;++iz) {
       double ez=exp(yy);
       double wr=1.+ez;
       double zr=wr/ez;
       double wz=1./wr;
       double zz=wz*ez;
       double coup = child!=22 ? 
 	_alphaS ->value(sqr(max(wz*q,_kinCutoff))) :
 	_alphaEM->value(sqr(max(wz*q,_kinCutoff)));
       double az=wz*zz*coup;
       // g -> q qbar
       if(ids[iflav]==ParticleID::g) {
 	// calculate splitting function   
 	// SP as q is always less than forcedSplitScale, the pdf scale is fixed 
 	// pdfval = _pdf->xfx(theBeamData,in,sqr(q),lastx*zr); 
 	double pdfval=_pdf->xfx(theBeamData,in,sqr(_forcedSplitScale),lastx*zr);
 	if(pdfval>0.) psum += pdfval*az*0.5*(sqr(zz)+sqr(wz));
       }
       // q -> q g
       else {
 	// calculate splitting function
 	// SP as q is always less than forcedSplitScale, the pdf scale is fixed 
  	// pdfval = _pdf->xfx(theBeamData,in,sqr(q),lastx*zr);
  	double pdfval=_pdf->xfx(theBeamData,in,sqr(_forcedSplitScale),lastx*zr);
 	if(pdfval>0.) psum += pdfval*az*4./3.*(1.+sqr(wz))*zr;
       }
       if(psum>0.) prob.push_back(psum);
       yy+=dely;
     }
     if(psum>0.) partonprob[ids[iflav]] = make_pair(psum,prob);
     ptotal+=psum;
   }
   // select the flavour
   if(ptotal==0.) throw ShowerHandler::ExtraScatterVeto();
   ptotal *= UseRandom::rnd();
   map<long,pair<double,vector<double> > >::const_iterator pit;
   for(pit=partonprob.begin();pit!=partonprob.end();++pit) {
     if(pit->second.first>=ptotal) break;
     else                          ptotal -= pit->second.first;
   }
   if(pit==partonprob.end()) 
     throw Exception() << "Can't select parton for forced backward evolution in "
 		      << "HwRemDecayer::forceSplit" << Exception::eventerror;
   // select z 
   unsigned int iz=0;
   for(;iz<pit->second.second.size();++iz) {
     if(pit->second.second[iz]>ptotal) break;
   }
   if(iz==pit->second.second.size()) --iz;
   double ey=exp(ymin+dely*(float(iz+1)-UseRandom::rnd()));
   double z=ey/(1.+ey);
   Energy2 pt2=sqr((1.-z)*q)- z*sqr(_kinCutoff);
   // create the particle
   if(pit->first!=ParticleID::g) child=pit->first;
   PPtr parton = getParticleData(child)->produceParticle();
   Energy2 emittedm2 = sqr(parton->dataPtr()->constituentMass());
   // Now boost pcm and pf to z only frame
   Lorentz5Momentum p_ref = Lorentz5Momentum(ZERO,  beam.vect());
   Lorentz5Momentum n_ref = Lorentz5Momentum(ZERO, -beam.vect());
   // generate phi and compute pt of branching
   double phi = Constants::twopi*UseRandom::rnd();
   Energy pt=sqrt(pt2);
   Lorentz5Momentum qt   = LorentzMomentum(pt*cos(phi), pt*sin(phi), ZERO, ZERO);
   Axis axis(p_ref.vect().unit());
   if(axis.perp2()>0.) {
     LorentzRotation rot;
     double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
     rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
     qt.transform(rot);
   }
   // compute alpha for previous particle
   Energy2 p_dot_n  = p_ref*n_ref;
   double lastalpha =    pf*n_ref/p_dot_n;
   Lorentz5Momentum qtout=qt;
   Energy2 qtout2=-qt*qt;
   double alphaout=(1.-z)/z*lastalpha;
   double betaout=0.5*(emittedm2+qtout2)/alphaout/p_dot_n;
   Lorentz5Momentum k=alphaout*p_ref+betaout*n_ref+qtout;
   k.rescaleMass();
   parton->set5Momentum(k);
   pf+=k;
   lastQ=q;
   lastx/=z;
   p += parton->momentum();
   thestep->addDecayProduct(rem,parton,false);
   return parton;
 }
 
 void HwRemDecayer::setRemMasses() const {
   // get the masses of the remnants
   Energy mrem[2];
   Lorentz5Momentum ptotal,pnew[2];
   vector<tRemPPtr> theprocessed;
   theprocessed.push_back(theRems.first);
   theprocessed.push_back(theRems.second);
   // one remnant in e.g. DIS
   if(!theprocessed[0]||!theprocessed[1]) {
     tRemPPtr rem = theprocessed[0] ? theprocessed[0] : theprocessed[1];
     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(!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;
  
     //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) - (pt*pt);
 
   if( pz2/GeV2 < 0.0 ){
     if(dbg)
       cerr << "EXCEPTION not enough energy...." << endl;
     throw ExtraSoftScatterVeto();
   }
   
   if(!multiPeriph_){
   	if(UseRandom::rndbool()){
-	//TODO:fix
     		pz = sqrt(pz2);
 
   	}else
                 pz = -sqrt(pz2);
   }else{
         pz = pz2 > ZERO ? sqrt(pz2) : ZERO;
 
   }
   
 
   if(dbg)
     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 -= 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:
+   //  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(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++) {
+  for(unsigned int i = 0; i<2*N; i++) {
       // choose constituents
       Energy newMass = ZERO;
       Energy toyMass;
       if(i<2){
         // u and d have the same mass so its enough to use u 
         toyMass = getParticleData(ParticleID::u)->constituentMass();
       }
       else{
         toyMass = getParticleData(ParticleID::g)->constituentMass();
       }
       Lorentz5Momentum cp(ZERO,ZERO,ZERO,newMass,newMass);
       // dummy container for the momentum that is used for momentum conservation
       Lorentz5Momentum dummy(ZERO,ZERO,ZERO,toyMass,toyMass);
       ladderMomenta.push_back(cp);
       sumMomenta+=dummy;
   }
 
   // Get the beam energy
   tcPPair beam(generator()->currentEventHandler()->currentCollision()->incoming());
   Lorentz5Momentum P1(beam.first->momentum()), P2(beam.second->momentum());
 
   // Calculate available energy for the partons
   double x_1(0.0), x_2(0.0);
   double x1max = (r1.e()+abs(r1.z()))/(P1.e() + abs(P1.z()));
   double x2max = (r2.e()+abs(r2.z()))/(P2.e() + abs(P2.z()));
   double x1;
   double param = (1/(2*valOfN_+1))*initTotRap_;
   do{
        // Need 1-x instead of x to get the proper final momenta
        x1 = UseRandom::rndGauss(gaussWidth_, 1 - (exp(param)-1)/exp(param));
   }while(x1 < 0 || x1>=1.0);
        x_1 = x1max*x1;
        x_2 = x2max*x1;
 
   // Remnants 1 and 2 need to be rescaled later by this amount
   Lorentz5Momentum ig1 = x1*r1;
   Lorentz5Momentum ig2 = x1*r2;
 
   // The available energy that is used to generate the ladder
   // sumMomenta is the the sum of rest masses of the ladder partons
   // the available energy goes all into the kinematics
   Energy availableEnergy = (ig1+ig2).m() - sumMomenta.m();
  
   // 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
+      // Different options for soft pt sampling
       //1) pT1>pT2...pTN
       //2) pT1>pT2>..>pTN
       //3) flat
       //4) y dependent
+      //5) Frist then flat
       int triesPt=0;
       Energy pt;
       Energy ptTest;
-//      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!)
+      Energy2 ptp = pt*pt;
       if(ptp <= ZERO) pt = - sqrt(-ptp);
       else pt = sqrt(ptp);
       // randomize azimuth
       Energy px,py;
       //randomize the azimuth, but the last one should cancel all others
       if(i<ncl-1){
        randAzm(pt,px,py);
        // set transverse momentum
        mom[i].setX(px);
        mom[i].setY(py);
-//       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::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_ << ladderMult_ << ladderbFactor_ << pomeronStructure_
      << ounit(mg_,GeV) << ounit(ptmin_,GeV) << ounit(beta_,sqr(InvGeV))
      << 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_ >> ladderMult_ >> ladderbFactor_ >> pomeronStructure_
      >> iunit(mg_,GeV) >> iunit(ptmin_,GeV) >> iunit(beta_,sqr(InvGeV))
      >> 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/UnderlyingEvent/MPIHandler.cc b/UnderlyingEvent/MPIHandler.cc
--- a/UnderlyingEvent/MPIHandler.cc
+++ b/UnderlyingEvent/MPIHandler.cc
@@ -1,857 +1,857 @@
 // -*- C++ -*-
 //
 // MPIHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2017 The Herwig Collaboration
 //
 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
 //
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the MPIHandler class.
 //
 
 #include "MPIHandler.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Handlers/StandardXComb.h"
 #include "ThePEG/Handlers/SubProcessHandler.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Parameter.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/Interface/RefVector.h"
 #include "ThePEG/Interface/ParVector.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/MatrixElement/MEBase.h"
 #include "ThePEG/Handlers/CascadeHandler.h"
 #include "ThePEG/Cuts/Cuts.h"
 #include "ThePEG/Cuts/JetCuts.h"
 #include "ThePEG/Cuts/SimpleKTCut.h"
 #include "ThePEG/PDF/PartonExtractor.h"
 
 #include "gsl/gsl_sf_bessel.h"
 
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 MPIHandler * MPIHandler::currentHandler_ = 0;
 
 bool MPIHandler::beamOK() const {
   return (HadronMatcher::Check(*eventHandler()->incoming().first)  &&
 	  HadronMatcher::Check(*eventHandler()->incoming().second) );
 }
 
 
 tStdXCombPtr MPIHandler::generate(unsigned int sel) {
   //generate a certain process
   if(sel+1 > processHandlers().size())
     throw Exception() << "MPIHandler::generate called with argument out of range"
                       << Exception::runerror;
 
   return processHandlers()[sel]->generate();
 }
 
 IBPtr MPIHandler::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr MPIHandler::fullclone() const {
   return new_ptr(*this);
 }
 
 void MPIHandler::finalize() {
   if( beamOK() ){
     statistics();
   }
 }
 
 void MPIHandler::initialize() {
   currentHandler_ = this;
   useMe();
   theHandler = generator()->currentEventHandler(); 
   //stop if the EventHandler is not present:
   assert(theHandler);
 
   //check if MPI is wanted
   if( !beamOK() ){
     throw Exception()  << "You have requested multiple parton-parton scattering,\n"
 		       << "but the model is not forseen for the beam setup you chose.\n" 
 		       << "You should therefore disable that by setting XXXGenerator:EventHandler:"
 		       << "CascadeHandler:MPIHandler to NULL"
                        << Exception::runerror;
   }
 
   numSubProcs_ = subProcesses().size();
 
   if( numSubProcs_ != cuts().size() ) 
     throw Exception() << "MPIHandler::each SubProcess needs a Cuts Object"
 		      << "ReferenceVectors are not equal in size"
 		      << Exception::runerror;
 
   if( additionalMultiplicities_.size()+1 != numSubProcs_ )
     throw Exception() << "MPIHandler: each additional SubProcess needs "
 		      << "a multiplicity assigned. This can be done in with "
 		      << "insert MPIHandler:additionalMultiplicities 0 1"
 		      << Exception::runerror;
 
   //identicalToUE_ = 0 hard process is identical to ue, -1 no one
   if( identicalToUE_ > (int)numSubProcs_ || identicalToUE_ < -1 )
     throw Exception() << "MPIHandler:identicalToUE has disallowed value"
 		      << Exception::runerror;
 
   // override the cuts for the additional scatters if energyExtrapolation_ is
   // set
   if (energyExtrapolation_ != 0 ) {
     overrideUECuts();
   }
 
   tcPDPtr gluon=getParticleData(ParticleID::g);
   //determine ptmin
   Ptmin_ = cuts()[0]->minKT(gluon);
 
   if(identicalToUE_ == -1){
     algorithm_ = 2;
   }else{
     if(identicalToUE_ == 0){
       //Need to work a bit, in case of LesHouches events for QCD2to2
       Ptr<StandardEventHandler>::pointer eH =
 	dynamic_ptr_cast<Ptr<StandardEventHandler>::pointer>(eventHandler());
 
       if( eH ) {
 	PtOfQCDProc_ = eH->cuts()->minKT(gluon);
 	// find the jet cut in the new style cuts
 	for(unsigned int ix=0;ix<eH->cuts()->multiCuts().size();++ix) {
 	  Ptr<JetCuts>::pointer jetCuts =
 	    dynamic_ptr_cast<Ptr<JetCuts>::pointer>(eH->cuts()->multiCuts()[ix]);
 	  if(jetCuts) {
 	    Energy ptMin=1e30*GeV;
 	    for(unsigned int iy=0;iy<jetCuts->jetRegions().size();++iy) {
 	      ptMin = min(ptMin,jetCuts->jetRegions()[iy]->ptMin());
 	    }
 	    if(ptMin<1e29*GeV&&ptMin>PtOfQCDProc_) PtOfQCDProc_ = ptMin;
 	  }
 	}
       }
       else {
 	if(PtOfQCDProc_ == -1.0*GeV)
 	  throw Exception() << "MPIHandler: You need to specify the pt cutoff "
 			    << "used to in the LesHouches file for QCD2to2 events"
 			    << Exception::runerror;
       }
     }
     else {
       PtOfQCDProc_ = cuts()[identicalToUE_]->minKT(gluon);
     }
 
     if(PtOfQCDProc_ > 2*Ptmin_)
       algorithm_ = 1;
     else
       algorithm_ = 0;
 
     if(PtOfQCDProc_ == ZERO)//pure MinBias mode
       algorithm_ = -1;
   }
 
   //Init all subprocesses
   for(unsigned int i=0; i<numSubProcs_; i++){
     theProcessHandlers.push_back(new_ptr(ProcessHandler()));
     processHandlers().back()->initialize(subProcesses()[i], 
 					 cuts()[i], eventHandler());
     processHandlers().back()->initrun();
   }
 
   //now calculate the individual Probabilities
   XSVector UEXSecs;
   UEXSecs.push_back(processHandlers()[0]->integratedXSec());
   //save the hard cross section
   hardXSec_ = UEXSecs.front();
 
   //determine sigma_soft and beta
   if(softInt_){//check that soft ints are requested
     GSLBisection rootFinder;
 
     if(twoComp_){
       //two component model
       /*
       GSLMultiRoot eqSolver;
       slopeAndTotalXSec eq(this);
       pair<CrossSection, Energy2> res = eqSolver.value(eq, 10*millibarn, 0.6*GeV2);
       softXSec_ = res.first;
       softMu2_ = res.second;
       */
       slopeBisection fs(this);
       try{
 	softMu2_ = rootFinder.value(fs, 0.3*GeV2, 1.*GeV2);
 	softXSec_ = fs.softXSec();
       }catch(GSLBisection::IntervalError){
         try{
 	// Very low energies (e.g. 200 GeV) need this.
 	// Very high energies (e.g. 100 TeV) produce issues with
 	// this choice. This is a temp. fix.  
           softMu2_ = rootFinder.value(fs, 0.25*GeV2, 1.3*GeV2);
           softXSec_ = fs.softXSec();
         }catch(GSLBisection::IntervalError){
 	throw Exception() <<
         "\n**********************************************************\n"	  "* Inconsistent MPI parameter choice for this beam energy *\n"
 	  "**********************************************************\n"
 	  "MPIHandler parameter choice is unable to reproduce\n"
 	  "the total cross section. Please check arXiv:0806.2949\n"
 	  "for the allowed parameter space."
 			  << Exception::runerror;
         }
       }
 
     }else{
       //single component model
       TotalXSecBisection fn(this);
       try{
 	softXSec_ = rootFinder.value(fn, 0*millibarn, 5000*millibarn);
       }catch(GSLBisection::IntervalError){
 	throw Exception() <<
 	"\n**********************************************************\n"
 	  "* Inconsistent MPI parameter choice for this beam energy *\n"
 	  "**********************************************************\n"
 	  "MPIHandler parameter choice is unable to reproduce\n"
 	  "the total cross section. Please check arXiv:0806.2949\n"
 	  "for the allowed parameter space."
 			  << Exception::runerror;      
       }
     }
 
     //now get the differential cross section at ptmin
     ProHdlPtr qcd = new_ptr(ProcessHandler());
     Energy eps = 0.1*GeV;
     Energy ptminPlus = Ptmin_ + eps;
     
     Ptr<SimpleKTCut>::pointer ktCut = new_ptr(SimpleKTCut(ptminPlus));
     ktCut->init();
     ktCut->initrun();
 
     CutsPtr qcdCut = new_ptr(Cuts(2*ptminPlus));
     qcdCut->add(dynamic_ptr_cast<tOneCutPtr>(ktCut));
     qcdCut->init();
     qcdCut->initrun();
       
     qcd->initialize(subProcesses()[0], qcdCut, eventHandler());
     qcd->initrun();
     
     // ds/dp_T^2 = 1/2/p_T ds/dp_T
     DiffXSec hardPlus = (hardXSec_-qcd->integratedXSec())/(2*Ptmin_*eps);
     
     betaBisection fn2(softXSec_, hardPlus, Ptmin_);
     try{
       beta_ = rootFinder.value(fn2, -10/GeV2, 2/GeV2);
     }catch(GSLBisection::IntervalError){
       try{
 	// Very low energies (e.g. 200 GeV) need this.
 	// Very high energies (e.g. 100 TeV) produce issues with
 	// this choice. This is a temp. fix.
         beta_ = rootFinder.value(fn2, -5/GeV2, 8./GeV2);
       }catch(GSLBisection::IntervalError){
       throw Exception() << "MPIHandler: slope of soft pt spectrum couldn't be "
 			<< "determined."
 			<< Exception::runerror;    
       }
     }
   }
 
   Probs(UEXSecs);
   //MultDistribution("probs.test");
 
   UEXSecs.clear();
 }
 
 void MPIHandler::MultDistribution(string filename) const {
   ofstream file;
   double p(0.0), pold(0.0);
   file.open(filename.c_str());
   //theMultiplicities  
   Selector<MPair>::const_iterator it = theMultiplicities.begin();
 
   while(it != theMultiplicities.end()){
     p = it->first;
     file << it->second.first << " " << it->second.second
 	 << " " << p-pold << '\n';
     it++;
     pold = p;
   }
 
   file << "sum of all probabilities: " << theMultiplicities.sum() 
        << endl;
 
   file.close();
 }
 
 void MPIHandler::statistics() const {
   ostream & file = generator()->misc();
   
   string line = "======================================="
     "=======================================\n";
 
   for(unsigned int i=0; i<cuts().size(); i++){
     Stat tot;
     if(i == 0)
       file << "Statistics for the UE process: \n";
     else
       file << "Statistics for additional hard Process " << i << ": \n";
 
     processHandlers()[i]->statistics(file, tot);
     file << "\n";
   }
 
   if(softInt_){
     file << line
 	 << "Eikonalized and soft cross sections:\n\n"
 	 << "Model parameters:                    "
 	 << "ptmin:   " << Ptmin_/GeV << " GeV"
       	 << ", mu2: " << invRadius_/sqr(1.*GeV) << " GeV2\n"
 	 << "                                     "
 	 << "DL mode: " << DLmode_
 	 << ", CMenergy: " << generator()->maximumCMEnergy()/GeV
 	 << " GeV" << '\n'
 	 << "hard inclusive cross section (mb):   "
 	 << hardXSec_/millibarn << '\n'
 	 << "soft inclusive cross section (mb):   "
 	 << softXSec_/millibarn << '\n'
 	 << "total cross section (mb):            "
 	 << totalXSecExp()/millibarn << '\n'
 	 << "inelastic cross section (mb):        "
 	 << inelXSec_/millibarn << '\n'
 	 << "soft inv radius (GeV2):              "
 	 << softMu2_/GeV2 << '\n'
 	 << "slope of soft pt spectrum (1/GeV2):  "
 	 << beta_*sqr(1.*GeV) << '\n'
 	 << "Average hard multiplicity:           "
 	 << avgNhard_ << '\n'
 	 << "Average soft multiplicity:           "
 	 << avgNsoft_ << '\n' << line << endl;
   }else{
     file << line
 	 << "Eikonalized and soft cross sections:\n\n"
 	 << "Model parameters:                    "
 	 << "ptmin:   " << Ptmin_/GeV << " GeV"
       	 << ", mu2: " << invRadius_/sqr(1.*GeV) << " GeV2\n"
 	 << "                                     "
 	 << ", CMenergy: " << generator()->maximumCMEnergy()/GeV
 	 << " GeV" << '\n'
 	 << "hard inclusive cross section (mb):   "
 	 << hardXSec_/millibarn << '\n'
 	 << "Average hard multiplicity:           "
 	 << avgNhard_ << '\n' << line << endl;
   }
 }
 
 unsigned int MPIHandler::multiplicity(unsigned int sel){
   if(sel==0){//draw from the pretabulated distribution
     MPair m = theMultiplicities.select(UseRandom::rnd());
     softMult_ = m.second;
     return m.first;
   } else{ //fixed multiplicities for the additional hard scatters
     if(additionalMultiplicities_.size() < sel)
       throw Exception() << "MPIHandler::multiplicity: process index "
 			<< "is out of range"
 			<< Exception::runerror;
 
     return additionalMultiplicities_[sel-1];
   }
 }
 
 
 void MPIHandler::Probs(XSVector UEXSecs) {
   GSLIntegrator integrator;
   unsigned int iH(1), iS(0);
   double P(0.0);
   double P0(0.0);//the probability for i hard and zero soft scatters
   Length bmax(500.0*sqrt(millibarn));
   //only one UE process will be drawn from a probability distribution,
   //so check that.
   assert(UEXSecs.size() == 1);
 
   for ( XSVector::const_iterator it = UEXSecs.begin();
         it != UEXSecs.end(); ++it ) {
 
     iH = 0; 
 
     //get the inel xsec
     Eikonalization inelint(this, -1, *it, softXSec_, softMu2_);
     inelXSec_ = integrator.value(inelint, ZERO, bmax);
 
     avgNhard_ = 0.0;
     avgNsoft_ = 0.0;
     bmax = 10.0*sqrt(millibarn);
     do{//loop over hard ints
       if(Algorithm()>-1 && iH==0){
 	iH++;
 	continue;
       }
       iS = 0;
       do{//loop over soft ints
 	if( ( Algorithm() == -1 && iS==0 && iH==0 ) ){
 	  iS++;
 	  continue;
 	}
 
 	Eikonalization integrand(this, iH*100+iS, *it, softXSec_, softMu2_);
       
 	if(Algorithm() > 0){
 	  P = integrator.value(integrand, ZERO, bmax) / *it;
 	}else{
 	  P = integrator.value(integrand, ZERO, bmax) / inelXSec_;
 	}
 	//store the probability
 	if(Algorithm()>-1){
 	  theMultiplicities.insert(P, make_pair(iH-1, iS));
 	  avgNhard_ += P*(iH-1);
 	}else{
 	  theMultiplicities.insert(P, make_pair(iH, iS));
 	  avgNhard_ += P*(iH);
 	}
 	avgNsoft_ += P*iS;
 	if(iS==0)
 	  P0 = P;
 
 	iS++;
       } while ( (iS < maxScatters_) && (iS < 5 || P > 1.e-15 ) && softInt_ );
       iH++;
     } while ( (iH < maxScatters_) && (iH < 5 || P0 > 1.e-15) );
   }
 }
 
 
 // calculate the integrand
 Length Eikonalization::operator() (Length b) const {
   unsigned int Nhard(0), Nsoft(0);
   //fac is just: d^2b=fac*db despite that large number
   Length fac(Constants::twopi*b);
   
   double chiTot(( theHandler->OverlapFunction(b)*hardXSec_ + 
 		  theHandler->OverlapFunction(b, softMu2_)*softXSec_ ) / 2.0);
 
   //total cross section wanted
   if(theoption == -2) return 2 * fac * ( 1 - exp(-chiTot) );
   //inelastic cross section
   if(theoption == -1) return fac * ( 1 - exp(- 2.0 * chiTot) );
 
   if(theoption >= 0){
     //encode multiplicities as: N_hard*100 + N_soft   
     Nhard = theoption/100;
     Nsoft = theoption%100;
 
     if(theHandler->Algorithm() > 0){
       //P_n*sigma_hard: n-1 extra scatters + 1 hard scatterer != hardXSec_
       return fac * Nhard * theHandler->poisson(b, hardXSec_, Nhard) *
 	theHandler->poisson(b, softXSec_, Nsoft, softMu2_);
     }else{
       //P_n*sigma_inel: n extra scatters
       return fac * theHandler->poisson(b, hardXSec_, Nhard) *
 	theHandler->poisson(b, softXSec_, Nsoft, softMu2_);
     }
 
   }else{
     throw Exception() << "Parameter theoption in Struct Eikonalization in " 
 		      << "MPIHandler.cc has not allowed value"
                       << Exception::runerror;
     return 0.0*meter;
   }
 }
 
 InvEnergy2 slopeBisection::operator() (Energy2 softMu2) const {
   GSLBisection root;
   //single component model
   TotalXSecBisection fn(handler_, softMu2);
   try{
     softXSec_ = root.value(fn, 0*millibarn, 5000*millibarn);
   }catch(GSLBisection::IntervalError){
     throw Exception() << "MPIHandler 2-Component model didn't work out."
 		      << Exception::runerror;
   }
   return handler_->slopeDiff(softXSec_, softMu2);
 }
 
 LengthDiff slopeInt::operator() (Length b) const {
   //fac is just: d^2b=fac*db
   Length fac(Constants::twopi*b);
   
   double chiTot(( handler_->OverlapFunction(b)*hardXSec_ + 
 		  handler_->OverlapFunction(b, softMu2_)*softXSec_ ) / 2.0);
 
   InvEnergy2 b2 = sqr(b/hbarc);
   //B*sigma_tot
   return fac * b2 * ( 1 - exp(-chiTot) );
 }
 
 double MPIHandler::factorial (unsigned int n) const {
 
   static double f[] = {1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880.,3.6288e6,
 		3.99168e7,4.790016e8,6.2270208e9,8.71782912e10,1.307674368e12,
 		2.0922789888e13,3.55687428096e14,6.402373705728e15,1.21645100408832e17,
 		2.43290200817664e18,5.10909421717094e19,1.12400072777761e21,
 		2.5852016738885e22,6.20448401733239e23,1.5511210043331e25,
 		4.03291461126606e26,1.08888694504184e28,3.04888344611714e29,
 		8.8417619937397e30,2.65252859812191e32,8.22283865417792e33,
 		2.63130836933694e35,8.68331761881189e36,2.95232799039604e38,
 		1.03331479663861e40,3.71993326789901e41,1.37637530912263e43,
 		5.23022617466601e44,2.03978820811974e46,8.15915283247898e47,
 		3.34525266131638e49,1.40500611775288e51,6.04152630633738e52,
 		2.65827157478845e54,1.1962222086548e56,5.50262215981209e57,
 		2.58623241511168e59,1.24139155925361e61,6.08281864034268e62,
 		3.04140932017134e64,1.55111875328738e66,8.06581751709439e67,
 		4.27488328406003e69,2.30843697339241e71,1.26964033536583e73,
 		7.10998587804863e74,4.05269195048772e76,2.35056133128288e78,
 		1.3868311854569e80,8.32098711274139e81,5.07580213877225e83,
 		3.14699732603879e85,1.98260831540444e87,1.26886932185884e89,
 		8.24765059208247e90,5.44344939077443e92,3.64711109181887e94,
 		2.48003554243683e96,1.71122452428141e98,1.19785716699699e100,
 		8.50478588567862e101,6.12344583768861e103,4.47011546151268e105,
 		3.30788544151939e107,2.48091408113954e109,1.88549470166605e111,
 		1.45183092028286e113,1.13242811782063e115,8.94618213078298e116,
 		7.15694570462638e118,5.79712602074737e120,4.75364333701284e122,
 		3.94552396972066e124,3.31424013456535e126,2.81710411438055e128,
 		2.42270953836727e130,2.10775729837953e132,1.85482642257398e134,
 		1.65079551609085e136,1.48571596448176e138,1.3520015276784e140,
 		1.24384140546413e142,1.15677250708164e144,1.08736615665674e146,
 		1.03299784882391e148,9.9167793487095e149,9.61927596824821e151,
 		9.42689044888325e153,9.33262154439442e155,9.33262154439442e157};
 
   if(n > maxScatters_) 
         throw Exception() << "MPIHandler::factorial called with too large argument"
                       << Exception::runerror;
   else
     return f[n];
 }
 
 InvArea MPIHandler::OverlapFunction(Length b, Energy2 mu2) const {
   if(mu2 == ZERO)
     mu2 = invRadius_;
 
   InvLength mu = sqrt(mu2)/hbarc;
   return (sqr(mu)/96/Constants::pi)*pow(mu*b, 3)*(gsl_sf_bessel_Kn(3, mu*b));
 }
 
 double MPIHandler::poisson(Length b, CrossSection sigma, unsigned int N, Energy2 mu2) const {
   if(sigma > 0*millibarn){
     return pow(OverlapFunction(b, mu2)*sigma, (double)N)/factorial(N)
       *exp(-OverlapFunction(b, mu2)*sigma);
   }else{
     return (N==0) ? 1.0 : 0.0;
   }
 }
 
 CrossSection MPIHandler::totalXSecDiff(CrossSection softXSec, 
 				       Energy2 softMu2) const {
   GSLIntegrator integrator;
   Eikonalization integrand(this, -2, hardXSec_, softXSec, softMu2);
   Length bmax = 500.0*sqrt(millibarn);
 
   CrossSection tot = integrator.value(integrand, ZERO, bmax);
   return (tot-totalXSecExp());
 }
 
 InvEnergy2 MPIHandler::slopeDiff(CrossSection softXSec, 
 				 Energy2 softMu2) const {
   GSLIntegrator integrator;
   Eikonalization integrand(this, -2, hardXSec_, softXSec, softMu2);
   Length bmax = 500.0*sqrt(millibarn);
 
   CrossSection tot = integrator.value(integrand, ZERO, bmax);
   
   slopeInt integrand2(this, hardXSec_, softXSec, softMu2);
   
   return integrator.value(integrand2, ZERO, bmax)/tot - slopeExp();
 }
 
 CrossSection MPIHandler::totalXSecExp() const {
   if(totalXSecExp_ != 0*millibarn)
     return totalXSecExp_;
 
   double pom_old = 0.0808;
   CrossSection coef_old = 21.7*millibarn;
   double pom_new_hard = 0.452;
   CrossSection coef_new_hard = 0.0139*millibarn;
   double pom_new_soft = 0.0667;
   CrossSection coef_new_soft = 24.22*millibarn;
 
   Energy energy(generator()->maximumCMEnergy());
 
   switch(DLmode_){
   case 1://old DL extrapolation
     return coef_old * pow(energy/GeV, 2*pom_old);
     break;
 
   case 2://old DL extrapolation fixed to CDF
     return 81.8*millibarn * pow(energy/1800.0/GeV, 2*pom_old);
     break;
     
   case 3://new DL extrapolation
-    return coef_new_hard * pow(energy/GeV, 2*pom_new_hard) + 
+    return 0.8*coef_new_hard * pow(energy/GeV, 2*pom_new_hard) + 
       coef_new_soft * pow(energy/GeV, 2*pom_new_soft);
     break;
     
   default:
     throw Exception() << "MPIHandler::totalXSecExp non-existing mode selected"
                       << Exception::runerror;   
   }
 }
 
 InvEnergy2 MPIHandler::slopeExp() const{
   //Currently return the slope as calculated in the pomeron fit by
   //Donnachie & Landshoff
   Energy energy(generator()->maximumCMEnergy());
   //slope at
   Energy e_0 = 1800*GeV;
   InvEnergy2 b_0 = 17/GeV2;
   return b_0 + log(energy/e_0)/GeV2;
 }
 
 void MPIHandler::overrideUECuts() {
   if(energyExtrapolation_==1)
     Ptmin_ = EEparamA_ * log(generator()->maximumCMEnergy() / EEparamB_);
   else if(energyExtrapolation_==2)
     Ptmin_ = pT0_*pow(double(generator()->maximumCMEnergy()/refScale_),b_);
   else
     assert(false);
   // create a new SimpleKTCut object with the calculated ptmin value
   Ptr<SimpleKTCut>::pointer newUEktCut = new_ptr(SimpleKTCut(Ptmin_));
   newUEktCut->init();
   newUEktCut->initrun();
 
   // create a new Cuts object with MHatMin = 2 * Ptmin_
   CutsPtr newUEcuts = new_ptr(Cuts(2*Ptmin_));
   newUEcuts->add(dynamic_ptr_cast<tOneCutPtr>(newUEktCut));
   newUEcuts->init();
   newUEcuts->initrun();
 
   // replace the old Cuts object
   cuts()[0] = newUEcuts;
 }
 
 void MPIHandler::persistentOutput(PersistentOStream & os) const {
   os << theMultiplicities << theHandler
      << theSubProcesses << theCuts << theProcessHandlers
      << additionalMultiplicities_ << identicalToUE_ 
      << ounit(PtOfQCDProc_, GeV) << ounit(Ptmin_, GeV) 
      << ounit(hardXSec_, millibarn) << ounit(softXSec_, millibarn)
      << ounit(beta_, 1/GeV2)
      << algorithm_ << ounit(invRadius_, GeV2)
      << numSubProcs_ << colourDisrupt_ << softInt_ << twoComp_ 
      << DLmode_ << ounit(totalXSecExp_, millibarn)
      << energyExtrapolation_ << ounit(EEparamA_, GeV) << ounit(EEparamB_, GeV)
      << ounit(refScale_,GeV) << ounit(pT0_,GeV) << b_
      << avgNhard_ << avgNsoft_ << softMult_ 
      << ounit(inelXSec_, millibarn) 
      << ounit(softMu2_, GeV2);
 }
 
 void MPIHandler::persistentInput(PersistentIStream & is, int) {
   is >> theMultiplicities >> theHandler
      >> theSubProcesses >> theCuts >> theProcessHandlers
      >> additionalMultiplicities_ >> identicalToUE_ 
      >> iunit(PtOfQCDProc_, GeV) >> iunit(Ptmin_, GeV)
      >> iunit(hardXSec_, millibarn) >> iunit(softXSec_, millibarn)
      >> iunit(beta_, 1/GeV2)
      >> algorithm_ >> iunit(invRadius_, GeV2)
      >> numSubProcs_ >> colourDisrupt_ >> softInt_ >> twoComp_ 
      >> DLmode_ >> iunit(totalXSecExp_, millibarn)
      >> energyExtrapolation_ >> iunit(EEparamA_, GeV) >> iunit(EEparamB_, GeV)
      >> iunit(refScale_,GeV) >> iunit(pT0_,GeV) >> b_
      >> avgNhard_ >> avgNsoft_ >> softMult_ 
      >> iunit(inelXSec_, millibarn) 
      >> iunit(softMu2_, GeV2);
   currentHandler_ = this;
 }
 
 void MPIHandler::clean() {
   // ThePEG's event handler's usual event cleanup doesn't reach these
   // XCombs. Need to do it by hand here.
   for ( size_t i = 0; i < theSubProcesses.size(); ++i ) {
     theSubProcesses[i]->pExtractor()->lastXCombPtr()->clean();
   }
 }
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<MPIHandler,Interfaced>
 describeHerwigMPIHandler("Herwig::MPIHandler",
 			 "JetCuts.so SimpleKTCut.so HwMPI.so");
 
 void MPIHandler::Init() {
 
   static ClassDocumentation<MPIHandler> documentation
     ("The MPIHandler class is the main administrator of the multiple interaction model", 
      "The underlying event was simulated with an eikonal model for multiple partonic interactions."
      "Details can be found in Ref.~\\cite{Bahr:2008dy,Bahr:2009ek}.", 
      "%\\cite{Bahr:2008dy}\n"
      "\\bibitem{Bahr:2008dy}\n"
      "  M.~Bahr, S.~Gieseke and M.~H.~Seymour,\n"
      "  ``Simulation of multiple partonic interactions in Herwig,''\n"
      "  JHEP {\\bf 0807}, 076 (2008)\n"
      "  [arXiv:0803.3633 [hep-ph]].\n"
      "  %%CITATION = JHEPA,0807,076;%%\n"
     "\\bibitem{Bahr:2009ek}\n"
      "  M.~Bahr, J.~M.~Butterworth, S.~Gieseke and M.~H.~Seymour,\n"
      "  ``Soft interactions in Herwig,''\n"
      "  arXiv:0905.4671 [hep-ph].\n"
      "  %%CITATION = ARXIV:0905.4671;%%\n"
      );
   
   static RefVector<MPIHandler,SubProcessHandler> interfaceSubhandlers
     ("SubProcessHandlers",
      "The list of sub-process handlers used in this EventHandler. ",
      &MPIHandler::theSubProcesses, -1, false, false, true, false, false);
 
   static RefVector<MPIHandler,Cuts> interfaceCuts
     ("Cuts",
      "List of cuts used for the corresponding list of subprocesses. These cuts "
      "should not be overidden in individual sub-process handlers.",
      &MPIHandler::theCuts, -1, false, false, true, false, false);
 
   static Parameter<MPIHandler,Energy2> interfaceInvRadius
     ("InvRadius",
      "The inverse hadron radius squared used in the overlap function",
      &MPIHandler::invRadius_, GeV2, 2.0*GeV2, 0.2*GeV2, 4.0*GeV2,
      true, false, Interface::limited);
 
   static ParVector<MPIHandler,int> interfaceadditionalMultiplicities
     ("additionalMultiplicities",
      "specify the multiplicities of secondary hard processes (multiple parton scattering)",
      &MPIHandler::additionalMultiplicities_, 
      -1, 0, 0, 3,
      false, false, true);
 
   static Parameter<MPIHandler,int> interfaceIdenticalToUE
     ("IdenticalToUE",
      "Specify which of the hard processes is identical to the UE one (QCD dijets)",
      &MPIHandler::identicalToUE_, -1, 0, 0,
      false, false, Interface::nolimits);
 
   static Parameter<MPIHandler,Energy> interfacePtOfQCDProc
     ("PtOfQCDProc",
      "Specify the value of the pt cutoff for the process that is identical to the UE one",
      &MPIHandler::PtOfQCDProc_, GeV, -1.0*GeV, ZERO, ZERO,
      false, false, Interface::nolimits);
 
   static Parameter<MPIHandler,double> interfacecolourDisrupt
     ("colourDisrupt",
      "Fraction of connections to additional subprocesses, which are colour disrupted.",
      &MPIHandler::colourDisrupt_, 
      0.0, 0.0, 1.0, 
      false, false, Interface::limited);
 
   
   static Switch<MPIHandler,bool> interfacesoftInt
     ("softInt",
      "Switch to enable soft interactions",
      &MPIHandler::softInt_, true, false, false);
 
   static SwitchOption interfacesoftIntYes
     (interfacesoftInt,
      "Yes",
      "enable the two component model",
      true);
   static SwitchOption interfacesoftIntNo
     (interfacesoftInt,
      "No",
      "disable the model",
      false);
 
 
   static Switch<MPIHandler,unsigned int> interEnergyExtrapolation
     ("EnergyExtrapolation",
      "Switch to ignore the cuts object at MPIHandler:Cuts[0]. "
      "Instead, extrapolate the pt cut.",
      &MPIHandler::energyExtrapolation_, 2, false, false);
   static SwitchOption interEnergyExtrapolationLog
     (interEnergyExtrapolation,
      "Log",
      "Use logarithmic dependence, ptmin = A * log (sqrt(s) / B).",
      1);
   static SwitchOption interEnergyExtrapolationPower
     (interEnergyExtrapolation,
      "Power",
      "Use power law, ptmin = pt_0 * (sqrt(s) / E_0)^b.",
      2);
   static SwitchOption interEnergyExtrapolationNo
     (interEnergyExtrapolation,
      "No",
      "Use manually set value for the minimal pt, "
      "specified in MPIHandler:Cuts[0]:OneCuts[0]:MinKT.",
      0);
 
   static Parameter<MPIHandler,Energy> interfaceEEparamA
     ("EEparamA",
      "Parameter A in the empirical parametrization "
      "ptmin = A * log (sqrt(s) / B)",
      &MPIHandler::EEparamA_, GeV, 0.6*GeV, ZERO, Constants::MaxEnergy,
      false, false, Interface::limited);
 
   static Parameter<MPIHandler,Energy> interfaceEEparamB
     ("EEparamB",
      "Parameter B in the empirical parametrization "
      "ptmin = A * log (sqrt(s) / B)",
      &MPIHandler::EEparamB_, GeV, 39.0*GeV, ZERO, Constants::MaxEnergy,
      false, false, Interface::limited);
 
   static Switch<MPIHandler,bool> interfacetwoComp
     ("twoComp",
      "switch to enable the model with a different radius for soft interactions",
      &MPIHandler::twoComp_, true, false, false);
 
   static SwitchOption interfacetwoCompYes
     (interfacetwoComp,
      "Yes",
      "enable the two component model",
      true);
   static SwitchOption interfacetwoCompNo
     (interfacetwoComp,
      "No",
      "disable the model",
      false);
 
 
   static Parameter<MPIHandler,CrossSection> interfaceMeasuredTotalXSec
     ("MeasuredTotalXSec",
      "Value for the total cross section (assuming rho=0). If non-zero, this "
      "overwrites the Donnachie-Landshoff parametrizations.",
      &MPIHandler::totalXSecExp_, millibarn, 0.0*millibarn, 0.0*millibarn, 0*millibarn,
      false, false, Interface::lowerlim);
 
   
   static Switch<MPIHandler,unsigned int> interfaceDLmode
     ("DLmode",
      "Choice of Donnachie-Landshoff parametrization for the total cross section.",
      &MPIHandler::DLmode_, 2, false, false);
   static SwitchOption interfaceDLmodeStandard
     (interfaceDLmode,
      "Standard",
      "Standard parametrization with s**0.08",
      1);
   static SwitchOption interfaceDLmodeCDF
     (interfaceDLmode,
      "CDF",
      "Standard parametrization but normalization fixed to CDF's measured value",
      2);
   static SwitchOption interfaceDLmodeNew
     (interfaceDLmode,
      "New",
      "Parametrization taking hard and soft pomeron contributions into account",
      3);
 
   static Parameter<MPIHandler,Energy> interfaceReferenceScale
     ("ReferenceScale",
      "The reference energy for power law energy extrapolation of pTmin",
      &MPIHandler::refScale_, GeV, 7000.0*GeV, 0.0*GeV, 20000.*GeV,
      false, false, Interface::limited);
 
   static Parameter<MPIHandler,Energy> interfacepTmin0
     ("pTmin0",
      "The pTmin at the reference scale for power law extrapolation of pTmin.",
      &MPIHandler::pT0_, GeV, 3.11*GeV, 0.0*GeV, 10.0*GeV,
      false, false, Interface::limited);
 
   static Parameter<MPIHandler,double> interfacePower
     ("Power",
      "The power for power law extrapolation of the pTmin cut-off.",
      &MPIHandler::b_, 0.21, 0.0, 10.0,
      false, false, Interface::limited);
 }
diff --git a/src/snippets/MB.in b/src/snippets/MB.in
--- a/src/snippets/MB.in
+++ b/src/snippets/MB.in
@@ -1,43 +1,47 @@
 ##################################################
 # MEMinBias Matrix Element
 ##################################################
 
 # MPI model settings
 set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0
 
 ## Report the correct cross section
 cd /Herwig/Generators
 create Herwig::MPIXSecReweighter MPIXSecReweighter
 insert EventGenerator:EventHandler:PostSubProcessHandlers 0 MPIXSecReweighter
 set EventGenerator:EventHandler:CascadeHandler NULL
 
 clear EventGenerator:EventHandler:SubProcessHandlers[0]
 
 ##################################################
 # Create separate SubProcessHandler for MinBias
 ##################################################
 cd /Herwig/MatrixElements/
 cp SubProcess QCDMinBias
 
 set QCDMinBias:CascadeHandler /Herwig/Shower/ShowerHandler
 set QCDMinBias:CascadeHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
 set QCDMinBias:DecayHandler /Herwig/Decays/DecayHandler
 
 # Due to numerics the pomeron could be seen as timelike.
 set /Herwig/Shower/ShowerHandler:SplitHardProcess No
 set /Herwig/DipoleShower/DipoleShowerHandler:SplitHardProcess No
 
 insert QCDMinBias:MatrixElements[0] MEMinBias
 
 cd /Herwig/Generators
 # MinBias parameters used for the new kinematics of soft MPI
-set /Herwig/Cuts/MinBiasCuts:X1Min 0.11
-set /Herwig/Cuts/MinBiasCuts:X2Min 0.11
+set /Herwig/Cuts/MinBiasCuts:X1Min 0.011
+set /Herwig/Cuts/MinBiasCuts:X2Min 0.011
+
+
 
 # Needed to get the correct fraction of diffractive events
 set /Herwig/MatrixElements/MEMinBias:csNorm 4.5584
 
+set /Herwig/MatrixElements/MEMinBias:Scale 2.0
+
 set EventGenerator:EventHandler:Cuts /Herwig/Cuts/MinBiasCuts
 
 cd /Herwig/MatrixElements/
 insert /Herwig/Generators/EventGenerator:EventHandler:SubProcessHandlers[0] QCDMinBias