No OneTemporary
Actions

Size

49 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/Hadron/MEDiffraction.cc b/MatrixElement/Hadron/MEDiffraction.cc
	--- a/MatrixElement/Hadron/MEDiffraction.cc
	+++ b/MatrixElement/Hadron/MEDiffraction.cc
	@@ -1,853 +1,873 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the MEDiffraction class.
	//

	#include "MEDiffraction.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Utilities/SimplePhaseSpace.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Handlers/StandardXComb.h"


	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace Herwig;

	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"

	#include "Herwig/Utilities/Kinematics.h"


	MEDiffraction::MEDiffraction()
	: HwMEBase(),
	deltaOnly(false),
	- isInRunPhase(false) {}
	+ isInRunPhase(false),
	+ theProtonMass(0.93827203*GeV)
	+{}


	void MEDiffraction::getDiagrams() const {
	//incoming particles
	cPDPair incomingHardons = generator()->eventHandler()->incoming();

	tcPDPtr pom = getParticleData(990);

	//get incoming particles
	tcPDPtr prt11 = getParticleData(incomingHardons.first->id());
	tcPDPtr prt12 = getParticleData(incomingHardons.second->id());

	//get sign of id
	int sign1=0, sign2=0;
	sign1 = (incomingHardons.first->id() > 0) ? 1 : -1;
	sign2 = (incomingHardons.second->id() > 0) ? 1 : -1;

	tcPDPtr prt21 = getParticleData(sign1*2214);//Delta+
	tcPDPtr prt22 = getParticleData(sign2*2214);//Delta+

	//for the left side
	tcPDPtr q11 = getParticleData(sign1*2); //u
	tcPDPtr q21 = getParticleData(sign1*1); //d
	//for the right side
	tcPDPtr q12 = getParticleData(sign2*2); //u
	tcPDPtr q22 = getParticleData(sign2*1); //d
	//for the left side
	tcPDPtr dq11 = getParticleData(sign1*2101); //ud_0
	tcPDPtr dq111 = getParticleData(sign1*2103); //ud_1
	tcPDPtr dq21 = getParticleData(sign1*2203); //uu_1
	//for the right side
	tcPDPtr dq12 = getParticleData(sign2*2101); //ud_0
	tcPDPtr dq112 = getParticleData(sign2*2103); //ud_1
	tcPDPtr dq22 = getParticleData(sign2*2203); //uu_1

	tcPDPtr gl = getParticleData(21);//gluon

	//switch between dissociation decays to different
	//number of clusters or dissociation into delta only
	//(Maybe can be automated???)
	//(Should be generalized to ppbar, for example!!!)
	switch(dissociationDecay){
	case 0: //one cluster or only delta in the final state
	if(deltaOnly) //only delta in the final state
	{

	switch (diffDirection){
	case 0:
	add(new_ptr((Tree2toNDiagram(3), prt11, pom, prt12, 1, prt21, 2, prt12, -1)));
	break;
	case 1:
	add(new_ptr((Tree2toNDiagram(3), prt11, pom, prt12, 1, prt11, 2, prt22, -1)));
	break;
	case 2:
	add(new_ptr((Tree2toNDiagram(3), prt11, pom, prt12, 1, prt21, 2, prt22, -1)));
	break;
	}


	}else
	{
	//switch between direction of dissociated proton for single diffraction or
	//double diffraction
	switch (diffDirection){
	case 0: //left
	//u -- ud_0
	add(new_ptr((Tree2toNDiagram(4), prt11, q11, pom, prt12, 3, prt12, 1, dq11, 2, q11, -1)));
	//d -- uu_1
	add(new_ptr((Tree2toNDiagram(4), prt11, q21, pom, prt12, 3, prt12, 1, dq21, 2, q21, -2)));
	break;
	case 1: //right
	//u -- ud_0
	add(new_ptr((Tree2toNDiagram(4), prt11, pom, q12, prt12, 1, prt11, 3, dq12, 2, q12, -1)));

	//d -- uu_1
	add(new_ptr((Tree2toNDiagram(4), prt11, pom, q22, prt12, 1, prt11, 3, dq22, 2, q22, -2)));
	break;
	case 2: //double
	//u -- ud_0 left u -- ud_0 right
	add(new_ptr((Tree2toNDiagram(5), prt11, q11, pom, q12, prt12, 1, dq11, 2, q11, 3, q12, 4, dq12, -1)));

	//u -- ud_0 left d -- uu_1 right
	add(new_ptr((Tree2toNDiagram(5), prt11, q11, pom, q22, prt12, 1, dq11, 2, q11, 3, q22, 4, dq22, -2)));

	//d -- uu_1 left u -- ud_0 right
	add(new_ptr((Tree2toNDiagram(5), prt11, q21, pom, q12, prt12, 1,dq21, 2, q21, 3, q12, 4, dq12, -3)));

	//d -- uu_1 left d -- uu_1 right
	add(new_ptr((Tree2toNDiagram(5), prt11, q21, pom, q22, prt12, 1, dq21, 2, q21, 3, q22, 4, dq22, -4)));
	break;
	}

	}
	break;
	case 1: //two clusters (cases with ud_1 not included)
	switch (diffDirection){
	case 0: //left
	//u -- ud_0
	add(new_ptr((Tree2toNDiagram(5), prt11, q11, gl, pom, prt12, 1, dq11, 2, q11, 3, gl, 4, prt12, -1)));
	//d -- uu_1
	add(new_ptr((Tree2toNDiagram(5), prt11, q21, gl, pom, prt12, 1, dq21, 2, q21, 3, gl, 4, prt12, -2)));
	break;
	case 1: //right
	//u -- ud_0
	add(new_ptr((Tree2toNDiagram(5), prt11, pom, gl, q12, prt12, 1, prt11, 2, gl, 3, q12, 4, dq12, -1)));
	//d -- ud_1
	add(new_ptr((Tree2toNDiagram(5), prt11, pom, gl, q22, prt12, 1, prt11, 2, gl, 3, q22, 4, dq22, -2)));
	break;
	case 2: //double
	//u -- ud_0 left u -- ud_0 right
	add(new_ptr((Tree2toNDiagram(7), prt11, q11, gl, pom, gl, q12, prt12, 1, dq11, 2, q11, 3, gl, 4,
	gl, 5, q12, 6, dq12, -1)));
	//u -- ud_0 left d -- uu_1 right
	add(new_ptr((Tree2toNDiagram(7), prt11, q11, gl, pom, gl, q22, prt12, 1, dq11, 2, q11, 3, gl, 4,
	gl, 5, q22, 6, dq22, -2)));
	//d -- uu_1 left u -- ud_0 right
	add(new_ptr((Tree2toNDiagram(7), prt11, q21, gl, pom, gl, q12, prt12, 1, dq21, 2, q21, 3, gl, 4,
	gl, 5, q12, 6, dq12, -3)));
	//d -- uu_1 left d -- uu_1 right
	add(new_ptr((Tree2toNDiagram(7), prt11, q21, gl, pom, gl, q22, prt12, 1, dq21, 2, q21, 3, gl, 4,
	gl, 5, q22, 6, dq22, -4)));
	break;
	}
	break;
	}
	}

	Energy2 MEDiffraction::scale() const {
	return sqr(10*GeV);
	}

	int MEDiffraction::nDim() const {
	return 0;
	}

	void MEDiffraction::setKinematics() {
	HwMEBase::setKinematics(); // Always call the base class method first
	}

	bool MEDiffraction::generateKinematics(const double * ) {
	// generate the masses of the particles
	for (size_t i = 2; i < meMomenta().size(); ++i)
	meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->generateMass());

	/* sample M12, M22 and t, characterizing the diffractive final state */
	const pair<pair<Energy2,Energy2>,Energy2> point = diffractiveMassAndMomentumTransfer();
	const Energy2 M12 (point.first.first);
	const Energy2 M22 (point.first.second);
	const Energy2 t(point.second);


	/* construct the hadronic momenta in the lab frame */
	const double phi = UseRandom::rnd() * Constants::twopi;
	const Energy cmEnergy = generator()->maximumCMEnergy();
	const Energy2 s = sqr(cmEnergy);

	//proton mass
	- const Energy2 m2 = sqr( getParticleData(2212)->mass() );
	+ const Energy2 m2 = sqr( theProtonMass );

	const Energy E3 = (s - M22 + M12) / (2.*cmEnergy);
	const Energy E4 = (s + M22 - M12) / (2.*cmEnergy);

	//Momentum of outgoing proton and dissociated proton
	const Energy pprime = sqrt(kallen(s, M12, M22)) / (2.*cmEnergy);

	//costheta of scattering angle
	double costheta = s(s + 2t - 2*m2 - M12 - M22)
	/ sqrt( kallen(s, M12, M22)*kallen(s, m2, m2) );

	assert(abs(costheta)<=1.);

	const Energy pzprime = pprime*costheta;
	const Energy pperp = pprime*sqrt(1 - sqr(costheta));

	/* momenta in the lab frame */
	const Lorentz5Momentum p3 = Lorentz5Momentum(pperpcos(phi), pperpsin(phi), pzprime, E3);
	const Lorentz5Momentum p4 = Lorentz5Momentum(-pperpcos(phi), -pperpsin(phi), -pzprime, E4);

	/* decay dissociated proton into quark-diquark */
	//squares of constituent masses of quark and diquark
	const Energy2 mqq2(sqr(mqq())), mq2(sqr(mq()));

	Energy2 Mx2;
	switch(diffDirection){
	case 0:
	Mx2=M12;
	break;
	case 1:
	Mx2=M22;
	break;
	}


	/* Select between left/right single diffraction and double diffraction */
	//check if we want only delta for the excited state


	//pair of momenta for double decay for a two cluster case
	pair<Lorentz5Momentum,Lorentz5Momentum> momPair, momPair1;
	//fraction of momenta
	double frac = UseRandom::rnd();

	switch(dissociationDecay){
	case 0:
	if(!deltaOnly)
	{

	pair<Lorentz5Momentum,Lorentz5Momentum> decayMomenta;
	pair<Lorentz5Momentum,Lorentz5Momentum> decayMomentaTwo;
	const double phiprime = UseRandom::rnd() * Constants::twopi;

	//aligned with outgoing dissociated proton
	const double costhetaprime = costheta;

	const double sinthetaprime=sqrt(1-sqr(costhetaprime));
	//axis along which diquark from associated proton is aligned
	Axis dir = Axis(sinthetaprimecos(phiprime), sinthetaprimesin(phiprime), costhetaprime);

	switch (diffDirection){
	case 0://Left single diffraction
	meMomenta()[4].setT(sqrt(mq2+sqr(meMomenta()[4].x())+sqr(meMomenta()[4].y())+sqr(meMomenta()[4].z())));
	////////////////////////////////////////////////////

	do{}
	while(!Kinematics::twoBodyDecay(p3,mqq(),mq(),-dir,decayMomenta.first,decayMomenta.second));
	///////////

	meMomenta()[2].setVect(p4.vect());
	meMomenta()[2].setT(p4.t());

	meMomenta()[3].setVect(decayMomenta.first.vect());
	meMomenta()[3].setT(decayMomenta.first.t());
	meMomenta()[4].setVect(decayMomenta.second.vect());
	meMomenta()[4].setT(decayMomenta.second.t());

	meMomenta()[2].rescaleEnergy();
	meMomenta()[3].rescaleEnergy();
	meMomenta()[4].rescaleEnergy();
	break;
	case 1://Right single diffraction
	meMomenta()[4].setT(sqrt(mq2+sqr(meMomenta()[4].x())+sqr(meMomenta()[4].y())+sqr(meMomenta()[4].z())));
	////////////////////////////////////////////////////

	do{}
	while(!Kinematics::twoBodyDecay(p4,mqq(),mq(),dir,decayMomenta.first,decayMomenta.second));

	meMomenta()[2].setVect(p3.vect());
	meMomenta()[2].setT(p3.t());

	meMomenta()[3].setVect(decayMomenta.first.vect());
	meMomenta()[3].setT(decayMomenta.first.t());
	meMomenta()[4].setVect(decayMomenta.second.vect());
	meMomenta()[4].setT(decayMomenta.second.t());

	meMomenta()[2].rescaleEnergy();
	meMomenta()[3].rescaleEnergy();
	meMomenta()[4].rescaleEnergy();
	break;
	case 2://double diffraction

	do{}
	while(!Kinematics::twoBodyDecay(p3,mqq(),mq(),-dir,decayMomenta.first,decayMomenta.second));

	do{}
	while(!Kinematics::twoBodyDecay(p4,mqq(),mq(),dir,decayMomentaTwo.first,decayMomentaTwo.second));

	meMomenta()[2].setVect(decayMomenta.first.vect());
	meMomenta()[2].setT(decayMomenta.first.t());
	meMomenta()[3].setVect(decayMomenta.second.vect());
	meMomenta()[3].setT(decayMomenta.second.t());

	meMomenta()[4].setVect(decayMomentaTwo.second.vect());
	meMomenta()[4].setT(decayMomentaTwo.second.t());
	meMomenta()[5].setVect(decayMomentaTwo.first.vect());
	meMomenta()[5].setT(decayMomentaTwo.first.t());


	meMomenta()[2].rescaleEnergy();
	meMomenta()[3].rescaleEnergy();
	meMomenta()[4].rescaleEnergy();

	meMomenta()[5].rescaleEnergy();

	break;
	}

	}else
	{
	meMomenta()[2+diffDirection].setVect(p3.vect());
	meMomenta()[2+diffDirection].setT(p3.t());
	meMomenta()[3-diffDirection].setVect(p4.vect());
	meMomenta()[3-diffDirection].setT(p4.t());

	meMomenta()[2].rescaleEnergy();
	meMomenta()[3].rescaleEnergy();

	}
	break;
	case 1:
	switch(diffDirection){
	case 0:
	//quark and diquark masses
	meMomenta()[2].setMass(mqq());
	meMomenta()[3].setMass(mq());

	//gluon constituent mass
	meMomenta()[4].setMass(getParticleData(21)->constituentMass());

	//outgoing proton
	meMomenta()[5].setVect(p4.vect());
	meMomenta()[5].setT(p4.t());

	//two body decay of the outgoing dissociation proton
	do{}
	while(!Kinematics::twoBodyDecay(p3,mqq()+mq(),getParticleData(21)->constituentMass(),
	p3.vect().unit(),momPair.first,momPair.second));
	//put gluon back-to-back with quark-diquark
	//set momenta of quark and diquark
	frac = mqq()/(mqq()+mq());
	meMomenta()[2].setVect(frac*momPair.first.vect());
	meMomenta()[2].setT(sqrt(sqr(frac)*momPair.first.vect().mag2()+sqr(mqq())));
	meMomenta()[3].setVect((1-frac)*momPair.first.vect());
	meMomenta()[3].setT(sqrt(sqr(1-frac)*momPair.first.vect().mag2()+sqr(mq())));
	//set momentum of gluon
	meMomenta()[4].setVect(momPair.second.vect());
	meMomenta()[4].setT(momPair.second.t());

	break;
	case 1:
	//quark and diquark masses
	meMomenta()[5].setMass(mqq());
	meMomenta()[4].setMass(mq());

	//gluon constituent mass
	meMomenta()[3].setMass(getParticleData(21)->constituentMass());

	//outgoing proton
	meMomenta()[2].setVect(p3.vect());
	meMomenta()[2].setT(p3.t());

	//two body decay of the outgoing dissociation proton
	do{}
	while(!Kinematics::twoBodyDecay(p4,mqq()+mq(),getParticleData(21)->constituentMass(),
	p4.vect().unit(),momPair.first,momPair.second));

	//put gluon back-to-back with quark-diquark
	//set momenta of quark and diquark
	frac = mqq()/(mqq()+mq());
	meMomenta()[5].setVect(frac*momPair.first.vect());
	meMomenta()[5].setT(sqrt(sqr(frac)*momPair.first.vect().mag2()+sqr(mqq())));
	meMomenta()[4].setVect((1-frac)*momPair.first.vect());
	meMomenta()[4].setT(sqrt(sqr(1-frac)*momPair.first.vect().mag2()+sqr(mq())));
	//set momentum of gluon
	meMomenta()[3].setVect(momPair.second.vect());
	meMomenta()[3].setT(momPair.second.t());



	break;
	case 2:
	//first dissociated proton constituents
	meMomenta()[2].setMass(mqq());
	meMomenta()[3].setMass(mq());
	meMomenta()[4].setMass(getParticleData(21)->constituentMass());
	//second dissociated proton constituents
	meMomenta()[5].setMass(getParticleData(21)->constituentMass());
	meMomenta()[6].setMass(mq());
	meMomenta()[7].setMass(mqq());


	//two body decay of the outgoing dissociation proton
	do{}
	while(!Kinematics::twoBodyDecay(p3,mqq()+mq(),getParticleData(21)->constituentMass(),
	p3.vect().unit(),momPair.first,momPair.second));

	do{}
	while(!Kinematics::twoBodyDecay(p4,mqq()+mq(),getParticleData(21)->constituentMass(),
	p4.vect().unit(),momPair1.first,momPair1.second));

	//put gluon back-to-back with quark-diquark
	frac = mqq()/(mqq()+mq());

	//first dissociated proton
	//set momenta of quark and diquark

	meMomenta()[2].setVect(frac*momPair.first.vect());
	meMomenta()[2].setT(sqrt(sqr(frac)*momPair.first.vect().mag2()+sqr(mqq())));
	meMomenta()[3].setVect((1-frac)*momPair.first.vect());
	meMomenta()[3].setT(sqrt(sqr(1-frac)*momPair.first.vect().mag2()+sqr(mq())));
	//set momentum of gluon
	meMomenta()[4].setVect(momPair.second.vect());
	meMomenta()[4].setT(momPair.second.t());

	//first dissociated proton
	//set momenta of quark and diquark

	meMomenta()[7].setVect(frac*momPair1.first.vect());
	meMomenta()[7].setT(sqrt(sqr(frac)*momPair1.first.vect().mag2()+sqr(mqq())));
	meMomenta()[6].setVect((1-frac)*momPair1.first.vect());
	meMomenta()[6].setT(sqrt(sqr(1-frac)*momPair1.first.vect().mag2()+sqr(mq())));
	//set momentum of gluon
	meMomenta()[5].setVect(momPair1.second.vect());
	meMomenta()[5].setT(momPair1.second.t());
	break;

	}
	meMomenta()[2].rescaleEnergy();
	meMomenta()[3].rescaleEnergy();
	meMomenta()[4].rescaleEnergy();
	meMomenta()[5].rescaleEnergy();
	if(diffDirection==2){
	meMomenta()[6].rescaleEnergy();
	meMomenta()[7].rescaleEnergy();
	}

	break;
	}

	jacobian(sqr(cmEnergy)/GeV2);
	return true;
	}
	//Generate masses of dissociated protons and momentum transfer from probability f(M2,t)
	//(for single diffraction). Sample according to f(M2,t)=f(M2)f(t\|M2).
	pair<pair<Energy2,Energy2>,Energy2> MEDiffraction::diffractiveMassAndMomentumTransfer() const {
	Energy2 theM12(ZERO),theM22(ZERO), thet(ZERO);
	int count = 0;
	//proton mass squared
	- const Energy2 m2 = sqr(getParticleData(2212)->mass());
	+ const Energy2 m2 = sqr(theProtonMass);
	//delta mass squared
	const Energy2 md2 = sqr(getParticleData(2214)->mass());
	Energy2 M2;
	bool condition = true;
	do {

	//check if we want only delta
	- if(deltaOnly)
	- {
	+ if(deltaOnly) {
	switch(diffDirection){
	case 0:
	theM12 = md2;
	theM22 = m2;
	M2 = md2;
	thet = randomt(md2);
	break;
	case 1:
	theM22 = md2;
	theM12 = m2;
	M2 = md2;
	thet = randomt(md2);
	break;
	case 2:
	theM12 = md2;
	theM22 = md2;
	M2 = md2;
	thet = doublediffrandomt(theM12,theM22);
	break;
	}

	- }else{
	+ }
	+ else {
	switch (diffDirection){
	case 0:
	M2=randomM2();
	thet = randomt(M2);
	theM12=M2;

	theM22=m2;

	break;
	case 1:

	theM12=m2;
	M2=randomM2();
	thet = randomt(M2);


	theM22=M2;
	break;
	case 2:
	theM12=randomM2();
	theM22=randomM2();
	M2=(theM12>theM22) ? theM12: theM22;

	thet = doublediffrandomt(theM12,theM22);

	break;
	}
	}
	count++;

	const Energy cmEnergy = generator()->maximumCMEnergy();
	const Energy2 s = sqr(cmEnergy);
	+ if(generator()->maximumCMEnergy()<sqrt(theM12)+sqrt(theM22)) {
	+ condition = true;
	+ }
	+ else {
	InvEnergy2 slope;
	if(diffDirection==2){
	slope = 2softPomeronSlope()log(.1+(sqr(cmEnergy)/softPomeronSlope())/(theM12*theM22));
	}else{
	slope = protonPomeronSlope()
	+ 2softPomeronSlope()log(sqr(cmEnergy)/M2);
	}


	const double expmax = exp(slope*tmaxfun(s,m2,M2));
	const double expmin = exp(slope*tminfun(s,m2,M2));


	//without (1-M2/s) constraint
	condition = (UseRandom::rnd()>(protonPomeronSlope()GeV2)(expmax-expmin)/(slope*GeV2))
	\|\|((theM12/GeV2)(theM22/GeV2)>=(sqr(cmEnergy)/GeV2)/(softPomeronSlope()GeV2));
	- }while(condition);
	+ }
	+ }
	+ while(condition);

	return make_pair (make_pair(theM12,theM22),thet);
	}

	//Decay of the excited proton to quark-diquark
	pair<Lorentz5Momentum,Lorentz5Momentum> MEDiffraction::twoBodyDecayMomenta(Lorentz5Momentum pp) const{
	//Decay of the excited proton
	const Energy2 Mx2(sqr(pp.mass())),mq2(sqr(mq())),mqq2(sqr(mqq()));

	const Energy2 psq = ((Mx2-sqr(mq()+mqq()))(Mx2-sqr(mq()-mqq())))/(4Mx2);

	assert(psq/GeV2>0);
	const Energy p(sqrt(psq));

	const double phi = UseRandom::rnd() * Constants::twopi;
	const double costheta =1-2*UseRandom::rnd();
	const double sintheta = sqrt(1-sqr(costheta));

	Lorentz5Momentum k1=Lorentz5Momentum(psinthetacos(phi), psinthetasin(phi), p*costheta, sqrt(mq2+psq));
	Lorentz5Momentum k2=Lorentz5Momentum(-psinthetacos(phi), -psinthetasin(phi), -p*costheta,sqrt(mqq2+psq));

	//find boost to pp center of mass
	const Boost betap3 = (pp).findBoostToCM();

	//k1 and k2 calculated at p3 center of mass, so boost back
	k1.boost(-betap3);
	k2.boost(-betap3);

	//first is quark, second diquark
	return make_pair(k1,k2);
	}



	Energy2 MEDiffraction::randomt(Energy2 M2) const {
	assert(protonPomeronSlope()*GeV2 > 0);
	//proton mass
	- const Energy2 m2 = sqr( getParticleData(2212)->mass() );
	+ const Energy2 m2 = sqr( theProtonMass );
	const Energy cmEnergy = generator()->maximumCMEnergy();
	const Energy2 ttmin = tminfun(sqr(cmEnergy),m2,M2);
	const Energy2 ttmax = tmaxfun(sqr(cmEnergy),m2,M2);

	const InvEnergy2 slope = protonPomeronSlope()
	+ 2softPomeronSlope()log(sqr(cmEnergy)/M2);
	return log( exp(slope*ttmin) +
	UseRandom::rnd()(exp(slopettmax) - exp(slope*ttmin)) ) / slope;
	}

	Energy2 MEDiffraction::doublediffrandomt(Energy2 M12, Energy2 M22) const {

	const Energy cmEnergy = generator()->maximumCMEnergy();

	const double shift = 0.1;
	const InvEnergy2 slope = 2softPomeronSlope()log(shift+(sqr(cmEnergy)/softPomeronSlope())/(M12*M22));

	const Energy2 ttmin = tminfun(sqr(cmEnergy),M12,M22);
	const Energy2 ttmax = tmaxfun(sqr(cmEnergy),M12,M22);
	-
	- return log( exp(slope*ttmin) +
	- UseRandom::rnd()(exp(slopettmax) - exp(slope*ttmin)) ) / slope;
	+ double r = UseRandom::rnd();
	+ Energy2 newVal;
	+ if(slopettmax>slopettmin) {
	+ newVal = ttmax + log( r + (1.-r)exp(slope(ttmin-ttmax)) ) / slope;
	+ }
	+ else {
	+ newVal = ttmin + log( 1. - r + rexp(slope(ttmax-ttmin))) / slope;
	+ }
	+ return newVal;
	}

	Energy2 MEDiffraction::randomM2() const {
	const double tmp = 1 - softPomeronIntercept();

	const Energy cmEnergy = generator()->maximumCMEnergy();
	return sqr(cmEnergy) * pow( pow(M2min()/sqr(cmEnergy),tmp) +
	UseRandom::rnd() * (pow(M2max()/sqr(cmEnergy),tmp) - pow(M2min()/sqr(cmEnergy),tmp)),
	1.0/tmp );
	}

	Energy2 MEDiffraction::tminfun(Energy2 s, Energy2 M12, Energy2 M22) const {
	- const Energy2 m2 = sqr( getParticleData(2212)->mass() );
	-
	- return (1/(2s))(-sqrt(kallen(s, m2, m2)kallen(s, M12, M22))-sqr(s)+2sm2+sM12+s*M22);
	+ const Energy2 m2 = sqr( theProtonMass );
	+ return 0.5/s(-sqrt(kallen(s, m2, m2)kallen(s, M12, M22))-sqr(s)+2sm2+sM12+sM22);
	}

	Energy2 MEDiffraction::tmaxfun(Energy2 s, Energy2 M12, Energy2 M22) const {
	- const Energy2 m2 = sqr( getParticleData(2212)->mass() );
	+ const Energy2 m2 = sqr( theProtonMass );

	- return (1/(2s))(sqrt(kallen(s, m2, m2)kallen(s, M12, M22))-sqr(s)+2sm2+sM12+s*M22);
	+ return 0.5/s(sqrt(kallen(s, m2, m2)kallen(s, M12, M22))-sqr(s)+2sm2+sM12+sM22);
	}


	double MEDiffraction::me2() const{
	return theme2;
	}

	CrossSection MEDiffraction::dSigHatDR() const {
	return me2()jacobian()/sHat()sqr(hbarc);
	}

	unsigned int MEDiffraction::orderInAlphaS() const {
	return 0;
	}

	unsigned int MEDiffraction::orderInAlphaEW() const {
	return 0;
	}

	Selector<MEBase::DiagramIndex>
	MEDiffraction::diagrams(const DiagramVector & diags) const {
	Selector<DiagramIndex> sel;
	if(!deltaOnly){
	if(diffDirection<2){

	for(unsigned int i = 0; i < diags.size(); i++){
	if(diags[0]->id()==-1)
	sel.insert(2./3.,i);
	else
	sel.insert(1./3.,i);
	}

	}else{
	for(unsigned int i = 0; i < diags.size(); i++){
	if(diags[0]->id()==-1)
	sel.insert(4./9.,i);
	else if(diags[0]->id()==-2)
	sel.insert(2./9.,i);
	else if(diags[0]->id()==-3)
	sel.insert(2./9.,i);
	else
	sel.insert(1./9.,i);
	}
	}
	}else{
	sel.insert(1.0,0);
	}

	return sel;
	}

	Selector<const ColourLines *>
	MEDiffraction::colourGeometries(tcDiagPtr ) const {
	Selector<const ColourLines *> sel;

	int sign1=0, sign2=0;
	sign1 = (generator()->eventHandler()->incoming().first->id() > 0) ? 1 : -1;
	sign2 = (generator()->eventHandler()->incoming().second->id() > 0) ? 1 : -1;

	switch(dissociationDecay){
	case 0:
	if(!deltaOnly)
	{
	if(diffDirection!=2){

	if (diffDirection == 0){
	if(sign1>0){
	static ColourLines dqq0=ColourLines("-6 2 7");
	sel.insert(1.0,&dqq0);
	}else{
	static ColourLines dqq0=ColourLines("6 -2 -7");
	sel.insert(1.0,&dqq0);
	}


	}
	else{
	if(sign2>0){
	static ColourLines dqq1=ColourLines("-6 3 7");
	sel.insert(1.0,&dqq1);
	}else{
	static ColourLines dqq1=ColourLines("6 -3 -7");
	sel.insert(1.0,&dqq1);
	}
	}

	}else{

	if(sign1>0 && sign2>0){
	static ColourLines ddqq0=ColourLines("-6 2 7, -9 4 8");
	sel.insert(1.0,&ddqq0);
	}else if(sign1<0 && sign2>0){
	static ColourLines ddqq0=ColourLines("6 -2 -7, -9 4 8");
	sel.insert(1.0,&ddqq0);
	}else if(sign1>0&& sign2<0){
	static ColourLines ddqq0=ColourLines("-6 2 7, 9 -4 -8");
	sel.insert(1.0,&ddqq0);
	}else{
	static ColourLines ddqq0=ColourLines("6 -2 -7, 9 -4 -8");
	sel.insert(1.0,&ddqq0);
	}


	}

	}else
	{
	static ColourLines cl("");

	sel.insert(1.0, &cl);
	}
	break;
	case 1:
	switch(diffDirection){
	case 0:
	static ColourLines clleft("-6 2 3 8, -8 -3 7");
	sel.insert(1.0, &clleft);
	break;
	case 1:
	static ColourLines clright("-9 4 3 7, -7 -3 8");
	sel.insert(1.0, &clright);
	break;
	case 2:
	static ColourLines cldouble("-8 2 3 10, -10 -3 9, -13 6 5 11, -11 -5 12");
	sel.insert(1.0, &cldouble);
	break;
	}
	break;
	}

	return sel;
	}

	+void MEDiffraction::doinit() {
	+ HwMEBase::doinit();
	+ theProtonMass = getParticleData(2212)->mass();
	+}
	+
	void MEDiffraction::doinitrun() {
	HwMEBase::doinitrun();
	isInRunPhase = true;
	}

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

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


	ClassDescription<MEDiffraction> MEDiffraction::initMEDiffraction;
	// Definition of the static class description member.

	void MEDiffraction::persistentOutput(PersistentOStream & os) const {
	os << theme2 << deltaOnly << diffDirection << theprotonPomeronSlope
	- << thesoftPomeronIntercept << thesoftPomeronSlope << dissociationDecay;
	+ << thesoftPomeronIntercept << thesoftPomeronSlope << dissociationDecay
	+ << ounit(theProtonMass,GeV);
	}

	void MEDiffraction::persistentInput(PersistentIStream & is, int) {
	is >> theme2 >> deltaOnly >> diffDirection >> theprotonPomeronSlope
	- >> thesoftPomeronIntercept >> thesoftPomeronSlope >> dissociationDecay;
	+ >> thesoftPomeronIntercept >> thesoftPomeronSlope >> dissociationDecay
	+ >> iunit(theProtonMass,GeV);
	}

	InvEnergy2 MEDiffraction::protonPomeronSlope() const{
	return theprotonPomeronSlope/GeV2;
	}

	double MEDiffraction::softPomeronIntercept() const {
	return thesoftPomeronIntercept;
	}

	InvEnergy2 MEDiffraction::softPomeronSlope() const {
	return thesoftPomeronSlope/GeV2;
	}

	void MEDiffraction::Init() {

	static ClassDocumentation<MEDiffraction> documentation
	("There is no documentation for the MEDiffraction class");

	static Parameter<MEDiffraction,double> interfaceme2
	("DiffractionAmplitude",
	"The square of the diffraction amplitude used to determine the "
	"cross section.",
	&MEDiffraction::theme2, 1.0, 0.00001, 100.0,
	false, false, Interface::limited);

	static Parameter<MEDiffraction,double> interfaceprotonPomeronSlope
	("ProtonPomeronSlope",
	"The proton-pomeron slope parameter.",
	&MEDiffraction::theprotonPomeronSlope, 10.1, 0.00001, 100.0,
	false, false, Interface::limited);

	static Parameter<MEDiffraction,double> interfacesoftPomeronIntercept
	("SoftPomeronIntercept",
	"The soft pomeron intercept.",
	&MEDiffraction::thesoftPomeronIntercept, 1.08, 0.00001, 100.0,
	false, false, Interface::limited);

	static Parameter<MEDiffraction,double> interfacesoftPomeronSlope
	("SoftPomeronSlope",
	"The soft pomeron slope parameter.",
	&MEDiffraction::thesoftPomeronSlope, 0.25, 0.00001, 100.0,
	false, false, Interface::limited);


	static Switch<MEDiffraction, bool> interfaceDeltaOnly
	("DeltaOnly",
	"proton-proton to proton-delta only",
	&MEDiffraction::deltaOnly, 0, false, false);
	static SwitchOption interfaceDeltaOnly0
	(interfaceDeltaOnly,"No","Final state with Delta only is OFF", 0);
	static SwitchOption interfaceDeltaOnly1
	(interfaceDeltaOnly,"Yes","Final state with Delta only is ON", 1);

	//Select if the left, right or both protons are excited
	static Switch<MEDiffraction, unsigned int> interfaceDiffDirection
	("DiffDirection",
	"Direction of the excited proton",
	&MEDiffraction::diffDirection, 0, false, false);
	static SwitchOption left
	(interfaceDiffDirection,"Left","Proton moving in the positive z direction", 0);
	static SwitchOption right
	(interfaceDiffDirection,"Right","Proton moving in the negative z direction", 1);
	static SwitchOption both
	(interfaceDiffDirection,"Both","Both protons", 2);

	//Select if two or three body decay
	static Switch<MEDiffraction, unsigned int> interfaceDissociationDecay
	("DissociationDecay",
	"Number of clusters the dissociated proton decays",
	&MEDiffraction::dissociationDecay, 0, false, false);
	static SwitchOption one
	(interfaceDissociationDecay,"One","Dissociated proton decays into one cluster", 0);
	static SwitchOption two
	(interfaceDissociationDecay,"Two","Dissociated proton decays into two clusters", 1);
	}
	diff --git a/MatrixElement/Hadron/MEDiffraction.h b/MatrixElement/Hadron/MEDiffraction.h
	--- a/MatrixElement/Hadron/MEDiffraction.h
	+++ b/MatrixElement/Hadron/MEDiffraction.h
	@@ -1,336 +1,344 @@
	// -- C++ --
	#ifndef HERWIG_MEDiffraction_H
	#define HERWIG_MEDiffraction_H
	//
	// This is the declaration of the MEDiffraction class.
	//

	#include "Herwig/MatrixElement/HwMEBase.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* The MEDiffraction class provides a simple colour singlet exchange matrix element
	* to be used in the soft component of the multiple scattering model of the
	* underlying event
	*
	* @see \ref MEDiffractionInterfaces "The interfaces"
	* defined for MEDiffraction.
	*/
	class MEDiffraction: public HwMEBase {

	public:

	MEDiffraction();

	/** @name Virtual functions required by the MEBase class. */
	//@{
	/**
	* Return the order in \f$\alpha_S\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInAlphaS() const;

	/**
	* Return the order in \f$\alpha_{EW}\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInAlphaEW() const;

	/**
	* The matrix element for the kinematical configuration
	* previously provided by the last call to setKinematics(), suitably
	* scaled by sHat() to give a dimension-less number.
	* @return the matrix element scaled with sHat() to give a
	* dimensionless number.
	*/
	virtual double me2() const;

	/**
	* Return the scale associated with the last set phase space point.
	*/
	virtual Energy2 scale() const;

	/**
	* Set the typed and momenta of the incoming and outgoing partons to
	* be used in subsequent calls to me() and colourGeometries()
	* according to the associated XComb object. If the function is
	* overridden in a sub class the new function must call the base
	* class one first.
	*/
	virtual void setKinematics();

	/**
	* The number of internal degrees of freedom used in the matrix
	* element.
	*/
	virtual int nDim() const;

	/**
	* Generate internal degrees of freedom given nDim() uniform
	* random numbers in the interval \f$ ]0,1[ \f$. To help the phase space
	* generator, the dSigHatDR should be a smooth function of these
	* numbers, although this is not strictly necessary.
	* @param r a pointer to the first of nDim() consecutive random numbers.
	* @return true if the generation succeeded, otherwise false.
	*/
	virtual bool generateKinematics(const double * r);

	/**
	* Return the matrix element squared differential in the variables
	* given by the last call to generateKinematics().
	*/
	virtual CrossSection dSigHatDR() const;

	/**
	* Add all possible diagrams with the add() function.
	*/
	virtual void getDiagrams() const;

	/**
	* Get diagram selector. With the information previously supplied with the
	* setKinematics method, a derived class may optionally
	* override this method to weight the given diagrams with their
	* (although certainly not physical) relative probabilities.
	* @param dv the diagrams to be weighted.
	* @return a Selector relating the given diagrams to their weights.
	*/
	virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;

	/**
	* Return a Selector with possible colour geometries for the selected
	* diagram weighted by their relative probabilities.
	* @param diag the diagram chosen.
	* @return the possible colour geometries weighted by their
	* relative probabilities.
	*/
	virtual Selector<const ColourLines *>
	colourGeometries(tcDiagPtr diag) const;
	//@}

	/**
	* Expect the incoming partons in the laboratory frame
	*/
	/* virtual bool wantCMS() const { return false; } */

	public:

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

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

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

	protected:

	+ /** @name Standard Interfaced functions. */
	+ //@{
	+ /**
	+ * Initialize this object after the setup phase before saving an
	+ * EventGenerator to disk.
	+ * @throws InitException if object could not be initialized properly.
	+ */
	+ virtual void doinit();
	+
	/**
	* Initialize this object. Called in the run phase just before a run begins.
	*/
	virtual void doinitrun();
	+ //@}

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	-
	-// If needed, insert declarations of virtual function defined in the
	-// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
	-
	-
	private:

	/* The matrix element squared */
	double theme2;

	/* Use only delta as excited state */
	bool deltaOnly;


	/* Direction of the excited proton */
	unsigned int diffDirection;

	/* Number of clusters the dissociated proton decays into */
	unsigned int dissociationDecay;

	/* The mass of the consitutent quark */
	Energy mq() const {return Energy(0.325*GeV);}

	/* The mass of the constituent diquark */
	Energy mqq() const {return Energy(0.650*GeV);}

	/* The proton-pomeron slope */
	double theprotonPomeronSlope;

	/* The soft pomeron intercept */
	double thesoftPomeronIntercept;

	/* The soft pomeron slope */
	double thesoftPomeronSlope;


	/**
	* Sample the diffractive mass squared M2 and the momentum transfer t
	*/
	pair<pair<Energy2,Energy2>,Energy2> diffractiveMassAndMomentumTransfer() const;


	/**
	* Random value for the diffractive mass squared M2 according to (M2/s0)^(-intercept)
	*/
	Energy2 randomM2() const;

	/**
	* Random value for t according to exp(diffSlope*t)
	*/
	Energy2 randomt(Energy2 M2) const;

	/**
	* Random value for t according to exp(diffSlope*t) for double diffraction
	*/

	Energy2 doublediffrandomt(Energy2 M12, Energy2 M22) const;


	/**
	* Returns the momenta of the two-body decay of momentum pp
	*/
	pair<Lorentz5Momentum,Lorentz5Momentum> twoBodyDecayMomenta(Lorentz5Momentum pp) const;

	/**
	* Returns the proton-pomeron slope
	*/
	InvEnergy2 protonPomeronSlope() const;

	/**
	* Returns the soft pomeron intercept
	*/
	double softPomeronIntercept() const;

	//M12 and M22 are masses squared of
	//outgoing particles

	/**
	* Returns the minimal possible value of momentum transfer t given the center
	* of mass energy and diffractive masses
	*/
	Energy2 tminfun(Energy2 s, Energy2 M12, Energy2 M22) const;

	/**
	* Returns the maximal possible value of momentum transfer t given the center
	* of mass energy and diffractive masses
	*/
	Energy2 tmaxfun(Energy2 s , Energy2 M12, Energy2 M22) const;

	/**
	* Returns the minimal possible value of diffractive mass
	*/
	//lowest possible mass given the constituent masses of quark and diquark
	Energy2 M2min() const{return sqr(getParticleData(2212)->mass()+mq()+mqq());}

	/**
	* Returns the maximal possible value of diffractive mass
	*/
	Energy2 M2max() const{
	return sqr(generator()->maximumCMEnergy()-getParticleData(2212)->mass());
	}//TODO:modify to get proper parameters

	InvEnergy2 softPomeronSlope() const;



	/* Kallen function */
	template<int L, int E, int Q, int DL, int DE, int DQ>
	Qty<2L,2E,2*Q,DL,DE,DQ> kallen(Qty<L,E,Q,DL,DE,DQ> a,
	Qty<L,E,Q,DL,DE,DQ> b,
	Qty<L,E,Q,DL,DE,DQ> c) const {
	return aa + bb + cc - 2.0(ab + bc + c*a);
	}



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

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

	+ bool isInRunPhase;

	- bool isInRunPhase;
	+
	+ /* The proton mass */
	+ Energy theProtonMass;

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the MEDiffraction class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::MEDiffraction>
	: public ClassTraitsBase<Herwig::MEDiffraction> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::MEDiffraction"; }
	/**
	* The name of a file containing the dynamic library where the class
	* MEDiffraction is implemented. It may also include several, space-separated,
	* libraries if the class MEDiffraction 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 "HwMEHadron.so";}
	};

	/** @endcond */

	}

	#endif /* HERWIG_MEDiffraction_H */
	diff --git a/src/LHC-MB-Soft.in b/src/LHC-MB-Soft.in
	--- a/src/LHC-MB-Soft.in
	+++ b/src/LHC-MB-Soft.in
	@@ -1,234 +1,234 @@
	# -- ThePEG-repository --

	################################################################################
	# This file contains our best tune to UE data from ATLAS at 7 TeV. More recent
	# tunes and tunes for other centre-of-mass energies as well as more usage
	# instructions can be obtained from this Herwig++ wiki page:
	# http://projects.hepforge.org/herwig/trac/wiki/MB_UE_tunes
	################################################################################

	read snippets/PPCollider.in

	##################################################
	# Technical parameters for this run
	##################################################
	cd /Herwig/Generators
	##################################################
	# LHC physics parameters (override defaults here)
	##################################################
	set EventGenerator:EventHandler:LuminosityFunction:Energy 7000.0

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


	##################################################
	# MEMinBias Matrix Element
	##################################################
	# MPI model settings
	set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0

	## Report the correct cross section
	cd /Herwig/Generators
	create Herwig::MPIXSecReweighter MPIXSecReweighter
	insert EventGenerator:EventHandler:PostSubProcessHandlers 0 MPIXSecReweighter
	set EventGenerator:EventHandler:CascadeHandler NULL

	clear EventGenerator:EventHandler:SubProcessHandlers[0]

	##################################################
	# Create separate SubProcessHandler for MinBias
	##################################################
	cd /Herwig/MatrixElements/
	cp SubProcess QCDMinBias

	set QCDMinBias:CascadeHandler /Herwig/Shower/ShowerHandler
	set QCDMinBias:CascadeHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
	set QCDMinBias:DecayHandler /Herwig/Decays/DecayHandler

	insert QCDMinBias:MatrixElements[0] MEMinBias

	cd /Herwig/Generators
	# MinBias parameters used for the new kinematics of soft MPI
	set /Herwig/Cuts/MinBiasCuts:X1Min 0.11
	set /Herwig/Cuts/MinBiasCuts:X2Min 0.11

	set EventGenerator:EventHandler:Cuts /Herwig/Cuts/MinBiasCuts

	set /Herwig/Partons/RemnantDecayer:colourDisrupt 0.0
	set /Herwig/Hadronization/ColourReconnector:ColourReconnection Yes
	# Use multiperipheral kinematics
	set /Herwig/Partons/RemnantDecayer:MultiPeriph Yes
	# Set ladder multiplicity factor
	set /Herwig/Partons/RemnantDecayer:ladderbFactor 0.0
	# Set gaussian width of longitudinal momentum fraction
	# flucutuation
	set /Herwig/Partons/RemnantDecayer:gaussWidth 0.03

	# Tuned to min-bias data
	set /Herwig/Hadronization/ColourReconnector:ReconnectionProbability 0.445557
	set /Herwig/UnderlyingEvent/MPIHandler:pTmin0 3.145333
	set /Herwig/UnderlyingEvent/MPIHandler:InvRadius 1.101560
	set /Herwig/UnderlyingEvent/MPIHandler:Power 0.709107
	set /Herwig/Partons/RemnantDecayer:ladderMult 0.700985

	cd /Herwig/MatrixElements/
	insert /Herwig/Generators/EventGenerator:EventHandler:SubProcessHandlers[0] QCDMinBias

	##################################################
	# Create separate SubProcessHandler for Diffraction
	##################################################

	cd /Herwig/MatrixElements
	create Herwig::MEDiffraction MEDiffractionLeft
	set MEDiffractionLeft:DiffDirection Left
	create Herwig::MEDiffraction MEDiffractionRight
	set MEDiffractionRight:DiffDirection Right
	create Herwig::MEDiffraction MEDiffractionDouble
	set MEDiffractionDouble:DiffDirection Both

	create Herwig::MEDiffraction MEDiffractionDeltaLeft
	set MEDiffractionDouble:DiffDirection Left
	create Herwig::MEDiffraction MEDiffractionDeltaRight
	set MEDiffractionDouble:DiffDirection Right
	create Herwig::MEDiffraction MEDiffractionDeltaDouble
	set MEDiffractionDouble:DiffDirection Both

	# Make a parton extractor for diffraction
	cd /Herwig/Partons
	-cp Extractor DiffPPExtractor
	+cp PPExtractor DiffPPExtractor
	set DiffPPExtractor:FirstPDF /Herwig/Partons/NoPDF
	set DiffPPExtractor:SecondPDF /Herwig/Partons/NoPDF

	cd /Herwig/MatrixElements/
	# Create Diffraction SubProcessHandler
	cp SubProcess QCDDiffraction

	# Assign the PartonExtractor to the SubProcessHandler
	set QCDDiffraction:PartonExtractor /Herwig/Partons/DiffPPExtractor

	# Use only Delta as final excited state (Yes/No)
	set MEDiffractionLeft:DeltaOnly No
	set MEDiffractionRight:DeltaOnly No
	set MEDiffractionDouble:DeltaOnly No

	set MEDiffractionDeltaLeft:DeltaOnly Yes
	set MEDiffractionDeltaRight:DeltaOnly Yes
	set MEDiffractionDeltaDouble:DeltaOnly Yes


	# Set weight for Diffraction
	set MEDiffractionLeft:DiffractionAmplitude 12
	set MEDiffractionRight:DiffractionAmplitude 12
	set MEDiffractionDouble:DiffractionAmplitude 8

	set MEDiffractionDeltaLeft:DiffractionAmplitude 4
	set MEDiffractionDeltaRight:DiffractionAmplitude 4
	set MEDiffractionDeltaDouble:DiffractionAmplitude 2

	# Set soft diffraction paramters
	set MEDiffractionLeft:ProtonPomeronSlope 10.1
	set MEDiffractionLeft:SoftPomeronIntercept 1.08
	set MEDiffractionLeft:SoftPomeronSlope 0.25

	set MEDiffractionRight:ProtonPomeronSlope 10.1
	set MEDiffractionRight:SoftPomeronIntercept 1.08
	set MEDiffractionRight:SoftPomeronSlope 0.25

	set MEDiffractionDouble:ProtonPomeronSlope 10.1
	set MEDiffractionDouble:SoftPomeronIntercept 1.08
	set MEDiffractionDouble:SoftPomeronSlope 0.25

	set MEDiffractionDeltaLeft:ProtonPomeronSlope 10.1
	set MEDiffractionDeltaLeft:SoftPomeronIntercept 1.08
	set MEDiffractionDeltaLeft:SoftPomeronSlope 0.25

	set MEDiffractionDeltaRight:ProtonPomeronSlope 10.1
	set MEDiffractionDeltaRight:SoftPomeronIntercept 1.08
	set MEDiffractionDeltaRight:SoftPomeronSlope 0.25

	set MEDiffractionDeltaDouble:ProtonPomeronSlope 10.1
	set MEDiffractionDeltaDouble:SoftPomeronIntercept 1.08
	set MEDiffractionDeltaDouble:SoftPomeronSlope 0.25

	# Set number of clusters for dissociation
	set MEDiffractionLeft:DissociationDecay One
	set MEDiffractionRight:DissociationDecay One
	set MEDiffractionDouble:DissociationDecay One

	set MEDiffractionDeltaLeft:DissociationDecay One
	set MEDiffractionDeltaRight:DissociationDecay One
	set MEDiffractionDeltaDouble:DissociationDecay One

	# Insert matrix elements
	insert QCDDiffraction:MatrixElements[0] MEDiffractionLeft
	insert QCDDiffraction:MatrixElements[0] MEDiffractionRight
	insert QCDDiffraction:MatrixElements[0] MEDiffractionDouble

	insert QCDDiffraction:MatrixElements[0] MEDiffractionDeltaLeft
	insert QCDDiffraction:MatrixElements[0] MEDiffractionDeltaRight
	insert QCDDiffraction:MatrixElements[0] MEDiffractionDeltaDouble
	# No cluster fission
	#set /Herwig/Hadronization/ClusterFissioner:ClMaxLight 3500

	set QCDDiffraction:CascadeHandler NULL
	set /Herwig/Generators/EventGenerator:EventHandler:CascadeHandler NULL

	insert /Herwig/Generators/EventGenerator:EventHandler:SubProcessHandlers[0] QCDDiffraction
	##########################################################################
	##########################################################################





	cd /Herwig/Analysis
	create ThePEG::RivetAnalysis RivetAnalysis RivetAnalysis.so

	cd /Herwig/Generators

	insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis

	#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 RivetATLAS_2012_I1084540
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2012_I1084540
	# need to comment for -d1
	###insert /Herwig/Analysis/RivetAnalysis:Analyses 0 CMS_2015_I1356998
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALICE_2012_I1181770

	#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2012_I1084540

	#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2010_S8591806
	#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 LHCB_2014_I1281685

	#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALICE_2010_S8624100
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2010_S8894728
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2011_S8994773
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2010_S8918562
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2011_I894867
	##insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2012_I1124167
	##insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2014_I1298811
	###insert /Herwig/Analysis/RivetAnalysis:Analyses 0 CMS_2010_S8656010
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 CMS_2011_S8884919
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 CMS_2011_S8978280

	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2012_I1091481
	#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 CMS_2012_I1193338
	##insert /Herwig/Analysis/RivetAnalysis:Analyses 0 CMS_2013_I1218372
	#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 CMS_QCD_10_024


	set /Herwig/Analysis/Plot:EventNumber 54
	cd /Herwig/Generators
	insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/Plot

	#insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
	#set /Herwig/Analysis/HepMCFile:PrintEvent 1000000
	#set /Herwig/Analysis/HepMCFile:Format GenEvent
	#set /Herwig/Analysis/HepMCFile:Units GeV_mm
	#set /Herwig/Analysis/HepMCFile:Filename events.fifo

	##################################################
	# Save run for later usage with 'Herwig++ run'
	##################################################
	saverun LHC-MB-Soft EventGenerator

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 6:21 PM (1 d, 21 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3790213
Default Alt Text: (49 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions