No OneTemporary
Actions

Size

193 KB

Subscribers

None

View Options

	diff --git a/Decay/Radiation/IFDipole.cc b/Decay/Radiation/IFDipole.cc
	--- a/Decay/Radiation/IFDipole.cc
	+++ b/Decay/Radiation/IFDipole.cc
	@@ -1,709 +1,706 @@
	// -- C++ --
	//
	// IFDipole.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the IFDipole class.
	//

	#include "IFDipole.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"

	using namespace ThePEG::Helicity;

	using namespace Herwig;

	void IFDipole::persistentOutput(PersistentOStream & os) const {
	os << _alpha << ounit(_emin,GeV) << _maxwgt
	<< _mode << _maxtry << _energyopt << _betaopt;
	}

	void IFDipole::persistentInput(PersistentIStream & is, int) {
	is >> _alpha >> iunit(_emin,GeV) >> _maxwgt
	>> _mode >> _maxtry >> _energyopt >> _betaopt;
	}

	ClassDescription<IFDipole> IFDipole::initIFDipole;
	// Definition of the static class description member.

	void IFDipole::Init() {
	static ClassDocumentation<IFDipole> documentation
	("The IFDipole class implements the initial-final dipole for the SOPTHY algorithm");

	static Switch<IFDipole,unsigned int> interfaceUnWeight
	("UnWeight",
	"Control the type of unweighting to perform, only one should be used the"
	" other options are for debugging purposes.",
	&IFDipole::_mode, 1, false, false);
	static SwitchOption interfaceUnWeightNoUnweighting
	(interfaceUnWeight,
	"NoUnweighting",
	"Perform no unweighting",
	0);
	static SwitchOption interfaceUnWeightAllWeights
	(interfaceUnWeight,
	"AllWeights",
	"Include all the weights",
	1);
	static SwitchOption interfaceUnWeightNoJacobian
	(interfaceUnWeight,
	"NoJacobian",
	"Only include the dipole and YFS weights",
	2);
	static SwitchOption interfaceUnWeightDipole
	(interfaceUnWeight,
	"Dipole",
	"Only include the dipole weight",
	3);
	static SwitchOption interfaceUnWeightYFS
	(interfaceUnWeight,
	"YFS",
	"Only include the YFS weight",
	4);

	static Parameter<IFDipole,unsigned int> interfaceMaximumTries
	("MaximumTries",
	"Maximum number of attempts to unweight",
	&IFDipole::_maxtry, 500, 10, 100000,
	false, false, Interface::limited);

	static Parameter<IFDipole,Energy> interfaceMinimumEnergyRest
	("MinimumEnergyRest",
	"The minimum energy of the photons in the rest frame of the decaying particle",
	&IFDipole::_emin, MeV, 1.MeV, ZERO, 10000.0MeV,
	false, false, Interface::limited);

	static Parameter<IFDipole,double> interfaceMaximumWeight
	("MaximumWeight",
	"The maximum weight for unweighting",
	&IFDipole::_maxwgt, 2.0, 0.0, 100.0,
	false, false, Interface::limited);

	static Switch<IFDipole,unsigned int> interfaceEnergyCutOff
	("EnergyCutOff",
	"The type of cut-off on the photon energy to apply",
	&IFDipole::_energyopt, 1, false, false);
	static SwitchOption interfaceEnergyCutOffRestFrame
	(interfaceEnergyCutOff,
	"RestFrame",
	"Apply cut-off in rest frame",
	1);
	static SwitchOption interfaceEnergyCutOff2
	(interfaceEnergyCutOff,
	"LabFrame",
	"Apply cut-off in lab frame",
	2);

	static Switch<IFDipole,unsigned int> interfaceBetaOption
	("BetaOption",
	"Option for the inclusive of the higher beta coefficients",
	&IFDipole::_betaopt, 1, false, false);
	static SwitchOption interfaceBetaOptionNone
	(interfaceBetaOption,
	"None",
	"No higher betas included",
	0);
	static SwitchOption interfaceBetaOptionCollinear
	(interfaceBetaOption,
	"Collinear",
	"Include the collinear approx",
	1);
	static SwitchOption interfaceBetaOptionCollinearVirtA
	(interfaceBetaOption,
	"CollinearVirtualA",
	"Include the collinear approx with virtual corrections",
	2);
	static SwitchOption interfaceBetaOptionCollinearVirtB
	(interfaceBetaOption,
	"CollinearVirtualB",
	"Include the collinear approx with virtual corrections",
	3);
	static SwitchOption interfaceBetaOptionExact
	(interfaceBetaOption,
	"Exact",
	"Include the exact higher order terms if available",
	4);

	}

	ParticleVector IFDipole::generatePhotons(const Particle & p,ParticleVector children) {
	// set parameters which won't change in the event loop
	// masses of the particles
	_m[0] = p.mass();
	_m[1] = children[0]->mass();
	_m[2] = children[1]->mass();
	// momenta before radiation in lab
	for(unsigned int ix=0;ix<2;++ix){_qlab[ix]=children[ix]->momentum();}
	// get the charges of the particles in units of the positron charge
	// chrg1 is the charge of the parent and chrg2 is the charge of the
	// charged child. Also we create a map between the arguments of
	// _q???[X] _m[X] etc so that
	// _q???[_map[0]] and _m[_map[0]] are the momenta and masses of
	// the charged child while
	// _q???[_map[1]] and _m[_map[1]] are the momenta and masses of
	// the neutral child.
	_chrg1 = p.dataPtr()->iCharge()/3.0;
	if(children[1]->dataPtr()->iCharge()/3.0==0.0) {
	_chrg2 = children[0]->dataPtr()->iCharge()/3.0;
	_map[0] = 0; _map[1] = 1;
	}
	else if(children[0]->dataPtr()->iCharge()/3.0==0.0) {
	_chrg2 = children[1]->dataPtr()->iCharge()/3.0;
	_map[0] = 1; _map[1] = 0;
	}
	// check the radiating particle is not massless
	// if(children[1]->mass()<
	if(children[_map[0]]->mass()<1e-4*GeV) {
	ostringstream message;
	message << "IFDipole::generatePhotons() trying to generate QED radiation from "
	<< children[_map[0]]->dataPtr()->PDGName() << "\n with mass " << children[_map[0]]->mass()/GeV
	<< "which is much smaller than the mass of the electron.\n"
	<< "This is probably due to reading events from a LHEF,\nskipping radiation in this case.\n";
	generator()->logWarning( Exception(message.str(), Exception::warning));
	return children;
	}
	- cerr << "testing mass "
	- << children[_map[0]]->mass()/GeV << " "
	- << children[_map[0]]->dataPtr()->mass()/GeV << "\n";
	// boost the momenta to the rest frame
	Boost boostv(p.momentum().boostVector());
	// boost the particles to the parent rest frame
	// and set the initial momenta of the charged particles
	// in the dipole rest frame: currently this is the same
	// as the boson rest frame...
	for(unsigned int ix=0;ix<2;++ix) {
	// KMH - 08/11/05 - This used to be boostv instead of -boostv
	// -boostv is the boost from the lab to the parent rest frame
	// whereas boostv goes the other way!!!
	children[ix]->deepBoost(-boostv);
	_qprf[ix]=children[ix]->momentum();
	}
	// perform the unweighting
	double wgt;
	unsigned int ntry(0);
	do {
	wgt =makePhotons(boostv,children);
	++ntry;
	// Record warnings about large and weird weights in the .log file.
	if(wgt>_maxwgt\|\|wgt<0.0\|\|isnan(wgt)) {
	generator()->log() << "IFDipole.cc:\n";
	if(wgt>_maxwgt) {
	generator()->log() << "Weight exceeds maximum for decay!\n";
	}
	if(wgt<0.0) {
	generator()->log() << "Weight is negative! \n";
	}
	if(isnan(wgt)) {
	generator()->log() << "Weight is NAN! \n";
	wgt = 0.;
	}
	generator()->log() << p.PDGName() << " "
	<< children[0]->PDGName() << " "
	<< children[1]->PDGName()
	<< endl
	<< " Current Maximum = " << _maxwgt
	<< endl
	<< " Current Weight = " << wgt
	<< endl;
	generator()->log() << "Photon Multiplicity : "
	<< _multiplicity << endl
	<< "Original Parent rest frame momenta: " << endl
	<< "charged child: " << ounit(_qprf[_map[0]],GeV) << endl
	<< "neutral child: " << ounit(_qprf[_map[1]],GeV) << endl
	<< "Parent rest frame momenta: " << endl
	<< "charged child: " << ounit(_qnewprf[_map[0]],GeV)<< endl
	<< "neutral child: " << ounit(_qnewprf[_map[1]],GeV)<< endl
	<< "photons : " << ounit(_bigLprf,GeV) << endl
	<< "Weights : " << endl
	<< "_dipolewgt : " << _dipolewgt << endl
	<< "_yfswgt : " << _yfswgt << endl
	<< "_jacobianwgt : " << _jacobianwgt << endl
	<< "_mewgt : " << _mewgt << endl;
	for(unsigned int ct=0;ct<_multiplicity;ct++) {
	generator()->log() << "_cosphot[" << ct << "]: " << _cosphot[ct] << endl;
	generator()->log() << "_sinphot[" << ct << "]: " << _sinphot[ct] << endl;
	}
	if(wgt>_maxwgt) {
	if(wgt<15.0) {
	generator()->log() << "Resetting maximum weight"
	<< endl << " New Maximum = " << wgt << endl;
	_maxwgt=wgt;
	} else {
	generator()->log() << "Maximum weight set to limit (15)" << endl;
	_maxwgt=15.0;
	}
	}
	}
	} while (wgt<(_maxwgt*UseRandom::rnd()) && ntry<_maxtry);
	if(ntry>=_maxtry) {
	generator()->log() << "IFDipole Failed to generate QED radiation for the decay "
	<< p.PDGName() << " -> "
	<< children[0]->PDGName() << " "
	<< children[1]->PDGName() << endl;
	return children;
	}
	// produce products after radiation if needed
	if(_multiplicity>0) {
	// change the momenta of the children, they are currently
	// in parent rest frame
	for(unsigned int ix=0;ix<2;++ix) {
	LorentzRotation boost(solveBoost(_qnewprf[ix],children[ix]->momentum()));
	children[ix]->deepTransform(boost);
	// boost back to the lab
	// KMH - 08/11/05 - This used to be -boostv instead of boostv
	// -boostv is the boost from the lab to the parent rest frame
	// whereas boostv goes the other way!!!
	children[ix]->deepBoost(boostv);
	}
	// add the photons to the event record
	tcPDPtr photon=getParticleData(ParticleID::gamma);
	for(unsigned int ix=0;ix<_multiplicity;++ix) {
	PPtr newphoton=new_ptr(Particle(photon));
	newphoton->set5Momentum(_llab[ix]);
	children.push_back(newphoton);
	}
	return children;
	}
	// otherwise just return the orginial particles
	// boosted back to lab
	else {
	for(unsigned int ix=0;ix<children.size();++ix)
	children[ix]->deepBoost(boostv);
	return children;
	}
	}

	// member which generates the photons
	double IFDipole::makePhotons(Boost boostv,ParticleVector children) {
	// set the initial parameters
	// number of photons (zero)
	_multiplicity=0;
	// zero size of photon vectors
	_lprf.clear();
	_llab.clear();
	// zero size of angle storage
	_sinphot.clear();
	_cosphot.clear();
	// zero total momenta of the photons
	_bigLprf=Lorentz5Momentum();
	// set the initial values of the reweighting factors to one
	_dipolewgt = 1.0;
	_yfswgt = 1.0;
	_jacobianwgt = 1.0;
	_mewgt = 1.0;
	// set the maximum photon energy (exact - no approximations here).
	double boost_factor = 1.0;
	_emax=(0.5(_m[0]-sqr(_m[1]+_m[2])/_m[0]))boost_factor;
	// calculate the velocities of the children (crude/overvalued)
	double beta1(sqrt( (_qprf[_map[0]].e()+_m[_map[0]+1])
	*(_qprf[_map[0]].e()-_m[_map[0]+1])
	)
	/_qprf[_map[0]].e());
	double beta2(sqrt( (_qprf[_map[1]].e()+_m[_map[1]+1])
	*(_qprf[_map[1]].e()-_m[_map[1]+1])
	)
	/_qprf[_map[1]].e());
	// calculate 1-beta to avoid numerical problems
	double ombeta1(sqr(_m[_map[0]+1]/_qprf[_map[0]].e())/(1.+beta1));
	double ombeta2(sqr(_m[_map[1]+1]/_qprf[_map[1]].e())/(1.+beta2));
	// calculate the average photon multiplicity
	double aver(nbar(beta1,ombeta1));
	// calculate the number of photons using the poisson
	_multiplicity = UseRandom::rndPoisson(aver);
	// calculate the first part of the YFS factor
	_yfswgt/=crudeYFSFormFactor(beta1,ombeta1);
	// generate the photon momenta with respect to q1
	// keeping track of the weight
	double dipoles(1.);
	for(unsigned int ix=0;ix<_multiplicity;++ix)
	{ dipoles *= photon(beta1,ombeta1); }
	// calculate contributions to the dipole weights so far
	_dipolewgt /=dipoles;

	// now do the momentum reshuffling
	Lorentz5Momentum pmom(ZERO,ZERO,ZERO,_m[0],_m[0]);
	if(_multiplicity>0) {
	// total energy and momentum of photons
	Energy L0(_bigLprf.e()),modL(_bigLprf.rho());
	// squared invariant mass of final state fermions...
	Energy2 m122 = sqr(_m[0]-L0)-sqr(modL);
	if(m122<sqr(_m[1]+_m[2])) return 0.;
	// 3-momenta of charged particles
	Energy modq(_qprf[_map[0]].rho());
	// total photon momentum perpendicular to charged child...
	Energy LT(_bigLprf.perp());
	// kallen function...
	Energy4 kallen = ( m122 - sqr(_m[1]+_m[2]) )
	* ( m122 - sqr(_m[1]-_m[2]) );
	// discriminant of rho...
	Energy4 droot = kallen-4.sqr(_m[_map[0]+1]LT);
	if(droot<ZERO) return 0.;
	double disc = (_m[0]-L0) * sqrt(droot) / (2.modq(m122+LT*LT));
	// calculate the energy rescaling factor
	double rho = disc-_bigLprf.z()
	* (m122+sqr(_m[_map[0]+1])-sqr(_m[_map[1]+1]))
	/ (2.modq(m122+LT*LT));
	// calculate the rescaled charged child momentum
	_qnewprf[_map[0]]=rho*_qprf[_map[0]];
	_qnewprf[_map[0]].setMass(_m[_map[0]+1]);
	_qnewprf[_map[0]].rescaleEnergy();
	// rotate the photons so in parent rest frame rather
	// than angle measured w.r.t q1 first work out the rotation
	SpinOneLorentzRotation rotation;
	rotation.setRotateZ(-_qprf[_map[0]].phi());
	rotation.rotateY(_qprf[_map[0]].theta());
	rotation.rotateZ(_qprf[_map[0]].phi());
	// rotate the total
	_bigLprf*=rotation;
	// rotate the photons
	for(unsigned int ix=0;ix<_multiplicity;++ix){_lprf[ix]*=rotation;}
	// calculate the rescaled neutral child momentum
	_qnewprf[_map[1]]=pmom-_qnewprf[_map[0]]-_bigLprf;
	_qnewprf[_map[1]].setMass(_m[_map[1]+1]);
	_qnewprf[_map[1]].rescaleEnergy();
	// calculate the new dipole weight
	// Note this (weight) is Lorentz invariant
	// calculate velocities and 1-velocites
	beta1=sqrt( (_qnewprf[_map[0]].e()+_m[_map[0]+1])
	*(_qnewprf[_map[0]].e()-_m[_map[0]+1]))
	/_qnewprf[_map[0]].e();
	beta2=sqrt( (_qnewprf[_map[1]].e()+_m[_map[1]+1])
	*(_qnewprf[_map[1]].e()-_m[_map[1]+1]))
	/_qnewprf[_map[1]].e();
	ombeta1=sqr(_m[_map[0]+1]/_qnewprf[_map[0]].e())/(1.+beta1);
	ombeta2=sqr(_m[_map[1]+1]/_qnewprf[_map[1]].e())/(1.+beta2);
	for(unsigned int ix=0;ix<_multiplicity;++ix)
	{_dipolewgt*=exactDipoleWeight(beta1,ombeta1,ix);}
	// calculate the second part of the yfs form factor
	_yfswgt*=exactYFSFormFactor(beta1,ombeta1,beta2,ombeta2);
	// Now boost from the parent rest frame to the lab frame
	SpinOneLorentzRotation boost(boostv);
	// Boosting charged particles
	for(unsigned int ix=0;ix<2;++ix){_qnewlab[ix]=boost*_qnewprf[ix];}
	// Boosting total photon momentum
	_bigLlab=boost*_bigLprf;
	// Boosting individual photon momenta
	for(unsigned int ix=0;ix<_multiplicity;++ix)
	{_llab.push_back(boost*_lprf[ix]);}
	// Calculating jacobian weight
	_jacobianwgt = jacobianWeight();
	// Calculating beta^1 weight
	_mewgt = meWeight(children);
	// Apply phase space vetos...
	if(kallen<(4.sqr(_m[_map[0]+1]LT))\|\|m122<sqr(_m[1]+_m[2])\|\|rho<0.0) {
	// generator()->log() << "Outside Phase Space" << endl;
	// generator()->log() << "Photon Multiplicity: "
	// << _multiplicity << endl
	// << "Original Parent rest frame momenta: " << endl
	// << "charged child: " << _qprf[_map[0]] << endl
	// << "neutral child: " << _qprf[_map[1]] << endl
	// << "rescaling : " << rho << endl
	// << "Parent rest frame momenta: " << endl
	// << "charged child: " << _qnewprf[_map[0]] << endl
	// << "neutral child: " << _qnewprf[_map[1]] << endl
	// << "photons : " << _bigLprf << endl
	// << endl;
	_dipolewgt = 0.0 ;
	_yfswgt = 0.0 ;
	_jacobianwgt = 0.0 ;
	_mewgt = 0.0 ;
	}
	_qprf[_map[0]].rescaleEnergy();
	_qprf[_map[1]].rescaleEnergy();
	_qnewprf[_map[0]].rescaleEnergy();
	_qnewprf[_map[1]].rescaleEnergy();
	if( ((abs(_m[0]-_bigLprf.e()-_qnewprf[0].e()-_qnewprf[1].e())>0.00001*MeV)\|\|
	(abs( _bigLprf.x()+_qnewprf[0].x()+_qnewprf[1].x())>0.00001*MeV)\|\|
	(abs( _bigLprf.y()+_qnewprf[0].y()+_qnewprf[1].y())>0.00001*MeV)\|\|
	(abs( _bigLprf.z()+_qnewprf[0].z()+_qnewprf[1].z())>0.00001*MeV))
	&&(_dipolewgt_jacobianwgt_yfswgt*_mewgt>0.0)) {
	Lorentz5Momentum ptotal = _bigLprf+_qnewprf[0]+_qnewprf[1];
	ptotal.setE(ptotal.e()-_m[0]);
	generator()->log()
	<< "Warning! Energy Not Conserved! tol = 0.00001 MeV"
	<< "\nwgt = " << _dipolewgt_yfswgt_jacobianwgt*_mewgt
	<< "\nrho = " << rho
	<< "\nmultiplicity = " << _multiplicity
	<< "\n_qprf[_map[0]] = " << _qprf[_map[0]]/GeV
	<< "\n_qprf[_map[1]] = " << _qprf[_map[1]]/GeV
	<< "\n_qnewprf[_map[0]] = " << _qnewprf[_map[0]]/GeV << " "
	<< _qnewprf[_map[0]].m()/GeV << " " << _m[_map[0]+1]/GeV
	<< "\n_qnewprf[_map[1]] = " << _qnewprf[_map[1]]/GeV << " "
	<< _qnewprf[_map[1]].m()/GeV << " " << _m[_map[1]+1]/GeV
	<< "\n_bigLprf = " << _bigLprf/GeV
	<< "\n_bigLprf.m2() = " << _bigLprf.m2()/GeV2
	<< "\n_total out -in = " << ptotal/GeV
	<< "\nRejecting Event. " << "\n";
	_dipolewgt = 0.0 ;
	_yfswgt = 0.0 ;
	_jacobianwgt = 0.0 ;
	_mewgt = 0.0 ;
	}
	}
	// otherwise copy momenta
	else
	{ for(unsigned int ix=0;ix<2;++ix) {
	_qnewprf[ix]=_qprf[ix];
	_qnewlab[ix]=_qlab[ix];
	}
	_jacobianwgt = 1.0;
	// calculate the second part of the yfs form factor
	_yfswgt*=exactYFSFormFactor(beta1,ombeta1,beta2,ombeta2);
	_dipolewgt = 1.0;
	}
	// Virtual corrections for beta_0:
	// These should be zero for the scalar case as there is no
	// collinear singularity going by the dipoles above...
	// Use mass of decaying particle...
	if(_betaopt==2) {
	if((children[_map[0]]->dataPtr()->iSpin())==2) {
	_mewgt += (0.5*_alpha/pi)
	* log(sqr(_m[0]
	/_m[_map[0]+1])
	);
	}
	}
	// OR Use invariant mass of final state children...
	if(_betaopt==3) {
	if((children[_map[0]]->dataPtr()->iSpin())==2) {
	_mewgt += (0.5*_alpha/pi)
	* log((_qnewprf[0]+_qnewprf[1]).m2()
	/sqr(_m[_map[0]+1])
	);
	}
	}
	// calculate the weight depending on the option
	double wgt;
	if(_mode==0){wgt=_maxwgt;}
	else if(_mode==1){wgt=_mewgt_jacobianwgt_yfswgt*_dipolewgt;}
	else if(_mode==2){wgt=_jacobianwgt_yfswgt_dipolewgt;}
	else if(_mode==3){wgt=_yfswgt*_dipolewgt;}
	else {wgt=_yfswgt;}
	return wgt;
	}

	double IFDipole::photon(double beta1,double ombeta1)
	{
	// generate the azimuthal angle randomly in -pi->+pi
	double phi(-pi+UseRandom::rnd()2.pi);
	// generate the polar angle
	double r(UseRandom::rnd());
	double costh,sinth,ombc;
	ombc = pow(1.+beta1,1.-r)*pow(ombeta1,r);
	costh = 1./beta1*(1.-ombc);
	sinth = sqrt(ombc(2.-ombc)-(1.+beta1)ombeta1*sqr(costh));
	// generate the ln(energy) uniformly in ln(_emin)->ln(_emax)
	Energy energy = pow(_emax/_emin,UseRandom::rnd())*_emin;
	// calculate the weight (omit the pre and energy factors
	// which would cancel later anyway)
	double wgt = 2./ombc;
	// store the angles
	_cosphot.push_back(costh);
	_sinphot.push_back(sinth);
	// store the four vector for the photon
	_lprf.push_back(Lorentz5Momentum(energysinthcos(phi),energysinthsin(phi),
	energy*costh,energy,ZERO));
	// add the photon momentum to the total
	_bigLprf+=_lprf.back();
	// return the weight
	return wgt;
	}

	double IFDipole::meWeight(ParticleVector children)
	{
	unsigned int spin = children[_map[0]]->dataPtr()->iSpin();
	double mewgt = 1.0;
	double beta1=sqrt( (_qnewprf[_map[0]].e()+_m[_map[0]+1])
	*(_qnewprf[_map[0]].e()-_m[_map[0]+1]))
	/_qnewprf[_map[0]].e();
	double ombeta1=sqr(_m[_map[0]+1]/_qnewprf[_map[0]].e())/(1.+beta1);
	// option which does nothing
	if(_betaopt==0){mewgt=1.;}
	// collinear approx
	else if(_betaopt==1\|\|_betaopt==2\|\|_betaopt==3)
	{
	double ombc;
	InvEnergy2 dipole;
	for(unsigned int i=0;i<_multiplicity;++i) {
	double opbc;
	if(_cosphot[i]<0.0)
	{ opbc=ombeta1+beta1*sqr(_sinphot[i])/(1.-_cosphot[i]); }
	// if cos is greater than zero use result accurate as cos->-1
	else
	{ opbc=1.+beta1*_cosphot[i]; }
	// if cos is greater than zero use result accurate as cos->1
	if(_cosphot[i]>0.0)
	{ ombc=ombeta1+beta1*sqr(_sinphot[i])/(1.+_cosphot[i]); }
	// if cos is less than zero use result accurate as cos->-1
	else
	{ ombc=1.-beta1*_cosphot[i]; }
	if(((_qnewprf[_map[0]].z()>ZERO)&&(_qprf[_map[0]].z()<ZERO))\|\|
	((_qnewprf[_map[0]].z()<ZERO)&&(_qprf[_map[0]].z()>ZERO))) {
	dipole = sqr(beta1_sinphot[i]/(opbc_lprf[i].e()));
	} else {
	dipole = sqr(beta1_sinphot[i]/(ombc_lprf[i].e()));
	}
	// here "dipole" is the exact dipole function divided by alpha/4pi^2.

	if(spin==2) {
	Energy magpi= sqrt( sqr(_qnewprf[_map[0]].x())
	+ sqr(_qnewprf[_map[0]].y())
	+ sqr(_qnewprf[_map[0]].z())
	);

	mewgt += sqr(_lprf[i].e())_qnewprf[_map[0]].e()ombc
	/ (sqr(magpi_sinphot[i])(_qnewprf[_map[0]].e()+_lprf[i].e()));
	}
	else if(spin==3) {
	Energy2 pik = _qnewprf[_map[0]].e()*_lprf[i].e()
	- _qnewprf[_map[0]].x()*_lprf[i].x()
	- _qnewprf[_map[0]].y()*_lprf[i].y()
	- _qnewprf[_map[0]].z()*_lprf[i].z();

	Energy2 pjk = _m[0]*_lprf[i].e();

	Energy2 pipj = _m[0]*_qnewprf[_map[0]].e();

	mewgt += (2.pjkpipj/(pik*sqr(pipj+pjk))
	+2.pjk/(pik(pipj+pik))
	)/dipole;
	}
	else {
	mewgt = 1.0;
	}
	}
	}
	return mewgt;
	}

	double IFDipole::exactYFSFormFactor(double beta1,double ombeta1,
	double beta2,double ombeta2) {
	double Y = 0.0 ;
	double b = beta1 ;
	double omb = ombeta1;
	double c = beta2 ;
	double omc = ombeta2;
	double arg1 = -omc/(2.*c);
	double arg2 = -ombomc/(2.(b+c));
	double arg3 = 2.*b/(1.+b);
	if(_m[_map[1]+1]!=ZERO) {
	Y = _chrg1_chrg2(_alpha/(2.pi))(
	log(_m[0]_m[_map[1]+1]/sqr(2._emin))
	+log(_m[_map[0]+1]_m[_map[1]+1]/sqr(2._emin))
	-(1./b )log((1.+b)/omb)log(sqr(_m[_map[1]+1]/(2.*_emin)))
	-(1./b )*log(omb/(1.+b))
	-(0.5/b )*sqr(log(omb/(1.+b)))
	+((b+c )/(bomc))log((b+c )/(b*omc))
	-((c+bc)/(bomc))log((c+bc)/(b*omc))
	+((b+c )/(b+bc))log((b+c )/(b+b*c))
	-((comb)/(b+bc))log((comb)/(b+b*c))
	+(0.5/b)( sqr(log( (b+c)/(bomc)))-sqr(log((c+bc)/(bomc)))
	+ sqr(log((comb)/(b+bc)))-sqr(log((b+ c)/(b+b*c)))
	)
	+(2./b )*( real(Math::Li2(arg1))
	- real(Math::Li2(arg2))
	- real(Math::Li2(arg3))
	)
	+(1./b )log((b+c)/(b+bc))log((1.+c)/(2.c))
	-(1./b )log((comb)/(b(1.+c)))log((1.+b)(1.+c)/(2.(b+c)))
	-(1./b )log((2.c/b)((b+c)/(omc(1.+c))))log((b+c)/(comb))
	);
	}
	else if(_m[_map[1]+1]==ZERO) {
	Y = _chrg1_chrg2(_alpha/(2.pi))(
	log(sqr(_m[0]/(2.*_emin)))
	+log(sqr(_m[_map[0]+1]/(2.*_emin)))
	-(1./b )*log((1.+b)/omb)
	log((sqr(_m[0])-sqr(_m[_map[0]+1]))/sqr(2._emin))
	-0.5log(omb(1.+b)/sqr(2.*b))
	+((1.+b)/(2.b))log((1.+b)/(2.*b))
	-( omb/(2.b))log( omb/(2.*b))
	-(1./b )*log((1.-b)/(1.+b))
	+1.
	+(0.5/b)sqr(log( omb/(2.b)))
	-(0.5/b)sqr(log((1.+b)/(2.b)))
	-(0.5/b)*sqr(log((1.-b)/(1.+b)))
	-(2. /b)*real(Math::Li2(arg3))
	);
	}
	return exp(Y);
	}

	double IFDipole::jacobianWeight() {
	// calculate the velocities of the children (crude/overvalued)
	Energy mag1old = sqrt( (_qprf[_map[0]].e() +_m[_map[0]+1])
	*(_qprf[_map[0]].e() -_m[_map[0]+1])
	);
	Energy mag1new = sqrt( (_qnewprf[_map[0]].e()+_m[_map[0]+1])
	*(_qnewprf[_map[0]].e()-_m[_map[0]+1])
	);
	Energy magL = sqrt( sqr(_bigLprf.x())
	+ sqr(_bigLprf.y())
	+ sqr(_bigLprf.z())
	);

	// 14/12/05 - KMH - This was another mistake. This is supposed to be
	// the angel between _qnewprf[_map[0]] and _bigLprf instead of
	// between _qnewprf[0] and _bigLprf. Stupid. Hopefully this weight
	// is correct now.
	// double cos1L = (_qnewprf[0].x()*_bigLprf.x()
	// +_qnewprf[0].y()*_bigLprf.y()
	// +_qnewprf[0].z()*_bigLprf.z()
	// )
	// /(mag1new*magL);
	double cos1L = (_qnewprf[_map[0]].x()*_bigLprf.x()
	+_qnewprf[_map[0]].y()*_bigLprf.y()
	+_qnewprf[_map[0]].z()*_bigLprf.z()
	)
	/(mag1new*magL);

	return abs( (_m[0]*sqr(mag1new)/mag1old)
	/ ( mag1new*(_m[0]-_bigLprf.e())
	+_qnewprf[_map[0]].e()magLcos1L
	)
	);
	}

	LorentzRotation IFDipole::solveBoost(const Lorentz5Momentum & q,
	const Lorentz5Momentum & p ) const {
	Energy modp = p.vect().mag();
	Energy modq = q.vect().mag();
	double betam = (p.e()modp-q.e()modq)/(sqr(modq)+sqr(modp)+p.mass2());
	Boost beta = -betam*q.vect().unit();
	ThreeVector<Energy2> ax = p.vect().cross( q.vect() );
	double delta = p.vect().angle( q.vect() );
	LorentzRotation R;
	using Constants::pi;
	if ( ax.mag2()/GeV2/MeV2 > 1e-16 ) {
	R.rotate( delta, unitVector(ax) ).boost( beta );
	}
	else {
	if(p.mass()>ZERO) {
	R.boost(p.findBoostToCM(),p.e()/p.mass());
	R.boost(q.boostVector(),q.e()/q.mass());
	}
	else {
	if(modp>modq) beta = -betam*p.vect().unit();
	R.boost( beta );
	}
	}
	return R;
	}

	void IFDipole::doinit() {
	Interfaced::doinit();
	// get the value fo alpha from the Standard Model object
	_alpha=generator()->standardModel()->alphaEM();
	}
	diff --git a/Hadronization/Cluster.cc b/Hadronization/Cluster.cc
	--- a/Hadronization/Cluster.cc
	+++ b/Hadronization/Cluster.cc
	@@ -1,258 +1,254 @@
	// -- C++ --
	//
	// Cluster.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the Cluster class.
	//

	#include "Cluster.h"
	#include <ThePEG/Repository/UseRandom.h>
	#include <ThePEG/Repository/CurrentGenerator.h>
	#include <ThePEG/PDT/ParticleData.h>
	#include "ClusterHadronizationHandler.h"

	using namespace Herwig;

	PPtr Cluster::clone() const {
	return new_ptr(*this);
	}

	PPtr Cluster::fullclone() const {
	return clone();
	}

	-cCluHadHdlPtr Cluster::_clusterHadHandler = cCluHadHdlPtr();
	-
	-Energy2 Cluster::_mg2 = ZERO;
	-
	ClassDescription<Cluster> Cluster::initCluster;

	Cluster::Cluster()
	: Particle(),
	_isAvailable(true),
	_hasReshuffled(false),
	_component(),
	_original(),
	_isBeamRemnant(),
	_isPerturbative(),
	_numComp(0),
	_id(0) {}

	void Cluster::persistentOutput(PersistentOStream & os) const {
	os << _isAvailable << _hasReshuffled << _component << _original
	- << _isBeamRemnant << _isPerturbative << _numComp << _id
	- << _clusterHadHandler << ounit(_mg2,GeV2);
	+ << _isBeamRemnant << _isPerturbative << _numComp << _id;
	}

	void Cluster::persistentInput(PersistentIStream & is, int) {
	is >> _isAvailable >> _hasReshuffled >> _component >> _original
	- >> _isBeamRemnant >> _isPerturbative >> _numComp >> _id
	- >> _clusterHadHandler >> iunit(_mg2,GeV2);
	+ >> _isBeamRemnant >> _isPerturbative >> _numComp >> _id;
	}



	Cluster::Cluster(tPPtr p1, tPPtr p2, tPPtr p3)
	: Particle(CurrentGenerator::current().
	getParticleData(long(ParticleID::Cluster))),
	_isAvailable(true), _hasReshuffled(false)
	{
	if(!dataPtr()) {
	cerr << "Cluster Particle Data not defined. Cannot complete Hadronization "
	<< "without ParticleData for id " << ParticleID::Cluster << '\n';
	}
	_component.push_back(new_ptr(Particle(*p1)));
	_component.push_back(new_ptr(Particle(*p2)));
	if(p3) _component.push_back(new_ptr(Particle(*p3)));
	_original.push_back(p1);
	_original.push_back(p2);
	if(p3) _original.push_back(p3);

	_isPerturbative.push_back(initPerturbative(p1));
	_isPerturbative.push_back(initPerturbative(p2));
	if(p3) _isPerturbative.push_back(initPerturbative(p3));
	else _isPerturbative.push_back(false);
	for(int i = 0; i<3; i++) _isBeamRemnant.push_back(false);

	if(p3) {
	_numComp = 3;
	_id = 100abs(p1->id()) + 10abs(p2->id()) + abs(p3->id());
	} else {
	_numComp = 2;
	if(p2->id() > 10)
	_id = 10*abs(p2->id()/100) + abs(p1->id());
	else if(p1->id() > 10)
	_id = 10*abs(p1->id()/100) + abs(p2->id());
	else
	_id = 10*abs(p1->id()) + abs(p2->id());
	}
	calculateP();
	calculateX();

	}

	Cluster::Cluster(tcEventPDPtr x)
	: Particle(x),
	_isAvailable(false),
	_hasReshuffled(false),
	_component(),
	_original(),
	_isBeamRemnant(),
	_isPerturbative(),
	_numComp(0),
	_id(0) {}

	Cluster::Cluster(const Particle &x)
	: Particle(x),
	_isAvailable(false),
	_hasReshuffled(false),
	_component(),
	_original(),
	_isBeamRemnant(),
	_isPerturbative(),
	_numComp(0),
	_id(0) {}

	Energy Cluster::sumConstituentMasses() const
	{
	if(_numComp == 3) {
	return _component[0]->mass() +
	_component[1]->mass() +
	_component[2]->mass();
	} else if(_numComp == 2)
	return _component[0]->mass() + _component[1]->mass();
	else return ZERO;
	}


	void Cluster::calculateP() {
	Lorentz5Momentum m;
	for(int i = 0; i<_numComp; i++)
	m += _component[i]->momentum();
	m.rescaleMass();
	set5Momentum(m);
	}


	void Cluster::calculateX() {
	if ( _numComp != 2 ) {
	// Only in the case of two components we have a definition of cluster
	// position in terms of the two components.
	setVertex(LorentzPoint());
	} else {
	// Get the needed parameters.
	- assert(_clusterHadHandler);
	- Energy2 vmin2 = _clusterHadHandler->minVirtuality2();
	- Length dmax = _clusterHadHandler->maxDisplacement();
	+ Energy2 vmin2
	+ = ClusterHadronizationHandler::currentHandler()->minVirtuality2();
	+ Length dmax
	+ = ClusterHadronizationHandler::currentHandler()->maxDisplacement();

	// Get the positions and displacements of the two components (Lab frame).
	LorentzPoint pos1 = _component[0]->vertex();
	Lorentz5Momentum p1 = _component[0]->momentum();
	LorentzDistance displace1 = -log( UseRandom::rnd() ) *
	hbarc * p1 * (1 / sqrt(sqr(p1.m2() - p1.mass2()) + sqr(vmin2)));
	if ( abs( displace1.m() ) > dmax ) {
	displace1 *= dmax / abs( displace1.m() );
	}
	LorentzPoint pos2 = _component[1]->vertex();
	Lorentz5Momentum p2 = _component[1]->momentum();
	LorentzDistance displace2 = -log( UseRandom::rnd() ) *
	hbarc * p2 * (1 / sqrt(sqr(p2.m2() - p2.mass2()) + sqr(vmin2)));
	if ( abs( displace2.m() ) > dmax ) {
	displace2 *= dmax / abs( displace2.m() );
	}
	double s1 = 0.0, s2 = 0.0;
	Lorentz5Momentum pcl = p1 + p2;
	if ( abs( pcl.vect().dot( displace1.vect() ) ) > 1.0e-20MeVmm &&
	abs( pcl.vect().dot( displace2.vect() ) ) > 1.0e-20MeVmm ) {
	// The displacement with the smallest projection along pcl.vect()
	// is scaled up such that both displacements have equal projections
	// along pcl.vect().
	double ratio = ( abs( pcl.vect().dot( displace1.vect() ) ) /
	abs( pcl.vect().dot( displace2.vect() ) ) );
	if ( pcl.vect().dot(displace1.vect()) *
	pcl.vect().dot(displace2.vect()) < 0.0sqr(MeVmm) ) {
	ratio *= -1;
	}
	if ( abs( ratio ) > 1.0 ) {
	displace2 *= ratio;
	} else {
	displace1 *= ratio;
	}
	// Now determine the s1 and s2 values.
	double s1minusS2 = ( pcl.vect().dot( pos2.vect() - pos1.vect() ) /
	pcl.vect().dot( displace1.vect() ) );
	if ( s1minusS2 < 0 ) {
	s1 = 1.0;
	s2 = s1 - s1minusS2;
	} else if ( s1minusS2 > 0 ) {
	s2 = 1;
	s1 = s2 + s1minusS2;
	}
	}
	// Now, finally, determine the cluster position
	setVertex(0.5 * (pos1 + pos2 + s1displace1 + s2displace2));
	} // end else part of if ( _collecCompPtr.size() != 2 )
	}

	bool Cluster::isBeamCluster() const {
	for(int i = 0; i<_numComp; i++)
	if(_isBeamRemnant[i]) return true;
	return false;
	}


	void Cluster::isBeamCluster(tPPtr part) {
	for(int i = 0; i<_numComp; i++) {
	if(_original[i] == part) {
	_isBeamRemnant[i] = true;
	break;
	}
	}
	}

	bool Cluster::isStatusFinal() const {
	int s = children().size();
	for(unsigned int i = 0; i<children().size(); i++)
	if(children()[i]->PDGName() == "Cluster") s--;
	return ( s > 0);
	}

	tPPtr Cluster::particle(int i) const {
	return (i < _numComp) ? _component[i] : PPtr();
	}

	tPPtr Cluster::colParticle(bool anti) const {
	if ( _numComp != 2 ) return PPtr();
	if ( _original[0]->hasColour(anti) ) return _original[0];
	else if ( _original[1]->hasColour(anti) ) return _original[1];
	else return PPtr();
	}

	tPPtr Cluster::antiColParticle() const {
	return colParticle(true);
	}

	bool Cluster::isPerturbative(int i) const {
	return _isPerturbative[i];
	}

	bool Cluster::isBeamRemnant(int i) const {
	return _isBeamRemnant[i];
	}

	void Cluster::setBeamRemnant(int i, bool b) {
	if(i < _numComp)
	_isBeamRemnant[i] = b;
	}

	-void Cluster::setPointerClusterHadHandler(tcCluHadHdlPtr gp) {
	- _clusterHadHandler = gp;
	- _mg2=sqr(_clusterHadHandler->
	- getParticleData(ParticleID::g)->constituentMass());
	+bool Cluster::initPerturbative(tPPtr p)
	+{
	+ Energy mg
	+ = CurrentGenerator::current().getParticleData(ParticleID::g)->constituentMass();
	+ return p->scale() > sqr(mg);
	}
	diff --git a/Hadronization/Cluster.h b/Hadronization/Cluster.h
	--- a/Hadronization/Cluster.h
	+++ b/Hadronization/Cluster.h
	@@ -1,356 +1,338 @@
	// -- C++ --
	//
	// Cluster.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_Cluster_H
	#define HERWIG_Cluster_H

	#include <ThePEG/EventRecord/Particle.h>
	#include "Herwig++/Utilities/EnumParticles.h"
	#include "CluHadConfig.h"
	#include "ClusterHadronizationHandler.fh"
	#include "Cluster.fh"
	#include "ThePEG/Utilities/ClassDescription.h"

	namespace Herwig {
	using namespace ThePEG;

	/** \ingroup Hadronization
	* \class Cluster
	* \brief This class describes a cluster object.
	* \author Philip Stephens
	* \author Alberto Ribon
	*
	* This class represents a cluster, which is a colour singlet made usually
	* of two components (quark-antiquark, quark-diquark, antiquark-antidiquark)
	* or rarely by three components (quark-quark-quark, antiquark-antiquark-
	* antiquark). A reference to the container with the pointers to its
	* Components is provided.
	*
	* The class provides access to the pointers which point to:
	*
	* - The cluster parent. In the case that the cluster it is a fission
	* product of a heavy cluster the parent is a cluster. If the cluster
	* is formed from the perturbative partons then the parents will be
	* the colour connected partons that formed the cluster.
	* - The children (usually two). In the case the cluster is a
	* heavy cluster that undergoes fission the children are clusters.
	* Occasionally the cluster has been "redefined" (re-interpreted). For
	* example in the case that three quark or anti-quark components
	* have been redefined as two components (quark+diquark, or antiquark+
	* antidiquark).
	* - The (eventual) reshuffling partner, necessary for energy-momentum
	* conservation when light clusters are decayed into single hadron. Not
	* all clusters will have a reshuffling partner.
	*
	* Notice that in order to determine the cluster position from the positions
	* of the components, the Cluster class needs some parameters.
	* Because the Cluster class is neither interfaced nor persistent,
	* a static pointer to the ClusterHadronizationHandler class instance,
	* where the parameters are, is used. This static pointer is
	* set via the method setPointerClusterHadHandler(), during the
	* run initialization, doinitrun() of ClusterHadronizationHandler.
	*
	* @see ClusterHadronizationHandler
	*/
	class Cluster : public Particle {

	protected:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Default constructor. Only used in PersistentIStream
	*/
	Cluster();

	/**
	* The ClassTraits<Cluster> class must be a friend to be able to
	* use the private default constructor.
	*/
	friend struct ClassTraits<Cluster>;

	public:
	/**
	* Constructor with a particleData pointer
	*/
	Cluster(tcEventPDPtr);

	/**
	* This creates a cluster from 2 (or 3) partons.
	*/
	Cluster(tPPtr part1, tPPtr part2, tPPtr part3 = tPPtr());

	/**
	* Also a constructor where a particle is given not a cluster.
	*/
	Cluster(const Particle &);
	//@}
	-
	- /**
	- * Set the static pointer to the ClusterHadronizationHandler object.
	- * The pointer is set in ClusterHadronizationHandler::doinitrun().
	- */
	- static void setPointerClusterHadHandler(tcCluHadHdlPtr gp);

	/**
	* Number of quark (diquark) constituents (normally two).
	*/
	int numComponents() const
	{ return _numComp; }

	/**
	* Sum of the constituent masses of the components of the cluster.
	*/
	Energy sumConstituentMasses() const;

	/**
	* Returns the ith constituent.
	*/
	tPPtr particle(int i) const;

	/**
	* Returns the original constituent carrying colour
	*/
	tPPtr colParticle(bool anti = false) const;

	/**
	* Returns the original constituent carrying anticolour
	*/
	tPPtr antiColParticle() const;

	/**
	* Returns whether the ith constituent is from a perturbative process.
	*/
	bool isPerturbative(int) const;

	/**
	* Indicates whether the ith constituent is a beam remnant.
	*/
	bool isBeamRemnant(int) const;

	/**
	* Sets whether the ith constituent is a beam remnant.
	*/
	void setBeamRemnant(int,bool);

	/**
	* Returns the clusters id, not the same as the PDG id.
	*/
	int clusterId() const { return _id; }

	public:

	/**
	* Returns true when a constituent is a beam remnant.
	*/
	bool isBeamCluster() const;

	/**
	* Set the pointer to the reshuffling partner cluster.
	*/
	void flagAsReshuffled()
	{ _hasReshuffled = true; }

	/**
	* Sets the component (if any) that points to "part" as a beam remnant.
	*/
	void isBeamCluster(tPPtr part);

	/**
	* Returns true if this cluster is to be handled by the hadronization.
	*/
	bool isAvailable() const
	{ return _isAvailable; }

	/**
	* Sets the value of availability.
	*/
	void isAvailable(bool inputAvailable)
	{ _isAvailable = inputAvailable; }

	/**
	* Return true if the cluster does not have cluster parent.
	*/
	bool isStatusInitial() const
	{ return parents().empty(); }

	/**
	* Return true if the cluster does not have cluster children and
	* it is not already decayed (i.e. it does not have hadron children)
	* (to be used only after the fission of heavy clusters).
	*/
	bool isReadyToDecay() const
	{ return children().empty(); }

	/**
	* Return true if the cluster has one and only one cluster children
	* and no hadron children: that means either that its three quarks or
	* anti-quarks components have been redefined as two components
	* (quark+diquark, or antiquark+antidiquark), or that the cluster
	* has been used as a partner for the momentum reshuffling necessary
	* to conserve energy-momentum when a light cluster is decayed into
	* a single hadron (notice that this latter light cluster has
	* isRedefined() false, because it has an hadron child).
	* In both cases, the unique cluster children is the new redefined
	* cluster. The two cases can be distinguish by the next method.
	*/
	bool isRedefined() const {
	return ( children().size() == 1
	&& children()[0]->id() == ParticleID::Cluster );
	}

	/**
	* Return true when it has a reshuffling partner.
	* Notice that a cluster can have hasBeenReshuffled() true but
	* isRedefined() false: this is the case of a light cluster
	* that decays into a single hadron.
	*/
	bool hasBeenReshuffled() const
	{ return _hasReshuffled; }

	/**
	* Return true if the cluster has hadron children.
	*/
	bool isStatusFinal() const;

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual PPtr 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 PPtr fullclone() const;
	//@}

	public:

	/** @name Input and output functions. */
	//@{
	/**
	* Standard function for writing to a persistent stream.
	*/
	void persistentOutput(PersistentOStream &) const;

	/**
	* Standard function for reading from a persistent stream.
	*/
	void persistentInput(PersistentIStream &, int);

	private:
	/**
	* Private and non-existent assignment operator.
	*/
	Cluster & operator=(const Cluster &);

	/**
	* Calculate the 5-momentum vector of the cluster
	* The 5-momentum of the cluster is given by
	* \f[ P = \sum_i p_i \f]
	* and the mass of the cluster is \f$m^2 = P^2\f$
	*/
	void calculateP();

	/**
	* Calculate the 4-position vector of the cluster
	* Displacement of the ith constituent given by momentum \f$p_i\f$
	* vertex \f$x_i\f$ and mass \f$m_i\f$ is
	* \f[ D_i = -C \log(r) \frac{p_i}{\sqrt{(p_i^2 - m_i^2)^2 + v^4}} \f]
	* where \f$r\f$ is a random number [0,1],
	* \f$v\f$ is the minimum virtuality and \f$C\f$ is
	* a conversion factor from GeV to millimeters. We can then find the
	* difference in \f$s\f$ factors as
	* \f[ (s_1-s_2) = \frac{(\vec{p}_1 + \vec{p}_2) \cdot (\vec{x}_2 -
	* \vec{x}_1)}{(\vec{p}_1 + \vec{p}_2) \cdot \vec{D}_1}.
	* \f]
	* if \f$s_2>s_1\f$ then \f$s_1 = 1.0\f$ otherwise \f$s_2 = 1.0\f$.
	* These are then used to determine the value of the clusters vertex as
	* \f[ X = \frac{1}{2} ( x_1 +x_2 + s_1 D_1 + s_2 D_2). \f]
	*/
	void calculateX();

	/**
	* Determines whether constituent p is perturbative or not.
	*/
	- bool initPerturbative(tPPtr p)
	- { return p->scale() > _mg2; }
	-
	- /**
	- * This is needed to determine if a cluster is from a perturbative quark.
	- */
	- static cCluHadHdlPtr _clusterHadHandler;
	-
	- /**
	- * The gluon mass is needed to determine if a cluster is from a perturbative quark
	- */
	- static Energy2 _mg2;
	-
	+ bool initPerturbative(tPPtr p);

	/**
	* Describe an abstract base class with persistent data.
	*/
	static ClassDescription<Cluster> initCluster;



	bool _isAvailable; //!< Whether the cluster is hadronizing
	bool _hasReshuffled; //!< Whether the cluster has been reshuffled
	ParticleVector _component; //!< The constituent partons
	- ParticleVector _original; //!< The original components
	+ tParticleVector _original; //!< The original components
	vector<bool> _isBeamRemnant; //!< Whether a parton is a beam remnant
	vector<bool> _isPerturbative; //!< Whether a parton is perturbative
	int _numComp; //!< The number of constituents
	long _id; //!< The id of this cluster
	};

	} // end namespace Herwig

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/**
	* The following template specialization informs ThePEG about the
	* base class of Cluster.
	*/
	template <>
	struct BaseClassTrait<Herwig::Cluster,1> {
	/** Typedef of the base class of Cluster. */
	typedef Particle NthBase;
	};

	/**
	* The following template specialization informs ThePEG about the
	* name of this class and the shared object where it is defined.
	*/
	template <>
	struct ClassTraits<Herwig::Cluster>:
	public ClassTraitsBase<Herwig::Cluster> {
	/** Return the class name. */
	static string className() { return "Herwig::Cluster"; }
	/** Create a Particle object. */
	static TPtr create() { return TPtr::Create(Herwig::Cluster()); }
	};

	/** @endcond */

	}




	#endif // HERWIG_Cluster_H
	diff --git a/Hadronization/ClusterHadronizationHandler.cc b/Hadronization/ClusterHadronizationHandler.cc
	--- a/Hadronization/ClusterHadronizationHandler.cc
	+++ b/Hadronization/ClusterHadronizationHandler.cc
	@@ -1,313 +1,312 @@
	// -- C++ --
	//
	// ClusterHadronizationHandler.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ClusterHadronizationHandler class.
	//

	#include "ClusterHadronizationHandler.h"
	#include <ThePEG/Interface/ClassDocumentation.h>
	#include <ThePEG/Persistency/PersistentOStream.h>
	#include <ThePEG/Persistency/PersistentIStream.h>
	#include <ThePEG/Interface/Parameter.h>
	#include <ThePEG/Interface/Reference.h>
	#include <ThePEG/Handlers/EventHandler.h>
	#include <ThePEG/Handlers/Hint.h>
	#include <ThePEG/PDT/ParticleData.h>
	#include <ThePEG/EventRecord/Particle.h>
	#include <ThePEG/EventRecord/Step.h>
	#include <ThePEG/PDT/PDT.h>
	#include <ThePEG/PDT/EnumParticles.h>
	#include <ThePEG/Utilities/Throw.h>
	#include "Herwig++/Utilities/EnumParticles.h"
	#include "CluHadConfig.h"
	#include "Cluster.h"
	#include <ThePEG/Utilities/DescribeClass.h>

	using namespace Herwig;

	+ClusterHadronizationHandler * ClusterHadronizationHandler::currentHandler_ = 0;
	+
	DescribeClass<ClusterHadronizationHandler,HadronizationHandler>
	describeClusterHadronizationHandler("Herwig::ClusterHadronizationHandler","");

	IBPtr ClusterHadronizationHandler::clone() const {
	return new_ptr(*this);
	}

	IBPtr ClusterHadronizationHandler::fullclone() const {
	return new_ptr(*this);
	}

	void ClusterHadronizationHandler::persistentOutput(PersistentOStream & os)
	const {
	os << _partonSplitter
	<< _clusterFinder
	<< _colourReconnector
	<< _clusterFissioner
	<< _lightClusterDecayer
	<< _clusterDecayer
	<< ounit(_minVirtuality2,GeV2)
	<< ounit(_maxDisplacement,mm)
	<< _underlyingEventHandler;
	}


	void ClusterHadronizationHandler::persistentInput(PersistentIStream & is, int) {
	is >> _partonSplitter
	>> _clusterFinder
	>> _colourReconnector
	>> _clusterFissioner
	>> _lightClusterDecayer
	>> _clusterDecayer
	>> iunit(_minVirtuality2,GeV2)
	>> iunit(_maxDisplacement,mm)
	>> _underlyingEventHandler;
	}


	void ClusterHadronizationHandler::Init() {

	static ClassDocumentation<ClusterHadronizationHandler> documentation
	("This is the main handler class for the Cluster Hadronization",
	"The hadronization was performed using the cluster model of \\cite{Webber:1983if}.",
	"%\\cite{Webber:1983if}\n"
	"\\bibitem{Webber:1983if}\n"
	" B.~R.~Webber,\n"
	" ``A QCD Model For Jet Fragmentation Including Soft Gluon Interference,''\n"
	" Nucl.\\ Phys.\\ B {\\bf 238}, 492 (1984).\n"
	" %%CITATION = NUPHA,B238,492;%%\n"
	// main manual
	);

	static Reference<ClusterHadronizationHandler,PartonSplitter>
	interfacePartonSplitter("PartonSplitter",
	"A reference to the PartonSplitter object",
	&Herwig::ClusterHadronizationHandler::_partonSplitter,
	false, false, true, false);

	static Reference<ClusterHadronizationHandler,ClusterFinder>
	interfaceClusterFinder("ClusterFinder",
	"A reference to the ClusterFinder object",
	&Herwig::ClusterHadronizationHandler::_clusterFinder,
	false, false, true, false);

	static Reference<ClusterHadronizationHandler,ColourReconnector>
	interfaceColourReconnector("ColourReconnector",
	"A reference to the ColourReconnector object",
	&Herwig::ClusterHadronizationHandler::_colourReconnector,
	false, false, true, false);

	static Reference<ClusterHadronizationHandler,ClusterFissioner>
	interfaceClusterFissioner("ClusterFissioner",
	"A reference to the ClusterFissioner object",
	&Herwig::ClusterHadronizationHandler::_clusterFissioner,
	false, false, true, false);

	static Reference<ClusterHadronizationHandler,LightClusterDecayer>
	interfaceLightClusterDecayer("LightClusterDecayer",
	"A reference to the LightClusterDecayer object",
	&Herwig::ClusterHadronizationHandler::_lightClusterDecayer,
	false, false, true, false);

	static Reference<ClusterHadronizationHandler,ClusterDecayer>
	interfaceClusterDecayer("ClusterDecayer",
	"A reference to the ClusterDecayer object",
	&Herwig::ClusterHadronizationHandler::_clusterDecayer,
	false, false, true, false);

	static Parameter<ClusterHadronizationHandler,Energy2> interfaceMinVirtuality2
	("MinVirtuality2",
	"Minimum virtuality^2 of partons to use in calculating distances (unit [GeV2]).",
	&ClusterHadronizationHandler::_minVirtuality2, GeV2, 0.1GeV2, ZERO, 10.0GeV2,false,false,false);

	static Parameter<ClusterHadronizationHandler,Length> interfaceMaxDisplacement
	("MaxDisplacement",
	"Maximum displacement that is allowed for a particle (unit [millimeter]).",
	&ClusterHadronizationHandler::_maxDisplacement, mm, 1.0e-10*mm,
	0.0mm, 1.0e-9mm,false,false,false);

	static Reference<ClusterHadronizationHandler,StepHandler> interfaceUnderlyingEventHandler
	("UnderlyingEventHandler",
	"Pointer to the handler for the Underlying Event. "
	"Set to NULL to disable.",
	&ClusterHadronizationHandler::_underlyingEventHandler, false, false, true, true, false);
	}

	-void ClusterHadronizationHandler::doinitrun() {
	- HadronizationHandler::doinitrun();
	- // The run initialization is used here to all Cluster to have access to the
	- // ClusterHadronizationHandler class instance, via a static pointer.
	- Cluster::setPointerClusterHadHandler(this);
	-}
	-
	namespace {
	void extractChildren(tPPtr p, set<PPtr> & all) {
	if (p->children().empty()) return;

	for (PVector::const_iterator child = p->children().begin();
	child != p->children().end(); ++child) {
	all.insert(*child);
	extractChildren(*child, all);
	}
	}
	}

	void ClusterHadronizationHandler::
	handle(EventHandler & ch, const tPVector & tagged,
	const Hint &) {
	useMe();
	+ currentHandler_ = this;
	PVector currentlist(tagged.begin(),tagged.end());
	// set the scale for coloured particles to just above the gluon mass squared
	// if less than this so they are classed as perturbative
	Energy2 Q02 = 1.01*sqr(getParticleData(ParticleID::g)->constituentMass());
	for(unsigned int ix=0;ix<currentlist.size();++ix) {
	if(currentlist[ix]->scale()<Q02) currentlist[ix]->scale(Q02);
	}
	// split the gluons
	_partonSplitter->split(currentlist);
	// form the clusters
	ClusterVector clusters =
	_clusterFinder->formClusters(currentlist);

	_clusterFinder->reduceToTwoComponents(clusters);

	// perform colour reconnection if needed and then
	// decay the clusters into one hadron
	bool lightOK = false;
	short tried = 0;
	const ClusterVector savedclusters = clusters;
	tPVector finalHadrons; // only needed for partonic decayer
	while (!lightOK && tried++ < 10) {

	// no colour reconnection with baryon-number-violating (BV) clusters
	ClusterVector CRclusters, BVclusters;
	CRclusters.reserve( clusters.size() );
	BVclusters.reserve( clusters.size() );
	for (size_t ic = 0; ic < clusters.size(); ++ic) {
	ClusterPtr cl = clusters.at(ic);
	bool hasClusterParent = false;
	for (unsigned int ix=0; ix < cl->parents().size(); ++ix) {
	if (cl->parents()[ix]->id() == ParticleID::Cluster) {
	hasClusterParent = true;
	break;
	}
	}
	if (cl->numComponents() > 2 \|\| hasClusterParent) BVclusters.push_back(cl);
	else CRclusters.push_back(cl);
	}

	// colour reconnection
	_colourReconnector->rearrange(CRclusters);

	// tag new clusters as children of the partons to hadronize
	_setChildren(CRclusters);

	// recombine vectors of (possibly) reconnected and BV clusters
	clusters.clear();
	clusters.insert( clusters.end(), CRclusters.begin(), CRclusters.end() );
	clusters.insert( clusters.end(), BVclusters.begin(), BVclusters.end() );

	// fission of heavy clusters
	// NB: during cluster fission, light hadrons might be produced straight away
	finalHadrons = _clusterFissioner->fission(clusters,isSoftUnderlyingEventON());

	lightOK = _lightClusterDecayer->decay(clusters,finalHadrons);

	// if the decay of the light clusters was not successful, undo the cluster
	// fission and decay steps and revert to the original state of the event
	// record
	if (!lightOK) {
	clusters = savedclusters;
	for_each(clusters.begin(),
	clusters.end(),
	mem_fun(&Particle::undecay));
	}
	}
	- if (!lightOK)
	+ if (!lightOK) {
	+ // currentHandler_ = 0;
	throw Exception("CluHad::handle(): tried LightClusterDecayer 10 times!",
	Exception::eventerror);
	+ }


	// decay the remaining clusters
	_clusterDecayer->decay(clusters,finalHadrons);

	// *****************************************
	// *****************************************
	// *****************************************

	StepPtr pstep = newStep();
	set<PPtr> allDecendants;
	for (tPVector::const_iterator it = tagged.begin();
	it != tagged.end(); ++it) {
	extractChildren(*it, allDecendants);
	}

	for(set<PPtr>::const_iterator it = allDecendants.begin();
	it != allDecendants.end(); ++it) {
	// this is a workaround because the set sometimes
	// re-orders parents after their children
	if ((*it)->children().empty())
	pstep->addDecayProduct(*it);
	else {
	pstep->addDecayProduct(*it);
	pstep->addIntermediate(*it);
	}
	}

	// *****************************************
	// *****************************************
	// *****************************************

	// soft underlying event if needed
	if (isSoftUnderlyingEventON()) {
	assert(_underlyingEventHandler);
	ch.performStep(_underlyingEventHandler,Hint::Default());
	}
	// zero all positions
	// extract all particles from the event
	tEventPtr event=ch.currentEvent();
	vector<tPPtr> particles;
	particles.reserve(256);
	event->select(back_inserter(particles), ThePEG::AllSelector());
	// and the final-state particles
	set<tPPtr> finalstate;
	event->selectFinalState(inserter(finalstate));
	for(vector<tPPtr>::const_iterator pit=particles.begin();
	pit!=particles.end();++pit) {
	// if a final-state particle just zero production
	if(finalstate.find(*pit)!=finalstate.end()) {
	(**pit).setVertex(LorentzPoint());
	}
	// if not zero the lot
	else {
	(**pit).setVertex(LorentzPoint());
	(**pit).setLifeLength(LorentzDistance());
	}
	}
	+ // currentHandler_ = 0;
	}

	void ClusterHadronizationHandler::_setChildren(ClusterVector clusters) const {

	// erase existing information about the partons' children
	tPVector partons;
	for (ClusterVector::const_iterator cl = clusters.begin();
	cl != clusters.end(); cl++) {
	partons.push_back( (*cl)->colParticle() );
	partons.push_back( (*cl)->antiColParticle() );
	}
	for_each(partons.begin(), partons.end(), mem_fun(&Particle::undecay));

	// give new parents to the clusters: their constituents
	for (ClusterVector::iterator cl = clusters.begin();
	cl != clusters.end(); cl++) {
	(cl)->colParticle()->addChild(cl);
	(cl)->antiColParticle()->addChild(cl);
	}

	}

	diff --git a/Hadronization/ClusterHadronizationHandler.h b/Hadronization/ClusterHadronizationHandler.h
	--- a/Hadronization/ClusterHadronizationHandler.h
	+++ b/Hadronization/ClusterHadronizationHandler.h
	@@ -1,213 +1,217 @@
	// -- C++ --
	//
	// ClusterHadronizationHandler.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_ClusterHadronizationHandler_H
	#define HERWIG_ClusterHadronizationHandler_H

	#include <ThePEG/Handlers/HadronizationHandler.h>
	#include "PartonSplitter.h"
	#include "ClusterFinder.h"
	#include "ColourReconnector.h"
	#include "ClusterFissioner.h"
	#include "LightClusterDecayer.h"
	#include "ClusterDecayer.h"
	#include "ClusterHadronizationHandler.fh"

	namespace Herwig {
	using namespace ThePEG;


	/** \ingroup Hadronization
	* \class ClusterHadronizationHandler
	* \brief Class that controls the cluster hadronization algorithm.
	* \author Philip Stephens // cerr << *ch.currentEvent() << '\n';
	cerr << finalHadrons.size() << '\n';

	cerr << "Finished hadronizing \n";

	* \author Alberto Ribon
	*
	* This class is the main driver of the cluster hadronization: it is
	* responsible for the proper handling of all other specific collaborating
	* classes PartonSplitter, ClusterFinder, ColourReconnector, ClusterFissioner,
	* LightClusterDecayer, ClusterDecayer;
	* and for the storing of the produced particles in the Event record.
	*
	* @see PartonSplitter
	* @see ClusterFinder
	* @see ColourReconnector
	* @see ClusterFissioner
	* @see LightClusterDecayer
	* @see ClusterDecayer
	* @see Cluster
	* @see \ref ClusterHadronizationHandlerInterfaces "The interfaces"
	* defined for ClusterHadronizationHandler.
	*/
	class ClusterHadronizationHandler: public HadronizationHandler {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	ClusterHadronizationHandler()
	: _minVirtuality2( 0.1GeV2 ), _maxDisplacement( 1.0e-10mm )
	{}

	//@}

	public:

	/**
	* The main method which manages the all cluster hadronization.
	*
	* This routine directs "traffic". It determines which function is called
	* and on which particles/clusters. This function also handles the
	* situation of vetos on the hadronization.
	*/
	virtual void handle(EventHandler & ch, const tPVector & tagged,
	const Hint & hint);

	/**
	* It returns minimum virtuality^2 of partons to use in calculating
	* distances. It is used both in the Showering and Hadronization.
	*/
	Energy2 minVirtuality2() const
	{ return _minVirtuality2; }

	/**
	* It returns the maximum displacement that is allowed for a particle
	* (used to determine the position of a cluster with two components).
	*/
	Length maxDisplacement() const
	{ return _maxDisplacement; }

	/**
	* It returns true/false according if the soft underlying model
	* is switched on/off.
	*/
	bool isSoftUnderlyingEventON() const
	{ return _underlyingEventHandler; }
	+
	+ /**
	+ * pointer to "this", the current HadronizationHandler.
	+ */
	+ static const ClusterHadronizationHandler * currentHandler() {
	+ assert(currentHandler_);
	+ return currentHandler_;
	+ }

	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);
	//@}

	/**
	* Standard Init function used to initialize the interfaces.
	*/
	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;
	//@}

	-protected:
	-
	- /** @name Standard Interfaced functions. */
	- //@{
	- /**
	- * Initialize this object at the begining of the run phase.
	- */
	- virtual void doinitrun();
	- //@}
	-
	private:

	/**
	* Private and non-existent assignment operator.
	*/
	ClusterHadronizationHandler & operator=(const ClusterHadronizationHandler &);

	/**
	* This is a pointer to a Herwig::PartonSplitter object.
	*/
	PartonSplitterPtr _partonSplitter;

	/**
	* This is a pointer to a Herwig::ClusterFinder object.
	*/
	ClusterFinderPtr _clusterFinder;

	/**
	* This is a pointer to a Herwig::ColourReconnector object.
	*/
	ColourReconnectorPtr _colourReconnector;

	/**
	* This is a pointer to a Herwig::ClusterFissioner object.
	*/
	ClusterFissionerPtr _clusterFissioner;

	/**
	* This is a pointer to a Herwig::LightClusterDecayer object.
	*/
	LightClusterDecayerPtr _lightClusterDecayer;

	/**
	* This is a pointer to a Herwig::ClusterDecayer object.
	*/
	ClusterDecayerPtr _clusterDecayer;

	/**
	* The minimum virtuality^2 of partons to use in calculating
	* distances.
	*/
	Energy2 _minVirtuality2;

	/**
	* The maximum displacement that is allowed for a particle
	* (used to determine the position of a cluster with two components).
	*/
	Length _maxDisplacement;

	/**
	* The pointer to the Underlying Event handler.
	*/
	StepHdlPtr _underlyingEventHandler;

	/**
	* Tag the constituents of the clusters as their parents
	*/
	void _setChildren(ClusterVector clusters) const;

	+ /**
	+ * pointer to "this", the current HadronizationHandler.
	+ */
	+ static ClusterHadronizationHandler * currentHandler_;
	+
	+
	};


	}

	#endif /* HERWIG_ClusterHadronizationHandler_H */
	diff --git a/MatrixElement/ProductionMatrixElement.cc b/MatrixElement/ProductionMatrixElement.cc
	--- a/MatrixElement/ProductionMatrixElement.cc
	+++ b/MatrixElement/ProductionMatrixElement.cc
	@@ -1,157 +1,445 @@
	// -- C++ --
	//
	// ProductionMatrixElement.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ProductionMatrixElement class.
	//
	// Author: Peter Richardson
	//
	#include "ProductionMatrixElement.h"
	#include "ThePEG/Interface/ClassDocumentation.h"

	using namespace Herwig;
	-
	+
	// calculate a decay matrix for one of the incoming particles
	RhoDMatrix ProductionMatrixElement::
	calculateDMatrix(int id, const RhoDMatrix & rhoin,
	const vector<RhoDMatrix> & rhoout) const {
	// vectors for the helicities
	vector<unsigned int> ihel1(_outspin.size()+2),ihel2(_outspin.size()+2);
	// rhomatrix to be returned
	RhoDMatrix output(_inspin[id], false);
	// loop over all helicity components of the matrix element
	// outer loop
	Complex temp;
	unsigned int ix,iy;
	int ixa,iya;
	for(ix=0;ix<_matrixelement.size();++ix)
	{
	// map the vector index to the helicities
	for(ixa=_outspin.size()+1;ixa>=0;--ixa)
	{ihel1[ixa]=(ix%_constants[ixa])/_constants[ixa+1];}
	// inner loop
	for(iy=0;iy<_matrixelement.size();++iy)
	{
	// map the vector index to the helicities
	for(iya=_outspin.size()+1;iya>=0;--iya)
	{ihel2[iya]=(iy%_constants[iya])/_constants[iya+1];}
	// matrix element piece
	temp=_matrixelement[ix]*conj(_matrixelement[iy]);
	// spin density matrices for the outgoing particles
	for(unsigned int iz=0;iz<_outspin.size();++iz)
	{temp*=rhoout[iz](ihel1[iz+2],ihel2[iz+2]);}
	// construct the spin density matrix
	if(id==0)
	{
	temp*=rhoin(ihel1[1],ihel2[1]);
	output(ihel1[0],ihel2[0])+=temp;
	}
	else
	{
	temp*=rhoin(ihel1[0],ihel2[0]);
	output(ihel1[1],ihel2[1])+=temp;
	}
	}
	}
	// return the answer
	return output;
	}

	// calculate the rho matrix for a given outgoing particle
	RhoDMatrix ProductionMatrixElement::
	calculateRhoMatrix(int id,const RhoDMatrix & rhoin0,
	const RhoDMatrix & rhoin1,
	const vector<RhoDMatrix> & rhoout) const {
	unsigned int ix,iy;
	int ixa,iya;
	// vectors for the helicities
	vector<unsigned int> ihel1(_outspin.size()+2),ihel2(_outspin.size()+2);
	// rhomatrix to be returned
	RhoDMatrix output(_outspin[id], false);
	// loop over all helicity components of the matrix element
	// outer loop
	Complex temp;
	for(ix=0;ix<_matrixelement.size();++ix) {
	// map the vector index to the helicities
	for(ixa=_outspin.size()+1;ixa>=0;--ixa)
	ihel1[ixa]=(ix%_constants[ixa])/_constants[ixa+1];
	// inner loop
	for(iy=0;iy<_matrixelement.size();++iy) {
	// map the vector index to the helicities
	for(iya=_outspin.size()+1;iya>=0;--iya)
	ihel2[iya] = (iy%_constants[iya])/_constants[iya+1];
	// matrix element piece
	temp = _matrixelement[ix]*conj(_matrixelement[iy]);
	// spin denisty matrix for the incoming particles
	temp *= rhoin0(ihel1[0],ihel2[0]);
	temp *= rhoin1(ihel1[1],ihel2[1]);
	// spin density matrix for the outgoing particles
	for(unsigned int iz=0;iz<_outspin.size()-1;++iz) {
	if(int(iz)<id) temp *= rhoout[iz](ihel1[iz+2],ihel2[iz+2]);
	else temp *= rhoout[iz](ihel1[iz+3],ihel2[iz+3]);
	}
	output(ihel1[id+2],ihel2[id+2])+=temp;
	}
	}
	// normalise the matrix so it has unit trace
	output.normalize();
	// return the answer
	return output;
	}

	double ProductionMatrixElement::average() const {
	double output(0.);
	for(unsigned int ix=0;ix<_matrixelement.size();++ix) {
	output += norm(_matrixelement[ix]);
	}
	return output;
	}

	double ProductionMatrixElement::average(const RhoDMatrix & in1,
	const RhoDMatrix & in2) const {
	Complex output(0.);
	for( int ihel1=0;ihel1<int(_inspin[0]);++ihel1) {
	for( int ihel2=0;ihel2<int(_inspin[1]);++ihel2) {
	int loc1 = ihel1_constants[1] + ihel2_constants[2];
	for( int jhel1=0;jhel1<int(_inspin[0]);++jhel1) {
	for( int jhel2=0;jhel2<int(_inspin[1]);++jhel2) {
	int loc2 = jhel1_constants[1] + jhel2_constants[2];
	Complex fact = in1(ihel1,jhel1)*in2(ihel2,jhel2);
	for(int ohel=0;ohel<_constants[2];++ohel) {
	output += fact
	_matrixelement[loc1+ohel]conj(_matrixelement[loc2+ohel]);
	}
	}
	}
	}
	}
	return real(output);
	}

	Complex ProductionMatrixElement::average(const ProductionMatrixElement & me2,
	const RhoDMatrix & in1,
	const RhoDMatrix & in2) const {
	Complex output(0.);
	for( int ihel1=0;ihel1<int(_inspin[0]);++ihel1) {
	for( int ihel2=0;ihel2<int(_inspin[1]);++ihel2) {
	int loc1 = ihel1_constants[1] + ihel2_constants[2];
	for( int jhel1=0;jhel1<int(_inspin[0]);++jhel1) {
	for( int jhel2=0;jhel2<int(_inspin[1]);++jhel2) {
	int loc2 = jhel1_constants[1] + jhel2_constants[2];
	Complex fact = in1(ihel1,jhel1)*in2(ihel2,jhel2);
	for(int ohel=0;ohel<_constants[2];++ohel) {
	output += fact
	_matrixelement[loc1+ohel]conj(me2._matrixelement[loc2+ohel]);
	}
	}
	}
	}
	}
	return output;
	}
	+
	+ProductionMatrixElement::ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out) {
	+ _nout=2;
	+ _inspin.resize(2);
	+ _inspin[0]=in1;
	+ _inspin[1]=in2;
	+ _outspin.push_back(out);
	+ setMESize();
	+}
	+
	+ProductionMatrixElement::ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out1,
	+ PDT::Spin out2) {
	+ _nout=2;
	+ _inspin.resize(2);
	+ _inspin[0]=in1;
	+ _inspin[1]=in2;
	+ _outspin.push_back(out1);
	+ _outspin.push_back(out2);
	+ setMESize();
	+}
	+
	+ProductionMatrixElement::ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out1,
	+ PDT::Spin out2,PDT::Spin out3) {
	+ _inspin.resize(2);
	+ _nout=3;
	+ _inspin[0]=in1;
	+ _inspin[1]=in2;
	+ _outspin.push_back(out1);
	+ _outspin.push_back(out2);
	+ _outspin.push_back(out3);
	+ setMESize();
	+}
	+
	+ProductionMatrixElement::ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out1,
	+ PDT::Spin out2,PDT::Spin out3, PDT::Spin out4) {
	+ _nout=4;
	+ _inspin.resize(2);
	+ _inspin[0]=in1;
	+ _inspin[1]=in2;
	+ _outspin.push_back(out1);
	+ _outspin.push_back(out2);
	+ _outspin.push_back(out3);
	+ _outspin.push_back(out4);
	+ setMESize();
	+}
	+
	+ProductionMatrixElement::ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out1,
	+ PDT::Spin out2,PDT::Spin out3, PDT::Spin out4,
	+ PDT::Spin out5) {
	+ _nout=5;
	+ _inspin.resize(2);
	+ _inspin[0]=in1;
	+ _inspin[1]=in2;
	+ _outspin.push_back(out1);
	+ _outspin.push_back(out2);
	+ _outspin.push_back(out3);
	+ _outspin.push_back(out4);
	+ _outspin.push_back(out5);
	+ setMESize();
	+}
	+
	+ProductionMatrixElement::ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out1,
	+ PDT::Spin out2,PDT::Spin out3, PDT::Spin out4,
	+ PDT::Spin out5, PDT::Spin out6) {
	+ _nout=6;
	+ _inspin.resize(2);
	+ _inspin[0]=in1;
	+ _inspin[1]=in2;
	+ _outspin.push_back(out1);
	+ _outspin.push_back(out2);
	+ _outspin.push_back(out3);
	+ _outspin.push_back(out4);
	+ _outspin.push_back(out5);
	+ _outspin.push_back(out6);
	+ setMESize();
	+}
	+
	+ProductionMatrixElement::ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,vector<PDT::Spin> out) {
	+ _inspin.resize(2);
	+ _nout=out.size();
	+ _inspin[0]=in1;
	+ _inspin[1]=in2;
	+ _outspin=out;
	+ setMESize();
	+}
	+
	+Complex ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out) const {
	+ assert(_outspin.size()==1);
	+ unsigned int iloc = in1_constants[1] + in2_constants[2] + out*_constants[3];
	+ assert(iloc<_matrixelement.size());
	+ return _matrixelement[iloc];
	+}
	+
	+Complex & ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out) {
	+ assert(_outspin.size()==1);
	+ unsigned int iloc = in1_constants[1] + in2_constants[2] + out*_constants[3];
	+ assert(iloc<_matrixelement.size());
	+ return _matrixelement[iloc];
	+}
	+
	+Complex ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out1,unsigned int out2) const {
	+ assert(_outspin.size()==2);
	+ unsigned int iloc = in1_constants[1] + in2_constants[2] +
	+ out1_constants[3] + out2_constants[4];
	+ assert(iloc<_matrixelement.size());
	+ return _matrixelement[iloc];
	+}
	+
	+Complex & ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out1,unsigned int out2) {
	+ assert(_outspin.size()==2);
	+ unsigned int iloc = in1_constants[1] + in2_constants[2] +
	+ out1_constants[3] + out2_constants[4];
	+ assert(iloc<_matrixelement.size());
	+ return _matrixelement[iloc];
	+}
	+
	+Complex ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out1,unsigned int out2,
	+ unsigned int out3) const {
	+ assert(_outspin.size()==3);
	+ vector<unsigned int> ivec(5);
	+ ivec[0]=in1;
	+ ivec[1]=in2;
	+ ivec[2]=out1;
	+ ivec[3]=out2;
	+ ivec[4]=out3;
	+ return (*this)(ivec);
	+}
	+
	+Complex & ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out1,unsigned int out2,
	+ unsigned int out3) {
	+ assert(_outspin.size()==3);
	+ vector<unsigned int> ivec(5);
	+ ivec[0]=in1;
	+ ivec[1]=in2;
	+ ivec[2]=out1;
	+ ivec[3]=out2;
	+ ivec[4]=out3;
	+ return (*this)(ivec);
	+}
	+
	+Complex ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out1,unsigned int out2,
	+ unsigned int out3,unsigned int out4) const {
	+ assert(_outspin.size()==4);
	+ vector<unsigned int> ivec(6);
	+ ivec[0]=in1;
	+ ivec[1]=in2;
	+ ivec[2]=out1;
	+ ivec[3]=out2;
	+ ivec[4]=out3;
	+ ivec[5]=out4;
	+ return (*this)(ivec);
	+}
	+
	+Complex & ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out1,unsigned int out2,
	+ unsigned int out3, unsigned int out4) {
	+ assert(_outspin.size()==4);
	+ vector<unsigned int> ivec(6);
	+ ivec[0]=in1;
	+ ivec[1]=in2;
	+ ivec[2]=out1;
	+ ivec[3]=out2;
	+ ivec[4]=out3;
	+ ivec[5]=out4;
	+ return (*this)(ivec);
	+}
	+
	+Complex ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out1,unsigned int out2,
	+ unsigned int out3,unsigned int out4,
	+ unsigned int out5) const {
	+ assert(_outspin.size()==5);
	+ vector<unsigned int> ivec(7);
	+ ivec[0]=in1;
	+ ivec[1]=in2;
	+ ivec[2]=out1;
	+ ivec[3]=out2;
	+ ivec[4]=out3;
	+ ivec[5]=out4;
	+ ivec[6]=out5;
	+ return (*this)(ivec);
	+}
	+
	+Complex & ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out1,unsigned int out2,
	+ unsigned int out3, unsigned int out4,
	+ unsigned int out5) {
	+ assert(_outspin.size()==5);
	+ vector<unsigned int> ivec(7);
	+ ivec[0]=in1;
	+ ivec[1]=in2;
	+ ivec[2]=out1;
	+ ivec[3]=out2;
	+ ivec[4]=out3;
	+ ivec[5]=out4;
	+ ivec[6]=out5;
	+ return (*this)(ivec);
	+}
	+
	+Complex ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out1,unsigned int out2,
	+ unsigned int out3,unsigned int out4,
	+ unsigned int out5,unsigned int out6) const {
	+ assert(_outspin.size()==6);
	+ vector<unsigned int> ivec(8);
	+ ivec[0]=in1;
	+ ivec[1]=in2;
	+ ivec[2]=out1;
	+ ivec[3]=out2;
	+ ivec[4]=out3;
	+ ivec[5]=out4;
	+ ivec[6]=out5;
	+ ivec[7]=out6;
	+ return (*this)(ivec);
	+}
	+
	+Complex & ProductionMatrixElement::operator () (unsigned int in1,unsigned int in2,
	+ unsigned int out1,unsigned int out2,
	+ unsigned int out3, unsigned int out4,
	+ unsigned int out5, unsigned int out6) {
	+ assert(_outspin.size()==6);
	+ vector<unsigned int> ivec(8);
	+ ivec[0]=in1;
	+ ivec[1]=in2;
	+ ivec[2]=out1;
	+ ivec[3]=out2;
	+ ivec[4]=out3;
	+ ivec[5]=out4;
	+ ivec[6]=out5;
	+ ivec[7]=out6;
	+ return (*this)(ivec);
	+}
	+
	+Complex ProductionMatrixElement::operator () (vector<unsigned int> hel) const {
	+ assert(_outspin.size() == hel.size()-2);
	+ unsigned int iloc(0),ix;
	+ // incoming and outgoing particles
	+ for(ix=0;ix<hel.size();++ix)
	+ iloc += hel[ix]*_constants[ix+1];
	+ assert(iloc<_matrixelement.size());
	+ return _matrixelement[iloc];
	+}
	+
	+Complex & ProductionMatrixElement::operator () (vector<unsigned int> hel) {
	+ assert(_outspin.size() == hel.size()-2);
	+ unsigned int iloc=0,ix;
	+ // incoming particles
	+ for(ix=0;ix<hel.size();++ix)
	+ iloc += hel[ix]*_constants[ix+1];
	+ assert(iloc<_matrixelement.size());
	+ return _matrixelement[iloc];
	+}
	+//@}
	+
	+void ProductionMatrixElement::reset(const ProductionMatrixElement & x) const {
	+ _nout = x._nout;
	+ _inspin = x._inspin;
	+ _outspin = x._outspin;
	+ _matrixelement = x._matrixelement;
	+ _constants = x._constants;
	+}
	+
	+void ProductionMatrixElement::setMESize() {
	+ unsigned int ix;
	+ int isize=_inspin[0]*_inspin[1];
	+ for(ix=0;ix<_outspin.size();++ix)
	+ isize*=_outspin[ix];
	+ // zero the matrix element
	+ _matrixelement.resize(isize,0.);
	+ // set up the constants for the mapping of helicity to vectro index
	+ _constants.resize(_outspin.size()+3);
	+ unsigned int temp=1;
	+ for(ix=_outspin.size()+1;ix>1;--ix) {
	+ temp*=_outspin[ix-2];
	+ _constants[ix]=temp;
	+ }
	+ temp*=_inspin[1];_constants[1]=temp;
	+ temp*=_inspin[0];_constants[0]=temp;
	+ _constants[_outspin.size()+2]=1;
	+}
	+
	diff --git a/MatrixElement/ProductionMatrixElement.h b/MatrixElement/ProductionMatrixElement.h
	--- a/MatrixElement/ProductionMatrixElement.h
	+++ b/MatrixElement/ProductionMatrixElement.h
	@@ -1,663 +1,454 @@
	// -- C++ --
	//
	// ProductionMatrixElement.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_ProductionMatrixElement_H
	#define HERWIG_ProductionMatrixElement_H
	//
	// This is the declaration of the ProductionMatrixElement class.

	#include <ThePEG/Config/ThePEG.h>
	#include <ThePEG/Utilities/ClassDescription.h>
	#include <ThePEG/EventRecord/RhoDMatrix.h>

	namespace Herwig {


	using namespace ThePEG;

	/** \ingroup Helicity
	*
	* The storage of the helicity amplitude expression for the matrix element
	* of a hard process. Two incoming particles and an arbitary number of
	* external particles are supported.
	*
	* @see DecayMatrixElement
	* @see RhoDMatrix
	* @see HardVertex
	*
	* \author Peter Richardson
	*/

	class ProductionMatrixElement {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Constructor for 2-1 scattering.
	* @param in1 \f$2S+1\f$ for the first incoming particle.
	* @param in2 \f$2S+1\f$ for the second incoming particle.
	* @param out \f$2S+1\f$ for the outgoing particle.
	*/
	- ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out) {
	- _nout=2;
	- _inspin.resize(2);
	- _inspin[0]=in1;
	- _inspin[1]=in2;
	- _outspin.push_back(out);
	- setMESize();
	- }
	+ ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out);

	/**
	* Constructor for 2-2 scattering.
	* @param in1 \f$2S+1\f$ for the first incoming particle.
	* @param in2 \f$2S+1\f$ for the second incoming particle.
	* @param out1 \f$2S+1\f$ for the first outgoing particle.
	* @param out2 \f$2S+1\f$ for the second outgoing particle.
	*/
	ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out1,
	- PDT::Spin out2) {
	- _nout=2;
	- _inspin.resize(2);
	- _inspin[0]=in1;
	- _inspin[1]=in2;
	- _outspin.push_back(out1);
	- _outspin.push_back(out2);
	- setMESize();
	- }
	+ PDT::Spin out2);

	/**
	* Constructor for 2-3 scattering.
	* @param in1 \f$2S+1\f$ for the first incoming particle.
	* @param in2 \f$2S+1\f$ for the second incoming particle.
	* @param out1 \f$2S+1\f$ for the first outgoing particle.
	* @param out2 \f$2S+1\f$ for the second outgoing particle.
	* @param out3 \f$2S+1\f$ for the third outgoing particle.
	*/
	ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out1,
	- PDT::Spin out2,PDT::Spin out3) {
	- _inspin.resize(2);
	- _nout=3;
	- _inspin[0]=in1;
	- _inspin[1]=in2;
	- _outspin.push_back(out1);
	- _outspin.push_back(out2);
	- _outspin.push_back(out3);
	- setMESize();
	- }
	+ PDT::Spin out2,PDT::Spin out3);

	/**
	* Constructor for 2-4 scattering.
	* @param in1 \f$2S+1\f$ for the first incoming particle.
	* @param in2 \f$2S+1\f$ for the second incoming particle.
	* @param out1 \f$2S+1\f$ for the first outgoing particle.
	* @param out2 \f$2S+1\f$ for the second outgoing particle.
	* @param out3 \f$2S+1\f$ for the third outgoing particle.
	* @param out4 \f$2S+1\f$ for the fourth outgoing particle.
	*/
	ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out1,
	- PDT::Spin out2,PDT::Spin out3, PDT::Spin out4) {
	- _nout=4;
	- _inspin.resize(2);
	- _inspin[0]=in1;
	- _inspin[1]=in2;
	- _outspin.push_back(out1);
	- _outspin.push_back(out2);
	- _outspin.push_back(out3);
	- _outspin.push_back(out4);
	- setMESize();
	- }
	+ PDT::Spin out2,PDT::Spin out3, PDT::Spin out4);

	/**
	* Constructor for 2-5 scattering.
	* @param in1 \f$2S+1\f$ for the first incoming particle.
	* @param in2 \f$2S+1\f$ for the second incoming particle.
	* @param out1 \f$2S+1\f$ for the first outgoing particle.
	* @param out2 \f$2S+1\f$ for the second outgoing particle.
	* @param out3 \f$2S+1\f$ for the third outgoing particle.
	* @param out4 \f$2S+1\f$ for the fourth outgoing particle.
	* @param out5 \f$2S+1\f$ for the fifth outgoing particle.
	*/
	ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out1,
	PDT::Spin out2,PDT::Spin out3, PDT::Spin out4,
	- PDT::Spin out5) {
	- _nout=5;
	- _inspin.resize(2);
	- _inspin[0]=in1;
	- _inspin[1]=in2;
	- _outspin.push_back(out1);
	- _outspin.push_back(out2);
	- _outspin.push_back(out3);
	- _outspin.push_back(out4);
	- _outspin.push_back(out5);
	- setMESize();
	- }
	+ PDT::Spin out5);

	/**
	* Constructor for 2-6 scattering.
	* @param in1 \f$2S+1\f$ for the first incoming particle.
	* @param in2 \f$2S+1\f$ for the second incoming particle.
	* @param out1 \f$2S+1\f$ for the first outgoing particle.
	* @param out2 \f$2S+1\f$ for the second outgoing particle.
	* @param out3 \f$2S+1\f$ for the third outgoing particle.
	* @param out4 \f$2S+1\f$ for the fourth outgoing particle.
	* @param out5 \f$2S+1\f$ for the fifth outgoing particle.
	* @param out6 \f$2S+1\f$ for the sixth outgoing particle.
	*/
	ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,PDT::Spin out1,
	PDT::Spin out2,PDT::Spin out3, PDT::Spin out4,
	- PDT::Spin out5, PDT::Spin out6) {
	- _nout=6;
	- _inspin.resize(2);
	- _inspin[0]=in1;
	- _inspin[1]=in2;
	- _outspin.push_back(out1);
	- _outspin.push_back(out2);
	- _outspin.push_back(out3);
	- _outspin.push_back(out4);
	- _outspin.push_back(out5);
	- _outspin.push_back(out6);
	- setMESize();
	- }
	+ PDT::Spin out5, PDT::Spin out6);

	/**
	* Constructor for 2-n scattering.
	* @param in1 \f$2S+1\f$ for the first incoming particle.
	* @param in2 \f$2S+1\f$ for the second incoming particle.
	* @param out A vector containing \f$2S+1\f$ for the outgoing particles.
	*/
	- ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,vector<PDT::Spin> out) {
	- _inspin.resize(2);
	- _nout=out.size();
	- _inspin[0]=in1;
	- _inspin[1]=in2;
	- _outspin=out;
	- setMESize();
	- }
	+ ProductionMatrixElement(PDT::Spin in1,PDT::Spin in2,vector<PDT::Spin> out);

	/**
	* Default constructor.
	*/
	- ProductionMatrixElement() {};
	+ ProductionMatrixElement() {}
	//@}

	public:

	/** @name Access to the spins of the particles. */
	//@{
	/**
	* Get the spins of the incoming particles particle
	* @return A vector containing \f$2S+1\f$ for the two incoming particles.
	*/
	vector<PDT::Spin> inspin() {return _inspin;}

	/**
	* Get the spins of the outgoing particles.
	* @return A vector containing \f$2S+1\f$ for the outgoing particles.
	*/
	vector<PDT::Spin> outspin() {return _outspin;}
	//@}

	public:

	/** @name Access to the individual helicity components. */
	//@{
	/**
	* Access the helicity components for a 2-1 scattering. This method supplies
	* the component but does not allow it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out The helicity of the outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex operator () (unsigned int in1,unsigned int in2,
	- unsigned int out) const {
	- assert(_outspin.size()==1);
	- unsigned int iloc = in1_constants[1] + in2_constants[2] + out*_constants[3];
	- assert(iloc<_matrixelement.size());
	- return _matrixelement[iloc];
	- }
	+ unsigned int out) const;

	/**
	* Access the helicity components for a 2-1 scattering. This method supplies
	* the component and allows it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out The helicity of the outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex & operator () (unsigned int in1,unsigned int in2,
	- unsigned int out) {
	- assert(_outspin.size()==1);
	- unsigned int iloc = in1_constants[1] + in2_constants[2] + out*_constants[3];
	- assert(iloc<_matrixelement.size());
	- return _matrixelement[iloc];
	- }
	+ unsigned int out);

	/**
	* Access the helicity components for a 2-2 scattering. This method supplies
	* the component but does not allow it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out1 The helicity of the first outgoing particle.
	* @param out2 The helicity of the second outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex operator () (unsigned int in1,unsigned int in2,
	- unsigned int out1,unsigned int out2) const {
	- assert(_outspin.size()==2);
	- unsigned int iloc = in1_constants[1] + in2_constants[2] +
	- out1_constants[3] + out2_constants[4];
	- assert(iloc<_matrixelement.size());
	- return _matrixelement[iloc];
	- }
	+ unsigned int out1,unsigned int out2) const;

	/**
	* Access the helicity components for a 2-2 scattering. This method supplies
	* the component and allows it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out1 The helicity of the first outgoing particle.
	* @param out2 The helicity of the second outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex & operator () (unsigned int in1,unsigned int in2,
	- unsigned int out1,unsigned int out2) {
	- assert(_outspin.size()==2);
	- unsigned int iloc = in1_constants[1] + in2_constants[2] +
	- out1_constants[3] + out2_constants[4];
	- assert(iloc<_matrixelement.size());
	- return _matrixelement[iloc];
	- }
	+ unsigned int out1,unsigned int out2);

	/**
	* Access the helicity components for a 2-3 scattering. This method supplies
	* the component but does not allow it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out1 The helicity of the first outgoing particle.
	* @param out2 The helicity of the second outgoing particle.
	* @param out3 The helicity of the third outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex operator () (unsigned int in1,unsigned int in2,
	unsigned int out1,unsigned int out2,
	- unsigned int out3) const {
	- assert(_outspin.size()==3);
	- vector<unsigned int> ivec(5);
	- ivec[0]=in1;
	- ivec[1]=in2;
	- ivec[2]=out1;
	- ivec[3]=out2;
	- ivec[4]=out3;
	- return (*this)(ivec);
	- }
	+ unsigned int out3) const;

	/**
	* Access the helicity components for a 2-3 scattering. This method supplies
	* the component and allows it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out1 The helicity of the first outgoing particle.
	* @param out2 The helicity of the second outgoing particle.
	* @param out3 The helicity of the third outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex & operator () (unsigned int in1,unsigned int in2,
	unsigned int out1,unsigned int out2,
	- unsigned int out3) {
	- assert(_outspin.size()==3);
	- vector<unsigned int> ivec(5);
	- ivec[0]=in1;
	- ivec[1]=in2;
	- ivec[2]=out1;
	- ivec[3]=out2;
	- ivec[4]=out3;
	- return (*this)(ivec);
	- }
	+ unsigned int out3);

	/**
	* Access the helicity components for a 2-4 scattering. This method supplies
	* the component but does not allow it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out1 The helicity of the first outgoing particle.
	* @param out2 The helicity of the second outgoing particle.
	* @param out3 The helicity of the third outgoing particle.
	* @param out4 The helicity of the fourth outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex operator () (unsigned int in1,unsigned int in2,
	unsigned int out1,unsigned int out2,
	- unsigned int out3,unsigned int out4) const {
	- assert(_outspin.size()==4);
	- vector<unsigned int> ivec(6);
	- ivec[0]=in1;
	- ivec[1]=in2;
	- ivec[2]=out1;
	- ivec[3]=out2;
	- ivec[4]=out3;
	- ivec[5]=out4;
	- return (*this)(ivec);
	- }
	+ unsigned int out3,unsigned int out4) const;

	/**
	* Access the helicity components for a 2-4 scattering. This method supplies
	* the component and allows it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out1 The helicity of the first outgoing particle.
	* @param out2 The helicity of the second outgoing particle.
	* @param out3 The helicity of the third outgoing particle.
	* @param out4 The helicity of the fourth outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex & operator () (unsigned int in1,unsigned int in2,
	unsigned int out1,unsigned int out2,
	- unsigned int out3, unsigned int out4) {
	- assert(_outspin.size()==4);
	- vector<unsigned int> ivec(6);
	- ivec[0]=in1;
	- ivec[1]=in2;
	- ivec[2]=out1;
	- ivec[3]=out2;
	- ivec[4]=out3;
	- ivec[5]=out4;
	- return (*this)(ivec);
	- }
	+ unsigned int out3, unsigned int out4);

	/**
	* Access the helicity components for a 2-5 scattering. This method supplies
	* the component but does not allow it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out1 The helicity of the first outgoing particle.
	* @param out2 The helicity of the second outgoing particle.
	* @param out3 The helicity of the third outgoing particle.
	* @param out4 The helicity of the fourth outgoing particle.
	* @param out5 The helicity of the fifth outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex operator () (unsigned int in1,unsigned int in2,
	unsigned int out1,unsigned int out2,
	unsigned int out3,unsigned int out4,
	- unsigned int out5) const {
	- assert(_outspin.size()==5);
	- vector<unsigned int> ivec(7);
	- ivec[0]=in1;
	- ivec[1]=in2;
	- ivec[2]=out1;
	- ivec[3]=out2;
	- ivec[4]=out3;
	- ivec[5]=out4;
	- ivec[6]=out5;
	- return (*this)(ivec);
	- }
	+ unsigned int out5) const;

	/**
	* Access the helicity components for a 2-5 scattering. This method supplies
	* the component and allows it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out1 The helicity of the first outgoing particle.
	* @param out2 The helicity of the second outgoing particle.
	* @param out3 The helicity of the third outgoing particle.
	* @param out4 The helicity of the fourth outgoing particle.
	* @param out5 The helicity of the fifth outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex & operator () (unsigned int in1,unsigned int in2,
	unsigned int out1,unsigned int out2,
	unsigned int out3, unsigned int out4,
	- unsigned int out5) {
	- assert(_outspin.size()==5);
	- vector<unsigned int> ivec(7);
	- ivec[0]=in1;
	- ivec[1]=in2;
	- ivec[2]=out1;
	- ivec[3]=out2;
	- ivec[4]=out3;
	- ivec[5]=out4;
	- ivec[6]=out5;
	- return (*this)(ivec);
	- }
	+ unsigned int out5);

	/**
	* Access the helicity components for a 2-6 scattering. This method supplies
	* the component but does not allow it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out1 The helicity of the first outgoing particle.
	* @param out2 The helicity of the second outgoing particle.
	* @param out3 The helicity of the third outgoing particle.
	* @param out4 The helicity of the fourth outgoing particle.
	* @param out5 The helicity of the fifth outgoing particle.
	* @param out6 The helicity of the sixth outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex operator () (unsigned int in1,unsigned int in2,
	unsigned int out1,unsigned int out2,
	unsigned int out3,unsigned int out4,
	- unsigned int out5,unsigned int out6) const {
	- assert(_outspin.size()==6);
	- vector<unsigned int> ivec(8);
	- ivec[0]=in1;
	- ivec[1]=in2;
	- ivec[2]=out1;
	- ivec[3]=out2;
	- ivec[4]=out3;
	- ivec[5]=out4;
	- ivec[6]=out5;
	- ivec[7]=out6;
	- return (*this)(ivec);
	- }
	+ unsigned int out5,unsigned int out6) const;

	/**
	* Access the helicity components for a 2-6 scattering. This method supplies
	* the component and allows it to be changed.
	* @param in1 The helicity of the first incoming particle.
	* @param in2 The helicity of the second incoming particle.
	* @param out1 The helicity of the first outgoing particle.
	* @param out2 The helicity of the second outgoing particle.
	* @param out3 The helicity of the third outgoing particle.
	* @param out4 The helicity of the fourth outgoing particle.
	* @param out5 The helicity of the fifth outgoing particle.
	* @param out6 The helicity of the sixth outgoing particle.
	* @return The matrix element for the given helicities.
	*/
	Complex & operator () (unsigned int in1,unsigned int in2,
	unsigned int out1,unsigned int out2,
	unsigned int out3, unsigned int out4,
	- unsigned int out5, unsigned int out6) {
	- assert(_outspin.size()==6);
	- vector<unsigned int> ivec(8);
	- ivec[0]=in1;
	- ivec[1]=in2;
	- ivec[2]=out1;
	- ivec[3]=out2;
	- ivec[4]=out3;
	- ivec[5]=out4;
	- ivec[6]=out5;
	- ivec[7]=out6;
	- return (*this)(ivec);
	- }
	+ unsigned int out5, unsigned int out6);

	/**
	* Access the helicity components for a 2-n scattering. This method supplies
	* the component but does not allow it to be changed.
	* @param hel The helicities of the incoming and outgoing particles
	* @return The matrix element for the given helicities.
	*/
	- Complex operator () (vector<unsigned int> hel) const {
	- assert(_outspin.size() == hel.size()-2);
	- unsigned int iloc(0),ix;
	- // incoming and outgoing particles
	- for(ix=0;ix<hel.size();++ix)
	- iloc += hel[ix]*_constants[ix+1];
	- assert(iloc<_matrixelement.size());
	- return _matrixelement[iloc];
	- }
	+ Complex operator () (vector<unsigned int> hel) const;

	/**
	* Access the helicity components for a 2-n scattering. This method supplies
	* the component and allows it to be changed.
	* @param hel The helicities of the incoming and outgoing particles
	* @return The matrix element for the given helicities.
	*/
	- Complex & operator () (vector<unsigned int> hel) {
	- assert(_outspin.size() == hel.size()-2);
	- unsigned int iloc=0,ix;
	- // incoming particles
	- for(ix=0;ix<hel.size();++ix)
	- iloc += hel[ix]*_constants[ix+1];
	- assert(iloc<_matrixelement.size());
	- return _matrixelement[iloc];
	- }
	+ Complex & operator () (vector<unsigned int> hel);
	//@}

	public:

	/**
	* Calculate the decay matrix for an incoming particle.
	*/
	RhoDMatrix calculateDMatrix(int,const RhoDMatrix &,
	const vector<RhoDMatrix> &) const;

	/**
	* Calculate the rho matrix for a given outgoing particle.
	*/
	RhoDMatrix calculateRhoMatrix(int,const RhoDMatrix &,
	const RhoDMatrix &,
	const vector<RhoDMatrix> &) const;

	/**
	* Compute the spin averaged matrix element
	*/
	double average() const;

	/**
	* Compute the spin average matrix element
	*/
	double average(const RhoDMatrix & in1,
	const RhoDMatrix & in2) const;

	/**
	* Compute the spin average matrix element
	*/
	Complex average(const ProductionMatrixElement & me2,
	const RhoDMatrix & in1,
	const RhoDMatrix & in2) const;

	public:

	/**
	* Reset the matrix element.
	*/
	- void reset(const ProductionMatrixElement & x) const {
	- _nout = x._nout;
	- _inspin = x._inspin;
	- _outspin = x._outspin;
	- _matrixelement = x._matrixelement;
	- _constants = x._constants;
	- }
	+ void reset(const ProductionMatrixElement & x) const;

	/**
	* Standard Init function used to initialize the interfaces.
	*/
	static void Init();

	private:

	/**
	* Set the size of the vector containing the matrix element.
	*/
	- void setMESize() {
	- unsigned int ix;
	- int isize=_inspin[0]*_inspin[1];
	- for(ix=0;ix<_outspin.size();++ix)
	- isize*=_outspin[ix];
	- // zero the matrix element
	- _matrixelement.resize(isize,0.);
	- // set up the constants for the mapping of helicity to vectro index
	- _constants.resize(_outspin.size()+3);
	- unsigned int temp=1;
	- for(ix=_outspin.size()+1;ix>1;--ix) {
	- temp*=_outspin[ix-2];
	- _constants[ix]=temp;
	- }
	- temp*=_inspin[1];_constants[1]=temp;
	- temp*=_inspin[0];_constants[0]=temp;
	- _constants[_outspin.size()+2]=1;
	- }
	+ void setMESize();

	private:

	/**
	* Number of outgoing particles.
	*/
	mutable unsigned int _nout;

	/**
	* Spin of the incoming particles as 2s+1.
	*/
	mutable vector<PDT::Spin> _inspin;

	/**
	* Spins of the outgoing particles.
	*/
	mutable vector<PDT::Spin> _outspin;

	/**
	* Storage of the matrix element, a vector is better for memory usage.
	*/
	mutable vector<Complex> _matrixelement;

	/**
	* Constants needed to map the index of the vector to a helicity structure.
	*/
	mutable vector<int> _constants;

	};

	}


	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/**
	* The following template specialization informs ThePEG about the
	* base class of ProductionMatrixElement.
	*/
	template <>
	struct BaseClassTrait<Herwig::ProductionMatrixElement,1> {
	/** Typedef of the base class of ProductionMatrixElement. */
	typedef Base NthBase;
	};

	/**
	* The following template specialization informs ThePEG about the
	* name of this class and the shared object where it is defined.
	*/
	template <>
	struct ClassTraits<Herwig::ProductionMatrixElement>
	: public ClassTraitsBase<Herwig::ProductionMatrixElement> {

	/**
	* Return the class name.
	*/
	static string className() { return "Herwig::ProductionMatrixElement"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_ProductionMatrixElement_H */
	diff --git a/Shower/ShowerHandler.cc b/Shower/ShowerHandler.cc
	--- a/Shower/ShowerHandler.cc
	+++ b/Shower/ShowerHandler.cc
	@@ -1,733 +1,734 @@
	// -- C++ --
	//
	// ShowerHandler.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ShowerHandler class.
	//

	#include "ShowerHandler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "ThePEG/PDF/PartonBinInstance.h"
	#include "Herwig++/PDT/StandardMatchers.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "Herwig++/Shower/Base/Evolver.h"
	#include "Herwig++/Shower/Base/ShowerParticle.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig++/Utilities/EnumParticles.h"
	#include "Herwig++/PDF/MPIPDF.h"
	#include "Herwig++/PDF/MinBiasPDF.h"
	#include "ThePEG/Handlers/EventHandler.h"
	#include "Herwig++/Shower/Base/ShowerTree.h"
	#include "Herwig++/Shower/Base/HardTree.h"
	#include "Herwig++/Shower/Base/KinematicsReconstructor.h"
	#include "Herwig++/Shower/Base/PartnerFinder.h"
	#include "Herwig++/PDF/HwRemDecayer.h"
	#include <cassert>

	using namespace Herwig;

	ShowerHandler::~ShowerHandler() {}

	ShowerHandler * ShowerHandler::currentHandler_ = 0;

	void ShowerHandler::doinit() {
	CascadeHandler::doinit();
	// copy particles to decay before showering from input vector to the
	// set used in the simulation
	particlesDecayInShower_.insert(inputparticlesDecayInShower_.begin(),
	inputparticlesDecayInShower_.end());

	// \todo DG: Disabled here because of momentum-violation problems
	// ShowerTree::_decayInShower = particlesDecayInShower_;
	}

	IBPtr ShowerHandler::clone() const {
	return new_ptr(*this);
	}

	IBPtr ShowerHandler::fullclone() const {
	return new_ptr(*this);
	}

	ShowerHandler::ShowerHandler() :
	pdfFreezingScale_(2.5*GeV),
	maxtry_(10),maxtryMPI_(10),maxtryDP_(10), subProcess_() {
	inputparticlesDecayInShower_.push_back( 6 ); // top
	inputparticlesDecayInShower_.push_back( 23 ); // Z0
	inputparticlesDecayInShower_.push_back( 24 ); // W+/-
	inputparticlesDecayInShower_.push_back( 25 ); // h0
	}

	void ShowerHandler::doinitrun(){
	CascadeHandler::doinitrun();
	//can't use isMPIOn here, because the EventHandler is not set at that stage
	if(MPIHandler_){
	MPIHandler_->initialize();
	if(MPIHandler_->softInt())
	remDec_->initSoftInteractions(MPIHandler_->Ptmin(), MPIHandler_->beta());
	}
	// \todo DG: Disabled here because of momentum-violation problems
	// Must set pre-cascade-handler to NewPhysics/DecayHandler instead
	// ShowerTree::_decayInShower = particlesDecayInShower_;
	}

	void ShowerHandler::dofinish(){
	CascadeHandler::dofinish();
	if(MPIHandler_) MPIHandler_->finalize();
	}

	void ShowerHandler::persistentOutput(PersistentOStream & os) const {
	os << evolver_ << remDec_ << ounit(pdfFreezingScale_,GeV) << maxtry_
	<< maxtryMPI_ << maxtryDP_ << inputparticlesDecayInShower_
	<< particlesDecayInShower_ << MPIHandler_ << PDFA_ << PDFB_;
	}

	void ShowerHandler::persistentInput(PersistentIStream & is, int) {
	is >> evolver_ >> remDec_ >> iunit(pdfFreezingScale_,GeV) >> maxtry_
	>> maxtryMPI_ >> maxtryDP_ >> inputparticlesDecayInShower_
	>> particlesDecayInShower_ >> MPIHandler_ >> PDFA_ >> PDFB_;
	}

	ClassDescription<ShowerHandler> ShowerHandler::initShowerHandler;
	// Definition of the static class description member.

	void ShowerHandler::Init() {

	static ClassDocumentation<ShowerHandler> documentation
	("Main driver class for the showering.",
	"The Shower evolution was performed using an algorithm described in "
	"\\cite{Marchesini:1983bm,Marchesini:1987cf,Gieseke:2003rz,Bahr:2008pv}.",
	"%\\cite{Marchesini:1983bm}\n"
	"\\bibitem{Marchesini:1983bm}\n"
	" G.~Marchesini and B.~R.~Webber,\n"
	" ``Simulation Of QCD Jets Including Soft Gluon Interference,''\n"
	" Nucl.\\ Phys.\\ B {\\bf 238}, 1 (1984).\n"
	" %%CITATION = NUPHA,B238,1;%%\n"
	"%\\cite{Marchesini:1987cf}\n"
	"\\bibitem{Marchesini:1987cf}\n"
	" G.~Marchesini and B.~R.~Webber,\n"
	" ``Monte Carlo Simulation of General Hard Processes with Coherent QCD\n"
	" Radiation,''\n"
	" Nucl.\\ Phys.\\ B {\\bf 310}, 461 (1988).\n"
	" %%CITATION = NUPHA,B310,461;%%\n"
	"%\\cite{Gieseke:2003rz}\n"
	"\\bibitem{Gieseke:2003rz}\n"
	" S.~Gieseke, P.~Stephens and B.~Webber,\n"
	" ``New formalism for QCD parton showers,''\n"
	" JHEP {\\bf 0312}, 045 (2003)\n"
	" [arXiv:hep-ph/0310083].\n"
	" %%CITATION = JHEPA,0312,045;%%\n"
	);

	static Reference<ShowerHandler,Evolver>
	interfaceEvolver("Evolver",
	"A reference to the Evolver object",
	&Herwig::ShowerHandler::evolver_,
	false, false, true, false);

	static Reference<ShowerHandler,HwRemDecayer>
	interfaceRemDecayer("RemDecayer",
	"A reference to the Remnant Decayer object",
	&Herwig::ShowerHandler::remDec_,
	false, false, true, false);

	static Parameter<ShowerHandler,Energy> interfacePDFFreezingScale
	("PDFFreezingScale",
	"The PDF freezing scale",
	&ShowerHandler::pdfFreezingScale_, GeV, 2.5GeV, 2.0GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<ShowerHandler,unsigned int> interfaceMaxTry
	("MaxTry",
	"The maximum number of attempts for the main showering loop",
	&ShowerHandler::maxtry_, 10, 1, 100,
	false, false, Interface::limited);

	static Parameter<ShowerHandler,unsigned int> interfaceMaxTryMPI
	("MaxTryMPI",
	"The maximum number of regeneration attempts for an additional scattering",
	&ShowerHandler::maxtryMPI_, 10, 0, 100,
	false, false, Interface::limited);

	static Parameter<ShowerHandler,unsigned int> interfaceMaxTryDP
	("MaxTryDP",
	"The maximum number of regeneration attempts for an additional hard scattering",
	&ShowerHandler::maxtryDP_, 10, 0, 100,
	false, false, Interface::limited);

	static ParVector<ShowerHandler,long> interfaceDecayInShower
	("DecayInShower",
	"PDG codes of the particles to be decayed in the shower",
	&ShowerHandler::inputparticlesDecayInShower_, -1, 0l, -10000000l, 10000000l,
	false, false, Interface::limited);

	static Reference<ShowerHandler,UEBase> interfaceMPIHandler
	("MPIHandler",
	"The object that administers all additional scatterings.",
	&ShowerHandler::MPIHandler_, false, false, true, true);

	static Reference<ShowerHandler,PDFBase> interfacePDFA
	("PDFA",
	"The PDF for beam particle A. Overrides the particle's own PDF setting.",
	&ShowerHandler::PDFA_, false, false, true, true, false);

	static Reference<ShowerHandler,PDFBase> interfacePDFB
	("PDFB",
	"The PDF for beam particle B. Overrides the particle's own PDF setting.",
	&ShowerHandler::PDFB_, false, false, true, true, false);

	}

	void ShowerHandler::cascade() {

	tcPDFPtr first = firstPDF().pdf();
	tcPDFPtr second = secondPDF().pdf();

	if ( PDFA_ ) first = PDFA_;
	if ( PDFB_ ) second = PDFB_;

	resetPDFs(make_pair(first,second));

	// get the incoming partons
	tPPair incomingPartons =
	eventHandler()->currentCollision()->primarySubProcess()->incoming();
	// and the parton bins
	PBIPair incomingBins =
	make_pair(lastExtractor()->partonBinInstance(incomingPartons.first),
	lastExtractor()->partonBinInstance(incomingPartons.second));
	// and the incoming hadrons
	tPPair incomingHadrons =
	eventHandler()->currentCollision()->incoming();
	remDec_->setHadronContent(incomingHadrons);
	// check if incoming hadron == incoming parton
	// and get the incoming hadron if exists or parton otherwise
	incoming_ = make_pair(incomingBins.first ?
	incomingBins.first ->particle() : incomingPartons.first,
	incomingBins.second ?
	incomingBins.second->particle() : incomingPartons.second);
	// check the collision is of the beam particles
	// and if not boost collision to the right frame
	// i.e. the hadron-hadron CMF of the collision
	bool btotal(false);
	LorentzRotation rtotal;
	if(incoming_.first != incomingHadrons.first \|\|
	incoming_.second != incomingHadrons.second ) {
	btotal = true;
	boostCollision(false);
	}
	// set the current ShowerHandler
	currentHandler_ = this;
	// first shower the hard process
	useMe();
	try {
	SubProPtr sub = eventHandler()->currentCollision()->primarySubProcess();
	incomingPartons = cascade(sub,lastXCombPtr());
	}
	catch(ShowerTriesVeto &veto){
	throw Exception() << "Failed to generate the shower after "
	<< veto.tries
	<< " attempts in ShowerHandler::cascade()"
	<< Exception::eventerror;
	}
	if(showerHardProcessVeto()) throw Veto();
	// if a non-hadron collision return (both incoming non-hadronic)
	if( ( !incomingBins.first\|\|
	!isResolvedHadron(incomingBins.first ->particle()))&&
	( !incomingBins.second\|\|
	!isResolvedHadron(incomingBins.second->particle()))) {
	// boost back to lab if needed
	if(btotal) boostCollision(true);
	// unset the current ShowerHandler
	currentHandler_ = 0;
	return;
	}
	// get the remnants for hadronic collision
	pair<tRemPPtr,tRemPPtr> remnants(getRemnants(incomingBins));
	// set the starting scale of the forced splitting to the PDF freezing scale
	remDec_->initialize(remnants, incoming_, *currentStep(), pdfFreezingScale());
	// do the first forcedSplitting
	try {
	remDec_->doSplit(incomingPartons, make_pair(firstPDF() .pdf(),
	secondPDF().pdf()), true);
	}
	catch (ExtraScatterVeto) {
	throw Exception() << "Remnant extraction failed in "
	<< "ShowerHandler::cascade() from primary interaction"
	<< Exception::eventerror;
	}
	// if no MPI return
	if( !isMPIOn() ) {
	remDec_->finalize();
	// boost back to lab if needed
	if(btotal) boostCollision(true);
	// unset the current ShowerHandler
	currentHandler_ = 0;
	return;
	}
	// generate the multiple scatters use modified pdf's now:
	setMPIPDFs();
	// additional "hard" processes
	unsigned int tries(0);
	// This is the loop over additional hard scatters (most of the time
	// only one, but who knows...)
	for(unsigned int i=1; i <= getMPIHandler()->additionalHardProcs(); i++){
	//counter for regeneration
	unsigned int multSecond = 0;
	// generate the additional scatters
	while( multSecond < getMPIHandler()->multiplicity(i) ) {
	// generate the hard scatter
	tStdXCombPtr lastXC = getMPIHandler()->generate(i);
	SubProPtr sub = lastXC->construct();
	// add to the Step
	newStep()->addSubProcess(sub);
	// increment the counters
	tries++;
	multSecond++;
	if(tries == maxtryDP_)
	throw Exception() << "Failed to establish the requested number "
	<< "of additional hard processes. If this error "
	<< "occurs often, your selection of additional "
	<< "scatter is probably unphysical"
	<< Exception::eventerror;
	// Generate the shower. If not possible veto the event
	try {
	incomingPartons = cascade(sub,lastXC);
	}
	catch(ShowerTriesVeto &veto){
	throw Exception() << "Failed to generate the shower of "
	<< "a secondary hard process after "
	<< veto.tries
	<< " attempts in Evolver::showerHardProcess()"
	<< Exception::eventerror;
	}
	try {
	// do the forcedSplitting
	remDec_->doSplit(incomingPartons, make_pair(firstPDF().pdf(),
	secondPDF().pdf()), false);
	}
	catch(ExtraScatterVeto){
	//remove all particles associated with the subprocess
	newStep()->removeParticle(incomingPartons.first);
	newStep()->removeParticle(incomingPartons.second);
	//remove the subprocess from the list
	newStep()->removeSubProcess(sub);
	//regenerate the scattering
	multSecond--;
	continue;
	}
	// connect with the remnants but don't set Remnant colour,
	// because that causes problems due to the multiple colour lines.
	if ( !remnants.first ->extract(incomingPartons.first , false) \|\|
	!remnants.second->extract(incomingPartons.second, false) )
	throw Exception() << "Remnant extraction failed in "
	<< "ShowerHandler::cascade() for additional scatter"
	<< Exception::runerror;
	}
	}
	// the underlying event processes
	unsigned int ptveto(1), veto(0);
	unsigned int max(getMPIHandler()->multiplicity());
	for(unsigned int i=0; i<max; i++) {
	// check how often this scattering has been regenerated
	if(veto > maxtryMPI_) break;
	//generate PSpoint
	tStdXCombPtr lastXC = getMPIHandler()->generate();
	SubProPtr sub = lastXC->construct();
	//If Algorithm=1 additional scatters of the signal type
	// with pt > ptmin have to be vetoed
	//with probability 1/(m+1), where m is the number of occurances in this event
	if( getMPIHandler()->Algorithm() == 1 ){
	//get the pT
	Energy pt = sub->outgoing().front()->momentum().perp();
	if(pt > getMPIHandler()->PtForVeto() && UseRandom::rnd() < 1./(ptveto+1) ){
	ptveto++;
	i--;
	continue;
	}
	}
	// add to the SubProcess to the step
	newStep()->addSubProcess(sub);
	// Run the Shower. If not possible veto the scattering
	try {
	incomingPartons = cascade(sub,lastXC);
	}
	// discard this extra scattering, but try the next one
	catch(ShowerTriesVeto) {
	newStep()->removeSubProcess(sub);
	//regenerate the scattering
	veto++;
	i--;
	continue;
	}
	try{
	//do the forcedSplitting
	remDec_->doSplit(incomingPartons, make_pair(firstPDF().pdf(),
	secondPDF().pdf()), false);
	}
	catch (ExtraScatterVeto) {
	//remove all particles associated with the subprocess
	newStep()->removeParticle(incomingPartons.first);
	newStep()->removeParticle(incomingPartons.second);
	//remove the subprocess from the list
	newStep()->removeSubProcess(sub);
	//regenerate the scattering
	veto++;
	i--;
	continue;
	}
	//connect with the remnants but don't set Remnant colour,
	//because that causes problems due to the multiple colour lines.
	if ( !remnants.first ->extract(incomingPartons.first , false) \|\|
	!remnants.second->extract(incomingPartons.second, false) )
	throw Exception() << "Remnant extraction failed in "
	<< "ShowerHandler::cascade() for MPI hard scattering"
	<< Exception::runerror;
	//reset veto counter
	veto = 0;
	}
	// finalize the remnants
	remDec_->finalize(getMPIHandler()->colourDisrupt(),
	getMPIHandler()->softMultiplicity());
	// boost back to lab if needed
	if(btotal) boostCollision(true);
	// unset the current ShowerHandler
	currentHandler_ = 0;
	+ getMPIHandler()->clean();
	}

	void ShowerHandler::fillEventRecord() {
	// create a new step
	StepPtr pstep = newStep();
	assert(!done_.empty());
	assert(done_[0]->isHard());
	// insert the steps
	for(unsigned int ix=0;ix<done_.size();++ix) {
	done_[ix]->fillEventRecord(pstep,
	evolver_->isISRadiationON(),
	evolver_->isFSRadiationON());
	}
	}

	void ShowerHandler::findShoweringParticles() {
	// clear the storage
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	// temporary storage of the particles
	set<PPtr> hardParticles;
	// outgoing particles from the hard process
	PVector outgoing = currentSubProcess()->outgoing();
	set<PPtr> outgoingset(outgoing.begin(),outgoing.end());
	// loop over the tagged particles
	tPVector thetagged;
	if( firstInteraction() ){
	thetagged = tagged();
	}
	else{
	thetagged.insert(thetagged.end(),
	outgoing.begin(),outgoing.end());
	}
	bool isHard=false;
	for (tParticleVector::const_iterator
	taggedP = thetagged.begin();
	taggedP != thetagged.end(); ++taggedP) {
	// if a remnant don't consider
	if(eventHandler()->currentCollision()->isRemnant(*taggedP))
	continue;
	// find the parent and whether its a colourless s-channel resonance
	bool isDecayProd=false;
	tPPtr parent;
	if(!(*taggedP)->parents().empty()) {
	parent = (*taggedP)->parents()[0];
	// check if from s channel decaying colourless particle
	isDecayProd = decayProduct(parent);
	}
	// add to list of outgoing hard particles if needed
	isHard \|=(outgoingset.find(*taggedP) != outgoingset.end());
	if(isDecayProd) hardParticles.insert(findParent(parent,isHard,outgoingset));
	else hardParticles.insert(*taggedP);
	}
	// there must be something to shower
	if(hardParticles.empty())
	throw Exception() << "No particles to shower in "
	<< "ShowerHandler::fillShoweringParticles"
	<< Exception::eventerror;
	if(!isHard)
	throw Exception() << "Starting on decay not yet implemented in "
	<< "ShowerHandler::findShoweringParticles()"
	<< Exception::runerror;
	// create the hard process ShowerTree
	ParticleVector out(hardParticles.begin(),hardParticles.end());
	hard_=new_ptr(ShowerTree(currentSubProcess()->incoming(),out, decay_));
	hard_->setParents();
	}

	void ShowerHandler::prepareCascade(tSubProPtr sub) {
	current_ = currentStep();
	subProcess_ = sub;
	}

	tPPair ShowerHandler::cascade(tSubProPtr sub,
	XCPtr xcomb) {
	prepareCascade(sub);
	// start of the try block for the whole showering process
	unsigned int countFailures=0;
	while (countFailures<maxtry_) {
	try {
	// find the particles in the hard process and the decayed particles to shower
	findShoweringParticles();
	// if no hard process
	if(!hard_) throw Exception() << "Shower starting with a decay"
	<< "is not implemented"
	<< Exception::runerror;
	// perform the shower for the hard process
	evolver_->showerHardProcess(hard_,xcomb);
	done_.push_back(hard_);
	hard_->updateAfterShower(decay_);
	// if no decaying particles to shower break out of the loop
	if(decay_.empty()) break;
	// shower the decay products
	while(!decay_.empty()) {
	// find particle whose production process has been showered
	ShowerDecayMap::iterator dit = decay_.begin();
	while(!dit->second->parent()->hasShowered() && dit!=decay_.end()) ++dit;
	assert(dit!=decay_.end());
	// get the particle
	ShowerTreePtr decayingTree = dit->second;
	// remove it from the multimap
	decay_.erase(dit);
	// make sure the particle has been decayed
	decayingTree->decay(decay_);
	// now shower the decay
	evolver_->showerDecay(decayingTree);
	done_.push_back(decayingTree);
	decayingTree->updateAfterShower(decay_);
	}
	// suceeded break out of the loop
	break;
	}
	catch (KinematicsReconstructionVeto) {
	++countFailures;
	}
	}
	// if loop exited because of too many tries, throw event away
	if (countFailures >= maxtry_) {
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	throw Exception() << "Too many tries for main while loop "
	<< "in ShowerHandler::cascade()."
	<< Exception::eventerror;
	}
	//enter the particles in the event record
	fillEventRecord();
	// clear storage
	hard_=ShowerTreePtr();
	decay_.clear();
	done_.clear();
	// non hadronic case return
	if (!isResolvedHadron(incoming_.first ) &&
	!isResolvedHadron(incoming_.second) )
	return incoming_;
	// remake the remnants (needs to be after the colours are sorted
	// out in the insertion into the event record)
	if ( firstInteraction() ) return remakeRemnant(sub->incoming());
	//Return the new pair of incoming partons. remakeRemnant is not
	//necessary here, because the secondary interactions are not yet
	//connected to the remnants.
	return make_pair(findFirstParton(sub->incoming().first ),
	findFirstParton(sub->incoming().second));
	}

	PPtr ShowerHandler::findParent(PPtr original, bool & isHard,
	set<PPtr> outgoingset) const {
	PPtr parent=original;
	isHard \|=(outgoingset.find(original) != outgoingset.end());
	if(!original->parents().empty()) {
	PPtr orig=original->parents()[0];
	if(current_->find(orig)&&decayProduct(orig)) {
	parent=findParent(orig,isHard,outgoingset);
	}
	}
	return parent;
	}

	ShowerHandler::RemPair
	ShowerHandler::getRemnants(PBIPair incomingBins) {
	RemPair remnants;
	// first beam particle
	if(incomingBins.first&&!incomingBins.first->remnants().empty()) {
	remnants.first =
	dynamic_ptr_cast<tRemPPtr>(incomingBins.first->remnants()[0] );
	if(remnants.first) {
	ParticleVector children=remnants.first->children();
	for(unsigned int ix=0;ix<children.size();++ix) {
	if(children[ix]->dataPtr()==remnants.first->dataPtr())
	remnants.first = dynamic_ptr_cast<RemPPtr>(children[ix]);
	}
	//remove existing colour lines from the remnants
	if(remnants.first->colourLine())
	remnants.first->colourLine()->removeColoured(remnants.first);
	if(remnants.first->antiColourLine())
	remnants.first->antiColourLine()->removeAntiColoured(remnants.first);
	}
	}
	// seconnd beam particle
	if(incomingBins.second&&!incomingBins. second->remnants().empty()) {
	remnants.second =
	dynamic_ptr_cast<tRemPPtr>(incomingBins.second->remnants()[0] );
	if(remnants.second) {
	ParticleVector children=remnants.second->children();
	for(unsigned int ix=0;ix<children.size();++ix) {
	if(children[ix]->dataPtr()==remnants.second->dataPtr())
	remnants.second = dynamic_ptr_cast<RemPPtr>(children[ix]);
	}
	//remove existing colour lines from the remnants
	if(remnants.second->colourLine())
	remnants.second->colourLine()->removeColoured(remnants.second);
	if(remnants.second->antiColourLine())
	remnants.second->antiColourLine()->removeAntiColoured(remnants.second);
	}
	}
	assert(remnants.first \|\| remnants.second);
	return remnants;
	}

	tPPair ShowerHandler::remakeRemnant(tPPair oldp){
	// get the parton extractor
	PartonExtractor & pex = *lastExtractor();
	// get the new partons
	tPPair newp = make_pair(findFirstParton(oldp.first ),
	findFirstParton(oldp.second));
	// if the same do nothing
	if(newp == oldp) return oldp;
	// Creates the new remnants and returns the new PartonBinInstances
	PBIPair newbins = pex.newRemnants(oldp, newp, newStep());
	newStep()->addIntermediate(newp.first);
	newStep()->addIntermediate(newp.second);
	// return the new partons
	return newp;
	}

	PPtr ShowerHandler::findFirstParton(tPPtr seed) const{
	if(seed->parents().empty()) return seed;
	tPPtr parent = seed->parents()[0];
	//if no parent there this is a loose end which will
	//be connected to the remnant soon.
	if(!parent \|\| parent == incoming_.first \|\|
	parent == incoming_.second ) return seed;
	else return findFirstParton(parent);
	}

	bool ShowerHandler::decayProduct(tPPtr particle) const {
	// must be time-like and not incoming
	if(particle->momentum().m2()<=ZERO\|\|
	particle == currentSubProcess()->incoming().first\|\|
	particle == currentSubProcess()->incoming().second) return false;
	// if only 1 outgoing and this is it
	if(currentSubProcess()->outgoing().size()==1 &&
	currentSubProcess()->outgoing()[0]==particle) return true;
	// must not be the s-channel intermediate otherwise
	if(find(currentSubProcess()->incoming().first->children().begin(),
	currentSubProcess()->incoming().first->children().end(),particle)!=
	currentSubProcess()->incoming().first->children().end()&&
	find(currentSubProcess()->incoming().second->children().begin(),
	currentSubProcess()->incoming().second->children().end(),particle)!=
	currentSubProcess()->incoming().second->children().end()&&
	currentSubProcess()->incoming().first ->children().size()==1&&
	currentSubProcess()->incoming().second->children().size()==1)
	return false;
	// if non-coloured this is enough
	if(!particle->dataPtr()->coloured()) return true;
	// if coloured must be unstable
	if(particle->dataPtr()->stable()) return false;
	// must not have same particle type as a child
	int id = particle->id();
	for(unsigned int ix=0;ix<particle->children().size();++ix)
	if(particle->children()[ix]->id()==id) return false;
	// otherwise its a decaying particle
	return true;
	}

	namespace {

	void addChildren(tPPtr in,set<tPPtr> & particles) {
	particles.insert(in);
	for(unsigned int ix=0;ix<in->children().size();++ix)
	addChildren(in->children()[ix],particles);
	}
	}

	void ShowerHandler::boostCollision(bool boost) {
	// calculate boost from lab to rest
	if(!boost) {
	Lorentz5Momentum ptotal=incoming_.first ->momentum()+incoming_.second->momentum();
	boost_ = LorentzRotation(-ptotal.boostVector());
	Axis axis((boost_*incoming_.first ->momentum()).vect().unit());
	if(axis.perp2()>0.) {
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	boost_.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	}
	}
	// first call performs the boost and second inverse
	// get the particles to be boosted
	set<tPPtr> particles;
	addChildren(incoming_.first,particles);
	addChildren(incoming_.second,particles);
	// apply the boost
	for(set<tPPtr>::const_iterator cit=particles.begin();
	cit!=particles.end();++cit) {
	(*cit)->transform(boost_);
	}
	if(!boost) boost_.invert();
	}

	void ShowerHandler::setMPIPDFs() {

	if ( !mpipdfs_.first ) {
	// first have to check for MinBiasPDF
	tcMinBiasPDFPtr first = dynamic_ptr_cast<tcMinBiasPDFPtr>(firstPDF().pdf());
	if(first)
	mpipdfs_.first = new_ptr(MPIPDF(first->originalPDF()));
	else
	mpipdfs_.first = new_ptr(MPIPDF(firstPDF().pdf()));
	}

	if ( !mpipdfs_.second ) {
	tcMinBiasPDFPtr second = dynamic_ptr_cast<tcMinBiasPDFPtr>(secondPDF().pdf());
	if(second)
	mpipdfs_.second = new_ptr(MPIPDF(second->originalPDF()));
	else
	mpipdfs_.second = new_ptr(MPIPDF(secondPDF().pdf()));
	}

	// reset the PDFs stored in the base class
	resetPDFs(mpipdfs_);

	}

	bool ShowerHandler::isResolvedHadron(tPPtr particle) {
	if(!HadronMatcher::Check(particle->data())) return false;
	for(unsigned int ix=0;ix<particle->children().size();++ix) {
	if(particle->children()[ix]->id()==ParticleID::Remnant) return true;
	}
	return false;
	}

	HardTreePtr ShowerHandler::generateCKKW(ShowerTreePtr ) const {
	return HardTreePtr();
	}
	diff --git a/Shower/UEBase.h b/Shower/UEBase.h
	--- a/Shower/UEBase.h
	+++ b/Shower/UEBase.h
	@@ -1,167 +1,172 @@
	// -- C++ --
	#ifndef HERWIG_UEBase_H
	#define HERWIG_UEBase_H
	//
	// This is the declaration of the UEBase class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "ThePEG/Handlers/StandardXComb.fh"
	#include "UEBase.fh"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* Abstract base class used to minimize the dependence between
	* MPIHandler and all Shower classes.
	*
	* \author Manuel B\"ahr
	*
	* @see \ref UEBaseInterfaces "The interfaces"
	* defined for UEBase.
	*/
	class UEBase: public Interfaced {

	public:

	/**
	* 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();

	/** @name Virtual functions used for the generation of additional
	interactions . */
	//@{
	/**
	* Some initialization code eventually.
	*/
	virtual void initialize() {}

	/**
	* Return true or false depending on the generator setup.
	*/
	virtual bool beamOK() const = 0;

	/**
	* Return true or false depending on whether soft interactions are enabled.
	*/
	virtual bool softInt() const {return false;}

	/**
	* Return the value of the pt cutoff.
	*/
	virtual Energy Ptmin() const = 0;

	/**
	* Return the slope of the soft pt spectrum. Only necessary when the
	* soft part is modelled.
	*/
	virtual InvEnergy2 beta() const {return ZERO;}

	/**
	* Some finalize code eventually.
	*/
	virtual void finalize() {}

	/**
	+ * Clean up method called after each event.
	+ */
	+ virtual void clean() {}
	+
	+ /**
	* Return the number of different hard processes. Use 0 as default to
	* not require implementation.
	*/
	virtual unsigned int additionalHardProcs() const {return 0;}

	/**
	* return the hard multiplicity of process i. Can't be constant in my
	* case because drawing from the probability distribution also
	* specifies the soft multiplicity that has to be stored....
	*/
	virtual unsigned int multiplicity(unsigned int i=0) = 0;

	/**
	* Generate a additional interaction for ProcessHandler sel. Method
	* can't be const because it saves the state of the underlying XComb
	* object on it's way.
	*/
	virtual tStdXCombPtr generate(unsigned int sel=0) = 0;

	/**
	* Return the type of algorithm.
	*/
	virtual int Algorithm() const = 0;

	/**
	* Return the value of the hard Process pt cutoff for vetoing.
	*/
	virtual Energy PtForVeto() const = 0;

	/**
	* Return the fraction of colour disrupted subprocesses. Use default 0
	* so that it is not required to implement.
	*/
	virtual double colourDisrupt() const {return 0.0;}

	/**
	* Return the soft multiplicity. Use 0 as default to not require
	* implementation.
	*/
	virtual unsigned int softMultiplicity() const {return 0;}
	//@}

	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is an abstract class with persistent data.
	*/
	static AbstractClassDescription<UEBase> initUEBase;

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	UEBase & operator=(const UEBase &);

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of UEBase. */
	template <>
	struct BaseClassTrait<Herwig::UEBase,1> {
	/** Typedef of the first base class of UEBase. */
	typedef Interfaced NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the UEBase class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::UEBase>
	: public ClassTraitsBase<Herwig::UEBase> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::UEBase"; }
	/**
	* The name of a file containing the dynamic library where the class
	* UEBase is implemented. It may also include several, space-separated,
	* libraries if the class UEBase depends on other classes (base classes
	* excepted). In this case the listed libraries will be dynamically
	* linked in the order they are specified.
	*/
	static string library() { return "HwShower.so"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_UEBase_H */
	diff --git a/UnderlyingEvent/MPIHandler.cc b/UnderlyingEvent/MPIHandler.cc
	--- a/UnderlyingEvent/MPIHandler.cc
	+++ b/UnderlyingEvent/MPIHandler.cc
	@@ -1,815 +1,825 @@
	// -- C++ --
	//
	// MPIHandler.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the MPIHandler class.
	//

	#include "MPIHandler.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/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
	if( dynamic_ptr_cast<Ptr<StandardEventHandler>::pointer>(eventHandler()) ){
	PtOfQCDProc_ = dynamic_ptr_cast
	<Ptr<StandardEventHandler>::pointer>(eventHandler())->cuts()->minKT(gluon);
	}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, 10millibarn, 0.6GeV2);
	softXSec_ = res.first;
	softMu2_ = res.second;
	*/
	slopeBisection fs(this);
	try{
	softMu2_ = rootFinder.value(fs, 0.3GeV2, 1.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, 0millibarn, 5000millibarn);
	}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())/(2Ptmin_eps);

	betaBisection fn2(softXSec_, hardPlus, Ptmin_);
	try{
	beta_ = rootFinder.value(fn2, -10/GeV2, 2/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, iH100+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, 0millibarn, 5000millibarn);
	}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 {

	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(mub, 3)(gsl_sf_bessel_Kn(3, mub));
	}

	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.8millibarn 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) +
	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
	+ 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
	+ 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();
	+ }
	+}
	+
	+
	ClassDescription<MPIHandler> MPIHandler::initMPIHandler;
	// Definition of the static class description member.

	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.0GeV2, 0.2GeV2, 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.0millibarn, 0.0millibarn, 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.0GeV, 0.0GeV, 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.11GeV, 0.0GeV, 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/UnderlyingEvent/MPIHandler.h b/UnderlyingEvent/MPIHandler.h
	--- a/UnderlyingEvent/MPIHandler.h
	+++ b/UnderlyingEvent/MPIHandler.h
	@@ -1,914 +1,920 @@
	// -- C++ --
	//
	// MPIHandler.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_MPIHandler_H
	#define HERWIG_MPIHandler_H
	//
	// This is the declaration of the MPIHandler class.
	//
	#include "ThePEG/Interface/Interfaced.h"
	#include "ThePEG/Handlers/StandardEventHandler.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig++/PDT/StandardMatchers.h"
	#include "Herwig++/Utilities/GSLBisection.h"
	//#include "Herwig++/Utilities/GSLMultiRoot.h"
	#include "Herwig++/Utilities/GSLIntegrator.h"
	#include "Herwig++/Shower/UEBase.h"

	#include <cassert>
	#include "ProcessHandler.h"
	#include "MPIHandler.fh"


	namespace Herwig {
	using namespace ThePEG;

	/** \ingroup UnderlyingEvent
	* \class MPIHandler
	* This class is responsible for generating additional
	* semi hard partonic interactions.
	*
	* \author Manuel B\"ahr
	*
	* @see \ref MPIHandlerInterfaces "The interfaces"
	* defined for MPIHandler.
	* @see ProcessHandler
	* @see ShowerHandler
	* @see HwRemDecayer
	*/

	class MPIHandler: public UEBase {

	/**
	* Maximum number of scatters
	*/
	static const unsigned int maxScatters_ = 99;

	/**
	* Class for the integration is a friend to access private members
	*/
	friend struct Eikonalization;
	friend struct TotalXSecBisection;
	friend struct slopeAndTotalXSec;
	friend struct slopeInt;
	friend struct slopeBisection;

	public:

	/** A vector of <code>SubProcessHandler</code>s. */
	typedef vector<SubHdlPtr> SubHandlerList;

	/** A vector of <code>Cut</code>s. */
	typedef vector<CutsPtr> CutsList;

	/** A vector of <code>ProcessHandler</code>s. */
	typedef vector<ProHdlPtr> ProcessHandlerList;

	/** A vector of cross sections. */
	typedef vector<CrossSection> XSVector;

	/** A pair of multiplicities: hard, soft. */
	typedef pair<unsigned int, unsigned int> MPair;

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	MPIHandler(): softMult_(0), identicalToUE_(-1),
	PtOfQCDProc_(-1.0GeV), Ptmin_(-1.0GeV),
	hardXSec_(0millibarn), softXSec_(0millibarn),
	totalXSecExp_(0*millibarn),
	softMu2_(ZERO), beta_(100.0/GeV2),
	algorithm_(2), numSubProcs_(0),
	colourDisrupt_(0.0), softInt_(true), twoComp_(true),
	DLmode_(2), avgNhard_(0.0), avgNsoft_(0.0),
	energyExtrapolation_(2), EEparamA_(0.6*GeV),
	EEparamB_(37.5GeV), refScale_(7000.GeV),
	pT0_(3.11*GeV), b_(0.21) {}

	/**
	* The destructor.
	*/
	virtual ~MPIHandler(){}
	//@}

	public:

	/** @name Methods for the MPI generation. */
	//@{

	/*
	* @return true if for this beam setup MPI can be generated
	*/
	virtual bool beamOK() const;

	/**
	* Return true or false depending on whether soft interactions are enabled.
	*/
	virtual bool softInt() const {return softInt_;}

	/**
	* Get the soft multiplicity from the pretabulated multiplicity
	* distribution. Generated in multiplicity in the first place.
	* @return the number of extra soft events in this collision
	*/
	virtual unsigned int softMultiplicity() const {return softMult_;}

	/**
	* Sample from the pretabulated multiplicity distribution.
	* @return the number of extra events in this collision
	*/
	virtual unsigned int multiplicity(unsigned int sel=0);

	/**
	* Select a StandardXComb according to it's weight
	* @return that StandardXComb Object
	* @param sel is the subprocess that should be returned,
	* if more than one is specified.
	*/
	virtual tStdXCombPtr generate(unsigned int sel=0);
	//@}


	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	/**
	* Initialize this Multiple Interaction handler and all related objects needed to
	* generate additional events.
	*/
	virtual void initialize();

	/**
	* Finalize this Multiple Interaction handler and all related objects needed to
	* generate additional events.
	*/
	virtual void finalize();

	/**
	+ * Clean up the XCombs from our subprocesses after each event.
	+ * ThePEG cannot see them, so the usual cleaning misses these.
	+ */
	+ virtual void clean();
	+
	+ /**
	* Write out accumulated statistics about integrated cross sections.
	*/
	void statistics() const;

	/**
	* The level of statistics. Controlls the amount of statistics
	* written out after each run to the <code>EventGenerator</code>s
	* <code>.out</code> file. Simply the EventHandler method is called here.
	*/
	int statLevel() const {return eventHandler()->statLevel();}

	/**
	* Return the hard cross section above ptmin
	*/
	CrossSection hardXSec() const { return hardXSec_; }

	/**
	* Return the soft cross section below ptmin
	*/
	CrossSection softXSec() const { return softXSec_; }

	/**
	* Return the inelastic cross section
	*/
	CrossSection inelasticXSec() const { return inelXSec_; }

	/** @name Simple access functions. */
	//@{

	/**
	* Return the ThePEG::EventHandler assigned to this handler.
	* This methods shadows ThePEG::StepHandler::eventHandler(), because
	* it is not virtual in ThePEG::StepHandler. This is ok, because this
	* method would give a null-pointer at some stages, whereas this method
	* gives access to the explicitely copied pointer (in initialize())
	* to the ThePEG::EventHandler.
	*/
	tEHPtr eventHandler() const {return theHandler;}

	/**
	* Return the current handler
	*/
	static const MPIHandler * currentHandler() {
	return currentHandler_;
	}

	/**
	* Return theAlgorithm.
	*/
	virtual int Algorithm() const {return algorithm_;}

	/**
	* Return the ptmin parameter of the model
	*/
	virtual Energy Ptmin() const {
	if(Ptmin_ > ZERO)
	return Ptmin_;
	else
	throw Exception() << "MPIHandler::Ptmin called without initialize before"
	<< Exception::runerror;
	}

	/**
	* Return the slope of the soft pt spectrum as calculated.
	*/
	virtual InvEnergy2 beta() const {
	if(beta_ != 100.0/GeV2)
	return beta_;
	else
	throw Exception() << "MPIHandler::beta called without initialization"
	<< Exception::runerror;
	}

	/**
	* Return the pt Cutoff of the Interaction that is identical to the UE
	* one.
	*/
	virtual Energy PtForVeto() const {return PtOfQCDProc_;}

	/**
	* Return the number of additional "hard" processes ( = multiple
	* parton scattering)
	*/
	virtual unsigned int additionalHardProcs() const {return numSubProcs_-1;}

	/**
	* Return the fraction of colour disrupted connections to the
	* suprocesses.
	*/
	virtual double colourDisrupt() const {return colourDisrupt_;}

	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;
	//@}

	private:

	/**
	* Access the list of sub-process handlers.
	*/
	const SubHandlerList & subProcesses()
	const {return theSubProcesses;}

	/**
	* Access the list of sub-process handlers.
	*/
	SubHandlerList & subProcesses() {return theSubProcesses;}

	/**
	* Access the list of cuts.
	*/
	const CutsList & cuts() const {return theCuts;}

	/**
	* Access the list of cuts.
	*/
	CutsList & cuts() {return theCuts;}

	/**
	* Access the list of sub-process handlers.
	*/
	const ProcessHandlerList & processHandlers()
	const {return theProcessHandlers;}

	/**
	* Access the list of sub-process handlers.
	*/
	ProcessHandlerList & processHandlers() {return theProcessHandlers;}


	/**
	* Method to calculate the individual probabilities for N scatters in the event.
	* @param UEXSecs is(are) the inclusiv cross section(s) for the UE process(es).
	*/
	void Probs(XSVector UEXSecs);

	/**
	* Debug method to check the individual probabilities.
	* @param filename is the file the output gets written to
	*/
	void MultDistribution(string filename) const;

	/**
	* Return the value of the Overlap function A(b) for a given impact
	* parameter \a b.
	* @param b impact parameter
	* @param mu2 = inv hadron radius squared. 0 will use the value of
	* invRadius_
	* @return inverse area.
	*/
	InvArea OverlapFunction(Length b, Energy2 mu2=ZERO) const;

	/**
	* Method to calculate the poisson probability for expectation value
	* \f$<n> = A(b)\sigma\f$, and multiplicity N.
	*/
	double poisson(Length b, CrossSection sigma,
	unsigned int N, Energy2 mu2=ZERO) const;

	/**
	* Return n!
	*/
	double factorial (unsigned int n) const;

	/**
	* Returns the total cross section for the current CMenergy. The
	* decision which parametrization will be used is steered by a
	* external parameter of this class.
	*/
	CrossSection totalXSecExp() const;

	/**
	* Difference of the calculated total cross section and the
	* experimental one from totalXSecExp.
	* @param softXSec = the soft cross section that is used
	* @param softMu2 = the soft radius, if 0 the hard radius will be used
	*/
	CrossSection totalXSecDiff(CrossSection softXSec,
	Energy2 softMu2=ZERO) const;

	/**
	* Difference of the calculated elastic slope and the
	* experimental one from slopeExp.
	* @param softXSec = the soft cross section that is used
	* @param softMu2 = the soft radius, if 0 the hard radius will be used
	*/
	InvEnergy2 slopeDiff(CrossSection softXSec,
	Energy2 softMu2=ZERO) const;

	/**
	* Returns the value of the elastic slope for the current CMenergy.
	* The decision which parametrization will be used is steered by a
	* external parameter of this class.
	*/
	InvEnergy2 slopeExp() const;


	/**
	* Calculate the minimal transverse momentum from the extrapolation
	*/
	void overrideUECuts();


	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<MPIHandler> initMPIHandler;

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	MPIHandler & operator=(const MPIHandler &);

	/**
	* A pointer to the EventHandler that calls us. Has to be saved, because the
	* method eventHandler() inherited from ThePEG::StepHandler returns a null-pointer
	* sometimes. Leif changed that in r1053 so that a valid pointer is present, when
	* calling doinitrun().
	*/
	tEHPtr theHandler;

	/**
	* The list of <code>SubProcessHandler</code>s.
	*/
	SubHandlerList theSubProcesses;

	/**
	* The kinematical cuts used for this collision handler.
	*/
	CutsList theCuts;

	/**
	* List of ProcessHandler used to sample different processes independently
	*/
	ProcessHandlerList theProcessHandlers;

	/**
	* A ThePEG::Selector where the individual Probabilities P_N are stored
	* and the actual Multiplicities can be selected.
	*/
	Selector<MPair> theMultiplicities;

	/**
	* Variable to store the soft multiplicity generated for a event. This
	* has to be stored as it is generated at the time of the hard
	* additional interactions but used later on.
	*/
	unsigned int softMult_;

	/**
	* Variable to store the multiplicity of the second hard process
	*/
	vector<int> additionalMultiplicities_;

	/**
	* Variable to store the information, which process is identical to
	* the UE one (QCD dijets).
	* 0 means "real" hard one
	* n>0 means the nth additional hard scatter
	* -1 means no one!
	*/
	int identicalToUE_;

	/**
	* Variable to store the minimal pt of the process that is identical
	* to the UE one. This only has to be set, if it can't be determined
	* automatically (i.e. when reading QCD LesHouches files in).
	*/
	Energy PtOfQCDProc_;

	/**
	* Variable to store the parameter ptmin
	*/
	Energy Ptmin_;

	/**
	* Variable to store the hard cross section above ptmin
	*/
	CrossSection hardXSec_;

	/**
	* Variable to store the final soft cross section below ptmin
	*/
	CrossSection softXSec_;

	/**
	* Variable to store the inelastic cross section
	*/
	CrossSection inelXSec_;

	/**
	* Variable to store the total pp cross section (assuming rho=0!) as
	* measured at LHC. If this variable is set, this value is used in the
	* subsequent run instead of any of the Donnachie-Landshoff
	* parametrizations.
	*/
	CrossSection totalXSecExp_;

	/**
	* Variable to store the soft radius, that is calculated during
	* initialization for the two-component model.
	*/
	Energy2 softMu2_;

	/**
	* slope to the non-perturbative pt spectrum: \f$d\sigma/dp_T^2 = A \exp
	* (- beta p_T^2)\f$. Its value is determined durint initialization.
	*/
	InvEnergy2 beta_;
	/**
	* Switch to be set from outside to determine the algorithm used for
	* UE activity.
	*/
	int algorithm_;

	/**
	* Inverse hadron Radius squared \f$ (\mu^2) \f$. Used inside the overlap function.
	*/
	Energy2 invRadius_;

	/**
	* Member variable to store the actual number of separate SubProcesses
	*/
	unsigned int numSubProcs_;

	/**
	* Variable to store the relative number of colour disrupted
	* connections to additional subprocesses. This variable is used in
	* Herwig::HwRemDecayer but store here, to have access to all
	* parameters through one Object.
	*/
	double colourDisrupt_;

	/**
	* Flag to store whether soft interactions, i.e. pt < ptmin should be
	* simulated.
	*/
	bool softInt_;

	/**
	* Flag to steer wheather the soft part has a different radius, that
	* will be dynamically fixed.
	*/
	bool twoComp_;

	/**
	* Switch to determine which Donnachie & Landshoff parametrization
	* should be used.
	*/
	unsigned int DLmode_;

	/**
	* Variable to store the average hard multiplicity.
	*/
	double avgNhard_;

	/**
	* Variable to store the average soft multiplicity.
	*/
	double avgNsoft_;

	/**
	* The current handler
	*/
	static MPIHandler * currentHandler_;

	/**
	* Flag to store whether to calculate the minimal UE pt according to an
	* extrapolation formula or whether to use MPIHandler:Cuts[0]:OneCuts[0]:MinKT
	*/
	unsigned int energyExtrapolation_;

	/**
	* Parameters for the energy extrapolation formula
	*/
	Energy EEparamA_;
	Energy EEparamB_;
	Energy refScale_;
	Energy pT0_;
	double b_;

	protected:

	/** @cond EXCEPTIONCLASSES */

	/**
	* Exception class used by the MultipleInteractionHandler, when something
	* during initialization went wrong.
	* \todo understand!!!
	*/
	class InitError: public Exception {};

	/** @endcond */

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of MPIHandler. */
	template <>
	struct BaseClassTrait<Herwig::MPIHandler,1> {
	/** Typedef of the first base class of MPIHandler. */
	typedef Interfaced NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the MPIHandler class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::MPIHandler>
	: public ClassTraitsBase<Herwig::MPIHandler> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::MPIHandler"; }
	/** Return the name(s) of the shared library (or libraries) be loaded to get
	* access to the MPIHandler class and any other class on which it depends
	* (except the base class). */
	static string library() { return "SimpleKTCut.so HwMPI.so"; }
	};

	/** @endcond */

	}

	namespace Herwig {

	/**
	* A struct for the 2D root finding that is necessary to determine the
	* soft cross section and the soft radius that is needed to describe
	* the total cross section correctly.
	* NOT IN USE CURRENTLY
	*/
	struct slopeAndTotalXSec : public GSLHelper<CrossSection, CrossSection> {

	public:

	/**
	* Constructor
	*/
	slopeAndTotalXSec(tcMPIHPtr handler): handler_(handler) {}

	/** second argument type */
	typedef Energy2 ArgType2;

	/** second value type */
	typedef InvEnergy2 ValType2;

	/** first element of the vector like function to find root for
	* @param softXSec soft cross-section
	* @param softMu2 \f$\mu^2\f$
	*/
	CrossSection f1(ArgType softXSec, ArgType2 softMu2) const {
	return handler_->totalXSecDiff(softXSec, softMu2);
	}

	/** second element of the vector like function to find root for
	* @param softXSec soft cross-section
	* @param softMu2 \f$\mu^2\f$
	*/
	InvEnergy2 f2(ArgType softXSec, ArgType2 softMu2) const {
	return handler_->slopeDiff(softXSec, softMu2);
	}

	/** provide the actual units of use */
	virtual ValType vUnit() const {return 1.0*millibarn;}

	/** otherwise rounding errors may get significant */
	virtual ArgType aUnit() const {return 1.0*millibarn;}

	/** provide the actual units of use */
	ValType2 vUnit2() const {return 1.0/GeV2;}

	/** otherwise rounding errors may get significant */
	ArgType2 aUnit2() const {return GeV2;}

	private:

	/**
	* Pointer to the handler
	*/
	tcMPIHPtr handler_;

	};

	/**
	* A struct for the root finding that is necessary to determine the
	* slope of the soft pt spectrum to match the soft cross section
	*/
	struct betaBisection : public GSLHelper<Energy2, InvEnergy2>{
	public:
	/**
	* Constructor.
	* @param soft = soft cross section, i.e. the integral of the soft
	* pt spectrum f(u=p_T^2) = dsig exp(-beta*u/u_min)
	* @param dsig = dsigma_hard/dp_T^2 at the p_T cutoff
	* @param ptmin = p_T cutoff
	*/
	betaBisection(CrossSection soft, DiffXSec dsig, Energy ptmin)
	: softXSec_(soft), dsig_(dsig), ptmin_(ptmin) {}

	/**
	* Operator that is used inside the GSLBisection class
	*/
	virtual Energy2 operator ()(InvEnergy2 beta) const
	{
	if( fabs(beta*GeV2) < 1.E-4 )
	beta = (beta > ZERO) ? 1.E-4/GeV2 : -1.E-4/GeV2;

	return (exp(beta*sqr(ptmin_)) - 1.0)/beta - softXSec_/dsig_;
	}

	/** provide the actual units of use */
	virtual ValType vUnit() const {return 1.0*GeV2;}

	/** provide the actual units of use */
	virtual ArgType aUnit() const {return 1.0/GeV2;}

	private:

	/** soft cross section */
	CrossSection softXSec_;

	/** dsigma/dp_T^2 at ptmin */
	DiffXSec dsig_;

	/** pt cutoff */
	Energy ptmin_;
	};

	/**
	* A struct for the root finding that is necessary to determine the
	* soft cross section and soft mu2 that are needed to describe the
	* total cross section AND elastic slope correctly.
	*/
	struct slopeBisection : public GSLHelper<InvEnergy2, Energy2> {
	public:
	/** Constructor */
	slopeBisection(tcMPIHPtr handler) : handler_(handler) {}

	/**
	* Return the difference of the calculated elastic slope to the
	* experimental one for a given value of the soft mu2. During that,
	* the soft cross section get fixed.
	*/
	InvEnergy2 operator ()(Energy2 arg) const;

	/** Return the soft cross section that has been calculated */
	CrossSection softXSec() const {return softXSec_;}

	private:
	/** const pointer to the MPIHandler to give access to member functions.*/
	tcMPIHPtr handler_;
	/** soft cross section that is determined on the fly.*/
	mutable CrossSection softXSec_;
	};

	/**
	* A struct for the root finding that is necessary to determine the
	* soft cross section that is needed to describe the total cross
	* section correctly.
	*/
	struct TotalXSecBisection : public GSLHelper<CrossSection, CrossSection> {
	public:

	/**
	* Constructor
	* @param handler The handler
	* @param softMu2 \f$\mu^2\f$
	*/
	TotalXSecBisection(tcMPIHPtr handler, Energy2 softMu2=ZERO):
	handler_(handler), softMu2_(softMu2) {}

	/**
	* operator to return the cross section
	* @param argument input cross section
	*/
	CrossSection operator ()(CrossSection argument) const {
	return handler_->totalXSecDiff(argument, softMu2_);
	}

	/** provide the actual units of use */
	virtual ValType vUnit() const {return 1.0*millibarn;}

	/** otherwise rounding errors may get significant */
	virtual ArgType aUnit() const {return 1.0*millibarn;}

	private:

	/**
	* The handler
	*/
	tcMPIHPtr handler_;

	/**
	* \f$\mu^2\f$
	*/
	Energy2 softMu2_;

	};

	/**
	* Typedef for derivative of the length
	*/
	typedef QTY<1,-2,0>::Type LengthDiff;

	/**
	* A struct for the integrand for the slope
	*/
	struct slopeInt : public GSLHelper<LengthDiff, Length>{

	public:
	/** Constructor
	* @param handler The handler
	* @param hard The hard cross section
	* @param soft The soft cross section
	* @param softMu2 \f$\mu^2\f$
	*/
	slopeInt(tcMPIHPtr handler, CrossSection hard,
	CrossSection soft=0*millibarn, Energy2 softMu2=ZERO)
	: handler_(handler), hardXSec_(hard),
	softXSec_(soft), softMu2_(softMu2) {}

	/**
	* Operator to return the answer
	* @param arg The argument
	*/
	ValType operator ()(ArgType arg) const;

	private:

	/**
	* Pointer to the Handler that calls this integrand
	*/
	tcMPIHPtr handler_;

	/**
	* The hard cross section to be eikonalized
	*/
	CrossSection hardXSec_;

	/**
	* The soft cross section to be eikonalized. Default is zero
	*/
	CrossSection softXSec_;

	/**
	* The inv radius^2 of the soft interactions.
	*/
	Energy2 softMu2_;

	};

	/**
	* A struct for the eikonalization of the inclusive cross section.
	*/
	struct Eikonalization : public GSLHelper<Length, Length>{

	/**
	* The constructor
	* @param handler is the pointer to the MPIHandler to get access to
	* MPIHandler::OverlapFunction and member variables of the MPIHandler.
	* @param option is a flag, whether the inelastic or the total
	* @param handler The handler
	* @param hard The hard cross section
	* @param soft The soft cross section
	* @param softMu2 \f$\mu^2\f$
	* cross section should be returned (-2 or -1). For option = N > 0 the integrand
	* is N(A(b)sigma)^N/N! exp(-A(b)sigma) this is the P_Nsigma where
	* P_N is the Probability of having exactly N interaction (including the hard one)
	* This is equation 14 from "Jimmy4: Multiparton Interactions in HERWIG for the LHC"
	*/
	Eikonalization(tcMPIHPtr handler, int option, CrossSection hard,
	CrossSection soft=0*millibarn, Energy2 softMu2=ZERO)
	: theHandler(handler), theoption(option), hardXSec_(hard),
	softXSec_(soft), softMu2_(softMu2) {}

	/**
	* Get the function value
	*/
	Length operator ()(Length argument) const;

	private:
	/**
	* Pointer to the Handler that calls this integrand
	*/
	tcMPIHPtr theHandler;

	/**
	* A flag to switch between the calculation of total and inelastic cross section
	* or calculations for the individual probabilities. See the constructor
	*/
	int theoption;

	/**
	* The hard cross section to be eikonalized
	*/
	CrossSection hardXSec_;

	/**
	* The soft cross section to be eikonalized. Default is zero
	*/
	CrossSection softXSec_;

	/**
	* The inv radius^2 of the soft interactions.
	*/
	Energy2 softMu2_;

	};
	}


	#ifndef ThePEG_TEMPLATES_IN_CC_FILE
	// #include "MPIHandler.tcc"
	#endif

	#endif /* HERWIG_MPIHandler_H */
	diff --git a/src/LHC-RPV.in b/src/LHC-RPV.in
	--- a/src/LHC-RPV.in
	+++ b/src/LHC-RPV.in
	@@ -1,69 +1,79 @@
	##################################################
	# Example generator for the MSSM
	# The best way to use this is to make your own
	# copy of this file and edit that as you require.
	#
	# The first section loads the model file which
	# does not contain anything that users need to touch.
	#
	# The second section contains the user settings.
	###################################################

	read MSSM.model
	cd /Herwig/NewPhysics

	##################################################
	#
	# This section contains the user defined settings
	#
	##################################################
	# --- Hard Process ----
	# The particle name can be found in the relevant model file
	# by searching for its PDG code and noting the text
	# '/Herwig/Particles/###' where the hashes denote the name

	# Switch to decide whether to include EW diagrams in the
	# hard process (On by default)
	set HPConstructor:IncludeEW No

	# Example hard process: Incoming proton, outgoing squarks
	insert HPConstructor:Incoming 0 /Herwig/Particles/g
	insert HPConstructor:Incoming 1 /Herwig/Particles/u
	insert HPConstructor:Incoming 2 /Herwig/Particles/ubar
	insert HPConstructor:Incoming 3 /Herwig/Particles/d
	insert HPConstructor:Incoming 4 /Herwig/Particles/dbar

	insert HPConstructor:Outgoing 0 /Herwig/Particles/~u_L
	insert HPConstructor:Outgoing 1 /Herwig/Particles/~u_Lbar
	insert HPConstructor:Outgoing 2 /Herwig/Particles/~d_L
	insert HPConstructor:Outgoing 3 /Herwig/Particles/~d_Lbar

	# --- Perturbative Decays ---
	# Read in the spectrum file and optional decay table.
	# If a decay table is in a separate file
	# then add another 'setup' line with that
	# file as the argument. The provided
	# spectrum file is an example using ISAJET
	setup MSSM/Model RPV3.1.slha

	# To disable a particular decay mode, add it's tag to the DisableModes
	# interface of the DecayConstructor object, i.e.
	#insert DecayConstructor:DisableModes 0 ~u_L->~chi_20,u;
	#insert DecayConstructor:DisableModes 1 ~chi_20->~e_R-,e+;
	# etc ...

	# To set a minimum allowed branching fraction (the default is shown)
	#set NewModel:MinimumBR 1e-6


	# Intrinsic pT tune extrapolated to LHC energy
	set /Herwig/Shower/Evolver:IntrinsicPtGaussian 2.2*GeV
	+cd /Herwig/Generators
	+set LHCGenerator:DumpPeriod 50
	+set LHCGenerator:KeepAllDumps Yes
	+
	+#set LHCGenerator:EventHandler:CascadeHandler NULL
	+#set LHCGenerator:EventHandler:CascadeHandler:MPIHandler NULL
	+#set LHCGenerator:EventHandler:HadronizationHandler NULL
	+#set LHCGenerator:EventHandler:DecayHandler NULL
	+#set /Herwig/Analysis/Basics:CheckQuark No
	+

	# Other parameters for run
	cd /Herwig/Generators
	set LHCGenerator:NumberOfEvents 10000000
	set LHCGenerator:RandomNumberGenerator:Seed 31122001
	-set LHCGenerator:PrintEvent 10
	+set LHCGenerator:PrintEvent 1000
	set LHCGenerator:MaxErrors 10000

	saverun LHC-RPV LHCGenerator
	diff --git a/src/LHC.in b/src/LHC.in
	--- a/src/LHC.in
	+++ b/src/LHC.in
	@@ -1,145 +1,157 @@
	##################################################
	# Example generator based on LHC parameters
	# usage: Herwig++ read LHC.in
	##################################################

	##################################################
	# Technical parameters for this run
	##################################################
	cd /Herwig/Generators
	set LHCGenerator:NumberOfEvents 10000000
	set LHCGenerator:RandomNumberGenerator:Seed 31122001
	-set LHCGenerator:PrintEvent 10
	+set LHCGenerator:PrintEvent 1000
	set LHCGenerator:MaxErrors 10000

	##################################################
	# LHC physics parameters (override defaults here)
	##################################################

	# Intrinsic pT tune extrapolated to LHC energy
	set /Herwig/Shower/Evolver:IntrinsicPtGaussian 2.2*GeV

	########################
	## sqrt(s) = 8000 GeV ##
	########################
	set LHCGenerator:EventHandler:LuminosityFunction:Energy 8000.0

	##################################################
	# Matrix Elements for hadron-hadron collisions
	# (by default only gamma/Z switched on)
	##################################################
	cd /Herwig/MatrixElements/

	#
	# Electroweak boson W/Z processes
	#
	# Drell-Yan Z/gamma
	insert SimpleQCD:MatrixElements[0] MEqq2gZ2ff
	#
	# Drell-Yan W
	# insert SimpleQCD:MatrixElements[0] MEqq2W2ff
	#
	# W+jet
	# insert SimpleQCD:MatrixElements[0] MEWJet
	#
	-# Z+jet
	+# Z+jet
	# insert SimpleQCD:MatrixElements[0] MEZJet
	#
	# WW/WZ/ZZ
	# insert SimpleQCD:MatrixElements[0] MEPP2VV
	#
	# Wgamma/Zgamma
	# insert SimpleQCD:MatrixElements[0] MEPP2VGamma

	#
	# QCD and gamma processes
	#
	# QCD 2-2 scattering
	# insert SimpleQCD:MatrixElements[0] MEQCD2to2
	#
	# top-antitop production
	# insert SimpleQCD:MatrixElements[0] MEHeavyQuark
	#
	# gamma+jet
	# insert SimpleQCD:MatrixElements[0] MEGammaJet
	#
	# gamma-gamma
	# insert SimpleQCD:MatrixElements[0] MEGammaGamma

	#
	# Higgs Processes
	#
	#
	# gg/qqbar -> Higgs (recommend including q qbar->Hg as not in ME correction)
	# insert SimpleQCD:MatrixElements[0] MEHiggs
	# insert SimpleQCD:MatrixElements[0] MEHiggsJet
	# set MEHiggsJet:Process qqbar
	# set /Herwig/Cuts/JetKtCut:MinKT 0.0*GeV
	#
	# higgs+jet
	# insert SimpleQCD:MatrixElements[0] MEHiggsJet
	#
	# higgs + W (N.B. if considering all W decay modes useful to set )
	# (jet pT cut to zero so no cut on W decay products )
	# insert SimpleQCD:MatrixElements[0] MEPP2WH
	# set /Herwig/Cuts/JetKtCut:MinKT 0.0*GeV
	#
	# higgs + Z (N.B. if considering all Z decay modes useful to set )
	# (jet pT cut to zero so no cut on Z decay products )
	# insert SimpleQCD:MatrixElements[0] MEPP2ZH
	# set /Herwig/Cuts/JetKtCut:MinKT 0.0*GeV
	#
	# VBF Higgs
	# insert SimpleQCD:MatrixElements[0] MEPP2HiggsVBF
	#
	# t tbar Higgs
	# insert SimpleQCD:MatrixElements[0] MEPP2ttbarH
	#
	# b bbar Higgs
	# insert SimpleQCD:MatrixElements[0] MEPP2bbbarH

	cd /Herwig/Generators

	##################################################
	# Useful analysis handlers for hadron-hadron physics
	##################################################
	# analysis of W/Z events
	# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/DrellYan
	# analysis of top-antitop events
	# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/TTbar
	# analysis of gamma+jet events
	# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/GammaJet
	# analysis of gamma-gamma events
	# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/GammaGamma
	# analysis of higgs-jet events
	# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/HiggsJet

	##################################################
	# Useful analysis handlers for HepMC related output
	##################################################
	# Schematic overview of an event (requires --with-hepmc to be set at configure time
	# and the graphviz program 'dot' to produce a plot)
	# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/Plot
	# A HepMC dump file (requires --with-hepmc to be set at configure time)
	# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
	# set /Herwig/Analysis/HepMCFile:PrintEvent 100
	# set /Herwig/Analysis/HepMCFile:Format GenEvent
	# set /Herwig/Analysis/HepMCFile:Units GeV_mm

	##################################################
	# Save run for later usage with 'Herwig++ run'
	##################################################
	+
	+
	+set LHCGenerator:DumpPeriod 50
	+set LHCGenerator:KeepAllDumps Yes
	+
	+#set LHCGenerator:EventHandler:CascadeHandler NULL
	+#set LHCGenerator:EventHandler:CascadeHandler:MPIHandler NULL
	+#set LHCGenerator:EventHandler:HadronizationHandler NULL
	+#set LHCGenerator:EventHandler:DecayHandler NULL
	+#set /Herwig/Analysis/Basics:CheckQuark No
	+
	+
	saverun LHC LHCGenerator

	##################################################
	# uncomment this section for an example batch run
	# of two repeats with different parameters
	#
	# Note that a separate call of 'Herwig run'
	# is not required in this case
	##################################################
	# set LHCGenerator:NumberOfEvents 10
	# run LHC-full LHCGenerator
	#
	# set LHCGenerator:EventHandler:LuminosityFunction:Energy 900.0
	# run LHC-initial LHCGenerator

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:54 PM (1 d, 7 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805937
Default Alt Text: (193 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions