No OneTemporary
Actions

Size

226 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/Gamma/MEGammaGamma2ff.cc b/MatrixElement/Gamma/MEGammaGamma2ff.cc
	--- a/MatrixElement/Gamma/MEGammaGamma2ff.cc
	+++ b/MatrixElement/Gamma/MEGammaGamma2ff.cc
	@@ -1,275 +1,275 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the MEGammaGamma2ff class.
	//

	#include "MEGammaGamma2ff.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
	#include "ThePEG/Interface/Switch.h"
	#include "Herwig/Models/StandardModel/StandardModel.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "Herwig/MatrixElement/HardVertex.h"

	using namespace Herwig;

	MEGammaGamma2ff::MEGammaGamma2ff() : process_(0) {
	massOption(vector<unsigned int>(2,1));
	}

	void MEGammaGamma2ff::doinit() {
	HwMEBase::doinit();
	// cast the SM pointer to the Herwig SM pointer
	tcHwSMPtr hwsm=ThePEG::dynamic_ptr_cast<tcHwSMPtr>(standardModel());
	// do the initialisation
	if(!hwsm)
	throw InitException() << "Must be the Herwig StandardModel class in "
	<< "MEGammaGamma2ff::doinit" << Exception::abortnow;
	vertex_ = hwsm->vertexFFP();
	}

	void MEGammaGamma2ff::getDiagrams() const {
	tcPDPtr gamma = getParticleData(ParticleID::gamma);
	for(int ix=1;ix<17;++ix) {
	// increment counter to switch between quarks and leptons
	if(ix==7) ix+=4;
	if(ix>11&&ix%2==0) ++ix;
	// is it a valid quark process
	bool quark = ix<=6 && (process_==0 \|\| process_==1 \|\| process_-10==ix);
	// is it a valid lepton process
	bool lepton= ix>=11 && ix<=16 &&
	(process_==0 \|\| process_==2 \|\| ( (ix-9)/2 ==process_-4) );
	// only add diagram if allowed
	if(!quark && !lepton) continue;
	tcPDPtr lm = getParticleData(ix);
	tcPDPtr lp = lm->CC();
	// first t-channel
	add(new_ptr((Tree2toNDiagram(3),gamma,lp,gamma,1,lm, 2,lp,-1)));
	// interchange
	add(new_ptr((Tree2toNDiagram(3),gamma,lp,gamma,2,lm, 1,lp,-2)));
	}
	}

	Energy2 MEGammaGamma2ff::scale() const {
	return 2.sHat()tHat()*uHat()/(sqr(sHat())+sqr(tHat())+sqr(uHat()));
	}

	double MEGammaGamma2ff::me2() const {
	VectorWaveFunction p1w(rescaledMomenta()[0],mePartonData()[0],incoming);
	VectorWaveFunction p2w(rescaledMomenta()[1],mePartonData()[1],incoming);
	SpinorBarWaveFunction fw(rescaledMomenta()[2],mePartonData()[2],outgoing);
	SpinorWaveFunction fbarw(rescaledMomenta()[3],mePartonData()[3],outgoing);
	vector<VectorWaveFunction> p1,p2;
	vector<SpinorBarWaveFunction> f;
	vector<SpinorWaveFunction> fbar;
	for(unsigned int ix=0;ix<2;++ix) {
	p1w.reset(2*ix);p1.push_back(p1w);
	p2w.reset(2*ix);p2.push_back(p2w);
	fw.reset(ix);f.push_back(fw);
	fbarw.reset(ix);fbar.push_back(fbarw);
	}
	// calculate the matrix element
	return helicityME(p1,p2,f,fbar,false);
	}

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

	unsigned int MEGammaGamma2ff::orderInAlphaEW() const {
	return 2;
	}

	double MEGammaGamma2ff::helicityME(vector<VectorWaveFunction> &p1,
	vector<VectorWaveFunction> &p2,
	vector<SpinorBarWaveFunction> & f,
	vector<SpinorWaveFunction> & fbar, bool calc) const {
	// scale (external photons so scale in couplings is 0)
	Energy2 mt(0.*GeV2);
	// matrix element to be stored
	if(calc) me_.reset(ProductionMatrixElement(PDT::Spin1,PDT::Spin1,
	PDT::Spin1Half,PDT::Spin1Half));
	// calculate the matrix element
	double output(0.),sumdiag[2]={0.,0.};
	Complex diag[2];
	SpinorWaveFunction inters;
	for(unsigned int ihel1=0;ihel1<2;++ihel1) {
	for(unsigned int ihel2=0;ihel2<2;++ihel2) {
	for(unsigned int ohel1=0;ohel1<2;++ohel1) {
	for(unsigned int ohel2=0;ohel2<2;++ohel2) {
	//first t-channel diagram
	inters =vertex_->evaluate(mt,1,fbar[ohel2].particle()->CC(),
	fbar[ohel2],p2[ihel2]);
	diag[0]=vertex_->evaluate(mt,inters,f[ohel1],p1[ihel1]);
	//second t-channel diagram
	inters =vertex_->evaluate(mt,1,fbar[ohel2].particle()->CC(),
	fbar[ohel2],p1[ihel1]);
	diag[1]=vertex_->evaluate(mt,inters,f[ohel1],p2[ihel2]);
	sumdiag[0] += norm(diag[0]);
	sumdiag[1] += norm(diag[1]);
	diag[0] += diag[1];
	output += norm(diag[0]);
	// store the me if needed
	if(calc) me_(2ihel1,2ihel2,ohel1,ohel2) = diag[0];
	}
	}
	}
	}
	// diagrams
	DVector save;
	save.push_back(sumdiag[0]);
	save.push_back(sumdiag[1]);
	meInfo(save);
	// colour factors if needed
	if(mePartonData()[2]->coloured()) output *= 3.;
	- // code to test vs the analytic result
	-// Energy2 m2 = sqr(f[0].particle()->mass());
	-// Energy2 tm = (p1[0].getMomentum()+f [0].getMomentum()).m2()-m2;
	-// Energy2 um = (p1[0].getMomentum()+fbar[0].getMomentum()).m2()-m2;
	-// double test = 8.um/tm+8.tm/um- 32m2/tm - 32m2/um
	-// -32sqr(double(m2/tm)) - 64sqr(m2)/tm/um - 32*sqr(double(m2/um));
	-// test = sqr(4.Constants::pi*SM().alphaEM());
	-// if(mePartonData()[2]->coloured()) test *= 3.;
	-// cerr << "testing ME " << (output-test)/(output+test) << "\n";
	+ // // code to test vs the analytic result
	+ // Energy2 m2 = sqr(f[0].particle()->mass());
	+ // Energy2 tm = (p1[0].momentum()+f [0].momentum()).m2()-m2;
	+ // Energy2 um = (p1[0].momentum()+fbar[0].momentum()).m2()-m2;
	+ // double test = 8.um/tm+8.tm/um- 32m2/tm - 32m2/um
	+ // -32sqr(double(m2/tm)) - 64sqr(m2)/tm/um - 32*sqr(double(m2/um));
	+ // test = sqr(4.Constants::pi*SM().alphaEM());
	+ // if(mePartonData()[2]->coloured()) test *= 3.;
	+ // cerr << "testing ME " << (output-test)/(output+test) << "\n";
	// spin factors
	return 0.25*output;
	}

	Selector<MEBase::DiagramIndex>
	MEGammaGamma2ff::diagrams(const DiagramVector & diags) const {
	Selector<DiagramIndex> sel;
	for ( DiagramIndex i = 0; i < diags.size(); ++i )
	if ( diags[i]->id() == -1 ) sel.insert(meInfo()[0], i);
	else if ( diags[i]->id() == -2 ) sel.insert(meInfo()[1], i);
	return sel;
	}

	Selector<const ColourLines *>
	MEGammaGamma2ff::colourGeometries(tcDiagPtr ) const {
	static ColourLines c1("");
	static ColourLines c2("4 -2 -5");
	Selector<const ColourLines *> sel;
	if(mePartonData()[2]->coloured()) sel.insert(1.0, &c2);
	else sel.insert(1.0, &c1);
	return sel;
	}

	void MEGammaGamma2ff::persistentOutput(PersistentOStream & os) const {
	os << process_ << vertex_;
	}

	void MEGammaGamma2ff::persistentInput(PersistentIStream & is, int) {
	is >> process_ >> vertex_;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<MEGammaGamma2ff,HwMEBase>
	describeHerwigMEGammaGamma2ff("Herwig::MEGammaGamma2ff", "HwMEGammaGamma.so");

	void MEGammaGamma2ff::Init() {

	static ClassDocumentation<MEGammaGamma2ff> documentation
	("The MEGammaGamma2ff class implements the matrix elemments for"
	" direct fermion-antifermion production in gamma-gamma collisions.");

	static Switch<MEGammaGamma2ff,int> interfaceProcess
	("Process",
	"Which process to included",
	&MEGammaGamma2ff::process_, 0, false, false);
	static SwitchOption interfaceProcessAll
	(interfaceProcess,
	"All",
	"Include all SM fermions as outgoing particles",
	0);
	static SwitchOption interfaceProcessQuarks
	(interfaceProcess,
	"Quarks",
	"All include the quarks as outgoing particles",
	1);
	static SwitchOption interfaceProcessLeptons
	(interfaceProcess,
	"Leptons",
	"Only include the leptons as outgoing particles",
	2);
	static SwitchOption interfaceProcessElectron
	(interfaceProcess,
	"Electron",
	"Only include e+e- as outgoing particles",
	5);
	static SwitchOption interfaceProcessMuon
	(interfaceProcess,
	"Muon",
	"Only include mu+mu- as outgoing particles",
	6);
	static SwitchOption interfaceProcessTau
	(interfaceProcess,
	"Tau",
	"Only include tau+tau- as outgoing particles",
	7);
	static SwitchOption interfaceProcessDown
	(interfaceProcess,
	"Down",
	"Only include d dbar as outgoing particles",
	11);
	static SwitchOption interfaceProcessUp
	(interfaceProcess,
	"Up",
	"Only include u ubar as outgoing particles",
	12);
	static SwitchOption interfaceProcessStrange
	(interfaceProcess,
	"Strange",
	"Only include s sbar as outgoing particles",
	13);
	static SwitchOption interfaceProcessCharm
	(interfaceProcess,
	"Charm",
	"Only include c cbar as outgoing particles",
	14);
	static SwitchOption interfaceProcessBottom
	(interfaceProcess,
	"Bottom",
	"Only include b bbar as outgoing particles",
	15);
	static SwitchOption interfaceProcessTop
	(interfaceProcess,
	"Top",
	"Only include t tbar as outgoing particles",
	16);
	}

	void MEGammaGamma2ff::constructVertex(tSubProPtr sub) {
	// extract the particles in the hard process
	ParticleVector hard;
	hard.push_back(sub->incoming().first);
	hard.push_back(sub->incoming().second);
	hard.push_back(sub->outgoing()[0]);
	hard.push_back(sub->outgoing()[1]);
	// order of particles
	unsigned int order[4]={0,1,2,3};
	// identify the process and calculate matrix element
	vector<VectorWaveFunction> p1,p2;
	if(hard[order[2]]->id()<0) swap(order[2],order[3]);
	vector<SpinorWaveFunction> fbar;
	vector<SpinorBarWaveFunction> f;
	VectorWaveFunction (p1 ,hard[order[0]],incoming,false,true);
	VectorWaveFunction (p2 ,hard[order[1]],incoming,false,true);
	SpinorBarWaveFunction(f ,hard[order[2]],outgoing,true );
	SpinorWaveFunction (fbar,hard[order[3]],outgoing,true );
	p1[1]=p1[2];
	p2[1]=p2[2];
	helicityME(p1,p2,f,fbar,true);
	// construct the vertex
	HardVertexPtr hardvertex=new_ptr(HardVertex());
	// set the matrix element for the vertex
	hardvertex->ME(me_);
	// set the pointers and to and from the vertex
	for(unsigned int ix=0;ix<4;++ix)
	hard[order[ix]]->spinInfo()->productionVertex(hardvertex);
	}
	diff --git a/PDF/HwRemDecayer.cc b/PDF/HwRemDecayer.cc
	--- a/PDF/HwRemDecayer.cc
	+++ b/PDF/HwRemDecayer.cc
	@@ -1,1536 +1,1537 @@
	// -- C++ --
	//
	// HwRemDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the HwRemDecayer class.
	//

	#include "HwRemDecayer.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Utilities/UtilityBase.h"
	#include "ThePEG/Utilities/SimplePhaseSpace.h"
	#include "ThePEG/Utilities/Throw.h"

	#include "Herwig/Shower/ShowerHandler.h"

	using namespace Herwig;

	namespace{

	const bool dbg = false;


	void reShuffle(Lorentz5Momentum &p1, Lorentz5Momentum &p2, Energy m1, Energy m2){

	Lorentz5Momentum ptotal(p1+p2);
	ptotal.rescaleMass();

	if( ptotal.m() < m1+m2 ) {
	if(dbg)
	cerr << "Not enough energy to perform reshuffling \n";
	throw HwRemDecayer::ExtraSoftScatterVeto();
	}

	Boost boostv = -ptotal.boostVector();
	ptotal.boost(boostv);

	p1.boost(boostv);
	// set the masses and energies,
	p1.setMass(m1);
	p1.setE(0.5/ptotal.m()*(ptotal.m2()+sqr(m1)-sqr(m2)));
	p1.rescaleRho();
	// boost back to the lab
	p1.boost(-boostv);

	p2.boost(boostv);
	// set the masses and energies,
	p2.setMass(m2);
	p2.setE(0.5/ptotal.m()*(ptotal.m2()+sqr(m2)-sqr(m1)));
	p2.rescaleRho();
	// boost back to the lab
	p2.boost(-boostv);
	}
	}

	void HwRemDecayer::initialize(pair<tRemPPtr, tRemPPtr> rems, tPPair beam, Step & step,
	Energy forcedSplitScale) {
	// the step
	thestep = &step;
	// valence content of the hadrons
	theContent.first = getHadronContent(beam.first);
	theContent.second = getHadronContent(beam.second);
	// momentum extracted from the hadrons
	theUsed.first = Lorentz5Momentum();
	theUsed.second = Lorentz5Momentum();
	theMaps.first.clear();
	theMaps.second.clear();
	theX.first = 0.0;
	theX.second = 0.0;
	theRems = rems;
	_forcedSplitScale = forcedSplitScale;
	// check remnants attached to the right hadrons
	if( (theRems.first && parent(theRems.first ) != beam.first ) \|\|
	(theRems.second && parent(theRems.second) != beam.second) )
	throw Exception() << "Remnant order wrong in "
	<< "HwRemDecayer::initialize(...)"
	<< Exception::runerror;
	return;
	}

	void HwRemDecayer::split(tPPtr parton, HadronContent & content,
	tRemPPtr rem, Lorentz5Momentum & used,
	PartnerMap &partners, tcPDFPtr pdf, bool first) {
	theBeam = parent(rem);
	theBeamData = dynamic_ptr_cast<Ptr<BeamParticleData>::const_pointer>
	(theBeam->dataPtr());
	double currentx = parton->momentum().rho()/theBeam->momentum().rho();
	double check = rem==theRems.first ? theX.first : theX.second;
	check += currentx;
	if(1.0-check < 1e-3) throw ShowerHandler::ExtraScatterVeto();
	bool anti;
	Lorentz5Momentum lastp(parton->momentum());
	int lastID(parton->id());
	Energy oldQ(_forcedSplitScale);
	_pdf = pdf;
	//do nothing if already valence quark
	if(first && content.isValenceQuark(parton)) {
	//set the extracted value, because otherwise no RemID could be generated.
	content.extract(lastID);
	// add the particle to the colour partners
	partners.push_back(make_pair(parton, tPPtr()));
	//set the sign
	anti = parton->hasAntiColour() && parton->id()!=ParticleID::g;
	if(rem==theRems.first) theanti.first = anti;
	else theanti.second = anti;
	// add the x and return
	if(rem==theRems.first) theX.first += currentx;
	else theX.second += currentx;
	return;
	}
	//or gluon for secondaries
	else if(!first && lastID == ParticleID::g) {
	partners.push_back(make_pair(parton, tPPtr()));
	// add the x and return
	if(rem==theRems.first) theX.first += currentx;
	else theX.second += currentx;
	return;
	}
	// if a sea quark.antiquark forced splitting to a gluon
	// Create the new parton with its momentum and parent/child relationship set
	PPtr newSea;
	if( !(lastID == ParticleID::g \|\|
	lastID == ParticleID::gamma) ) {
	newSea = forceSplit(rem, -lastID, oldQ, currentx, lastp, used,content);
	ColinePtr cl = new_ptr(ColourLine());
	if(newSea->id() > 0) cl-> addColoured(newSea);
	else cl->addAntiColoured(newSea);
	// if a secondard scatter finished so return
	if(!first \|\| content.isValenceQuark(ParticleID::g) ){
	partners.push_back(make_pair(parton, newSea));
	// add the x and return
	if(rem==theRems.first) theX.first += currentx;
	else theX.second += currentx;
	if(first) content.extract(ParticleID::g);
	return;
	}
	}
	// otherwise evolve back to valence
	// final valence splitting
	PPtr newValence = forceSplit(rem,
	lastID!=ParticleID::gamma ?
	ParticleID::g : ParticleID::gamma,
	oldQ, currentx , lastp, used, content);
	// extract from the hadron to allow remnant to be determined
	content.extract(newValence->id());
	// case of a gluon going into the hard subprocess
	if( lastID == ParticleID::g ) {
	partners.push_back(make_pair(parton, tPPtr()));
	anti = newValence->hasAntiColour();
	if(rem==theRems.first) theanti.first = anti;
	else theanti.second = anti;
	parton->colourLine(!anti)->addColoured(newValence, anti);
	return;
	}
	else if( lastID == ParticleID::gamma) {
	partners.push_back(make_pair(parton, newValence));
	anti = newValence->hasAntiColour();
	ColinePtr newLine(new_ptr(ColourLine()));
	newLine->addColoured(newValence, anti);
	if(rem==theRems.first) theanti.first = anti;
	else theanti.second = anti;
	// add the x and return
	if(rem==theRems.first) theX.first += currentx;
	else theX.second += currentx;
	return;
	}
	//The valence quark will always be connected to the sea quark with opposite sign
	tcPPtr particle;
	if(lastID*newValence->id() < 0){
	particle = parton;
	partners.push_back(make_pair(newSea, tPPtr()));
	}
	else {
	particle = newSea;
	partners.push_back(make_pair(parton, tPPtr()));
	}
	anti = newValence->hasAntiColour();
	if(rem==theRems.first) theanti.first = anti;
	else theanti.second = anti;

	if(particle->colourLine())
	particle->colourLine()->addAntiColoured(newValence);
	if(particle->antiColourLine())
	particle->antiColourLine()->addColoured(newValence);
	// add the x and return
	if(rem==theRems.first) theX.first += currentx;
	else theX.second += currentx;
	return;
	}

	void HwRemDecayer::doSplit(pair<tPPtr, tPPtr> partons,
	pair<tcPDFPtr, tcPDFPtr> pdfs,
	bool first) {
	if(theRems.first) {
	ParticleVector children=theRems.first->children();
	for(unsigned int ix=0;ix<children.size();++ix) {
	if(children[ix]->dataPtr()==theRems.first->dataPtr())
	theRems.first = dynamic_ptr_cast<RemPPtr>(children[ix]);
	}
	}
	if(theRems.second) {
	ParticleVector children=theRems.second->children();
	for(unsigned int ix=0;ix<children.size();++ix) {
	if(children[ix]->dataPtr()==theRems.second->dataPtr())
	theRems.second = dynamic_ptr_cast<RemPPtr>(children[ix]);
	}
	}
	// forced splitting for first parton
	if(isPartonic(partons.first )) {
	try {
	split(partons.first, theContent.first, theRems.first,
	theUsed.first, theMaps.first, pdfs.first, first);
	}
	catch(ShowerHandler::ExtraScatterVeto) {
	throw ShowerHandler::ExtraScatterVeto();
	}
	}
	// forced splitting for second parton
	if(isPartonic(partons.second)) {
	try {
	split(partons.second, theContent.second, theRems.second,
	theUsed.second, theMaps.second, pdfs.second, first);
	// additional check for the remnants
	// if can't do the rescale veto the emission
	if(!first&&partons.first->data().coloured()&&
	partons.second->data().coloured()) {
	Lorentz5Momentum pnew[2]=
	{theRems.first->momentum() - theUsed.first - partons.first->momentum(),
	theRems.second->momentum() - theUsed.second - partons.second->momentum()};

	pnew[0].setMass(getParticleData(theContent.first.RemID())->constituentMass());
	pnew[0].rescaleEnergy();
	pnew[1].setMass(getParticleData(theContent.second.RemID())->constituentMass());
	pnew[1].rescaleEnergy();

	for(unsigned int iy=0; iy<theRems.first->children().size(); ++iy)
	pnew[0] += theRems.first->children()[iy]->momentum();

	for(unsigned int iy=0; iy<theRems.second->children().size(); ++iy)
	pnew[1] += theRems.second->children()[iy]->momentum();

	Lorentz5Momentum ptotal=
	theRems.first ->momentum()-partons.first ->momentum()+
	theRems.second->momentum()-partons.second->momentum();
	// add x limits
	if(ptotal.m() < (pnew[0].m() + pnew[1].m()) ) {
	if(partons.second->id() != ParticleID::g){
	if(partons.second==theMaps.second.back().first)
	theUsed.second -= theMaps.second.back().second->momentum();
	else
	theUsed.second -= theMaps.second.back().first->momentum();

	thestep->removeParticle(theMaps.second.back().first);
	thestep->removeParticle(theMaps.second.back().second);
	}
	theMaps.second.pop_back();
	theX.second -= partons.second->momentum().rho()/
	parent(theRems.second)->momentum().rho();
	throw ShowerHandler::ExtraScatterVeto();
	}
	}
	}
	catch(ShowerHandler::ExtraScatterVeto){
	if(!partons.first\|\|!partons.second\|\|
	!theRems.first\|\|!theRems.second)
	throw ShowerHandler::ExtraScatterVeto();
	//case of the first forcedSplitting worked fine
	theX.first -= partons.first->momentum().rho()/
	parent(theRems.first)->momentum().rho();
	//case of the first interaction
	//throw veto immediately, because event get rejected anyway.
	if(first) throw ShowerHandler::ExtraScatterVeto();
	//secondary interactions have to end on a gluon, if parton
	//was NOT a gluon, the forced splitting particles must be removed
	if(partons.first->id() != ParticleID::g) {
	if(partons.first==theMaps.first.back().first)
	theUsed.first -= theMaps.first.back().second->momentum();
	else
	theUsed.first -= theMaps.first.back().first->momentum();

	thestep->removeParticle(theMaps.first.back().first);
	thestep->removeParticle(theMaps.first.back().second);
	}
	theMaps.first.pop_back();
	throw ShowerHandler::ExtraScatterVeto();
	}
	}
	// veto if not enough energy for extraction
	if( !first &&(theRems.first ->momentum().e() -
	partons.first ->momentum().e() < 1.0e-3*MeV \|\|
	theRems.second->momentum().e() -
	partons.second->momentum().e() < 1.0e-3*MeV )) {
	if(partons.first->id() != ParticleID::g) {
	if(partons.first==theMaps.first.back().first)
	theUsed.first -= theMaps.first.back().second->momentum();
	else
	theUsed.first -= theMaps.first.back().first->momentum();
	thestep->removeParticle(theMaps.first.back().first);
	thestep->removeParticle(theMaps.first.back().second);
	}
	theMaps.first.pop_back();
	if(partons.second->id() != ParticleID::g) {
	if(partons.second==theMaps.second.back().first)
	theUsed.second -= theMaps.second.back().second->momentum();
	else
	theUsed.second -= theMaps.second.back().first->momentum();
	thestep->removeParticle(theMaps.second.back().first);
	thestep->removeParticle(theMaps.second.back().second);
	}
	theMaps.second.pop_back();
	throw ShowerHandler::ExtraScatterVeto();
	}
	}

	void HwRemDecayer::mergeColour(tPPtr pold, tPPtr pnew, bool anti) const {
	ColinePtr clnew, clold;

	//save the corresponding colour lines
	clold = pold->colourLine(anti);
	clnew = pnew->colourLine(!anti);

	assert(clold);

	// There is already a colour line (not the final diquark)
	if(clnew){

	if( (clnew->coloured().size() + clnew->antiColoured().size()) > 1 ){
	if( (clold->coloured().size() + clold->antiColoured().size()) > 1 ){
	//join the colour lines
	//I don't use the join method, because potentially only (anti)coloured
	//particles belong to one colour line
	if(clold!=clnew){//procs are not already connected
	while ( !clnew->coloured().empty() ) {
	tPPtr p = clnew->coloured()[0];
	clnew->removeColoured(p);
	clold->addColoured(p);
	}
	while ( !clnew->antiColoured().empty() ) {
	tPPtr p = clnew->antiColoured()[0];
	clnew->removeAntiColoured(p);
	clold->addAntiColoured(p);
	}
	}

	}else{
	//if pold is the only member on it's
	//colour line, remove it.
	clold->removeColoured(pold, anti);
	//and add it to clnew
	clnew->addColoured(pold, anti);
	}
	} else{//pnnew is the only member on it's colour line.
	clnew->removeColoured(pnew, !anti);
	clold->addColoured(pnew, !anti);
	}
	} else {//there is no coline at all for pnew
	clold->addColoured(pnew, !anti);
	}
	}

	void HwRemDecayer::fixColours(PartnerMap partners, bool anti,
	double colourDisrupt) const {
	PartnerMap::iterator prev;
	tPPtr pnew, pold;
	assert(partners.size()>=2);

	PartnerMap::iterator it=partners.begin();
	while(it != partners.end()) {
	//skip the first one to have a partner
	if(it==partners.begin()){
	it++;
	continue;
	}

	prev = it - 1;
	//determine the particles to work with
	pold = prev->first;
	if(prev->second) {
	if(!pold->coloured())
	pold = prev->second;
	else if(pold->hasAntiColour() != anti)
	pold = prev->second;
	}
	assert(pold);

	pnew = it->first;
	if(it->second) {
	if(it->second->colourLine(!anti)) //look for the opposite colour
	pnew = it->second;
	}
	assert(pnew);

	// Implement the disruption of colour connections
	if( it != partners.end()-1 ) {//last one is diquark-has to be connected

	//has to be inside the if statement, so that the probability is
	//correctly counted:
	if( UseRandom::rnd() < colourDisrupt ){

	if(!it->second){//check, whether we have a gluon
	mergeColour(pnew, pnew, anti);
	}else{
	if(pnew==it->first)//be careful about the order
	mergeColour(it->second, it->first, anti);
	else
	mergeColour(it->first, it->second, anti);
	}

	it = partners.erase(it);
	continue;
	}

	}
	// regular merging
	mergeColour(pold, pnew, anti);
	//end of loop
	it++;
	}
	return;
	}

	PPtr HwRemDecayer::forceSplit(const tRemPPtr rem, long child, Energy &lastQ,
	double &lastx, Lorentz5Momentum &pf,
	Lorentz5Momentum &p,
	HadronContent & content) const {
	static const double eps=1e-6;
	// beam momentum
	Lorentz5Momentum beam = theBeam->momentum();
	// the last scale is minimum of last value and upper limit
	Energy minQ=_range_kinCutoffsqrt(lastx)/(1-lastx);
	if(minQ>lastQ) lastQ=minQ;
	// generate the new value of qtilde
	// weighted towards the lower value: dP/dQ = 1/Q -> Q(R) =
	// Q0 (Qmax/Q0)^R
	Energy q;
	unsigned int ntry=0,maxtry=100;
	double xExtracted = rem==theRems.first ? theX.first : theX.second;
	double zmin= lastx/(1.-xExtracted) ,zmax,yy;

	if(1-lastx<eps) throw ShowerHandler::ExtraScatterVeto();
	do {
	q = minQ*pow(lastQ/minQ,UseRandom::rnd());
	yy = 1.+0.5*sqr(_kinCutoff/q);
	zmax = yy - sqrt(sqr(yy)-1.);
	++ntry;
	}
	while(zmax<zmin&&ntry<maxtry);
	if(ntry==maxtry) throw ShowerHandler::ExtraScatterVeto();
	if(zmax-zmin<eps) throw ShowerHandler::ExtraScatterVeto();
	// now generate z as in FORTRAN HERWIG
	// use y = ln(z/(1-z)) as integration variable
	double ymin=log(zmin/(1.-zmin));
	double ymax=log(zmax/(1.-zmax));
	double dely=ymax-ymin;
	unsigned int nz=_nbinmax;
	dely/=nz;
	yy=ymin+0.5*dely;
	vector<int> ids;
	if(child==21\|\|child==22) {
	ids=content.flav;
	for(unsigned int ix=0;ix<ids.size();++ix) ids[ix] *= content.sign;
	}
	else {
	ids.push_back(ParticleID::g);
	}
	// probabilities of the different types of possible splitting
	map<long,pair<double,vector<double> > > partonprob;
	double ptotal(0.);
	for(unsigned int iflav=0;iflav<ids.size();++iflav) {
	// only do each parton once
	if(partonprob.find(ids[iflav])!=partonprob.end()) continue;
	// particle data object
	tcPDPtr in = getParticleData(ids[iflav]);
	double psum(0.);
	vector<double> prob;
	for(unsigned int iz=0;iz<nz;++iz) {
	double ez=exp(yy);
	double wr=1.+ez;
	double zr=wr/ez;
	double wz=1./wr;
	double zz=wz*ez;
	double coup = child!=22 ?
	_alphaS ->value(sqr(max(wz*q,_kinCutoff))) :
	_alphaEM->value(sqr(max(wz*q,_kinCutoff)));
	double az=wzzzcoup;
	// g -> q qbar
	if(ids[iflav]==ParticleID::g) {
	// calculate splitting function
	// SP as q is always less than forcedSplitScale, the pdf scale is fixed
	// pdfval = _pdf->xfx(theBeamData,in,sqr(q),lastx*zr);
	double pdfval=_pdf->xfx(theBeamData,in,sqr(_forcedSplitScale),lastx*zr);
	if(pdfval>0.) psum += pdfvalaz0.5*(sqr(zz)+sqr(wz));
	}
	// q -> q g
	else {
	// calculate splitting function
	// SP as q is always less than forcedSplitScale, the pdf scale is fixed
	// pdfval = _pdf->xfx(theBeamData,in,sqr(q),lastx*zr);
	double pdfval=_pdf->xfx(theBeamData,in,sqr(_forcedSplitScale),lastx*zr);
	if(pdfval>0.) psum += pdfvalaz4./3.(1.+sqr(wz))zr;
	}
	if(psum>0.) prob.push_back(psum);
	yy+=dely;
	}
	if(psum>0.) partonprob[ids[iflav]] = make_pair(psum,prob);
	ptotal+=psum;
	}
	// select the flavour
	if(ptotal==0.) throw ShowerHandler::ExtraScatterVeto();
	ptotal *= UseRandom::rnd();
	map<long,pair<double,vector<double> > >::const_iterator pit;
	for(pit=partonprob.begin();pit!=partonprob.end();++pit) {
	if(pit->second.first>=ptotal) break;
	else ptotal -= pit->second.first;
	}
	if(pit==partonprob.end())
	throw Exception() << "Can't select parton for forced backward evolution in "
	<< "HwRemDecayer::forceSplit" << Exception::eventerror;
	// select z
	unsigned int iz=0;
	for(;iz<pit->second.second.size();++iz) {
	if(pit->second.second[iz]>ptotal) break;
	}
	if(iz==pit->second.second.size()) --iz;
	double ey=exp(ymin+dely*(float(iz+1)-UseRandom::rnd()));
	double z=ey/(1.+ey);
	Energy2 pt2=sqr((1.-z)q)- zsqr(_kinCutoff);
	// create the particle
	if(pit->first!=ParticleID::g) child=pit->first;
	PPtr parton = getParticleData(child)->produceParticle();
	Energy2 emittedm2 = sqr(parton->dataPtr()->constituentMass());
	// Now boost pcm and pf to z only frame
	Lorentz5Momentum p_ref = Lorentz5Momentum(ZERO, beam.vect());
	Lorentz5Momentum n_ref = Lorentz5Momentum(ZERO, -beam.vect());
	// generate phi and compute pt of branching
	double phi = Constants::twopi*UseRandom::rnd();
	Energy pt=sqrt(pt2);
	Lorentz5Momentum qt = LorentzMomentum(ptcos(phi), ptsin(phi), ZERO, ZERO);
	Axis axis(p_ref.vect().unit());
	if(axis.perp2()>0.) {
	LorentzRotation rot;
	double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
	rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
	qt.transform(rot);
	}
	// compute alpha for previous particle
	Energy2 p_dot_n = p_ref*n_ref;
	double lastalpha = pf*n_ref/p_dot_n;
	Lorentz5Momentum qtout=qt;
	Energy2 qtout2=-qt*qt;
	double alphaout=(1.-z)/z*lastalpha;
	double betaout=0.5*(emittedm2+qtout2)/alphaout/p_dot_n;
	Lorentz5Momentum k=alphaoutp_ref+betaoutn_ref+qtout;
	k.rescaleMass();
	parton->set5Momentum(k);
	pf+=k;
	lastQ=q;
	lastx/=z;
	p += parton->momentum();
	thestep->addDecayProduct(rem,parton,false);
	return parton;
	}

	void HwRemDecayer::setRemMasses() const {
	// get the masses of the remnants
	Energy mrem[2];
	Lorentz5Momentum ptotal,pnew[2];
	vector<tRemPPtr> theprocessed;
	theprocessed.push_back(theRems.first);
	theprocessed.push_back(theRems.second);
	// one remnant in e.g. DIS
	if(!theprocessed[0]\|\|!theprocessed[1]) {
	tRemPPtr rem = theprocessed[0] ? theprocessed[0] : theprocessed[1];
	+ int remid = theprocessed[0] ? theContent.first.RemID() : theContent.second.RemID();
	Lorentz5Momentum deltap(rem->momentum());
	// find the diquark and momentum we still need in the energy
	tPPtr diquark;
	vector<PPtr> progenitors;
	for(unsigned int ix=0;ix<rem->children().size();++ix) {
	- if(!DiquarkMatcher::Check(rem->children()[ix]->data())) {
	+ if(rem->children()[ix]->id()!=remid) {
	progenitors.push_back(rem->children()[ix]);
	deltap -= rem->children()[ix]->momentum();
	}
	else
	diquark = rem->children()[ix];
	}
	// now find the total momentum of the hadronic final-state to
	// reshuffle against
	// find the hadron for this remnant
	tPPtr hadron=rem;
	do hadron=hadron->parents()[0];
	while(!hadron->parents().empty());
	// find incoming parton to hard process from this hadron
	tPPtr hardin =
	generator()->currentEvent()->primaryCollision()->incoming().first==hadron ?
	generator()->currentEvent()->primarySubProcess()->incoming().first :
	generator()->currentEvent()->primarySubProcess()->incoming().second;
	tPPtr parent=hardin;
	vector<PPtr> tempprog;
	// find the outgoing particles emitted from the backward shower
	do {
	assert(!parent->parents().empty());
	tPPtr newparent=parent->parents()[0];
	if(newparent==hadron) break;
	for(unsigned int ix=0;ix<newparent->children().size();++ix) {
	if(newparent->children()[ix]!=parent)
	findChildren(newparent->children()[ix],tempprog);
	}
	parent=newparent;
	}
	while(parent!=hadron);
	// add to list of potential particles to reshuffle against in right order
	for(unsigned int ix=tempprog.size();ix>0;--ix) progenitors.push_back(tempprog[ix-1]);
	// final-state particles which are colour connected
	tColinePair lines = make_pair(hardin->colourLine(),hardin->antiColourLine());
	vector<PPtr> others;
	for(ParticleVector::const_iterator
	cit = generator()->currentEvent()->primarySubProcess()->outgoing().begin();
	cit!= generator()->currentEvent()->primarySubProcess()->outgoing().end();++cit) {
	// colour connected
	if(lines.first&&lines.first==(**cit).colourLine()) {
	findChildren(*cit,progenitors);
	continue;
	}
	// anticolour connected
	if(lines.second&&lines.second==(**cit).antiColourLine()) {
	findChildren(*cit,progenitors);
	continue;
	}
	// not connected
	for(unsigned int ix=0;ix<(**cit).children().size();++ix)
	others.push_back((**cit).children()[ix]);
	}
	// work out how much of the system needed for rescaling
	unsigned int iloc=0;
	Lorentz5Momentum psystem,ptotal;
	do {
	psystem+=progenitors[iloc]->momentum();
	ptotal = psystem + deltap;
	ptotal.rescaleMass();
	psystem.rescaleMass();
	++iloc;
	if(ptotal.mass() > psystem.mass() + diquark->mass() &&
	psystem.mass()>1MeV && DISRemnantOpt_<2 && ptotal.e() > 0.GeV ) break;
	}
	while(iloc<progenitors.size());
	if(ptotal.mass() > psystem.mass() + diquark->mass()) --iloc;
	if(iloc==progenitors.size()) {
	// try touching the lepton in dis as a last restort
	for(unsigned int ix=0;ix<others.size();++ix) {
	progenitors.push_back(others[ix]);
	psystem+=progenitors[iloc]->momentum();
	ptotal = psystem + deltap;
	ptotal.rescaleMass();
	psystem.rescaleMass();
	++iloc;
	}
	--iloc;
	if(ptotal.mass() > psystem.mass() + diquark->mass()) {
	if(DISRemnantOpt_==0\|\|DISRemnantOpt_==2)
	Throw<Exception>() << "Warning had to adjust the momentum of the"
	<< " non-colour connected"
	<< " final-state, e.g. the scattered lepton in DIS"
	<< Exception::warning;
	else
	throw Exception() << "Can't set remnant momentum without adjusting "
	<< "the momentum of the"
	<< " non-colour connected"
	<< " final-state, e.g. the scattered lepton in DIS"
	<< " vetoing event"
	<< Exception::eventerror;
	}
	else {
	throw Exception() << "Can't put the remnant on-shell in HwRemDecayer::setRemMasses()"
	<< Exception::eventerror;
	}
	}
	psystem.rescaleMass();
	LorentzRotation R = Utilities::getBoostToCM(make_pair(psystem, deltap));
	Energy pz = SimplePhaseSpace::getMagnitude(sqr(ptotal.mass()),
	psystem.mass(), diquark->mass());
	LorentzRotation Rs(-(R*psystem).boostVector());
	Rs.boost(0.0, 0.0, pz/sqrt(sqr(pz) + sqr(psystem.mass())));
	Rs = Rs*R;
	// put remnant on shell
	deltap.transform(R);
	deltap.setMass(diquark->mass());
	deltap.setE(sqrt(sqr(diquark->mass())+sqr(pz)));
	deltap.rescaleRho();
	R.invert();
	deltap.transform(R);
	Rs = R*Rs;
	// apply transformation to required particles to absorb recoil
	for(unsigned int ix=0;ix<=iloc;++ix) {
	progenitors[ix]->deepTransform(Rs);
	}
	diquark->set5Momentum(deltap);
	}
	// two remnants
	else {
	for(unsigned int ix=0;ix<2;++ix) {
	if(!theprocessed[ix]) continue;
	pnew[ix]=Lorentz5Momentum();
	for(unsigned int iy=0;iy<theprocessed[ix]->children().size();++iy) {
	pnew[ix]+=theprocessed[ix]->children()[iy]->momentum();
	}
	mrem[ix]=sqrt(pnew[ix].m2());
	}
	// now find the remnant remnant cmf frame
	Lorentz5Momentum prem[2]={theprocessed[0]->momentum(),
	theprocessed[1]->momentum()};
	ptotal=prem[0]+prem[1];
	ptotal.rescaleMass();
	// boost momenta to this frame
	if(ptotal.m()< (pnew[0].m()+pnew[1].m()))
	throw Exception() << "Not enough energy in both remnants in "
	<< "HwRemDecayer::setRemMasses() "
	<< Exception::eventerror;

	Boost boostv(-ptotal.boostVector());
	ptotal.boost(boostv);
	for(unsigned int ix=0;ix<2;++ix) {
	prem[ix].boost(boostv);
	// set the masses and energies,
	prem[ix].setMass(mrem[ix]);
	prem[ix].setE(0.5/ptotal.m()*(sqr(ptotal.m())+sqr(mrem[ix])-sqr(mrem[1-ix])));
	prem[ix].rescaleRho();
	// boost back to the lab
	prem[ix].boost(-boostv);
	// set the momenta of the remnants
	theprocessed[ix]->set5Momentum(prem[ix]);
	}
	// boost the decay products
	Lorentz5Momentum ptemp;
	for(unsigned int ix=0;ix<2;++ix) {
	Boost btorest(-pnew[ix].boostVector());
	Boost bfmrest( prem[ix].boostVector());
	for(unsigned int iy=0;iy<theprocessed[ix]->children().size();++iy) {
	ptemp=theprocessed[ix]->children()[iy]->momentum();
	ptemp.boost(btorest);
	ptemp.boost(bfmrest);
	theprocessed[ix]->children()[iy]->set5Momentum(ptemp);
	}
	}
	}
	}

	void HwRemDecayer::initSoftInteractions(Energy ptmin, InvEnergy2 beta){
	ptmin_ = ptmin;
	beta_ = beta;
	}

	Energy HwRemDecayer::softPt() const {
	Energy2 pt2(ZERO);
	double xmin(0.0), xmax(1.0), x(0);

	if(beta_ == ZERO){
	return UseRandom::rnd(0.0,(double)(ptmin_/GeV))*GeV;
	}

	if(beta_ < ZERO){
	xmin = 1.0;
	xmax = exp( -beta_*sqr(ptmin_) );
	}else{
	xmin = exp( -beta_*sqr(ptmin_) );
	xmax = 1.0;
	}
	x = UseRandom::rnd(xmin, xmax);
	pt2 = 1.0/beta_ * log(1/x);

	if( pt2 < ZERO \|\| pt2 > sqr(ptmin_) )
	throw Exception() << "HwRemDecayer::softPt generation of pt "
	<< "outside allowed range [0," << ptmin_/GeV << "]."
	<< Exception::runerror;

	return sqrt(pt2);
	}

	void HwRemDecayer::softKinematics(Lorentz5Momentum &r1, Lorentz5Momentum &r2,
	Lorentz5Momentum &g1, Lorentz5Momentum &g2) const {

	g1 = Lorentz5Momentum();
	g2 = Lorentz5Momentum();
	//All necessary variables for the two soft gluons
	Energy pt(softPt()), pz(ZERO);
	Energy2 pz2(ZERO);
	double phi(UseRandom::rnd(2.*Constants::pi));
	double x_g1(0.0), x_g2(0.0);

	//Get the external momenta
	tcPPair beam(generator()->currentEventHandler()->currentCollision()->incoming());
	Lorentz5Momentum P1(beam.first->momentum()), P2(beam.second->momentum());

	if(dbg){
	cerr << "new event --------------------\n"
	<< (beam.first) << (softRems_.first)
	<< "-------------------\n"
	<< (beam.second) << (softRems_.second) << endl;
	}

	//parton mass
	Energy mp;
	if(quarkPair_){
	mp = getParticleData(ParticleID::u)->constituentMass();
	}else{
	mp = mg_;
	}

	//Get x_g1 and x_g2
	//first limits
	double xmin = sqr(ptmin_)/4.0/(P1+P2).m2();
	double x1max = (r1.e()+abs(r1.z()))/(P1.e() + abs(P1.z()));
	double x2max = (r2.e()+abs(r2.z()))/(P2.e() + abs(P2.z()));
	double x1;

	if(!multiPeriph_){
	//now generate according to 1/x
	x_g1 = xmin * exp(UseRandom::rnd(log(x1max/xmin)));
	x_g2 = xmin * exp(UseRandom::rnd(log(x2max/xmin)));
	}else{
	if(valOfN_==0) return;

	double param = (1/(2valOfN_+1))initTotRap_;
	do{
	// need 1-x instead of x to get the proper final momenta
	x1 = UseRandom::rndGauss(gaussWidth_, 1 - (exp(param)-1)/exp(param));
	}while(x1 < 0 \|\| x1>=1.0);
	x_g1 = x1max*x1;
	x_g2 = x2max*x1;
	}


	if(dbg)
	cerr << x1max << " " << x_g1 << endl << x2max << " " << x_g2 << endl;

	Lorentz5Momentum ig1, ig2, cmf;

	ig1 = x_g1*P1;
	ig2 = x_g2*P2;

	ig1.setMass(mp);
	ig2.setMass(mp);
	ig1.rescaleEnergy();
	ig2.rescaleEnergy();

	cmf = ig1 + ig2;
	//boost vector from cmf to lab
	Boost boostv(cmf.boostVector());

	//outgoing gluons in cmf
	g1.setMass(mp);
	g2.setMass(mp);

	g1.setX(pt*cos(phi));
	g2.setX(-pt*cos(phi));
	g1.setY(pt*sin(phi));
	g2.setY(-pt*sin(phi));

	pz2 = cmf.m2()/4 - sqr(mp) - sqr(pt);

	if(pz2/GeV2 < 0.0){
	if(dbg)
	cerr << "EXCEPTION not enough energy...." << endl;
	throw ExtraSoftScatterVeto();
	}

	if(!multiPeriph_){
	if(UseRandom::rndbool())
	pz = sqrt(pz2);
	else
	pz = -sqrt(pz2);
	}else{
	pz = pz2 > ZERO ? sqrt(pz2) : ZERO;
	}


	if(dbg)
	cerr << "pz has been calculated to: " << pz/GeV << endl;


	g1.setZ(pz);
	g2.setZ(-pz);
	g1.rescaleEnergy();
	g2.rescaleEnergy();

	if(dbg){
	cerr << "check inv mass in cmf frame: " << (g1+g2).m()/GeV
	<< " vs. lab frame: " << (ig1+ig2).m()/GeV << endl;
	}

	g1.boost(boostv);
	g2.boost(boostv);

	//recalc the remnant momenta
	Lorentz5Momentum r1old(r1), r2old(r2);

	r1 -= ig1;
	r2 -= ig2;

	try{
	reShuffle(r1, r2, r1old.m(), r2old.m());
	}catch(ExtraSoftScatterVeto){
	r1 = r1old;
	r2 = r2old;
	throw ExtraSoftScatterVeto();
	}

	if(dbg){
	cerr << "remnant 1,2 momenta: " << r1/GeV << "--" << r2/GeV << endl;
	cerr << "remnant 1,2 masses: " << r1.m()/GeV << " " << r2.m()/GeV << endl;
	cerr << "check momenta in the lab..." << (-r1old-r2old+r1+r2+g1+g2)/GeV << endl;
	}
	}

	void HwRemDecayer::doSoftInteractions_old(unsigned int N) {
	if(N == 0) return;
	if(!softRems_.first \|\| !softRems_.second)
	throw Exception() << "HwRemDecayer::doSoftInteractions: no "
	<< "Remnants available."
	<< Exception::runerror;

	if( ptmin_ == -1.*GeV )
	throw Exception() << "HwRemDecayer::doSoftInteractions: init "
	<< "code has not been called! call initSoftInteractions."
	<< Exception::runerror;

	Lorentz5Momentum g1, g2;
	Lorentz5Momentum r1(softRems_.first->momentum()), r2(softRems_.second->momentum());
	unsigned int tries(1), i(0);

	for(i=0; i<N; i++){
	//check how often this scattering has been regenerated
	if(tries > maxtrySoft_) break;

	if(dbg){
	cerr << "new try \n" << softRems_.first << softRems_.second << endl;
	}

	try{
	softKinematics(r1, r2, g1, g2);

	}catch(ExtraSoftScatterVeto){
	tries++;
	i--;
	continue;
	}

	PPair oldrems = softRems_;
	PPair gluons = make_pair(addParticle(softRems_.first, ParticleID::g, g1),
	addParticle(softRems_.second, ParticleID::g, g2));
	//now reset the remnants with the new ones
	softRems_.first = addParticle(softRems_.first, softRems_.first->id(), r1);
	softRems_.second = addParticle(softRems_.second, softRems_.second->id(), r2);

	//do the colour connections
	pair<bool, bool> anti = make_pair(oldrems.first->hasAntiColour(),
	oldrems.second->hasAntiColour());

	ColinePtr cl1 = new_ptr(ColourLine());
	ColinePtr cl2 = new_ptr(ColourLine());
	if( UseRandom::rnd() < colourDisrupt_ ){//this is the member variable, i.e. SOFT colour disruption
	//connect the remnants independent of the gluons
	oldrems.first->colourLine(anti.first)->addColoured(softRems_.first, anti.first);
	oldrems.second->colourLine(anti.second)->addColoured(softRems_.second, anti.second);
	//connect the gluons to each other
	cl1->addColoured(gluons.first);
	cl1->addAntiColoured(gluons.second);
	cl2->addColoured(gluons.second);
	cl2->addAntiColoured(gluons.first);
	}else{
	//connect the remnants to the gluons
	oldrems.first->colourLine(anti.first)->addColoured(gluons.first, anti.first);
	oldrems.second->colourLine(anti.second)->addColoured(gluons.second, anti.second);
	//and the remaining colour line to the final remnant
	cl1->addColoured(softRems_.first, anti.first);
	cl1->addColoured(gluons.first, !anti.first);
	cl2->addColoured(softRems_.second, anti.second);
	cl2->addColoured(gluons.second, !anti.second);
	}

	//reset counter
	tries = 1;
	}

	if(dbg)
	cerr << "generated " << i << "th soft scatters\n";

	}

	void HwRemDecayer::doSoftInteractions_multiPeriph(unsigned int N) {
	if(N == 0) return;
	int Nmpi = N;
	for(int j=0;j<Nmpi;j++){
	///////////////////////
	// Average number of partons in the ladder

	// Parameterize the ladder multiplicity to: ladderMult_ = A_0 * (s/1TeV^2)^alpha
	// with the two tunable parameters A_0 =ladderNorm_ and alpha = ladderPower_

	// Get the collision energy
	Energy energy(generator()->maximumCMEnergy());

	double reference = sqr(energy/TeV);
	double ladderMult_;

	// Parametrization of the ladder multiplicity
	ladderMult_ = ladderNorm_ * pow( ( reference ) , ladderPower_ );


	int avgN = floor(ladderMult_*log((softRems_.first->momentum()
	+softRems_.second->momentum()).m()/mg_) + ladderbFactor_);
	initTotRap_ = abs(softRems_.first->momentum().rapidity())
	+abs(softRems_.second->momentum().rapidity());

	//generate the poisson distribution with mean avgN
	double L = exp(-double(avgN));
	int k = 0;
	double p = 1;
	do {
	k++;
	p *= UseRandom::rnd();
	} while( p > L);
	N=k-1;

	valOfN_=N;
	if(N == 0) return;

	if(!softRems_.first \|\| !softRems_.second)
	throw Exception() << "HwRemDecayer::doSoftInteractions: no "
	<< "Remnants available."
	<< Exception::runerror;

	if( ptmin_ == -1.*GeV )
	throw Exception() << "HwRemDecayer::doSoftInteractions: init "
	<< "code has not been called! call initSoftInteractions."
	<< Exception::runerror;

	Lorentz5Momentum g1, g2, q1, q2;
	//intermediate gluons; erased in if there is
	//another step
	Lorentz5Momentum gint1, gint2;

	//momenta of the gluon pair generated in
	//each step
	vector< pair<Lorentz5Momentum,Lorentz5Momentum> > gluonMomPairs;
	//momenta of remnants
	Lorentz5Momentum r1(softRems_.first->momentum()), r2(softRems_.second->momentum());

	unsigned int tries(1), i(0);
	//generate the momenta of particles in the ladder
	for(i = 0; i < N; i++){
	//check how often this scattering has been regenerated
	//and break if exeeded maximal number
	if(tries > maxtrySoft_) break;

	if(dbg){
	cerr << "new try \n" << softRems_.first << softRems_.second << endl;
	}

	try{
	if(i==0){
	//generated partons in the ladder are quark-antiquark
	quarkPair_ = true;
	//first splitting: remnant -> remnant + quark
	softKinematics(r1, r2, q1, q2);
	}else if(i==1){

	//generated pair is gluon-gluon
	quarkPair_ = false;
	//second splitting; quark -> quark + gluon
	softKinematics(q1, q2, g1, g2);
	//record generated gluon pair
	gluonMomPairs.push_back(make_pair(g1,g2));

	}else{
	//consequent splittings gluon -> gluon+gluon
	//but, the previous gluon gets deleted
	quarkPair_ = false;
	//save first the previous gluon momentum
	g1=gluonMomPairs.back().first;
	g2=gluonMomPairs.back().second;
	//split gluon momentum
	softKinematics(g1, g2, gint1, gint2);

	//erase the last entry
	gluonMomPairs.pop_back();

	//add new gluons
	gluonMomPairs.push_back(make_pair(g1,g2));
	gluonMomPairs.push_back(make_pair(gint1,gint2));
	}

	}catch(ExtraSoftScatterVeto){
	tries++;
	i--;
	continue;
	}
	//reset counter
	tries = 1;
	}
	//if no gluons discard the ladder
	if(gluonMomPairs.size()==0) return;

	if(dbg)
	cerr << "generated " << i << "th soft scatters\n";

	//gluons from previous step
	PPair oldgluons;
	//quark-antiquark pair
	PPair quarks;
	//colour direction of quark-antiquark (true if anticolour)
	pair<bool, bool> anti_q;
	//generate particles (quark-antiquark and gluons) in the
	//ladder from momenta generated above
	for(i = 0; i <= gluonMomPairs.size(); i++){
	//remnants before splitting
	PPair oldrems = softRems_;
	//current gluon pair
	PPair gluons;

	//add quarks
	if(i==0){
	quarks = make_pair(addParticle(softRems_.first, ParticleID::u, q1),
	addParticle(softRems_.second, ParticleID::ubar, q2));
	anti_q = make_pair(quarks.first->hasAntiColour(),
	quarks.second->hasAntiColour());
	}



	if(i>0){

	gluons = make_pair(addParticle(softRems_.first, ParticleID::g,
	gluonMomPairs[i-1].first),
	addParticle(softRems_.second, ParticleID::g,
	gluonMomPairs[i-1].second));
	}

	//now reset the remnants with the new ones
	softRems_.first = addParticle(softRems_.first, softRems_.first->id(), r1);
	softRems_.second = addParticle(softRems_.second, softRems_.second->id(), r2);


	//colour direction of remnats
	pair<bool, bool> anti = make_pair(oldrems.first->hasAntiColour(),
	oldrems.second->hasAntiColour());



	ColinePtr cl1 = new_ptr(ColourLine());
	ColinePtr cl2 = new_ptr(ColourLine());
	//first colour connect remnants and if no
	//gluons the quark-antiquark pair
	if(i==0){
	oldrems.first->colourLine(anti.first)->addColoured(softRems_.first,
	anti.first);
	oldrems.second->colourLine(anti.second)->addColoured(softRems_.second,
	anti.second);
	if(gluonMomPairs.size()==0){
	ColinePtr clq = new_ptr(ColourLine());
	clq->addColoured(quarks.first, anti_q.first);
	clq->addColoured(quarks.second, anti_q.second);
	}
	}//colour connect remnants from previous step and gluons to quarks or gluons
	if(i!=0){
	oldrems.first->colourLine(anti.first)->addColoured(softRems_.first,
	anti.first);
	oldrems.second->colourLine(anti.second)->addColoured(softRems_.second,
	anti.second);
	if(i==1){
	cl1->addColoured(quarks.first, anti_q.first);
	cl1->addColoured(gluons.first, !anti_q.first);
	cl2->addColoured(quarks.second, anti_q.second);
	cl2->addColoured(gluons.second, !anti_q.second);
	}else{
	cl1->addColoured(oldgluons.first, anti_q.first);
	cl1->addColoured(gluons.first, !anti_q.first);
	cl2->addColoured(oldgluons.second, anti_q.second);
	cl2->addColoured(gluons.second, !anti_q.second);
	}

	} //last step; connect last gluons
	if(i > 0 && i == gluonMomPairs.size()){
	ColinePtr clg = new_ptr(ColourLine());
	clg->addColoured(gluons.first, anti_q.first);
	clg->addColoured(gluons.second, anti_q.second);
	}


	//save gluons for the next step
	if(i < gluonMomPairs.size()) oldgluons = gluons;
	}

	////////
	}
	////////
	}

	void HwRemDecayer::finalize(double colourDisrupt, unsigned int softInt){
	PPair diquarks;
	//Do the final Rem->Diquark or Rem->quark "decay"
	if(theRems.first) {
	diquarks.first = finalSplit(theRems.first, theContent.first.RemID(),
	theUsed.first);
	theMaps.first.push_back(make_pair(diquarks.first, tPPtr()));
	}
	if(theRems.second) {
	diquarks.second = finalSplit(theRems.second, theContent.second.RemID(),
	theUsed.second);
	theMaps.second.push_back(make_pair(diquarks.second, tPPtr()));
	}
	setRemMasses();
	if(theRems.first) {
	fixColours(theMaps.first, theanti.first, colourDisrupt);
	if(theContent.first.hadron->id()==ParticleID::pomeron&&
	pomeronStructure_==0) fixColours(theMaps.first, !theanti.first, colourDisrupt);
	}
	if(theRems.second) {
	fixColours(theMaps.second, theanti.second, colourDisrupt);
	if(theContent.second.hadron->id()==ParticleID::pomeron&&
	pomeronStructure_==0) fixColours(theMaps.second, !theanti.second, colourDisrupt);
	}

	if( !theRems.first \|\| !theRems.second ) return;
	//stop here if we don't have two remnants
	softRems_ = diquarks;
	doSoftInteractions(softInt);
	}


	HwRemDecayer::HadronContent
	HwRemDecayer::getHadronContent(tcPPtr hadron) const {
	HadronContent hc;
	hc.hadron = hadron->dataPtr();
	long id(hadron->id());
	// baryon
	if(BaryonMatcher::Check(hadron->data())) {
	hc.sign = id < 0? -1: 1;
	hc.flav.push_back((id = abs(id)/10)%10);
	hc.flav.push_back((id /= 10)%10);
	hc.flav.push_back((id /= 10)%10);
	hc.extracted = -1;
	}
	else if(hadron->data().id()==ParticleID::gamma \|\|
	(hadron->data().id()==ParticleID::pomeron && pomeronStructure_==1)) {
	hc.sign = 1;
	for(int ix=1;ix<6;++ix) {
	hc.flav.push_back( ix);
	hc.flav.push_back(-ix);
	}
	}
	else if(hadron->data().id()==ParticleID::pomeron ) {
	hc.sign = 1;
	hc.flav.push_back(ParticleID::g);
	hc.flav.push_back(ParticleID::g);
	}
	else if(hadron->data().id()==ParticleID::reggeon ) {
	hc.sign = 1;
	for(int ix=1;ix<3;++ix) {
	hc.flav.push_back( ix);
	hc.flav.push_back(-ix);
	}
	}
	hc.pomeronStructure = pomeronStructure_;
	return hc;
	}

	long HwRemDecayer::HadronContent::RemID() const{
	if(extracted == -1)
	throw Exception() << "Try to build a Diquark id without "
	<< "having extracted something in "
	<< "HwRemDecayer::RemID(...)"
	<< Exception::runerror;
	//the hadron was a meson or photon
	if(flav.size()==2) return sign*flav[(extracted+1)%2];

	long remId;
	int id1(sign*flav[(extracted+1)%3]),
	id2(sign*flav[(extracted+2)%3]),
	sign(0), spin(0);

	if (abs(id1) > abs(id2)) swap(id1, id2);
	sign = (id1 < 0) ? -1 : 1; // Needed for the spin 0/1 part
	remId = id21000+id1100;
	// Now decide if we have spin 0 diquark or spin 1 diquark
	if(id1 == id2) spin = 3; // spin 1
	else spin = 1; // otherwise spin 0
	remId += sign*spin;
	return remId;
	}


	tPPtr HwRemDecayer::addParticle(tcPPtr parent, long id, Lorentz5Momentum p) const {

	PPtr newp = new_ptr(Particle(getParticleData(id)));
	newp->set5Momentum(p);
	// Add the new remnant to the step, but don't do colour connections
	thestep->addDecayProduct(parent,newp,false);
	return newp;
	}

	void HwRemDecayer::findChildren(tPPtr part,vector<PPtr> & particles) const {
	if(part->children().empty()) particles.push_back(part);
	else {
	for(unsigned int ix=0;ix<part->children().size();++ix)
	findChildren(part->children()[ix],particles);
	}
	}

	ParticleVector HwRemDecayer::decay(const DecayMode &,
	const Particle &, Step &) const {
	throw Exception() << "HwRemDecayer::decay(...) "
	<< "must not be called explicitely."
	<< Exception::runerror;
	}

	void HwRemDecayer::persistentOutput(PersistentOStream & os) const {
	os << ounit(_kinCutoff, GeV) << _range << _zbin << _ybin
	<< _nbinmax << _alphaS << _alphaEM << DISRemnantOpt_
	<< maxtrySoft_ << colourDisrupt_ << ladderPower_<< ladderNorm_ << ladderbFactor_ << pomeronStructure_
	<< ounit(mg_,GeV) << ounit(ptmin_,GeV) << ounit(beta_,sqr(InvGeV))
	<< allowTop_ << multiPeriph_ << valOfN_ << initTotRap_;
	}

	void HwRemDecayer::persistentInput(PersistentIStream & is, int) {
	is >> iunit(_kinCutoff, GeV) >> _range >> _zbin >> _ybin
	>> _nbinmax >> _alphaS >> _alphaEM >> DISRemnantOpt_
	>> maxtrySoft_ >> colourDisrupt_ >> ladderPower_ >> ladderNorm_ >> ladderbFactor_ >> pomeronStructure_
	>> iunit(mg_,GeV) >> iunit(ptmin_,GeV) >> iunit(beta_,sqr(InvGeV))
	>> allowTop_ >> multiPeriph_ >> valOfN_ >> initTotRap_;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<HwRemDecayer,RemnantDecayer>
	describeHerwigHwRemDecayer("Herwig::HwRemDecayer", "HwShower.so");

	void HwRemDecayer::Init() {

	static ClassDocumentation<HwRemDecayer> documentation
	("The HwRemDecayer class decays the remnant for Herwig");

	static Parameter<HwRemDecayer,double> interfaceZBinSize
	("ZBinSize",
	"The size of the vbins in z for the interpolation of the splitting function.",
	&HwRemDecayer::_zbin, 0.05, 0.001, 0.1,
	false, false, Interface::limited);

	static Parameter<HwRemDecayer,int> interfaceMaxBin
	("MaxBin",
	"Maximum number of z bins",
	&HwRemDecayer::_nbinmax, 100, 10, 1000,
	false, false, Interface::limited);

	static Reference<HwRemDecayer,ShowerAlpha> interfaceAlphaS
	("AlphaS",
	"Pointer to object to calculate the strong coupling",
	&HwRemDecayer::_alphaS, false, false, true, false, false);

	static Reference<HwRemDecayer,ShowerAlpha> interfaceAlphaEM
	("AlphaEM",
	"Pointer to object to calculate the electromagnetic coupling",
	&HwRemDecayer::_alphaEM, false, false, true, false, false);

	static Parameter<HwRemDecayer,Energy> interfaceKinCutoff
	("KinCutoff",
	"Parameter kinCutoff used to constrain qtilde",
	&HwRemDecayer::_kinCutoff, GeV, 0.75GeV, 0.5GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<HwRemDecayer,double> interfaceEmissionRange
	("EmissionRange",
	"Factor above the minimum possible value in which the forced splitting is allowed.",
	&HwRemDecayer::_range, 1.1, 1.0, 10.0,
	false, false, Interface::limited);


	static Switch<HwRemDecayer,unsigned int> interfaceDISRemnantOption
	("DISRemnantOption",
	"Options for the treatment of the remnant in DIS",
	&HwRemDecayer::DISRemnantOpt_, 0, false, false);
	static SwitchOption interfaceDISRemnantOptionDefault
	(interfaceDISRemnantOption,
	"Default",
	"Use the minimum number of particles needed to take the recoil"
	" and allow the lepton to be used if needed",
	0);
	static SwitchOption interfaceDISRemnantOptionNoLepton
	(interfaceDISRemnantOption,
	"NoLepton",
	"Use the minimum number of particles needed to take the recoil but"
	" veto events where the lepton kinematics would need to be altered",
	1);
	static SwitchOption interfaceDISRemnantOptionAllParticles
	(interfaceDISRemnantOption,
	"AllParticles",
	"Use all particles in the colour connected system to take the recoil"
	" and use the lepton if needed.",
	2);
	static SwitchOption interfaceDISRemnantOptionAllParticlesNoLepton
	(interfaceDISRemnantOption,
	"AllParticlesNoLepton",
	"Use all the particles in the colour connected system to take the"
	" recoil but don't use the lepton.",
	3);

	static Parameter<HwRemDecayer,unsigned int> interfaceMaxTrySoft
	("MaxTrySoft",
	"The maximum number of regeneration attempts for an additional soft scattering",
	&HwRemDecayer::maxtrySoft_, 10, 0, 100,
	false, false, Interface::limited);

	static Parameter<HwRemDecayer,double> interfacecolourDisrupt
	("colourDisrupt",
	"Fraction of connections to additional soft subprocesses, which are colour disrupted.",
	&HwRemDecayer::colourDisrupt_,
	1.0, 0.0, 1.0,
	false, false, Interface::limited);

	static Parameter<HwRemDecayer,double> interaceladderPower
	("ladderPower",
	"The power factor in the ladder parameterization.",
	&HwRemDecayer::ladderPower_,
	1.0, -5.0, 10.0,
	false, false, Interface::limited);

	static Parameter<HwRemDecayer,double> interfaceladderNorm
	("ladderNorm",
	"The normalization factor in the ladder parameterization",
	&HwRemDecayer::ladderNorm_,
	1.0, 0.0, 10.0,
	false, false, Interface::limited);


	static Parameter<HwRemDecayer,double> interfaceladderbFactor
	("ladderbFactor",
	"The additive factor in the multiperipheral ladder multiplicity.",
	&HwRemDecayer::ladderbFactor_,
	1.0, 0.0, 10.0,
	false, false, Interface::limited);

	static Parameter<HwRemDecayer,double> interfacegaussWidth
	("gaussWidth",
	"The gaussian width of the fluctuation of longitudinal momentum fraction.",
	&HwRemDecayer::gaussWidth_,
	0.1, 0.0, 1.0,
	false, false, Interface::limited);


	static Switch<HwRemDecayer,unsigned int> interfacePomeronStructure
	("PomeronStructure",
	"Option for the treatment of the valance structure of the pomeron",
	&HwRemDecayer::pomeronStructure_, 0, false, false);
	static SwitchOption interfacePomeronStructureGluon
	(interfacePomeronStructure,
	"Gluon",
	"Assume the pomeron is a two gluon state",
	0);
	static SwitchOption interfacePomeronStructureQQBar
	(interfacePomeronStructure,
	"QQBar",
	"Assumne the pomeron is q qbar as for the photon,"
	" this option is not recommended and is provide for compatiblity with POMWIG",
	1);

	static Switch<HwRemDecayer,bool> interfaceAllowTop
	("AllowTop",
	"Allow top quarks in the hadron",
	&HwRemDecayer::allowTop_, false, false, false);
	static SwitchOption interfaceAllowTopNo
	(interfaceAllowTop,
	"No",
	"Don't allow them",
	false);
	static SwitchOption interfaceAllowTopYes
	(interfaceAllowTop,
	"Yes",
	"Allow them",
	true);

	static Switch<HwRemDecayer,bool> interfaceMultiPeriph
	("MultiPeriph",
	"Use multiperipheral kinematics",
	&HwRemDecayer::multiPeriph_, false, false, false);
	static SwitchOption interfaceMultiPeriphNo
	(interfaceMultiPeriph,
	"No",
	"Don't use multiperipheral",
	false);
	static SwitchOption interfaceMultiPeriphYes
	(interfaceMultiPeriph,
	"Yes",
	"Use multiperipheral kinematics",
	true);

	}

	bool HwRemDecayer::canHandle(tcPDPtr particle, tcPDPtr parton) const {
	if(! (StandardQCDPartonMatcher::Check(*parton) \|\| parton->id()==ParticleID::gamma) ) {
	if(abs(parton->id())==ParticleID::t) {
	if(!allowTop_)
	throw Exception() << "Top is not allow as a parton in the remant handling, please "
	<< "use a PDF which does not contain top for the remnant"
	<< " handling (preferred) or allow top in the remnant using\n"
	<< " set " << fullName() << ":AllowTop Yes\n"
	<< Exception::runerror;
	}
	else
	return false;
	}
	return HadronMatcher::Check(*particle) \|\| particle->id()==ParticleID::gamma
	\|\| particle->id()==ParticleID::pomeron \|\| particle->id()==ParticleID::reggeon;
	}

	bool HwRemDecayer::isPartonic(tPPtr parton) const {
	if(parton->parents().empty()) return false;
	tPPtr parent = parton->parents()[0];
	bool partonic = false;
	for(unsigned int ix=0;ix<parent->children().size();++ix) {
	if(dynamic_ptr_cast<tRemPPtr>(parent->children()[ix])) {
	partonic = true;
	break;
	}
	}
	return partonic;
	}
	diff --git a/PDF/Makefile.am b/PDF/Makefile.am
	--- a/PDF/Makefile.am
	+++ b/PDF/Makefile.am
	@@ -1,58 +1,64 @@
	EXTRA_DIST = diffraction

	pkglib_LTLIBRARIES = HwPomeronPDF.la
	HwPomeronPDF_la_SOURCES = \
	PomeronPDF.cc PomeronPDF.h
	## add this to produce tests of the PDFs
	## HwDIFFRACTIVEPDF_la_CPPFLAGS=$(AM_CPPFLAGS) -DDIFFRACTIVEPDF_TESTING
	HwPomeronPDF_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 4:0:0

	pkglib_LTLIBRARIES += HwReggeonPDF.la
	HwReggeonPDF_la_SOURCES = \
	ReggeonPDF.cc ReggeonPDF.h
	## add this to produce tests of the PDFs
	## HwDIFFRACTIVEPDF_la_CPPFLAGS=$(AM_CPPFLAGS) -DDIFFRACTIVEPDF_TESTING
	HwReggeonPDF_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 4:0:0

	pkglib_LTLIBRARIES += HwPomeronFlux.la
	HwPomeronFlux_la_SOURCES = \
	PomeronFlux.h PomeronFlux.cc
	HwPomeronFlux_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 5:0:0

	## bound into HwShower.la
	noinst_LTLIBRARIES = libHwRemDecayer.la
	libHwRemDecayer_la_SOURCES = \
	HwRemDecayer.h HwRemDecayer.cc HwRemDecayer.fh

	## bound into HwShower.la
	noinst_LTLIBRARIES += libHwMPIPDF.la
	libHwMPIPDF_la_SOURCES = \
	MPIPDF.h MPIPDF.cc MPIPDF.fh \
	MinBiasPDF.h MinBiasPDF.cc MinBiasPDF.fh

	pkglib_LTLIBRARIES += HwSatPDF.la
	HwSatPDF_la_SOURCES = \
	SatPDF.h SatPDF.cc

	HwSatPDF_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 7:0:0

	+pkglib_LTLIBRARIES += HwPartonExtractor.la
	+HwPartonExtractor_la_SOURCES = \
	+MultiPartonExtractor.h MultiPartonExtractor.cc
	+
	+HwPartonExtractor_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 7:0:0
	+
	pkglib_LTLIBRARIES += HwIncomingPhotonEvolver.la
	HwIncomingPhotonEvolver_la_SOURCES = \
	IncomingPhotonEvolver.h IncomingPhotonEvolver.cc

	HwIncomingPhotonEvolver_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 7:0:0

	pkglib_LTLIBRARIES += HwSaSPhotonPDF.la
	HwSaSPhotonPDF_la_SOURCES = \
	SaSGamma.f SaSPhotonPDF.cc SaSPhotonPDF.h

	HwSaSPhotonPDF_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 2:0:0

	install-data-local:
	for i in `find $(srcdir)/diffraction -name '*.data'`; \
	do \
	$(install_sh_DATA) $$i $(DESTDIR)$(pkgdatadir)/PDF/$${i#$(srcdir)/}; \
	done

	uninstall-local:
	rm -rf $(DESTDIR)$(pkgdatadir)/PDF
	diff --git a/PDF/MultiPartonExtractor.cc b/PDF/MultiPartonExtractor.cc
	new file mode 100644
	--- /dev/null
	+++ b/PDF/MultiPartonExtractor.cc
	@@ -0,0 +1,109 @@
	+// -- C++ --
	+//
	+// This is the implementation of the non-inlined, non-templated member
	+// functions of the MultiPartonExtractor class.
	+//
	+
	+#include "MultiPartonExtractor.h"
	+#include "ThePEG/Interface/ClassDocumentation.h"
	+#include "ThePEG/Interface/RefVector.h"
	+#include "ThePEG/EventRecord/Particle.h"
	+#include "ThePEG/Repository/UseRandom.h"
	+#include "ThePEG/Repository/EventGenerator.h"
	+#include "ThePEG/Utilities/DescribeClass.h"
	+#include "ThePEG/Persistency/PersistentOStream.h"
	+#include "ThePEG/Persistency/PersistentIStream.h"
	+#include "ThePEG/PDF/RemnantHandler.h"
	+#include "ThePEG/PDF/NoPDF.h"
	+
	+using namespace Herwig;
	+
	+IBPtr MultiPartonExtractor::clone() const {
	+ return new_ptr(*this);
	+}
	+
	+IBPtr MultiPartonExtractor::fullclone() const {
	+ return new_ptr(*this);
	+}
	+
	+void MultiPartonExtractor::persistentOutput(PersistentOStream & os) const {
	+ os << firstPDF_ << secondPDF_;
	+}
	+
	+void MultiPartonExtractor::persistentInput(PersistentIStream & is, int) {
	+ is >> firstPDF_ >> secondPDF_;
	+}
	+
	+// The following static variable is needed for the type
	+// description system in ThePEG.
	+DescribeClass<MultiPartonExtractor,PartonExtractor>
	+describeHerwigMultiPartonExtractor("Herwig::MultiPartonExtractor", "HwPartonExtractor.so");
	+
	+void MultiPartonExtractor::Init() {
	+
	+ static ClassDocumentation<MultiPartonExtractor> documentation
	+ ("The MultiPartonExtractor class allows more control over the PDFs used"
	+ " than the default PartonExtractor of ThePEG");
	+
	+ static RefVector<MultiPartonExtractor,PDFBase> interfaceFirstPDFs
	+ ("FirstPDFs",
	+ "The PDFs for the first beam",
	+ &MultiPartonExtractor::firstPDF_, -1, false, false, true, false, false);
	+
	+ static RefVector<MultiPartonExtractor,PDFBase> interfaceSecondPDFs
	+ ("SecondPDFs",
	+ "The PDFs for the second beam",
	+ &MultiPartonExtractor::secondPDF_, -1, false, false, true, false, false);
	+
	+}
	+
	+PartonPairVec MultiPartonExtractor::
	+getPartons(Energy maxEnergy, const cPDPair & incoming,
	+ const Cuts & kc) const {
	+ PartonPairVec result;
	+ PartonVector first;
	+ PDFCuts cuts1(kc, true, maxEnergy);
	+ PBPtr p1 =
	+ new_ptr(PartonBin(PDPtr(), PBPtr(), incoming.first, PDFPtr(), cuts1));
	+ std::deque<tcPDFPtr> pdfs(firstPDF_.begin(),firstPDF_.end());
	+ addPartons(p1, cuts1, pdfs, first);
	+ PartonVector second;
	+ PDFCuts cuts2(kc, false, maxEnergy);
	+ pdfs = std::deque<tcPDFPtr>(secondPDF_.begin(),secondPDF_.end());
	+ PBPtr p2 =
	+ new_ptr(PartonBin(PDPtr(), PBPtr(), incoming.second, PDFPtr(), cuts2));
	+ addPartons(p2, cuts2, pdfs, second);
	+ for ( PartonVector::iterator it1 = first.begin();
	+ it1 != first.end(); ++it1 )
	+ for ( PartonVector::iterator it2 = second.begin();
	+ it2 != second.end(); ++it2 )
	+ result.push_back(PBPair(it1, it2));
	+
	+ // We add the original parton bins as well to avoid them being
	+ // deleted.
	+ result.push_back(PBPair(p1, p2));
	+ return result;
	+}
	+
	+void MultiPartonExtractor::
	+addPartons(tPBPtr incoming, const PDFCuts & cuts, std::deque<tcPDFPtr> pdfs,
	+ PartonVector & pbins) const {
	+ tcPDFPtr pdf;
	+ if(!pdfs.empty()) {
	+ pdf = *pdfs.begin();
	+ pdfs.pop_front();
	+ }
	+ if(!pdf) pdf = getPDF(incoming->parton());
	+ if ( dynamic_ptr_cast<Ptr<NoPDF>::tcp>(pdf) \|\|
	+ incoming->parton() == incoming->particle() ) {
	+ pbins.push_back(incoming);
	+ return;
	+ }
	+ cPDVector partons = pdf->partons(incoming->parton());
	+ for ( int i = 0, N = partons.size(); i < N; ++i ) {
	+ PBPtr pb =
	+ new_ptr(PartonBin(incoming->parton(), incoming, partons[i], pdf, cuts));
	+ incoming->addOutgoing(pb);
	+ addPartons(pb, cuts, pdfs, pbins);
	+ }
	+}
	diff --git a/PDF/MultiPartonExtractor.h b/PDF/MultiPartonExtractor.h
	new file mode 100644
	--- /dev/null
	+++ b/PDF/MultiPartonExtractor.h
	@@ -0,0 +1,114 @@
	+// -- C++ --
	+#ifndef Herwig_MultiPartonExtractor_H
	+#define Herwig_MultiPartonExtractor_H
	+//
	+// This is the declaration of the MultiPartonExtractor class.
	+//
	+
	+#include "ThePEG/PDF/PartonExtractor.h"
	+#include <deque>
	+
	+namespace Herwig {
	+
	+using namespace ThePEG;
	+
	+/**
	+ * The MultiPartonExtractor class inherits from the PartonExtractor of ThePEG
	+ * but allows more control over the PDFs used in the case that there are multiple
	+ * stages of parton extraction
	+ *
	+ * @see \ref MultiPartonExtractorInterfaces "The interfaces"
	+ * defined for MultiPartonExtractor.
	+ */
	+class MultiPartonExtractor: public PartonExtractor {
	+
	+public:
	+
	+ /**
	+ * The default constructor.
	+ */
	+ MultiPartonExtractor() {};
	+
	+ /**
	+ * Return a vector of possible pairs of parton bins which can be
	+ * produced within a given maximum total particle-particle
	+ * invariant mass squared, \a maxEnergy sBin.
	+ */
	+ virtual PartonPairVec getPartons(Energy maxEnergy, const cPDPair &,
	+ const Cuts &) const;
	+
	+public:
	+
	+ /** @name Functions used by the persistent I/O system. */
	+ //@{
	+ /**
	+ * Function used to write out object persistently.
	+ * @param os the persistent output stream written to.
	+ */
	+ void persistentOutput(PersistentOStream & os) const;
	+
	+ /**
	+ * Function used to read in object persistently.
	+ * @param is the persistent input stream read from.
	+ * @param version the version number of the object when written.
	+ */
	+ void persistentInput(PersistentIStream & is, int version);
	+ //@}
	+
	+ /**
	+ * The standard Init function used to initialize the interfaces.
	+ * Called exactly once for each class by the class description system
	+ * before the main function starts or
	+ * when this class is dynamically loaded.
	+ */
	+ static void Init();
	+
	+protected:
	+
	+ /**
	+ * Add parton bins to pbins for the given incoming particle and the
	+ * specified cuts.
	+ */
	+ virtual void addPartons(tPBPtr incoming ,const PDFCuts & cuts,
	+ std::deque<tcPDFPtr> pdf ,PartonVector & pbins) const;
	+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:
	+
	+ /**
	+ * The assignment operator is private and must never be called.
	+ * In fact, it should not even be implemented.
	+ */
	+ MultiPartonExtractor & operator=(const MultiPartonExtractor &);
	+
	+
	+ /**
	+ * PDFBase object to override first PDF
	+ */
	+ vector<PDFPtr> firstPDF_;
	+
	+ /**
	+ * PDFBase object to override second PDF
	+ */
	+ vector<PDFPtr> secondPDF_;
	+
	+};
	+
	+}
	+
	+#endif /* Herwig_MultiPartonExtractor_H */
	diff --git a/Tests/Makefile.am b/Tests/Makefile.am
	--- a/Tests/Makefile.am
	+++ b/Tests/Makefile.am
	@@ -1,379 +1,397 @@
	AM_LDFLAGS += -module -avoid-version -rpath /dummy/path/not/used

	EXTRA_DIST = Inputs python Rivet

	EXTRA_LTLIBRARIES = LeptonTest.la GammaTest.la HadronTest.la DISTest.la

	if WANT_LIBFASTJET
	EXTRA_LTLIBRARIES += HadronJetTest.la LeptonJetTest.la
	HadronJetTest_la_SOURCES = \
	Hadron/VHTest.h Hadron/VHTest.cc\
	Hadron/VTest.h Hadron/VTest.cc\
	Hadron/HTest.h Hadron/HTest.cc
	HadronJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
	-I$(FASTJETPATH)
	HadronJetTest_la_LIBADD = $(FASTJETLIBS)
	LeptonJetTest_la_SOURCES = \
	Lepton/TopDecay.h Lepton/TopDecay.cc
	LeptonJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
	-I$(FASTJETPATH)
	LeptonJetTest_la_LIBADD = $(FASTJETLIBS)
	endif

	LeptonTest_la_SOURCES = \
	Lepton/VVTest.h Lepton/VVTest.cc \
	Lepton/VBFTest.h Lepton/VBFTest.cc \
	Lepton/VHTest.h Lepton/VHTest.cc \
	Lepton/FermionTest.h Lepton/FermionTest.cc

	GammaTest_la_SOURCES = \
	Gamma/GammaMETest.h Gamma/GammaMETest.cc \
	Gamma/GammaPMETest.h Gamma/GammaPMETest.cc

	DISTest_la_SOURCES = \
	DIS/DISTest.h DIS/DISTest.cc

	HadronTest_la_SOURCES = \
	Hadron/HadronVVTest.h Hadron/HadronVVTest.cc\
	Hadron/HadronVBFTest.h Hadron/HadronVBFTest.cc\
	Hadron/WHTest.h Hadron/WHTest.cc\
	Hadron/ZHTest.h Hadron/ZHTest.cc\
	Hadron/VGammaTest.h Hadron/VGammaTest.cc\
	Hadron/ZJetTest.h Hadron/ZJetTest.cc\
	Hadron/WJetTest.h Hadron/WJetTest.cc\
	Hadron/QQHTest.h Hadron/QQHTest.cc


	REPO = $(top_builddir)/src/HerwigDefaults.rpo
	HERWIG = $(top_builddir)/src/Herwig
	HWREAD = $(HERWIG) read --repo $(REPO) -L $(builddir)/.libs -i $(top_builddir)/src
	HWBUILD = $(HERWIG) build --repo $(REPO) -L $(builddir)/.libs -i $(top_builddir)/src
	HWINTEGRATE = $(HERWIG) integrate
	HWRUN = $(HERWIG) run -N $${NUMEVENTS:-10000}


	tests : tests-LEP tests-DIS tests-LHC tests-Gamma


	LEPDEPS = \
	test-LEP-VV \
	test-LEP-VH \
	test-LEP-VBF \
	test-LEP-BB \
	test-LEP-Quarks \
	test-LEP-Leptons

	if WANT_LIBFASTJET
	LEPDEPS += test-LEP-TopDecay
	endif

	tests-LEP : $(LEPDEPS)

	tests-DIS : test-DIS-Charged test-DIS-Neutral


	LHCDEPS = \
	test-LHC-WW test-LHC-WZ test-LHC-ZZ \
	test-LHC-ZGamma test-LHC-WGamma \
	test-LHC-ZH test-LHC-WH \
	test-LHC-ZJet test-LHC-WJet \
	test-LHC-Z test-LHC-W \
	test-LHC-ZZVBF test-LHC-VBF \
	test-LHC-WWVBF \
	test-LHC-bbH test-LHC-ttH \
	test-LHC-GammaGamma test-LHC-GammaJet \
	test-LHC-Higgs test-LHC-HiggsJet \
	test-LHC-QCDFast test-LHC-QCD \
	test-LHC-Top


	if WANT_LIBFASTJET
	LHCDEPS += \
	test-LHC-Bottom \
	test-LHC-WHJet test-LHC-ZHJet test-LHC-HJet \
	test-LHC-ZShower test-LHC-WShower \
	test-LHC-WHJet-Powheg test-LHC-ZHJet-Powheg test-LHC-HJet-Powheg \
	test-LHC-ZShower-Powheg test-LHC-WShower-Powheg
	endif

	tests-LHC : $(LHCDEPS)

	tests-Gamma : test-Gamma-FF test-Gamma-WW test-Gamma-P



	LEPLIBS = LeptonTest.la
	HADLIBS = HadronTest.la

	if WANT_LIBFASTJET
	LEPLIBS += LeptonJetTest.la
	HADLIBS += HadronJetTest.la
	endif


	test-LEP-% : Inputs/LEP-%.in $(LEPLIBS)
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<))

	test-Gamma-% : Inputs/Gamma-%.in GammaTest.la
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<))

	test-DIS-% : Inputs/DIS-%.in DISTest.la
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<))

	test-LHC-% : Inputs/LHC-%.in GammaTest.la $(HADLIBS)
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<))



	tests-Rivet : Rivet-LEP Rivet-BFactory Rivet-DIS Rivet-Star Rivet-SppS \
	Rivet-TVT-WZ Rivet-TVT-Photon Rivet-TVT-Jets \
	Rivet-LHC-Jets Rivet-LHC-EW Rivet-LHC-Photon Rivet-LHC-Higgs

	Rivet-%.run : Rivet/%.in
	$(HWBUILD) -c .cache/$(subst .run,,$@) $<

	Rivet-Matchbox-%.yoda : Rivet-Matchbox-%.run
	$(HWINTEGRATE) -c .cache/$(subst .run,,$<) $<
	$(HWRUN) -c .cache/$(subst .run,,$<) $<

	Rivet-%.yoda : Rivet-%.run
	$(HWRUN) $<

	Rivet/%.in :
	python/make_input_files.py $(notdir $(subst .in,,$@))



	Rivet-inputfiles: $(shell echo Rivet/LEP{,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg,-Merging}-{9.4,12,13,17,27.6,29,30.2,30.7,30,31.3,34.8,41,43.6,45,50,52,55,56,57,60.8,60,61.4,66,76,82,8510,12.8,22,26.8,35,44,48.0,91,93.0,130,133,136,161,172,177,183,189,192,196,197,200,202,206,91-nopi}.in) \
	$(shell echo Rivet/LEP{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Powheg,-Matchbox-Powheg}-14.in) \
	$(shell echo Rivet/LEP{,-Dipole}-{10.5,11.96,12.8,13.96,16.86,21.84,26.8,28.48,35.44,48.0,97.0}-gg.in) \
	$(shell echo Rivet/BFactory{,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg}-{10.52,10.52-sym,10.54,10.45,10.47}.in) \
	$(shell echo Rivet/BFactory{,-Dipole}-{Upsilon,Upsilon2,Upsilon4,Tau,Phi,10.58-res}.in) \
	$(shell echo Rivet/DIS{,-NoME,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg,-Merging}-{e--LowQ2,e+-LowQ2,e+-HighQ2}.in) \
	$(shell echo Rivet/TVT{,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg,-Merging}-{Run-I-Z,Run-I-W,Run-I-WZ,Run-II-Z-e,Run-II-Z-{,LowMass-,HighMass-}mu,Run-II-W}.in) \
	$(shell echo Rivet/TVT{,-Dipole}-Run-II-{DiPhoton-GammaGamma,DiPhoton-GammaJet,PromptPhoton}.in) \
	$(shell echo Rivet/TVT-Powheg-Run-II-{DiPhoton-GammaGamma,DiPhoton-GammaJet}.in) \
	$(shell echo Rivet/TVT{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{Run-II-Jets-{0..11},Run-I-Jets-{1..8}}.in ) \
	$(shell echo Rivet/TVT{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{630-Jets-{1..3},300-Jets-1,900-Jets-1}.in ) \
	$(shell echo Rivet/TVT{,-Dipole}-{Run-I,Run-II,300,630,900}-UE.in) \
	$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-DiJets-{1..7}-{A,B,C}.in ) \
	$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-13-DiJets-{6..11}-A.in ) \
	$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{7,8,13}-Jets-{0..10}.in ) \
	$(shell echo Rivet/LHC{,-Dipole}-{900,2360,2760,7,8,13}-UE.in ) \
	$(shell echo Rivet/LHC{,-Dipole}-{900,7,13}-UE-Long.in ) \
	$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-Charm-{1..5}.in) \
	$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-Bottom-{0..9}.in) \
	$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-Top-{L,SL}.in) \
	$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{8,13}-Top-{All,L,SL}.in) \
	$(shell echo Rivet/Star{,-Dipole}-{UE,Jets-{1..4}}.in ) \
	$(shell echo Rivet/SppS{,-Dipole}-{200,500,900,546}-UE.in ) \
	$(shell echo Rivet/LHC{,-Matchbox,-Matchbox-Powheg,-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{W-{e,mu},13-Z-{e,mu},Z-HighMass{1,2}-e,{8,13}-W-mu,8-Z-Mass{1..4}-{e,mu},Z-{e,mu,mu-SOPHTY},Z-LowMass-{e,mu},Z-MedMass-e,WZ,WW-{emu,ll},ZZ-{ll,lv},{8,13}-WZ,8-ZZ-lv,8-WW-ll,Z-mu-Short}.in) \
	$(shell echo Rivet/LHC{,-Dipole}-7-{W,Z}Gamma-{e,mu}.in) \
	$(shell echo Rivet/LHC{,-Dipole}-8-ZGamma-{e,mu}.in) \
	$(shell echo Rivet/LHC{,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{7-W-Jet-{1..3}-e,7-Z-Jet-{0..3}-e,7-Z-Jet-0-mu}.in) \
	$(shell echo Rivet/LHC{-Matchbox,-Matchbox-Powheg,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{Z-b,Z-bb,8-Z-b,8-Z-bb,W-b,8-Z-jj}.in) \
	$(shell echo Rivet/LHC{,-Dipole}-{7,8,13}-PromptPhoton-{1..4}.in) Rivet/LHC-GammaGamma-7.in \
	$(shell echo Rivet/LHC{,-Powheg}-{7,8}-{DiPhoton-GammaGamma,DiPhoton-GammaJet}.in) \
	$(shell echo Rivet/LHC{,-Powheg,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{ggH,VBF,WH,ZH}.in) \
	$(shell echo Rivet/LHC{,-Powheg,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-8-{{ggH,VBF,WH,ZH}{,-GammaGamma},ggH-WW}.in) \
	$(shell echo Rivet/LHC{,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-ggHJet.in)
	# $(shell echo Rivet/ISR-{30,44,53,62}-UE.in ) $(shell echo Rivet/SppS-{53,63}-UE.in )

	+Rivet-GammaGamma: Rivet-GammaGamma/done
	+ touch $@

	+Rivet-GammaGamma/done: $(shell echo Rivet-GammaGamma-mumu-{3.5,4.5,5.5,6.5,7.5,9.0,12.5,17.5,30.0}.yoda )
	+ rm -rf Rivet-GammaGamma
	+ python/merge-GammaGamma GammaGamma
	+ rivet-mkhtml -o Rivet-GammaGamma GammaGamma.yoda:Hw
	+ touch $@
	+
	+Rivet-LEP-Gamma: Rivet-LEP-Gamma/done
	+ touch $@
	+
	+Rivet-LEP-Gamma/done: $(shell echo Rivet-LEP-Gamma-Direct-mumu-{161,172,183,189,196,206}.yoda ) \
	+ $(shell echo Rivet-LEP-Gamma-Direct-tautau-{189,196,206}.yoda ) \
	+ $(shell echo Rivet-LEP-Gamma-{Direct,Single-Resolved,Double-Resolved}-Jets-{198,206}.yoda )
	+ rm -rf Rivet-LEP-Gamma
	+ python/merge-LEP-Gamma LEP-Gamma
	+ rivet-mkhtml -o Rivet-LEP-Gamma LEP-Gamma.yoda:Hw
	+ touch $@

	Rivet-LEP : Rivet-LEP/done
	touch $@

	Rivet-LEP/done : $(shell echo Rivet{,-Powheg}-LEP-{9.4,12,13,14,17,27.6,29,30.2,30.7,30,31.3,34.8,43.6,45,50,52,55,56,57,60.8,60,61.4,66,76,10,12.8,22,26.8,35,41,44,48.0,82,85,91,93.0,130,133,136,161,172,177,183,189,192,196,197,200,202,206,91-nopi}.yoda) \
	$(shell echo Rivet-LEP-{10.5,11.96,12.8,13.96,16.86,21.84,26.8,28.48,35.44,48.0,97.0}-gg.yoda)
	rm -rf Rivet-LEP
	python/merge-LEP --with-gg LEP
	python/merge-LEP Powheg-LEP
	rivet-mkhtml -o Rivet-LEP LEP.yoda:Hw Powheg-LEP.yoda:Hw-Powheg
	touch $@

	Rivet-LowEnergy-%.yoda:
	$(HWBUILD) -c .cache/$(subst .yoda,,$@) Rivet/$(subst .yoda,.in,$@)
	$(HWRUN) $(subst .yoda,.run,$@)

	Rivet-LowEnergy-%:
	OUTPUT=`python/LowEnergy.py --process $(subst Rivet-LowEnergy-,,$@) --non-perturbative --perturbative`; make -j12 $$OUTPUT NUMEVENTS=$${NUMEVENTS:-10000};
	python/mergeLowEnergy.py $(subst Rivet-LowEnergy-,,$@)
	rivet-mkhtml -o Rivet-LowEnergy-$(subst Rivet-LowEnergy-,,$@) LowEnergy-EE-NonPerturbative-$(subst Rivet-LowEnergy-,,$@).yoda:"Non-Pert" LowEnergy-EE-Perturbative-$(subst Rivet-LowEnergy-,,$@).yoda:"Pert"

	Rivet-R:
	OUTPUT=`python/R.py --perturbative`; make -j12 $$OUTPUT NUMEVENTS=$${NUMEVENTS:-10000};
	python/mergeLowEnergy.py R
	rivet-mkhtml -o Rivet-R LowEnergy-EE-Perturbative-R.yoda:"Pert"
	#LowEnergy-EE-NonPerturbative-R.yoda:"Non-Pert"

	Rivet-BFactory : Rivet-BFactory/done
	touch $@

	Rivet-BFactory/done: $(shell echo Rivet-BFactory-{10.52,10.52-sym,10.54,10.45,10.47,Upsilon,Upsilon2,Upsilon4,Tau,Phi,10.58-res,10.58}.yoda)
	rm -rf Rivet-BFactory
	python/merge-BFactory BFactory
	rivet-mkhtml -o Rivet-BFactory BFactory.yoda:Hw
	touch $@



	Rivet-DIS : Rivet-DIS/done
	touch $@

	Rivet-DIS/done: $(shell echo Rivet{-DIS,-DIS-NoME,-Powheg-DIS}-{e--LowQ2,e+-LowQ2,e+-HighQ2}.yoda)
	rm -rf Rivet-DIS
	python/merge-DIS DIS
	python/merge-DIS Powheg-DIS
	python/merge-DIS DIS-NoME
	rivet-mkhtml -o Rivet-DIS DIS.yoda:Hw Powheg-DIS.yoda:Hw-Powheg DIS-NoME.yoda:Hw-NoME
	touch $@



	Rivet-TVT-EW : Rivet-TVT-EW/done
	touch $@

	Rivet-TVT-EW/done: $(shell echo Rivet{,-Powheg}-TVT-{Run-I-Z,Run-I-W,Run-I-WZ,Run-II-Z-{e,{,LowMass-,HighMass-}mu},Run-II-W}.yoda)
	rm -rf Rivet-TVT-EW
	python/merge-TVT-EW TVT
	python/merge-TVT-EW Powheg-TVT
	rivet-mkhtml -o Rivet-TVT-EW TVT-EW.yoda:Hw Powheg-TVT-EW.yoda:Hw-Powheg
	touch $@

	Rivet-TVT-Photon : Rivet-TVT-Photon/done
	touch $@

	Rivet-TVT-Photon/done: $(shell echo Rivet{,-Powheg}-TVT-Run-II-{DiPhoton-GammaGamma,DiPhoton-GammaJet}.yoda Rivet-TVT-Run-II-PromptPhoton.yoda)
	rm -rf Rivet-TVT-Photon
	python/merge-TVT-Photon TVT
	python/merge-TVT-Photon Powheg-TVT
	rivet-mkhtml -o Rivet-TVT-Photon TVT-Photon.yoda:Hw Powheg-TVT-Photon.yoda:Hw-Powheg
	touch $@



	Rivet-TVT-Jets : Rivet-TVT-Jets/done
	touch $@

	Rivet-TVT-Jets/done: $(shell echo Rivet-TVT-{Run-II-Jets-{0..11},Run-I-Jets-{1..8}}.yoda ) \
	$(shell echo Rivet-TVT-{630-Jets-{1..3},300-Jets-1,900-Jets-1}.yoda ) \
	$(shell echo Rivet-TVT-{Run-I,Run-II,300,630,900}-UE.yoda)
	rm -rf Rivet-TVT-Jets
	python/merge-TVT-Jets TVT
	rivet-mkhtml -o Rivet-TVT-Jets TVT-Jets.yoda:Hw
	touch $@

	Rivet-Star : Rivet-Star/done
	touch $@

	Rivet-Star/done : $(shell echo Rivet-Star-{UE,Jets-{1..4}}.yoda )
	rm -rf Rivet-Star
	python/merge-Star Star
	rivet-mkhtml -o Rivet-Star Star.yoda:Hw
	touch $@



	Rivet-SppS : Rivet-SppS/done
	touch $@

	## $(shell echo Rivet-ISR-{30,44,53,62}-UE.yoda ) \
	## {53,63,200,500,900,546} EHS-UE

	Rivet-SppS/done : $(shell echo Rivet-SppS-{200,500,900,546}-UE.yoda )
	rm -rf Rivet-SppS
	python/merge-SppS SppS
	rivet-mkhtml -o Rivet-SppS SppS.yoda:Hw
	touch $@






	Rivet-LHC-Jets : Rivet-LHC-Jets/done
	touch $@

	Rivet-LHC-Jets/done : \
	$(shell echo Rivet-LHC-7-DiJets-{1..7}-{A,B,C}.yoda ) \
	$(shell echo Rivet-LHC-13-DiJets-{6..11}-A.yoda ) \
	$(shell echo Rivet-LHC-{7,8,13}-Jets-{0..10}.yoda ) \
	$(shell echo Rivet-LHC-2760-Jets-{1..3}.yoda ) \
	$(shell echo Rivet-LHC-{900,2360,2760,7,8,13}-UE.yoda ) \
	$(shell echo Rivet-LHC-{900,7,13}-UE-Long.yoda ) \
	$(shell echo Rivet-LHC-7-Charm-{1..5}.yoda ) \
	$(shell echo Rivet-LHC-7-Bottom-{0..9}.yoda ) \
	$(shell echo Rivet-LHC-{7,8,13}-Top-{L,SL}.yoda ) \
	$(shell echo Rivet-LHC-{8,13}-Top-All.yoda )
	rm -rf Rivet-LHC-Jets
	python/merge-LHC-Jets LHC
	rivet-mkhtml -o Rivet-LHC-Jets LHC-Jets.yoda:Hw
	touch $@



	Rivet-LHC-EW : Rivet-LHC-EW/done
	touch $@

	Rivet-LHC-EW/done: \
	$(shell echo Rivet{,-Powheg}-LHC-{13-Z-{e,mu},{8,13}-W-mu,Z-HighMass{1,2}-e,8-Z-Mass{1..4}-{e,mu},W-{e,mu},Z-{e,mu,mu-SOPHTY},Z-LowMass-{e,mu},Z-MedMass-e,WZ,WW-{emu,ll},ZZ-{ll,lv},{8,13}-WZ,8-ZZ-lv,8-WW-ll,Z-mu-Short}.yoda) \
	$(shell echo Rivet-LHC-{7-W-Jet-{1..3}-e,7-Z-Jet-{0..3}-e,7-Z-Jet-0-mu}.yoda) \
	$(shell echo Rivet-LHC-7-{W,Z}Gamma-{e,mu}.yoda) \
	$(shell echo Rivet-LHC-8-ZGamma-{e,mu}.yoda)
	rm -rf Rivet-LHC-EW;
	python/merge-LHC-EW LHC
	python/merge-LHC-EW Powheg-LHC
	rivet-mkhtml -o Rivet-LHC-EW LHC-EW.yoda:Hw Powheg-LHC-EW.yoda:Hw-Powheg \
	Rivet-LHC-Z-mu-SOPHTY.yoda:Hw Rivet-Powheg-LHC-Z-mu-SOPHTY.yoda:Hw-Powheg
	touch $@




	Rivet-LHC-Photon : Rivet-LHC-Photon/done
	touch $@

	Rivet-LHC-Photon/done: \
	$(shell echo Rivet-LHC-{7,8,13}-PromptPhoton-{1..4}.yoda) \
	Rivet-LHC-GammaGamma-7.yoda \
	$(shell echo Rivet{,-Powheg}-LHC-{7,8}-{DiPhoton-GammaGamma,DiPhoton-GammaJet}.yoda)
	rm -rf Rivet-LHC-Photon
	python/merge-LHC-Photon LHC
	python/merge-LHC-Photon Powheg-LHC
	rivet-mkhtml -o Rivet-LHC-Photon LHC-Photon.yoda:Hw Powheg-LHC-Photon.yoda:Hw-Powheg
	touch $@




	Rivet-LHC-Higgs : Rivet-LHC-Higgs/done
	touch $@

	Rivet-LHC-Higgs/done: \
	$(shell echo Rivet{,-Powheg}-LHC-{ggH,VBF,WH,ZH}.yoda) \
	$(shell echo Rivet{,-Powheg}-LHC-8-{{ggH,VBF,WH,ZH}{,-GammaGamma},ggH-WW}.yoda) \
	Rivet-LHC-ggHJet.yoda
	yodamerge --add Rivet-Powheg-LHC-8-{ggH{-GammaGamma,-WW,},{VBF,ZH,WH}{,-GammaGamma}}.yoda -o Powheg-LHC-Higgs.yoda
	yodamerge --add Rivet-LHC-8-{ggH{-GammaGamma,-WW,},{VBF,ZH,WH}{,-GammaGamma}}.yoda -o LHC-Higgs.yoda
	rm -rf Rivet-LHC-Higgs
	rivet-mkhtml -o Rivet-LHC-Higgs Powheg-LHC-Higgs.yoda:Hw-Powheg LHC-Higgs.yoda:Hw\
	Rivet-Powheg-LHC-ggH.yoda:gg-Powheg Rivet-LHC-ggH.yoda:gg Rivet-LHC-ggHJet.yoda:HJet \
	Rivet-Powheg-LHC-VBF.yoda:VBF-Powheg Rivet-LHC-VBF.yoda:VBF Rivet-LHC-WH.yoda:WH Rivet-LHC-ZH.yoda:ZH \
	Rivet-Powheg-LHC-WH.yoda:WH-Powheg Rivet-Powheg-LHC-ZH.yoda:ZH-Powheg
	touch $@





	clean-local:
	rm -f .out .log .tex .top .run .dump .mult .Bmult .yoda Rivet/.in
	rm -rf Rivet-*

	distclean-local:
	rm -rf .cache
	diff --git a/Tests/Rivet/BFactory/BFactory-10.47.in b/Tests/Rivet/BFactory/BFactory-10.47.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/BFactory/BFactory-10.47.in
	@@ -0,0 +1,10 @@
	+create ThePEG::LuminosityFunction /Herwig/EventHandlers/BFactoryLuminosity
	+set /Herwig/EventHandlers/BFactoryLuminosity:BeamEMaxA 5.235
	+set /Herwig/EventHandlers/BFactoryLuminosity:BeamEMaxB 5.235
	+set /Herwig/Generators/EventGenerator:EventHandler:Cuts:MHatMin 9.99
	+set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF
	+set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction /Herwig/EventHandlers/BFactoryLuminosity
	+
	+# BELLE charm hadron production
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ARGUS_1992_I319102
	diff --git a/Tests/Rivet/GammaGamma/GammaGamma-mumu-12.5.in b/Tests/Rivet/GammaGamma/GammaGamma-mumu-12.5.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/GammaGamma/GammaGamma-mumu-12.5.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 12.5
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127:PROCESS=GG
	diff --git a/Tests/Rivet/GammaGamma/GammaGamma-mumu-17.5.in b/Tests/Rivet/GammaGamma/GammaGamma-mumu-17.5.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/GammaGamma/GammaGamma-mumu-17.5.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 17.5
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127:PROCESS=GG
	diff --git a/Tests/Rivet/GammaGamma/GammaGamma-mumu-3.5.in b/Tests/Rivet/GammaGamma/GammaGamma-mumu-3.5.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/GammaGamma/GammaGamma-mumu-3.5.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 3.5
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127:PROCESS=GG
	diff --git a/Tests/Rivet/GammaGamma/GammaGamma-mumu-30.0.in b/Tests/Rivet/GammaGamma/GammaGamma-mumu-30.0.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/GammaGamma/GammaGamma-mumu-30.0.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 30.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127:PROCESS=GG
	diff --git a/Tests/Rivet/GammaGamma/GammaGamma-mumu-4.5.in b/Tests/Rivet/GammaGamma/GammaGamma-mumu-4.5.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/GammaGamma/GammaGamma-mumu-4.5.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 4.5
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127:PROCESS=GG
	diff --git a/Tests/Rivet/GammaGamma/GammaGamma-mumu-5.5.in b/Tests/Rivet/GammaGamma/GammaGamma-mumu-5.5.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/GammaGamma/GammaGamma-mumu-5.5.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 5.5
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127:PROCESS=GG
	diff --git a/Tests/Rivet/GammaGamma/GammaGamma-mumu-6.5.in b/Tests/Rivet/GammaGamma/GammaGamma-mumu-6.5.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/GammaGamma/GammaGamma-mumu-6.5.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 6.5
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127:PROCESS=GG
	diff --git a/Tests/Rivet/GammaGamma/GammaGamma-mumu-7.5.in b/Tests/Rivet/GammaGamma/GammaGamma-mumu-7.5.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/GammaGamma/GammaGamma-mumu-7.5.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 7.5
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127:PROCESS=GG
	diff --git a/Tests/Rivet/GammaGamma/GammaGamma-mumu-9.0.in b/Tests/Rivet/GammaGamma/GammaGamma-mumu-9.0.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/GammaGamma/GammaGamma-mumu-9.0.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 9.0
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127:PROCESS=GG
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-Jets-198.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-Jets-198.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-Jets-198.in
	@@ -0,0 +1,12 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 198.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 jets
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I661114
	+# OPAL jets
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2003_I611415
	\ No newline at end of file
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-Jets-206.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-Jets-206.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-Jets-206.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 206.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2008_I754316
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-161.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-161.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-161.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 161.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-172.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-172.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-172.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 172.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-183.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-183.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-183.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 183.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-189.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-189.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-189.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 189.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-196.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-196.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-196.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 196.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-206.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-206.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-mumu-206.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 206.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-tautau-189.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-tautau-189.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-tautau-189.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 189.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-tautau-196.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-tautau-196.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-tautau-196.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 196.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-tautau-206.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-tautau-206.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Direct-tautau-206.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 206.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I645127
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Double-Resolved-Jets-198.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Double-Resolved-Jets-198.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Double-Resolved-Jets-198.in
	@@ -0,0 +1,12 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 198.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 jets
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I661114
	+# OPAL jets
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2003_I611415
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Double-Resolved-Jets-206.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Double-Resolved-Jets-206.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Double-Resolved-Jets-206.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 206.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2008_I754316
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Single-Resolved-Jets-198.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Single-Resolved-Jets-198.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Single-Resolved-Jets-198.in
	@@ -0,0 +1,12 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 198.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 jets
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I661114
	+# OPAL jets
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2003_I611415
	diff --git a/Tests/Rivet/LEP-Gamma/LEP-Gamma-Single-Resolved-Jets-206.in b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Single-Resolved-Jets-206.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/LEP-Gamma/LEP-Gamma-Single-Resolved-Jets-206.in
	@@ -0,0 +1,10 @@
	+# -- ThePEG-repository --
	+##################################################
	+# LEP physics parameters (override defaults)
	+##################################################
	+set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 206.
	+##################################################
	+# select the analyses
	+##################################################
	+# L3 mu+mu-
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2008_I754316
	diff --git a/Tests/Rivet/LEP/LEP-91.in b/Tests/Rivet/LEP/LEP-91.in
	--- a/Tests/Rivet/LEP/LEP-91.in
	+++ b/Tests/Rivet/LEP/LEP-91.in
	@@ -1,87 +1,89 @@
	# -- ThePEG-repository --
	##################################################
	# LEP physics parameters (override defaults)
	##################################################
	set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 91.2
	##################################################
	# select the analyses
	##################################################
	# Validated
	##################################################
	# ALEPH charged particle multiplicity
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1991_S2435284
	# OPAL charged particle multiplicity
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1992_I321190
	# DELPHI charged particle multiplicity
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1991_I301657
	# ALEPH main LEP I QCD summary paper
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1996_S3486095
	# ALEPH D*
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1999_S4193598
	+# OPAL D*
	+insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1995_I382219
	# OPAL charged hadron analysis
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1994_S2927284
	# OPAL Delta++ analysis
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1995_S3198391
	# OPAL J/Psi analysis analysis
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1996_S3257789
	# ALEPH eta/omega analysis
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_2002_S4823664
	# OPAL K*0 analysis
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1997_S3608263
	# OPAL flavour specific charged multiplicities etc
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1998_S3780481
	# OPAL f_0,f_2 and phi production
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1998_S3702294
	# OPAL gamma,pi0,eta,eta',rho+/-,a0+/-
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1998_S3749908
	# OPAL K0
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2000_S4418603
	# OPAL K* +/-
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1993_I342766
	# SLD flavour specific charged multiplicities etc
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_1996_S3398250
	# SLD flavour specific charged multiplicities etc
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_1999_S3743934
	# SLD flavour specific charged multiplicities etc
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_2004_S5693039
	# OPAL event shapes and multiplicities at different energies
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2004_S6132243
	# ALEPH jet and event shapes at many energies
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_2004_S5765862
	# OPAL/JADE jet rates at many energies
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 JADE_OPAL_2000_S4300807
	# DELPHI strange baryon production
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1995_S3137023
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_2000_I524694
	# DELPHI f_0, rho_0 and f_2 production
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1999_S3960137
	# OPAL strange baryon production
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1997_S3396100
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1997_I421977
	# DELPHI tuning paper
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1996_S3430090
	# DELPHI b quark
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_2011_I890503
	# ALEPH b quark
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_2001_S4656318
	# SLD b quark
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_2002_S4869273
	# OPAL b quark
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2003_I599181
	# PDG hadron multiplicities and ratios
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 PDG_HADRON_MULTIPLICITIES
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 PDG_HADRON_MULTIPLICITIES_RATIOS
	# OPAL from gluon paper
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2004_I648738
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2004_I631361:PROCESS=QQ
	# L3 eta
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_1992_I336180
	# L3 jet rates
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I652683
	# ALEPH pi+-, K+- and (p, anti-p)
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1995_I382179
	##################################################
	# unvalidated
	##################################################
	# OPAL 4 jet angles
	insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2001_S4553896
	diff --git a/Tests/Rivet/Templates/GammaGamma.in b/Tests/Rivet/Templates/GammaGamma.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/Templates/GammaGamma.in
	@@ -0,0 +1,44 @@
	+# -- ThePEG-repository --
	+#
	+# DO NOT EDIT - autogenerated by make_input_files.py
	+#
	+##################################################
	+# base parameters for LEP analyses
	+##################################################
	+read snippets/EECollider.in
	+set /Herwig/EventHandlers/EventHandler:BeamA /Herwig/Particles/gamma
	+set /Herwig/EventHandlers/EventHandler:BeamB /Herwig/Particles/gamma
	+set /Herwig/Particles/gamma:PDF /Herwig/Partons/NoPDF
	+##################################################
	+# Technical parameters for this run
	+##################################################
	+cd /Herwig/Generators
	+set EventGenerator:MaxErrors 1000000
	+
	+##################################################
	+## Shower and scheme selection.
	+## Should be empty unless using the Dipole Shower.
	+##################################################
	+${shower}
	+
	+##################################################
	+# Create the Herwig analysis
	+##################################################
	+cd /Herwig/Generators
	+create ThePEG::RivetAnalysis /Herwig/Analysis/RivetAnalysis RivetAnalysis.so
	+insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
	+
	+##################################################
	+# Use the q qbar matrix element
	+##################################################
	+# default e+e- > q qbar (5 flavours d,u,s,c,b)
	+${process}
	+set /Herwig/Shower/PartnerFinder:QEDPartner IIandFF
	+
	+read ${parameterFile}
	+
	+##################################################
	+# Save run for later usage with 'Herwig run'
	+##################################################
	+cd /Herwig/Generators
	+saverun ${runname} EventGenerator
	diff --git a/Tests/Rivet/Templates/LEP-Gamma-Direct.in b/Tests/Rivet/Templates/LEP-Gamma-Direct.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/Templates/LEP-Gamma-Direct.in
	@@ -0,0 +1,46 @@
	+# -- ThePEG-repository --
	+#
	+# DO NOT EDIT - autogenerated by make_input_files.py
	+#
	+##################################################
	+# base parameters for LEP analyses
	+##################################################
	+read snippets/EECollider.in
	+##################################################
	+# Technical parameters for this run
	+##################################################
	+cd /Herwig/Generators
	+set EventGenerator:MaxErrors 1000000
	+
	+##################################################
	+# Switch off ISR
	+##################################################
	+set /Herwig/Particles/e+:PDF /Herwig/Partons/WWPDF
	+set /Herwig/Particles/e-:PDF /Herwig/Partons/WWPDF
	+##################################################
	+## Shower and scheme selection.
	+## Should be empty unless using the Dipole Shower.
	+##################################################
	+${shower}
	+
	+##################################################
	+# Create the Herwig analysis
	+##################################################
	+cd /Herwig/Generators
	+create ThePEG::RivetAnalysis /Herwig/Analysis/RivetAnalysis RivetAnalysis.so
	+insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
	+
	+##################################################
	+# Use the q qbar matrix element
	+##################################################
	+# default e+e- > q qbar (5 flavours d,u,s,c,b)
	+${process}
	+set /Herwig/Shower/PartnerFinder:QEDPartner IIandFF
	+
	+read ${parameterFile}
	+
	+##################################################
	+# Save run for later usage with 'Herwig run'
	+##################################################
	+cd /Herwig/Generators
	+saverun ${runname} EventGenerator
	diff --git a/Tests/Rivet/Templates/LEP-Gamma-Double-Resolved.in b/Tests/Rivet/Templates/LEP-Gamma-Double-Resolved.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/Templates/LEP-Gamma-Double-Resolved.in
	@@ -0,0 +1,47 @@
	+# -- ThePEG-repository --
	+#
	+# DO NOT EDIT - autogenerated by make_input_files.py
	+#
	+##################################################
	+# base parameters for LEP analyses
	+##################################################
	+read snippets/EECollider.in
	+##################################################
	+# Technical parameters for this run
	+##################################################
	+cd /Herwig/Generators
	+set EventGenerator:MaxErrors 1000000
	+
	+##################################################
	+# Switch off ISR
	+##################################################
	+set /Herwig/Particles/e+:PDF /Herwig/Partons/WWPDF
	+set /Herwig/Particles/e-:PDF /Herwig/Partons/WWPDF
	+set /Herwig/Particles/gamma:PDF /Herwig/Partons/SaSPDF
	+##################################################
	+## Shower and scheme selection.
	+## Should be empty unless using the Dipole Shower.
	+##################################################
	+${shower}
	+
	+##################################################
	+# Create the Herwig analysis
	+##################################################
	+cd /Herwig/Generators
	+create ThePEG::RivetAnalysis /Herwig/Analysis/RivetAnalysis RivetAnalysis.so
	+insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
	+
	+##################################################
	+# Use the q qbar matrix element
	+##################################################
	+# default e+e- > q qbar (5 flavours d,u,s,c,b)
	+${process}
	+set /Herwig/Shower/PartnerFinder:QEDPartner All
	+
	+read ${parameterFile}
	+
	+##################################################
	+# Save run for later usage with 'Herwig run'
	+##################################################
	+cd /Herwig/Generators
	+saverun ${runname} EventGenerator
	diff --git a/Tests/Rivet/Templates/LEP-Gamma-Single-Resolved.in b/Tests/Rivet/Templates/LEP-Gamma-Single-Resolved.in
	new file mode 100644
	--- /dev/null
	+++ b/Tests/Rivet/Templates/LEP-Gamma-Single-Resolved.in
	@@ -0,0 +1,72 @@
	+# -- ThePEG-repository --
	+#
	+# DO NOT EDIT - autogenerated by make_input_files.py
	+#
	+##################################################
	+# base parameters for LEP analyses
	+##################################################
	+read snippets/EECollider.in
	+##################################################
	+# Technical parameters for this run
	+##################################################
	+cd /Herwig/Generators
	+set EventGenerator:MaxErrors 1000000
	+
	+##################################################
	+# Select the PDFs
	+##################################################
	+set /Herwig/Particles/e+:PDF /Herwig/Partons/WWPDF
	+set /Herwig/Particles/e-:PDF /Herwig/Partons/WWPDF
	+set /Herwig/Particles/gamma:PDF /Herwig/Partons/SaSPDF
	+
	+
	+DISABLEREADONLY
	+create Herwig::MultiPartonExtractor /Herwig/Partons/GGExtractor HwPartonExtractor.so
	+newdef /Herwig/Partons/GGExtractor:NoPDF /Herwig/Partons/NoPDF
	+insert /Herwig/Partons/GGExtractor:FirstPDFs 0 /Herwig/Partons/WWPDF
	+insert /Herwig/Partons/GGExtractor:FirstPDFs 1 /Herwig/Partons/SaSPDF
	+insert /Herwig/Partons/GGExtractor:SecondPDFs 0 /Herwig/Partons/WWPDF
	+insert /Herwig/Partons/GGExtractor:SecondPDFs 1 /Herwig/Partons/NoPDF
	+
	+create Herwig::MultiPartonExtractor /Herwig/Partons/GG2Extractor HwPartonExtractor.so
	+newdef /Herwig/Partons/GG2Extractor:NoPDF /Herwig/Partons/NoPDF
	+insert /Herwig/Partons/GG2Extractor:FirstPDFs 0 /Herwig/Partons/WWPDF
	+insert /Herwig/Partons/GG2Extractor:FirstPDFs 1 /Herwig/Partons/NoPDF
	+insert /Herwig/Partons/GG2Extractor:SecondPDFs 0 /Herwig/Partons/WWPDF
	+insert /Herwig/Partons/GG2Extractor:SecondPDFs 1 /Herwig/Partons/SaSPDF
	+
	+
	+set /Herwig/MatrixElements/SubProcess:PartonExtractor /Herwig/Partons/GGExtractor
	+cp /Herwig/MatrixElements/SubProcess /Herwig/MatrixElements/SubProcess2
	+set /Herwig/MatrixElements/SubProcess2:PartonExtractor /Herwig/Partons/GG2Extractor
	+
	+insert /Herwig/Generators/EventGenerator:EventHandler:SubProcessHandlers 1 /Herwig/MatrixElements/SubProcess2
	+
	+
	+##################################################
	+## Shower and scheme selection.
	+## Should be empty unless using the Dipole Shower.
	+##################################################
	+${shower}
	+
	+##################################################
	+# Create the Herwig analysis
	+##################################################
	+cd /Herwig/Generators
	+create ThePEG::RivetAnalysis /Herwig/Analysis/RivetAnalysis RivetAnalysis.so
	+insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis
	+
	+##################################################
	+# Use the q qbar matrix element
	+##################################################
	+# default e+e- > q qbar (5 flavours d,u,s,c,b)
	+${process}
	+set /Herwig/Shower/PartnerFinder:QEDPartner All
	+
	+read ${parameterFile}
	+
	+##################################################
	+# Save run for later usage with 'Herwig run'
	+##################################################
	+cd /Herwig/Generators
	+saverun ${runname} EventGenerator
	diff --git a/Tests/python/make_input_files.py b/Tests/python/make_input_files.py
	--- a/Tests/python/make_input_files.py
	+++ b/Tests/python/make_input_files.py
	@@ -1,1876 +1,1927 @@
	#! /usr/bin/env python
	import logging,sys,os
	from string import strip, Template

	import sys
	if sys.version_info[:3] < (2,4,0):
	print "rivet scripts require Python version >= 2.4.0... exiting"
	sys.exit(1)

	if __name__ == "__main__":
	import logging
	from optparse import OptionParser, OptionGroup
	parser = OptionParser(usage="%prog name [...]")


	simulation=""

	numberOfAddedProcesses=0
	def addProcess(thefactory,theProcess,Oas,Oew,scale,mergedlegs,NLOprocesses):
	global numberOfAddedProcesses
	global simulation
	numberOfAddedProcesses+=1
	res ="set "+thefactory+":OrderInAlphaS "+Oas+"\n"
	res+="set "+thefactory+":OrderInAlphaEW "+Oew+"\n"
	res+="do "+thefactory+":Process "+theProcess+" "
	if ( mergedlegs != 0 ):
	if simulation!="Merging":
	print "simulation is not Merging, trying to add merged legs."
	sys.exit(1)
	res+="["
	for j in range(mergedlegs):
	res+=" j "
	res+="]"
	res+="\n"
	if (NLOprocesses!=0):
	if simulation!="Merging":
	print "simulation is not Merging, trying to add NLOProcesses."
	sys.exit(1)
	res+="set MergingFactory:NLOProcesses %s \n" % NLOprocesses
	if ( scale != "" ):
	res+="set "+thefactory+":ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/"+scale+"\n"
	return res


	def addLeptonPairCut(minmass,maxmass):
	return "set /Herwig/Cuts/LeptonPairMassCut:MinMass "+minmass+"GeV\nset /Herwig/Cuts/LeptonPairMassCut:MaxMass "+maxmass+"GeV\n"

	didaddfirstjet=False
	def addFirstJet(ptcut):
	global didaddfirstjet
	if(didaddfirstjet):
	logging.error("Can only add jetcut once.")
	sys.exit(1)

	res="set /Herwig/Cuts/Cuts:JetFinder /Herwig/Cuts/JetFinder\n"
	res+="insert /Herwig/Cuts/Cuts:MultiCuts 0 /Herwig/Cuts/JetCuts\n"
	res+="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/FirstJet\n"
	if(ptcut!=""):
	res+="set /Herwig/Cuts/FirstJet:PtMin "+ptcut+".*GeV\n"
	didaddfirstjet=True
	return res



	didaddsecondjet=False
	def addSecondJet(ptcut):
	global didaddsecondjet
	if(didaddsecondjet):
	logging.error("Can only add second jetcut once.")
	sys.exit(1)
	res="insert /Herwig/Cuts/JetCuts:JetRegions 0 /Herwig/Cuts/SecondJet\n"
	res+="set /Herwig/Cuts/SecondJet:PtMin "+ptcut+".*GeV\n"
	didaddsecondjet=True
	return res


	didaddjetpair=False
	def addJetPairCut(minmass):
	global didaddjetpair
	if(didaddjetpair):
	logging.error("Can only add second jetcut once.")
	sys.exit(1)
	res="""\
	create ThePEG::JetPairRegion /Herwig/Cuts/JetPairMass JetCuts.so
	set /Herwig/Cuts/JetPairMass:FirstRegion /Herwig/Cuts/FirstJet
	set /Herwig/Cuts/JetPairMass:SecondRegion /Herwig/Cuts/SecondJet
	insert /Herwig/Cuts/JetCuts:JetPairRegions 0 /Herwig/Cuts/JetPairMass
	set /Herwig/Cuts/JetPairMass:MassMin {mm}.*GeV
	""".format(mm=minmass)
	didaddjetpair=True
	return res

	addedBRReweighter=False
	def addBRReweighter():
	global addedBRReweighter
	if(addedBRReweighter):
	logging.error("Can only add BRReweighter once.")
	sys.exit(1)
	res="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
	res+="insert /Herwig/Generators/EventGenerator:EventHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
	addedBRReweighter=True
	return res

	def setHardProcessWidthToZero(list1):
	res=""
	for i in list1:
	res+="set /Herwig/Particles/"+i+":HardProcessWidth 0.\n"
	return res

	selecteddecaymode=False
	def selectDecayMode(particle,decaymodes):
	global selecteddecaymode
	res="do /Herwig/Particles/"+particle+":SelectDecayModes"
	for decay in decaymodes:
	res+=" /Herwig/Particles/"+particle+"/"+decay
	res+="\n"
	selecteddecaymode=True
	return res

	def jet_kt_cut(energy):
	return "set /Herwig/Cuts/JetKtCut:MinKT {E}*GeV\n".format(E=energy)

	def mhatmin_cut(energy):
	return "set /Herwig/Cuts/Cuts:MHatMin {E}*GeV\n".format(E=energy)

	def mhat_minm_maxm(e1,e2,e3):
	return """\
	set /Herwig/Cuts/Cuts:MHatMin {e1}*GeV
	set /Herwig/Cuts/MassCut:MinM {e2}*GeV
	set /Herwig/Cuts/MassCut:MaxM {e3}*GeV
	""".format(**locals())

	def collider_lumi(energy):
	return "set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy {E}*GeV\n".format(E=energy)

	-def insert_ME(me,process=None,ifname='Process'):
	- result = "insert /Herwig/MatrixElements/SubProcess:MatrixElements 0 /Herwig/MatrixElements/{me}\n".format(**locals())
	+def insert_ME(me,process=None,ifname='Process',subprocess="SubProcess"):
	+ result = "insert /Herwig/MatrixElements/{subprocess}:MatrixElements 0 /Herwig/MatrixElements/{me}\n".format(**locals())
	if process is not None:
	result += "set /Herwig/MatrixElements/{me}:{ifname} {process}".format(**locals())
	return result

	def particlegroup(name,*particles):
	result = ["do /Herwig/MatrixElements/Matchbox/Factory:StartParticleGroup {n}".format(n=name)]
	for p in particles:
	result.append(
	"insert /Herwig/MatrixElements/Matchbox/Factory:ParticleGroup 0 /Herwig/Particles/{p}".format(p=p)
	)
	result.append("do /Herwig/MatrixElements/Matchbox/Factory:EndParticleGroup")
	return '\n'.join(result)

	# settings for four flavour scheme
	fourFlavour="""
	read Matchbox/FourFlavourScheme.in
	{bjetgroup}
	set /Herwig/Cuts/MatchboxJetMatcher:Group bjet
	""".format(bjetgroup=particlegroup('bjet','b','bbar','c', 'cbar',
	's','sbar','d','dbar','u','ubar','g'))

	ME_Upsilon = """\
	create Herwig::MEee2VectorMeson /Herwig/MatrixElements/MEUpsilon HwMELepton.so
	set /Herwig/MatrixElements/MEUpsilon:VectorMeson /Herwig/Particles/Upsilon(4S)
	set /Herwig/MatrixElements/MEUpsilon:Coupling 96.72794
	""" + insert_ME("MEUpsilon")


	(opts, args) = parser.parse_args()
	## Check args
	if len(args) != 1:
	logging.error("Must specify at least input file")
	sys.exit(1)

	name = args[0]
	print name


	# select the template to load
	# collider
	KNOWN_COLLIDERS = [
	+ "GammaGamma",
	+ "LEP-Gamma",
	"BFactory",
	"LEP",
	"DIS",
	"TVT",
	"LHC-GammaGamma",
	"LHC",
	"ISR",
	"SppS",
	"Star",
	"EHS",
	]
	collider = ""
	for cand_collider in KNOWN_COLLIDERS:
	if cand_collider in name:
	collider = cand_collider
	break
	del cand_collider
	assert collider
	have_hadronic_collider = collider in ["TVT","LHC","ISR","SppS","Star","EHS"]

	-
	thefactory="Factory"

	parameters = {
	'shower' : '',
	'bscheme' : '',
	}
	# istart determines how many name parts need to be skipped
	istart = 1
	# Dipole shower with Matchbox Powheg
	if "Dipole-Matchbox-Powheg" in name :
	istart = 4
	simulation="Matchbox"
	parameters["shower"] = "read Matchbox/Powheg-DipoleShower.in\n"

	# Dipole shower with internal Powheg - Todo: Finish modifying template files.
	'''
	elif "Dipole-Powheg" in name :
	istart = 3
	simulation="Powheg"
	parameters["shower"] = "set /Herwig/EventHandlers/EventHandler:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler\nread snippets/Dipole_AutoTune_prel.in\n"
	'''

	# Dipole shower with MCatNLO
	elif "Dipole-MCatNLO" in name :
	istart = 3
	simulation="Matchbox"
	parameters["shower"] = "read Matchbox/MCatNLO-DipoleShower.in\n"

	# Dipole shower with Matchbox LO
	elif "Dipole-Matchbox-LO" in name :
	istart = 4
	simulation="Matchbox"
	parameters["shower"] = "read Matchbox/LO-DipoleShower.in\n"

	# Dipole shower with internal LO
	elif "Dipole" in name :
	istart = 2
	simulation=""
	parameters["shower"] = "set /Herwig/EventHandlers/EventHandler:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler\nread snippets/Dipole_AutoTune_prel.in\n"

	# AO shower with Matchbox Powheg
	elif "Matchbox-Powheg" in name :
	istart = 3
	simulation="Matchbox"
	parameters["shower"] = "read Matchbox/Powheg-DefaultShower.in\n"

	# AO shower with MCatNLO
	elif "Matchbox" in name :
	istart = 2
	simulation="Matchbox"
	parameters["shower"] = "read Matchbox/MCatNLO-DefaultShower.in\n"

	# AO shower with inernal Powheg
	elif "Powheg" in name :
	istart = 2
	simulation="Powheg"

	# Dipole shower with merging
	elif "Merging" in name :
	istart = 2
	simulation="Merging"
	thefactory="MergingFactory"

	# Flavour settings for Matchbox
	if simulation=="Matchbox" :
	parameters["bscheme"] = "read Matchbox/FiveFlavourScheme.in\n"

	if "Dipole" in parameters["shower"] :
	parameters["bscheme"] += "read Matchbox/FiveFlavourNoBMassScheme.in\n"

	if collider not in ['DIS','LEP'] :
	parameters["nlo"] = "read Matchbox/MadGraph-OpenLoops.in\n"

	# Flavour settings for dipole shower with internal ME
	if simulation=="" and "Dipole" in parameters["shower"] :
	parameters["bscheme"] = "read snippets/DipoleShowerFiveFlavours.in"



	# find the template
	if simulation=="" :
	if collider=="LHC-GammaGamma" :
	istart += 1
	templateName="Hadron-Gamma.in"
	elif have_hadronic_collider :
	templateName="Hadron.in"
	+ elif collider=="LEP-Gamma" :
	+ istart+=1
	+ if("Direct" in name) :
	+ templateName="LEP-Gamma-Direct.in"
	+ elif("Single-Resolved" in name) :
	+ templateName="LEP-Gamma-Single-Resolved.in"
	+ elif("Double-Resolved" in name) :
	+ templateName="LEP-Gamma-Double-Resolved.in"
	+ else :
	+ print "Unknown type of LEP-Gamma event ",name
	+ quit()
	+ elif collider=="GammaGamma" :
	+ templateName="GammaGamma.in"
	elif collider != "BFactory" :
	templateName= "%s.in" % collider
	else :
	templateName= "LEP.in"
	else :
	if have_hadronic_collider :
	templateName= "Hadron-%s.in" % simulation
	elif collider != "BFactory" :
	templateName= "%s-%s.in" % (collider,simulation)
	else :
	templateName= "LEP-%s.in" % simulation

	# work out the name of the parameter file
	parameterName="-".join(name.split("-")[istart:])
	del istart

	class StringBuilder(object):
	"""
	Avoid expensive string additions until the end
	by building up a list first.

	This helper class avoids rewriting all the += lower down
	to list operations.
	"""
	def __init__(self, init = None):
	self.lines = [] if init is None else [init]

	def __iadd__(self, line):
	self.lines.append(line)
	return self

	def __str__(self):
	return '\n'.join(self.lines)

	# work out the process and parameters
	process=StringBuilder()
	# Bfactory
	if(collider=="BFactory") :
	if(simulation=="") :
	if(parameterName=="10.58-res") :
	process += ME_Upsilon
	elif(parameterName=="10.58") :
	process += ME_Upsilon
	process += "set /Herwig/MatrixElements/MEee2gZ2qq:MaximumFlavour 4\n"
	else :
	process+=insert_ME("MEee2gZ2qq")
	process+= "set /Herwig/MatrixElements/MEee2gZ2qq:MaximumFlavour 4\n"
	elif(simulation=="Powheg") :
	process = StringBuilder("set /Herwig/MatrixElements/PowhegMEee2gZ2qq:MaximumFlavour 4\n")
	elif(simulation=="Matchbox" ) :
	process = StringBuilder(addProcess(thefactory,"e- e+ -> u ubar","0","2","",0,0))
	process+=addProcess(thefactory,"e- e+ -> d dbar","0","2","",0,0)
	process+=addProcess(thefactory,"e- e+ -> c cbar","0","2","",0,0)
	process+=addProcess(thefactory,"e- e+ -> s sbar","0","2","",0,0)
	elif(simulation=="Merging" ) :
	logging.warning("BFactory not explicitly tested for %s " % simulation)
	sys.exit(0)

	# DIS
	elif(collider=="DIS") :
	if(simulation=="") :
	if "NoME" in parameterName :
	process = StringBuilder("set /Herwig/Shower/ShowerHandler:HardEmission None")
	parameterName=parameterName.replace("NoME-","")
	else :
	process = StringBuilder("")
	elif(simulation=="Powheg") :
	process = StringBuilder("")
	elif(simulation=="Matchbox" ) :
	if "e-" in parameterName :
	process = StringBuilder(addProcess(thefactory,"e- p -> e- j","0","2","",0,0))
	else :
	process = StringBuilder(addProcess(thefactory,"e+ p -> e+ j","0","2","",0,0))
	elif(simulation=="Merging" ) :
	if "e-" in parameterName :
	process = StringBuilder(addProcess(thefactory,"e- p -> e- j","0","2","",2,2))
	else :
	process = StringBuilder(addProcess(thefactory,"e+ p -> e+ j","0","2","",2,2))

	# LEP
	elif(collider=="LEP") :
	if(simulation=="") :
	if "gg" in parameterName :
	process = StringBuilder("create Herwig::MEee2Higgs2SM /Herwig/MatrixElements/MEee2Higgs2SM\n")
	process+=insert_ME("MEee2Higgs2SM","Gluon","Allowed")
	else :
	process = StringBuilder(insert_ME("MEee2gZ2qq"))
	if(parameterName=="10") :
	process+="set /Herwig/MatrixElements/MEee2gZ2qq:MaximumFlavour 4"
	elif(simulation=="Powheg") :
	process = StringBuilder()
	if(parameterName=="10") :
	process = StringBuilder("set /Herwig/MatrixElements/PowhegMEee2gZ2qq:MaximumFlavour 4")
	elif(simulation=="Matchbox" ) :
	if(parameterName=="10") :
	process = StringBuilder(addProcess(thefactory,"e- e+ -> u ubar","0","2","",0,0))
	process+=addProcess(thefactory,"e- e+ -> d dbar","0","2","",0,0)
	process+=addProcess(thefactory,"e- e+ -> c cbar","0","2","",0,0)
	process+=addProcess(thefactory,"e- e+ -> s sbar","0","2","",0,0)
	else :
	process = StringBuilder(addProcess(thefactory,"e- e+ -> j j","0","2","",0,0))

	elif(simulation=="Merging" ) :
	if(parameterName=="10") :
	process = StringBuilder(addProcess(thefactory,"e- e+ -> j j","0","2","",2,2))
	process+="read Matchbox/FourFlavourScheme.in"
	else :
	process = StringBuilder(addProcess(thefactory,"e- e+ -> j j","0","2","",2,2))
	+# LEP-Gamma
	+elif(collider=="LEP-Gamma") :
	+ if(simulation=="") :
	+ if("mumu" in parameterName) :
	+ process = StringBuilder(insert_ME("MEgg2ff","Muon"))
	+ process +="set /Herwig/Cuts/Cuts:MHatMin 3.\n"
	+ elif( "tautau" in parameterName) :
	+ process = StringBuilder(insert_ME("MEgg2ff","Tau"))
	+ process +="set /Herwig/Cuts/Cuts:MHatMin 3.\n"
	+ elif( "Jets" in parameterName) :
	+ if("Direct" in parameterName ) :
	+ process = StringBuilder(insert_ME("MEgg2ff","Quarks"))
	+ elif("Single-Resolved" in parameterName ) :
	+ process = StringBuilder(insert_ME("MEGammaP2Jets",None,"Process","SubProcess"))
	+ process+= insert_ME("MEGammaP2Jets",None,"Process","SubProcess2")
	+ else :
	+ process = StringBuilder(insert_ME("MEQCD2to2"))
	+ process+="insert /Herwig/Cuts/Cuts:OneCuts[0] /Herwig/Cuts/JetKtCut"
	+ process+="set /Herwig/Cuts/JetKtCut:MinKT 3."
	+ else :
	+ print "process not supported for Gamma Gamma processes at LEP"
	+ quit()
	+ else :
	+ print "Only internal matrix elements currently supported for Gamma Gamma processes at LEP"
	+ quit()
	+elif(collider=="GammaGamma") :
	+ if(simulation=="") :
	+ if("mumu" in parameterName) :
	+ process = StringBuilder(insert_ME("MEgg2ff"))
	+ process +="set /Herwig/MatrixElements/MEgg2ff:Process Muon\n"
	+ process +="set /Herwig/Cuts/Cuts:MHatMin 3.\n"
	+ else :
	+ print "process not supported for Gamma Gamma processes at LEP"
	+ quit()
	+ else :
	+ print "Only internal matrix elements currently supported for Gamma Gamma processes at LEP"
	+ quit()
	# TVT
	elif(collider=="TVT") :
	process = StringBuilder("set /Herwig/Generators/EventGenerator:EventHandler:BeamB /Herwig/Particles/pbar-\n")
	if "Run-II" in parameterName : process+=collider_lumi(1960.0)
	elif "Run-I" in parameterName : process+=collider_lumi(1800.0)
	elif "900" in parameterName : process+=collider_lumi(900.0)
	elif "630" in parameterName : process+=collider_lumi(630.0)
	elif "300" in parameterName : process+=collider_lumi(300.0)

	if(simulation=="") :
	if "PromptPhoton" in parameterName :
	process+=insert_ME("MEGammaJet")
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 15.\n"
	elif "DiPhoton-GammaGamma" in parameterName :
	process+=insert_ME("MEGammaGamma")
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
	parameterName=parameterName.replace("-GammaGamma","")
	elif "DiPhoton-GammaJet" in parameterName :
	process+=insert_ME("MEGammaJet")
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
	parameterName=parameterName.replace("-GammaJet","")
	elif "UE" in parameterName :
	if "Dipole" in parameters["shower"]:
	process+="read snippets/MB-DipoleShower.in\n"
	else:
	process+="read snippets/MB.in\n"
	process+="read snippets/Diffraction.in\n"

	process += "set /Herwig/Decays/DecayHandler:LifeTimeOption 0\n"
	process += "set /Herwig/Decays/DecayHandler:MaxLifeTime 10*mm\n"
	elif "Jets" in parameterName :
	process+=insert_ME("MEQCD2to2")
	process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
	if "Run-II-Jets-10" in parameterName : process+=jet_kt_cut( 30.)+mhatmin_cut(500.)
	elif "Run-II-Jets-11" in parameterName: process+=jet_kt_cut( 30.)+mhatmin_cut(900.)
	elif "Run-I-Jets-1" in parameterName : process+=jet_kt_cut( 20.)
	elif "Run-I-Jets-2" in parameterName : process+=jet_kt_cut( 40.)
	elif "Run-I-Jets-3" in parameterName : process+=jet_kt_cut( 65.)
	elif "Run-I-Jets-4" in parameterName : process+=jet_kt_cut( 90.)
	elif "Run-I-Jets-5" in parameterName : process+=jet_kt_cut(160.)
	elif "Run-I-Jets-6" in parameterName : process+=jet_kt_cut( 30.)+mhatmin_cut(100.)
	elif "Run-I-Jets-7" in parameterName : process+=jet_kt_cut( 30.)+mhatmin_cut(400.)
	elif "Run-I-Jets-8" in parameterName : process+=jet_kt_cut( 30.)+mhatmin_cut(700.)
	elif "Run-II-Jets-0" in parameterName : process+=jet_kt_cut( 15.)
	elif "Run-II-Jets-1" in parameterName : process+=jet_kt_cut( 25.)
	elif "Run-II-Jets-2" in parameterName : process+=jet_kt_cut( 40.)
	elif "Run-II-Jets-3" in parameterName : process+=jet_kt_cut( 60.)
	elif "Run-II-Jets-4" in parameterName : process+=jet_kt_cut( 85.)
	elif "Run-II-Jets-5" in parameterName : process+=jet_kt_cut(110.)
	elif "Run-II-Jets-6" in parameterName : process+=jet_kt_cut(160.)
	elif "Run-II-Jets-7" in parameterName : process+=jet_kt_cut(250.)
	elif "Run-II-Jets-8" in parameterName : process+=jet_kt_cut( 30.)+mhatmin_cut(100.)
	elif "Run-II-Jets-9" in parameterName : process+=jet_kt_cut( 30.)+mhatmin_cut(300.)
	elif "900-Jets-1" in parameterName : process+=jet_kt_cut( 10.)
	elif "300-Jets-1" in parameterName : process+=jet_kt_cut( 6.)
	elif "630-Jets-1" in parameterName : process+=jet_kt_cut( 20.)
	elif "630-Jets-2" in parameterName : process+=jet_kt_cut( 40.)
	elif "630-Jets-3" in parameterName : process+=jet_kt_cut( 75.)
	elif "900-Jets-1" in parameterName : process+=jet_kt_cut( 10.)
	elif "Run-I-WZ" in parameterName :
	process+=insert_ME("MEqq2W2ff","Electron")
	process+=insert_ME("MEqq2gZ2ff","Electron")
	elif "Run-II-W" in parameterName or "Run-I-W" in parameterName :
	process+=insert_ME("MEqq2W2ff","Electron")
	elif "Run-II-Z-e" in parameterName or "Run-I-Z" in parameterName :
	process +=insert_ME("MEqq2gZ2ff","Electron")
	elif "Run-II-Z-LowMass-mu" in parameterName :
	process +=insert_ME("MEqq2gZ2ff","Muon")
	process+=addLeptonPairCut("25","70")
	elif "Run-II-Z-HighMass-mu" in parameterName :
	process +=insert_ME("MEqq2gZ2ff","Muon")
	process+=addLeptonPairCut("150","600")
	elif "Run-II-Z-mu" in parameterName :
	process +=insert_ME("MEqq2gZ2ff","Muon")
	elif(simulation=="Powheg") :
	if "Run-I-WZ" in parameterName :
	process+=insert_ME("PowhegMEqq2W2ff","Electron")
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	elif "Run-II-W" in parameterName or "Run-I-W" in parameterName :
	process+=insert_ME("PowhegMEqq2W2ff","Electron")
	elif "Run-II-Z-e" in parameterName or "Run-I-Z" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	elif "Run-II-Z-LowMass-mu" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
	process+=addLeptonPairCut("25","70")
	elif "Run-II-Z-HighMass-mu" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
	process+=addLeptonPairCut("150","600")
	elif "Run-II-Z-mu" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
	elif "DiPhoton-GammaGamma" in parameterName :
	process+=insert_ME("MEGammaGammaPowheg","GammaGamma")
	process+=insert_ME("MEGammaGamma","gg")
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
	process+=jet_kt_cut(5.)
	parameterName=parameterName.replace("-GammaGamma","")
	elif "DiPhoton-GammaJet" in parameterName :
	process+=insert_ME("MEGammaGammaPowheg","VJet")
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
	process+=jet_kt_cut(5.)
	parameterName=parameterName.replace("-GammaJet","")
	elif(simulation=="Matchbox" or simulation=="Merging" ) :
	if "Jets" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p -> j j","2","0","MaxJetPtScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p -> j j","2","0","MaxJetPtScale",1,0)
	process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
	if "Run-II-Jets-10" in parameterName :
	process+=addFirstJet("30")
	process+=addSecondJet("25")
	process+=addJetPairCut("500")
	elif "Run-II-Jets-11" in parameterName :
	process+=addFirstJet("30")
	process+=addSecondJet("25")
	process+=addJetPairCut("900")
	elif "Run-II-Jets-12" in parameterName :
	process+=addFirstJet("30")
	process+=addSecondJet("25")
	process+=addJetPairCut("300")
	elif "Run-I-Jets-1" in parameterName : process+=addFirstJet("20")
	elif "Run-I-Jets-2" in parameterName : process+=addFirstJet("40")
	elif "Run-I-Jets-3" in parameterName : process+=addFirstJet("65")
	elif "Run-I-Jets-4" in parameterName : process+=addFirstJet("90")
	elif "Run-I-Jets-5" in parameterName : process+=addFirstJet("160")
	elif "Run-I-Jets-6" in parameterName : process+=addFirstJet("30")+addSecondJet("25")+addJetPairCut("100")
	elif "Run-I-Jets-7" in parameterName : process+=addFirstJet("30")+addSecondJet("25")+addJetPairCut("400")
	elif "Run-I-Jets-8" in parameterName : process+=addFirstJet("30")+addSecondJet("25")+addJetPairCut("700")
	elif "Run-II-Jets-0" in parameterName : process+=addFirstJet("15")
	elif "Run-II-Jets-1" in parameterName : process+=addFirstJet("25")
	elif "Run-II-Jets-2" in parameterName : process+=addFirstJet("40")
	elif "Run-II-Jets-3" in parameterName : process+=addFirstJet("60")
	elif "Run-II-Jets-4" in parameterName : process+=addFirstJet("85")
	elif "Run-II-Jets-5" in parameterName : process+=addFirstJet("110")
	elif "Run-II-Jets-6" in parameterName : process+=addFirstJet("160")
	elif "Run-II-Jets-7" in parameterName : process+=addFirstJet("250")
	elif "Run-II-Jets-8" in parameterName : process+=addFirstJet("30")+addSecondJet("25")+addJetPairCut("100")
	elif "Run-II-Jets-9" in parameterName : process+=addFirstJet("30")+addSecondJet("25")+addJetPairCut("300")
	elif "900-Jets-1" in parameterName : process+=addFirstJet("10")
	elif "300-Jets-1" in parameterName : process+=addFirstJet("6")
	elif "630-Jets-1" in parameterName : process+=addFirstJet("20")
	elif "630-Jets-2" in parameterName : process+=addFirstJet("40")
	elif "630-Jets-3" in parameterName : process+=addFirstJet("75")
	elif "900-Jets-1" in parameterName : process+=addFirstJet("10")
	else :
	logging.error("Exit 00007")
	sys.exit(1)
	elif "Run-I-WZ" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p pbar e+ e-","0","2","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p pbar e+ nu","0","2","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p pbar e- nu","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=particlegroup('epm','e+','e-')
	process+=particlegroup('epmnu','e+','e-','nu_e','nu_ebar')
	process+=addProcess(thefactory,"p pbar epm epmnu","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("60","120")
	elif "Run-II-W" in parameterName or "Run-I-W" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p pbar e+ nu","0","2","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p pbar e- nu","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=particlegroup('epm','e+','e-')
	process+=addProcess(thefactory,"p pbar epm nu","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("60","120")
	elif "Run-II-Z-e" in parameterName or "Run-I-Z" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p pbar e+ e-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p pbar e+ e-","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("60","120")
	elif "Run-II-Z-LowMass-mu" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("25","70")
	elif "Run-II-Z-HighMass-mu" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("150","600")
	elif "Run-II-Z-mu" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("60","120")
	# Star
	elif(collider=="Star" ) :
	process = StringBuilder("set /Herwig/Decays/DecayHandler:LifeTimeOption 0\n")
	process+= "set /Herwig/Decays/DecayHandler:MaxLifeTime 10*mm\n"
	process+= "set /Herwig/Generators/EventGenerator:EventHandler:BeamB /Herwig/Particles/p+\n"
	process+= collider_lumi(200.0)
	process+= "set /Herwig/Cuts/Cuts:X2Min 0.01\n"
	if(simulation=="") :
	if "UE" in parameterName :
	if "Dipole" in parameters["shower"]:
	process+="read snippets/MB-DipoleShower.in\n"
	else:
	process+="read snippets/MB.in\n"
	process+="read snippets/Diffraction.in\n"


	else :
	process+=insert_ME("MEQCD2to2")
	process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
	if "Jets-1" in parameterName : process+=jet_kt_cut(2.)
	elif "Jets-2" in parameterName : process+=jet_kt_cut(5.)
	elif "Jets-3" in parameterName : process+=jet_kt_cut(20.)
	elif "Jets-4" in parameterName : process+=jet_kt_cut(25.)
	else :
	logging.error("Star not supported for %s " % simulation)
	sys.exit(1)
	# ISR and SppS
	elif(collider=="ISR" or collider =="SppS" or collider == "EHS" ) :
	process = StringBuilder("set /Herwig/Decays/DecayHandler:LifeTimeOption 0\n")
	process+="set /Herwig/Decays/DecayHandler:MaxLifeTime 10*mm\n"
	if(collider=="SppS") :
	process = StringBuilder("set /Herwig/Generators/EventGenerator:EventHandler:BeamB /Herwig/Particles/pbar-\n")
	if "30" in parameterName : process+=collider_lumi( 30.4)
	elif "44" in parameterName : process+=collider_lumi( 44.4)
	elif "53" in parameterName : process+=collider_lumi( 53.0)
	elif "62" in parameterName : process+=collider_lumi( 62.2)
	elif "63" in parameterName : process+=collider_lumi( 63.0)
	elif "200" in parameterName : process+=collider_lumi(200.0)
	elif "500" in parameterName : process+=collider_lumi(500.0)
	elif "546" in parameterName : process+=collider_lumi(546.0)
	elif "900" in parameterName : process+=collider_lumi(900.0)
	if(simulation=="") :
	if "Dipole" in parameters["shower"]:
	process+="read snippets/MB-DipoleShower.in\n"
	else:
	process+="read snippets/MB.in\n"
	process+="read snippets/Diffraction.in\n"
	else :
	logging.error(" SppS and ISR not supported for %s " % simulation)
	sys.exit(1)
	# LHC
	elif(collider=="LHC") :
	if parameterName.startswith("7-") : process = StringBuilder(collider_lumi(7000.0))
	elif parameterName.startswith("8-") : process = StringBuilder(collider_lumi(8000.0))
	elif parameterName.startswith("13-") : process = StringBuilder(collider_lumi(13000.0))
	elif parameterName.startswith("900") : process = StringBuilder(collider_lumi(900.0))
	elif parameterName.startswith("2360") : process = StringBuilder(collider_lumi(2360.0))
	elif parameterName.startswith("2760") : process = StringBuilder(collider_lumi(2760.0))
	else : process = StringBuilder(collider_lumi(7000.0))
	if(simulation=="") :
	if "VBF" in parameterName :
	process+=insert_ME("MEPP2HiggsVBF")
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	addedBRReweighter = True
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	addedBRReweighter = True
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	addedBRReweighter = True
	elif "8-" not in parameterName :
	process+=selectDecayMode("h0",["h0->tau-,tau+;"])
	addedBRReweighter = True
	process+="set /Herwig/Particles/tau-:Stable Stable\n"

	elif "ggHJet" in parameterName :
	process+=selectDecayMode("h0",["h0->tau-,tau+;"])
	addedBRReweighter = True
	process+="set /Herwig/Particles/tau-:Stable Stable\n"
	process+=insert_ME("MEHiggsJet")
	process+=jet_kt_cut(20.)
	elif "ggH" in parameterName :
	process+=insert_ME("MEHiggs")
	process+=insert_ME("MEHiggsJet","qqbar")
	process+=jet_kt_cut(0.0)
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	addedBRReweighter = True
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	addedBRReweighter = True
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	addedBRReweighter = True
	elif "8-" not in parameterName :
	process+=selectDecayMode("h0",["h0->tau-,tau+;"])
	addedBRReweighter = True
	process+="set /Herwig/Particles/tau-:Stable Stable\n"

	elif "PromptPhoton" in parameterName :
	process+=insert_ME("MEGammaJet")
	if "PromptPhoton-1" in parameterName :
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
	elif "PromptPhoton-2" in parameterName :
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 25.\n"
	elif "PromptPhoton-3" in parameterName :
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 80.\n"
	elif "PromptPhoton-4" in parameterName :
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 150.\n"
	elif "DiPhoton-GammaGamma" in parameterName :
	process+=insert_ME("MEGammaGamma")
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
	parameterName=parameterName.replace("-GammaGamma","")
	elif "DiPhoton-GammaJet" in parameterName :
	process+=insert_ME("MEGammaJet")
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
	parameterName=parameterName.replace("-GammaJet","")
	elif "8-WH" in parameterName :
	process+=insert_ME("MEPP2WH")
	process+=jet_kt_cut(0.0)
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	addedBRReweighter = True
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	addedBRReweighter = True
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	addedBRReweighter = True
	elif "8-ZH" in parameterName :
	process+=insert_ME("MEPP2ZH")
	process+=jet_kt_cut(0.0)
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	addedBRReweighter = True
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	addedBRReweighter = True
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	addedBRReweighter = True
	elif "WH" in parameterName :
	process+=selectDecayMode("h0",["h0->b,bbar;"])
	process+=selectDecayMode("W+",["W+->nu_e,e+;",
	"W+->nu_mu,mu+;"])
	addedBRReweighter = True
	process+=insert_ME("MEPP2WH")
	process+=jet_kt_cut(0.0)
	elif "ZH" in parameterName :
	process+=selectDecayMode("h0",["h0->b,bbar;"])
	process+=selectDecayMode("Z0",["Z0->e-,e+;",
	"Z0->mu-,mu+;"])
	addedBRReweighter = True
	process+=insert_ME("MEPP2ZH")
	process+=jet_kt_cut(0.0)
	elif "UE" in parameterName :
	if "Dipole" in parameters["shower"]:
	process+="read snippets/MB-DipoleShower.in\n"
	else:
	process+="set /Herwig/Shower/ShowerHandler:IntrinsicPtGaussian 2.2*GeV\n"
	process+="read snippets/MB.in\n"
	process+="read snippets/Diffraction.in\n"
	if "Long" in parameterName :
	process += "set /Herwig/Decays/DecayHandler:MaxLifeTime 100*mm\n"
	elif "8-DiJets" in parameterName or "7-DiJets" in parameterName or "13-DiJets" in parameterName :
	process+=insert_ME("MEQCD2to2")
	process+="set MEQCD2to2:MaximumFlavour 5\n"
	process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
	if "13-DiJets" not in parameterName :
	if "-A" in parameterName :
	process+=jet_kt_cut(45.)
	process+="set /Herwig/Cuts/JetKtCut:MinEta -3.\n"
	process+="set /Herwig/Cuts/JetKtCut:MaxEta 3.\n"
	elif "-B" in parameterName :
	process+=jet_kt_cut(20.)
	process+="set /Herwig/Cuts/JetKtCut:MinEta -2.7\n"
	process+="set /Herwig/Cuts/JetKtCut:MaxEta 2.7\n"
	elif "-C" in parameterName :
	process+=jet_kt_cut(20.)
	process+="set /Herwig/Cuts/JetKtCut:MinEta -4.8\n"
	process+="set /Herwig/Cuts/JetKtCut:MaxEta 4.8\n"
	else :
	if "-A" in parameterName :
	process+=jet_kt_cut(180.)
	process+="set /Herwig/Cuts/JetKtCut:MinEta -3.\n"
	process+="set /Herwig/Cuts/JetKtCut:MaxEta 3.\n"

	if "DiJets-1" in parameterName : process+=mhatmin_cut(90.)
	elif "DiJets-2" in parameterName : process+=mhatmin_cut(200.)
	elif "DiJets-3" in parameterName : process+=mhatmin_cut(450.)
	elif "DiJets-4" in parameterName : process+=mhatmin_cut(750.)
	elif "DiJets-5" in parameterName : process+=mhatmin_cut(950.)
	elif "DiJets-6" in parameterName : process+=mhatmin_cut(1550.)
	elif "DiJets-7" in parameterName : process+=mhatmin_cut(2150.)
	elif "DiJets-8" in parameterName : process+=mhatmin_cut(2750.)
	elif "DiJets-9" in parameterName : process+=mhatmin_cut(3750.)
	elif "DiJets-10" in parameterName : process+=mhatmin_cut(4750.)
	elif "DiJets-11" in parameterName : process+=mhatmin_cut(5750.)
	elif( "7-Jets" in parameterName
	or "8-Jets" in parameterName
	or "13-Jets" in parameterName
	or "2760-Jets" in parameterName
	) :
	process+=insert_ME("MEQCD2to2")
	process+="set MEQCD2to2:MaximumFlavour 5\n"
	process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
	if "Jets-10" in parameterName : process+=jet_kt_cut(1800.)
	elif "Jets-0" in parameterName : process+=jet_kt_cut(5.)
	elif "Jets-1" in parameterName : process+=jet_kt_cut(10.)
	elif "Jets-2" in parameterName : process+=jet_kt_cut(20.)
	elif "Jets-3" in parameterName : process+=jet_kt_cut(40.)
	elif "Jets-4" in parameterName : process+=jet_kt_cut(70.)
	elif "Jets-5" in parameterName : process+=jet_kt_cut(150.)
	elif "Jets-6" in parameterName : process+=jet_kt_cut(200.)
	elif "Jets-7" in parameterName : process+=jet_kt_cut(300.)
	elif "Jets-8" in parameterName : process+=jet_kt_cut(500.)
	elif "Jets-9" in parameterName : process+=jet_kt_cut(800.)
	elif( "7-Charm" in parameterName or "7-Bottom" in parameterName
	or "8-Bottom" in parameterName) :

	if("8-Bottom" in parameterName) :
	addBRReweighter()
	process+=selectDecayMode("Jpsi",["Jpsi->mu-,mu+;"])

	if "Bottom" in parameterName :
	process+="cp MEHeavyQuark MEBottom\n"
	process+="set MEBottom:QuarkType Bottom\n"
	process+=insert_ME("MEBottom")
	else :
	process+="cp MEHeavyQuark MECharm\n"
	process+="set MECharm:QuarkType Charm\n"
	process+=insert_ME("MECharm")
	process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
	if "-0" in parameterName :
	if "Bottom" in parameterName :
	process+="set MEBottom:Process Pair\n"
	process+=jet_kt_cut(0.)
	elif "-1" in parameterName : process+=jet_kt_cut(5.)
	elif "-2" in parameterName : process+=jet_kt_cut(15.)
	elif "-3" in parameterName : process+=jet_kt_cut(20.)
	elif "-4" in parameterName : process+=jet_kt_cut(50.)
	elif "-5" in parameterName : process+=jet_kt_cut(80.)
	elif "-6" in parameterName : process+=jet_kt_cut(110.)
	elif "-7" in parameterName : process+=jet_kt_cut(30.)+mhatmin_cut(90.)
	elif "-8" in parameterName : process+=jet_kt_cut(30.)+mhatmin_cut(340.)
	elif "-9" in parameterName : process+=jet_kt_cut(30.)+mhatmin_cut(500.)
	elif "Top-L" in parameterName :
	process+="set MEHeavyQuark:QuarkType Top\n"
	process+=insert_ME("MEHeavyQuark")
	process+=selectDecayMode("t",["t->nu_e,e+,b;",
	"t->nu_mu,mu+,b;"])
	process+=addBRReweighter()

	elif "Top-SL" in parameterName :
	process+="set MEHeavyQuark:QuarkType Top\n"
	process+=insert_ME("MEHeavyQuark")
	process+="set /Herwig/Particles/t:Synchronized Not_synchronized\n"
	process+="set /Herwig/Particles/tbar:Synchronized Not_synchronized\n"
	process+=selectDecayMode("t",["t->nu_e,e+,b;","t->nu_mu,mu+,b;"])
	process+=selectDecayMode("tbar",["tbar->b,bbar,cbar;",
	"tbar->bbar,cbar,d;",
	"tbar->bbar,cbar,s;",
	"tbar->bbar,s,ubar;",
	"tbar->bbar,ubar,d;"])
	process+=addBRReweighter()

	elif "Top-All" in parameterName :
	process+="set MEHeavyQuark:QuarkType Top\n"
	process+=insert_ME("MEHeavyQuark")
	elif "WZ" in parameterName :
	process+=insert_ME("MEPP2VV","WZ")
	process+=selectDecayMode("W+",["W+->nu_e,e+;",
	"W+->nu_mu,mu+;"])
	process+=selectDecayMode("W-",["W-->nu_ebar,e-;",
	"W-->nu_mubar,mu-;"])
	process+=selectDecayMode("Z0",["Z0->e-,e+;",
	"Z0->mu-,mu+;"])
	addedBRReweighter = True

	elif "WW-emu" in parameterName :
	process+=insert_ME("MEPP2VV","WW")
	process+="set /Herwig/Particles/W+:Synchronized 0\n"
	process+="set /Herwig/Particles/W-:Synchronized 0\n"
	process+=selectDecayMode("W+",["W+->nu_e,e+;"])
	process+=selectDecayMode("W-",["W-->nu_mubar,mu-;"])
	addedBRReweighter = True

	elif "WW-ll" in parameterName :
	process+=insert_ME("MEPP2VV","WW")
	process+=selectDecayMode("W+",["W+->nu_e,e+;","W+->nu_mu,mu+;","W+->nu_tau,tau+;"])
	addedBRReweighter = True

	elif "ZZ-ll" in parameterName :
	process+=insert_ME("MEPP2VV","ZZ")
	process+=selectDecayMode("Z0",["Z0->e-,e+;",
	"Z0->mu-,mu+;",
	"Z0->tau-,tau+;"])
	addedBRReweighter = True

	elif "ZZ-lv" in parameterName :
	process+=insert_ME("MEPP2VV","ZZ")
	process+=selectDecayMode("Z0",["Z0->e-,e+;",
	"Z0->mu-,mu+;",
	"Z0->tau-,tau+;",
	"Z0->nu_e,nu_ebar;",
	"Z0->nu_mu,nu_mubar;",
	"Z0->nu_tau,nu_taubar;"])
	addedBRReweighter = True

	elif "W-Z-e" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Electron")
	process+=insert_ME("MEqq2W2ff","Electron")
	elif "W-Z-mu" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Muon")
	process+=insert_ME("MEqq2W2ff","Muon")
	elif "W-e" in parameterName :
	process+=insert_ME("MEqq2W2ff","Electron")
	elif "W-mu" in parameterName :
	process+=insert_ME("MEqq2W2ff","Muon")
	elif "Z-e" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Electron")
	elif "Z-mu" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Muon")
	elif "Z-LowMass-e" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Electron")
	process+=mhat_minm_maxm(20,20,70)
	elif "Z-MedMass-e" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Electron")
	process+=mhat_minm_maxm(40,40,130)
	elif "Z-LowMass-mu" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Muon")
	process+=mhat_minm_maxm(10,10,70)
	elif "Z-Mass1" in parameterName :
	process+=mhat_minm_maxm(10,10,35)
	if "-e" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("MEqq2gZ2ff","Muon")
	elif "Z-Mass2" in parameterName :
	process+=mhat_minm_maxm(25,25,70)
	if "-e" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("MEqq2gZ2ff","Muon")
	elif "Z-Mass3" in parameterName :
	process+=mhat_minm_maxm(60,60,120)
	if "-e" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("MEqq2gZ2ff","Muon")
	elif "Z-Mass4" in parameterName :
	process+=mhat_minm_maxm(110,110,8000)
	if "-e" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("MEqq2gZ2ff","Muon")
	elif "Z-HighMass1" in parameterName :
	process+=mhat_minm_maxm(116,116,400)
	if "-e" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("MEqq2gZ2ff","Muon")
	elif "Z-HighMass2" in parameterName :
	process+=mhat_minm_maxm(400,400,7000)
	if "-e" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("MEqq2gZ2ff","Muon")
	elif "W-Jet" in parameterName :
	process+=insert_ME("MEWJet","Electron","WDecay")
	if "W-Jet-1-e" in parameterName :
	process+="set /Herwig/Cuts/WBosonKtCut:MinKT 100.0*GeV\n"
	parameterName=parameterName.replace("W-Jet-1-e","W-Jet-e")
	elif "W-Jet-2-e" in parameterName :
	process+="set /Herwig/Cuts/WBosonKtCut:MinKT 190.0*GeV\n"
	parameterName=parameterName.replace("W-Jet-2-e","W-Jet-e")
	elif "W-Jet-3-e" in parameterName :
	process+="set /Herwig/Cuts/WBosonKtCut:MinKT 270.0*GeV\n"
	parameterName=parameterName.replace("W-Jet-3-e","W-Jet-e")
	elif "Z-Jet" in parameterName :
	if "-e" in parameterName :
	process+=insert_ME("MEZJet","Electron","ZDecay")
	if "Z-Jet-0-e" in parameterName :
	process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 35.0*GeV\n"
	parameterName=parameterName.replace("Z-Jet-0-e","Z-Jet-e")
	elif "Z-Jet-1-e" in parameterName :
	process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 100.0*GeV\n"
	parameterName=parameterName.replace("Z-Jet-1-e","Z-Jet-e")
	elif "Z-Jet-2-e" in parameterName :
	process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 190.0*GeV\n"
	parameterName=parameterName.replace("Z-Jet-2-e","Z-Jet-e")
	elif "Z-Jet-3-e" in parameterName :
	process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 270.0*GeV\n"
	parameterName=parameterName.replace("Z-Jet-3-e","Z-Jet-e")
	else :
	process+=insert_ME("MEZJet","Muon","ZDecay")
	process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 35.0*GeV\n"
	parameterName=parameterName.replace("Z-Jet-0-mu","Z-Jet-mu")
	elif "WGamma" in parameterName :
	process+=insert_ME("MEPP2VGamma","1")
	process+="set MEPP2VGamma:MassOption 1"
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 10.\n"


	if "-e" in parameterName :
	process+=selectDecayMode("W+",["W+->nu_e,e+;"])
	addedBRReweighter=True
	else :
	process+=selectDecayMode("W+",["W+->nu_mu,mu+;"])
	addedBRReweighter=True
	elif "ZGamma" in parameterName :
	process+=insert_ME("MEPP2VGamma","2")
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 10.\n"
	if "-e" in parameterName :
	process+=selectDecayMode("Z0",["Z0->e-,e+;"])
	addedBRReweighter=True
	else :
	process+=selectDecayMode("Z0",["Z0->mu-,mu+;"])
	addedBRReweighter=True
	else :
	logging.error(" Process %s not supported for internal matrix elements" % name)
	sys.exit(1)
	elif(simulation=="Powheg") :
	if "VBF" in parameterName :
	process+=insert_ME("PowhegMEPP2HiggsVBF")
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	addedBRReweighter = True
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	addedBRReweighter = True
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	addedBRReweighter = True
	elif "8-" not in parameterName :
	process+=selectDecayMode("h0",["h0->tau-,tau+;"])
	addedBRReweighter = True
	process+="set /Herwig/Particles/tau-:Stable Stable\n"

	elif "ggHJet" in parameterName :
	logging.error(" Process %s not supported for POWHEG matrix elements" % name)
	sys.exit(1)
	elif "ggH" in parameterName :
	process+=insert_ME("PowhegMEHiggs")
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	addedBRReweighter = True
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	addedBRReweighter = True
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	addedBRReweighter = True
	elif "8-" not in parameterName :
	process+=selectDecayMode("h0",["h0->tau-,tau+;"])
	addedBRReweighter = True
	process+="set /Herwig/Particles/tau-:Stable Stable\n"
	elif "8-WH" in parameterName :
	process+=insert_ME("PowhegMEPP2WH")
	process+=jet_kt_cut(0.0)
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	addedBRReweighter = True
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	addedBRReweighter = True
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	addedBRReweighter = True
	elif "8-ZH" in parameterName :
	process+=insert_ME("PowhegMEPP2ZH")
	process+=jet_kt_cut(0.0)
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	addedBRReweighter = True
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	addedBRReweighter = True
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	addedBRReweighter = True
	elif "WH" in parameterName :
	process+=selectDecayMode("h0",["h0->b,bbar;"])
	process+=selectDecayMode("W+",["W+->nu_e,e+;",
	"W+->nu_mu,mu+;"])
	addedBRReweighter = True
	process+=insert_ME("PowhegMEPP2WH")
	process+=jet_kt_cut(0.0)
	elif "ZH" in parameterName :
	process+=selectDecayMode("h0",["h0->b,bbar;"])
	process+=selectDecayMode("Z0",["Z0->e-,e+;",
	"Z0->mu-,mu+;"])
	addedBRReweighter = True
	process+=insert_ME("PowhegMEPP2ZH")
	process+=jet_kt_cut(0.0)
	elif "UE" in parameterName :
	logging.error(" Process %s not supported for powheg matrix elements" % name)
	sys.exit(1)
	elif "WZ" in parameterName :
	process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
	process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
	process+="set /Herwig/Shower/ShowerHandler:SplitHardProcess No\n";
	process+="set /Herwig/Decays/ZDecayer:PhotonGenerator NULL\n";
	process+="set /Herwig/Decays/WDecayer:PhotonGenerator NULL\n";
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
	process+=insert_ME("PowhegMEPP2VV","WZ")
	process+=selectDecayMode("W+",["W+->nu_e,e+;",
	"W+->nu_mu,mu+;"])
	process+=selectDecayMode("W-",["W-->nu_ebar,e-;",
	"W-->nu_mubar,mu-;"])
	process+=selectDecayMode("Z0",["Z0->e-,e+;",
	"Z0->mu-,mu+;"])
	addedBRReweighter = True

	elif "WW-emu" in parameterName :
	process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
	process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
	process+="set /Herwig/Shower/ShowerHandler:SplitHardProcess No\n";
	process+="set /Herwig/Decays/ZDecayer:PhotonGenerator NULL\n";
	process+="set /Herwig/Decays/WDecayer:PhotonGenerator NULL\n";
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
	process+=insert_ME("PowhegMEPP2VV","WW")
	process+="set /Herwig/Particles/W+:Synchronized 0\n"
	process+="set /Herwig/Particles/W-:Synchronized 0\n"
	process+=selectDecayMode("W+",["W+->nu_e,e+;"])
	process+=selectDecayMode("W-",["W-->nu_mubar,mu-;"])
	addedBRReweighter = True

	elif "WW-ll" in parameterName :
	process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
	process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
	process+="set /Herwig/Shower/ShowerHandler:SplitHardProcess No\n";
	process+="set /Herwig/Decays/ZDecayer:PhotonGenerator NULL\n";
	process+="set /Herwig/Decays/WDecayer:PhotonGenerator NULL\n";
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
	process+=insert_ME("PowhegMEPP2VV","WW")
	process+=selectDecayMode("W+",["W+->nu_e,e+;",
	"W+->nu_mu,mu+;",
	"W+->nu_tau,tau+;"])
	addedBRReweighter = True

	elif "ZZ-ll" in parameterName :
	process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
	process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
	process+="set /Herwig/Shower/ShowerHandler:SplitHardProcess No\n";
	process+="set /Herwig/Decays/ZDecayer:PhotonGenerator NULL\n";
	process+="set /Herwig/Decays/WDecayer:PhotonGenerator NULL\n";
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
	process+=insert_ME("PowhegMEPP2VV","ZZ")
	process+=selectDecayMode("Z0",["Z0->e-,e+;",
	"Z0->mu-,mu+;",
	"Z0->tau-,tau+;"])
	addedBRReweighter = True

	elif "ZZ-lv" in parameterName :
	process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
	process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
	process+="set /Herwig/Shower/ShowerHandler:SplitHardProcess No\n";
	process+="set /Herwig/Decays/ZDecayer:PhotonGenerator NULL\n";
	process+="set /Herwig/Decays/WDecayer:PhotonGenerator NULL\n";
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
	process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
	process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
	process+=insert_ME("PowhegMEPP2VV","ZZ")
	process+=selectDecayMode("Z0",["Z0->e-,e+;",
	"Z0->mu-,mu+;",
	"Z0->tau-,tau+;",
	"Z0->nu_e,nu_ebar;",
	"Z0->nu_mu,nu_mubar;",
	"Z0->nu_tau,nu_taubar;"])
	addedBRReweighter = True

	elif "W-Z-e" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	process+=insert_ME("PowhegMEqq2W2ff","Electron")
	elif "W-Z-mu" in parameterName :
	process+=insert_ME("MEqq2gZ2ff","Muon")
	process+=insert_ME("MEqq2W2ff","Muon")
	elif "W-e" in parameterName :
	process+=insert_ME("PowhegMEqq2W2ff","Electron")
	elif "W-mu" in parameterName :
	process+=insert_ME("PowhegMEqq2W2ff","Muon")
	elif "Z-e" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	elif "Z-mu" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
	elif "Z-LowMass-e" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	process+=mhat_minm_maxm(20,20,70)
	elif "Z-MedMass-e" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	process+=mhat_minm_maxm(40,40,130)
	elif "Z-LowMass-mu" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
	process+=mhat_minm_maxm(10,10,70)
	elif "Z-Mass1" in parameterName :
	process+=mhat_minm_maxm(10,10,35)
	if "-e" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
	elif "Z-Mass2" in parameterName :
	process+=mhat_minm_maxm(25,25,70)
	if "-e" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
	elif "Z-Mass3" in parameterName :
	process+=mhat_minm_maxm(60,60,120)
	if "-e" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
	elif "Z-Mass4" in parameterName :
	process+=mhat_minm_maxm(110,110,8000)
	if "-e" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
	elif "Z-HighMass1" in parameterName :
	process+=mhat_minm_maxm(116,116,400)
	if "-e" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
	elif "Z-HighMass2" in parameterName :
	process+=mhat_minm_maxm(400,400,7000)
	if "-e" in parameterName :
	process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
	else :
	process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
	elif "DiPhoton-GammaGamma" in parameterName :
	process+=insert_ME("MEGammaGammaPowheg","GammaGamma")
	process+=insert_ME("MEGammaGamma","gg")
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
	process+=jet_kt_cut(5.)
	parameterName=parameterName.replace("-GammaGamma","")
	elif "DiPhoton-GammaJet" in parameterName :
	process+=insert_ME("MEGammaGammaPowheg","VJet")
	process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
	process+=jet_kt_cut(5.)
	parameterName=parameterName.replace("-GammaJet","")
	else :
	logging.error(" Process %s not supported for internal POWHEG matrix elements" % name)
	sys.exit(1)

	elif( simulation=="Matchbox" or simulation=="Merging" ) :
	if "VBF" in parameterName :
	parameters["nlo"] = "read Matchbox/VBFNLO.in\n"
	if(simulation=="Merging"):
	process+="cd /Herwig/Merging/\n"
	process+="insert "+thefactory+":DiagramGenerator:RestrictLines 0 /Herwig/Particles/Z0\n"
	process+="insert "+thefactory+":DiagramGenerator:RestrictLines 0 /Herwig/Particles/W+\n"
	process+="insert "+thefactory+":DiagramGenerator:RestrictLines 0 /Herwig/Particles/W-\n"
	process+="insert "+thefactory+":DiagramGenerator:RestrictLines 0 /Herwig/Particles/gamma\n"
	process+="do "+thefactory+":DiagramGenerator:TimeLikeRange 0 0\n"
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p h0 j j","0","3","FixedScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p h0 j j","0","3","FixedScale",1,1)
	process+=setHardProcessWidthToZero(["h0"])
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 125.7\n"
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	process+=addBRReweighter()
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	process+=addBRReweighter()
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	process+=addBRReweighter()
	elif "8-" not in parameterName :
	process+=selectDecayMode("h0",["h0->tau-,tau+;"])
	process+=addBRReweighter()
	process+="set /Herwig/Particles/tau-:Stable Stable\n"
	elif "ggHJet" in parameterName :
	if(simulation=="Merging"):
	logging.warning("ggHJet not explicitly tested for %s " % simulation)
	sys.exit(0)
	parameters["nlo"] = "read Matchbox/MadGraph-GoSam.in\nread Matchbox/HiggsEffective.in\n"
	process+=selectDecayMode("h0",["h0->tau-,tau+;"])
	process+=addBRReweighter()
	process+="set /Herwig/Particles/tau-:Stable Stable\n"
	process+=setHardProcessWidthToZero(["h0"])
	process+=addProcess(thefactory,"p p h0 j","3","1","FixedScale",0,0)
	process+=addFirstJet("20")
	process+="set "+thefactory+":ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 125.7\n"
	elif "ggH" in parameterName :
	parameters["nlo"] = "read Matchbox/MadGraph-GoSam.in\nread Matchbox/HiggsEffective.in\n"
	if(simulation=="Merging"):
	process+= "cd /Herwig/MatrixElements/Matchbox/Amplitudes\nset OpenLoops:HiggsEff On\nset MadGraph:Model heft\n"
	process+="cd /Herwig/Merging/\n"
	process+=setHardProcessWidthToZero(["h0"])
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p h0","2","1","FixedScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p h0","2","1","FixedScale",2,2)
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 125.7\n"
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	process+=addBRReweighter()
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	process+=addBRReweighter()
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	process+=addBRReweighter()
	elif "8-" not in parameterName :
	process+=selectDecayMode("h0",["h0->tau-,tau+;"])
	process+=addBRReweighter()
	process+="set /Herwig/Particles/tau-:Stable Stable\n"
	elif "8-WH" in parameterName :
	if(simulation=="Merging"):
	logging.warning("8-WH not explicitly tested for %s " % simulation)
	sys.exit(0)
	process+=setHardProcessWidthToZero(["h0","W+","W-"])
	process+=addProcess(thefactory,"p p W+ h0","0","2","FixedScale",0,0)
	process+=addProcess(thefactory,"p p W- h0","0","2","FixedScale",0,0)
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 125.7\n"
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	process+=addBRReweighter()
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	process+=addBRReweighter()
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	process+=addBRReweighter()

	elif "8-ZH" in parameterName :
	if(simulation=="Merging"):
	logging.warning("8-ZH not explicitly tested for %s " % simulation)
	sys.exit(0)
	process+=setHardProcessWidthToZero(["h0","Z0"])
	process+=addProcess(thefactory,"p p Z0 h0","0","2","FixedScale",0,0)
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 125.7\n"
	if "GammaGamma" in parameterName :
	process+=selectDecayMode("h0",["h0->gamma,gamma;"])
	process+=addBRReweighter()
	elif "WW" in parameterName :
	process+=selectDecayMode("h0",["h0->W+,W-;"])
	process+=addBRReweighter()
	elif "ZZ" in parameterName :
	process+=selectDecayMode("h0",["h0->Z0,Z0;"])
	process+=addBRReweighter()

	elif "WH" in parameterName :
	if(simulation=="Merging"):
	logging.warning("WH not explicitly tested for %s " % simulation)
	sys.exit(0)
	process+=selectDecayMode("h0",["h0->b,bbar;"])
	process+=addBRReweighter()
	process+=setHardProcessWidthToZero(["h0"])
	process+=addProcess(thefactory,"p p e+ nu h0","0","3","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p p e- nu h0","0","3","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p p mu+ nu h0","0","3","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p p mu- nu h0","0","3","LeptonPairMassScale",0,0)
	process+=addLeptonPairCut("60","120")
	elif "ZH" in parameterName :
	if(simulation=="Merging"):
	logging.warning("ZH not explicitly tested for %s " % simulation)
	sys.exit(0)
	process+=selectDecayMode("h0",["h0->b,bbar;"])
	process+=addBRReweighter()
	process+=setHardProcessWidthToZero(["h0"])
	process+=addProcess(thefactory,"p p e+ e- h0","0","3","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p p mu+ mu- h0","0","3","LeptonPairMassScale",0,0)
	process+=addLeptonPairCut("60","120")
	elif "UE" in parameterName :
	logging.error(" Process %s not supported for Matchbox matrix elements" % name)
	sys.exit(1)
	elif "8-DiJets" in parameterName or "7-DiJets" in parameterName or "13-DiJets" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p j j","2","0","MaxJetPtScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p j j","2","0","MaxJetPtScale",1,1)
	process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
	if "13-DiJets" not in parameterName :
	if "-A" in parameterName :
	process+=addFirstJet("45")
	process+=addSecondJet("25")
	process+="set /Herwig/Cuts/FirstJet:YRange -3. 3.\n"
	process+="set /Herwig/Cuts/SecondJet:YRange -3. 3.\n"
	elif "-B" in parameterName :
	process+=addFirstJet("20")
	process+=addSecondJet("15")
	process+="set /Herwig/Cuts/FirstJet:YRange -2.7 2.7\n"
	process+="set /Herwig/Cuts/SecondJet:YRange -2.7 2.7\n"
	elif "-C" in parameterName :
	process+=addFirstJet("20")
	process+=addSecondJet("15")
	process+="set /Herwig/Cuts/FirstJet:YRange -4.8 4.8\n"
	process+="set /Herwig/Cuts/SecondJet:YRange -4.8 4.8\n"
	else :
	logging.error("Exit 00001")
	sys.exit(1)
	else :
	if "-A" in parameterName :
	process+=addFirstJet("220")
	process+=addSecondJet("180.")
	process+="set /Herwig/Cuts/JetKtCut:MinEta -3.\n"
	process+="set /Herwig/Cuts/JetKtCut:MaxEta 3.\n"
	else :
	logging.error("Exit 00001")
	sys.exit(1)


	if "DiJets-1" in parameterName : process+=addJetPairCut("90")
	elif "DiJets-2" in parameterName : process+=addJetPairCut("200")
	elif "DiJets-3" in parameterName : process+=addJetPairCut("450")
	elif "DiJets-4" in parameterName : process+=addJetPairCut("750")
	elif "DiJets-5" in parameterName : process+=addJetPairCut("950")
	elif "DiJets-6" in parameterName : process+=addJetPairCut("1550")
	elif "DiJets-7" in parameterName : process+=addJetPairCut("2150")
	elif "DiJets-8" in parameterName : process+=addJetPairCut("2750")
	elif "DiJets-9" in parameterName : process+=mhatmin_cut(3750.)
	elif "DiJets-10" in parameterName : process+=mhatmin_cut(4750.)
	elif "DiJets-11" in parameterName : process+=mhatmin_cut(5750.)
	else :
	logging.error("Exit 00002")
	sys.exit(1)


	elif( "7-Jets" in parameterName
	or "8-Jets" in parameterName
	or "13-Jets" in parameterName
	or "2760-Jets" in parameterName
	) :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p j j","2","0","MaxJetPtScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p j j","2","0","MaxJetPtScale",1,1)
	process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
	if "Jets-10" in parameterName : process+=addFirstJet("1800")
	elif "Jets-0" in parameterName : process+=addFirstJet("5")
	elif "Jets-1" in parameterName : process+=addFirstJet("10")
	elif "Jets-2" in parameterName : process+=addFirstJet("20")
	elif "Jets-3" in parameterName : process+=addFirstJet("40")
	elif "Jets-4" in parameterName : process+=addFirstJet("70")
	elif "Jets-5" in parameterName : process+=addFirstJet("150")
	elif "Jets-6" in parameterName : process+=addFirstJet("200")
	elif "Jets-7" in parameterName : process+=addFirstJet("300")
	elif "Jets-8" in parameterName : process+=addFirstJet("500")
	elif "Jets-9" in parameterName : process+=addFirstJet("800")
	else :
	logging.error("Exit 00003")
	sys.exit(1)
	elif( "7-Charm" in parameterName or "7-Bottom" in parameterName
	or "8-Bottom" in parameterName) :
	parameters["bscheme"]=fourFlavour
	process+="set /Herwig/Particles/b:HardProcessMass 4.2*GeV\n"
	process+="set /Herwig/Particles/bbar:HardProcessMass 4.2*GeV\n"

	if("8-Bottom" in parameterName) :
	addBRReweighter()
	process+=selectDecayMode("Jpsi",["Jpsi->mu-,mu+;"])

	if "Bottom" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p b bbar","2","0","MaxJetPtScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p b bbar","2","0","MaxJetPtScale",1,0)
	else:
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p c cbar","2","0","MaxJetPtScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p c cbar","2","0","MaxJetPtScale",1,0)

	process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
	if "-0" in parameterName : process+=addFirstJet("0")
	elif "-1" in parameterName : process+=addFirstJet("5")
	elif "-2" in parameterName : process+=addFirstJet("15")
	elif "-3" in parameterName : process+=addFirstJet("20")
	elif "-4" in parameterName : process+=addFirstJet("50")
	elif "-5" in parameterName : process+=addFirstJet("80")
	elif "-6" in parameterName : process+=addFirstJet("110")
	elif "-7" in parameterName :
	process+=addFirstJet("30")
	process+=addSecondJet("25")
	process+=addJetPairCut("90")
	elif "-8" in parameterName :
	process+=addFirstJet("30")
	process+=addSecondJet("25")
	process+=addJetPairCut("340")
	elif "-9" in parameterName :
	process+=addFirstJet("30")
	process+=addSecondJet("25")
	process+=addJetPairCut("500")
	else :
	logging.error("Exit 00004")
	sys.exit(1)

	elif "Top-L" in parameterName :
	process+=setHardProcessWidthToZero(["t","tbar"])
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",2,2)
	process+=selectDecayMode("t",["t->nu_e,e+,b;",
	"t->nu_mu,mu+,b;"])
	process+=addBRReweighter()

	elif "Top-SL" in parameterName :
	process+=setHardProcessWidthToZero(["t","tbar"])
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",2,2)
	process+="set /Herwig/Particles/t:Synchronized Not_synchronized\n"
	process+="set /Herwig/Particles/tbar:Synchronized Not_synchronized\n"
	process+=selectDecayMode("t",["t->nu_e,e+,b;",
	"t->nu_mu,mu+,b;"])
	process+=selectDecayMode("tbar",["tbar->b,bbar,cbar;",
	"tbar->bbar,cbar,d;",
	"tbar->bbar,cbar,s;",
	"tbar->bbar,s,ubar;",
	"tbar->bbar,ubar,d;"])
	process+=addBRReweighter()

	elif "Top-All" in parameterName :
	process+=setHardProcessWidthToZero(["t","tbar"])
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",2,2)
	elif "WZ" in parameterName :
	if(simulation=="Merging"):
	logging.warning("WZ not explicitly tested for %s " % simulation)
	sys.exit(0)
	process+=setHardProcessWidthToZero(["W+","W-","Z0"])
	process+=addProcess(thefactory,"p p W+ Z0","0","2","FixedScale",0,0)
	process+=addProcess(thefactory,"p p W- Z0","0","2","FixedScale",0,0)
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 171.6*GeV\n\n"
	process+=selectDecayMode("W+",["W+->nu_e,e+;",
	"W+->nu_mu,mu+;"])
	process+=selectDecayMode("W-",["W-->nu_ebar,e-;",
	"W-->nu_mubar,mu-;"])
	process+=selectDecayMode("Z0",["Z0->e-,e+;",
	"Z0->mu-,mu+;"])
	process+=addBRReweighter()
	process+=addLeptonPairCut("60","120")
	elif "WW-emu" in parameterName :
	if(simulation=="Merging"):
	logging.warning("WW-emu not explicitly tested for %s " % simulation)
	sys.exit(0)

	process+=setHardProcessWidthToZero(["W+","W-","Z0"])
	process+=addProcess(thefactory,"p p W+ W-","0","2","FixedScale",0,0)
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 160.8*GeV\n"
	process+="set /Herwig/Particles/W+:Synchronized 0\n"
	process+="set /Herwig/Particles/W-:Synchronized 0\n"
	process+=selectDecayMode("W+",["W+->nu_e,e+;"])
	process+=selectDecayMode("W-",["W-->nu_mubar,mu-;"])
	process+=addBRReweighter()
	parameters["bscheme"] = "read Matchbox/FourFlavourScheme.in\n"

	process+=addLeptonPairCut("60","120")
	elif "WW-ll" in parameterName :
	if(simulation=="Merging"):
	logging.warning("WW-ll not explicitly tested for %s " % simulation)
	sys.exit(0)
	process+=setHardProcessWidthToZero(["W+","W-","Z0"])
	process+=addProcess(thefactory,"p p W+ W-","0","2","FixedScale",0,0)
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 160.8*GeV\n"
	process+=selectDecayMode("W+",["W+->nu_e,e+;",
	"W+->nu_mu,mu+;",
	"W+->nu_tau,tau+;"])
	process+=addBRReweighter()
	process+=addLeptonPairCut("60","120")
	parameters["bscheme"] = "read Matchbox/FourFlavourScheme.in\n"

	elif "ZZ-ll" in parameterName :
	if(simulation=="Merging"):
	logging.warning("ZZ-ll not explicitly tested for %s " % simulation)
	sys.exit(0)
	process+=setHardProcessWidthToZero(["W+","W-","Z0"])
	process+=addProcess(thefactory,"p p Z0 Z0","0","2","FixedScale",0,0)
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 182.2*GeV\n"
	process+=selectDecayMode("Z0",["Z0->e-,e+;",
	"Z0->mu-,mu+;",
	"Z0->tau-,tau+;"])
	process+=addBRReweighter()
	process+=addLeptonPairCut("60","120")
	elif "ZZ-lv" in parameterName :
	if(simulation=="Merging"):
	logging.warning("ZZ-lv not explicitly tested for %s " % simulation)
	sys.exit(0)
	process+=setHardProcessWidthToZero(["W+","W-","Z0"])
	process+=addProcess(thefactory,"p p Z0 Z0","0","2","FixedScale",0,0)
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 182.2*GeV\n"
	process+=selectDecayMode("Z0",["Z0->e-,e+;",
	"Z0->mu-,mu+;",
	"Z0->tau-,tau+;",
	"Z0->nu_e,nu_ebar;",
	"Z0->nu_mu,nu_mubar;",
	"Z0->nu_tau,nu_taubar;"])
	process+=addBRReweighter()
	process+=addLeptonPairCut("60","120")
	elif "W-Z-e" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p p e+ nu","0","2","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p p e- nu","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=particlegroup('epm','e+','e-')
	process+=particlegroup('epmnu','e+','e-','nu_e','nu_ebar')
	process+=addProcess(thefactory,"p p epm epmnu","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("60","120")
	elif "W-Z-mu" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p p mu+ nu","0","2","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p p mu- nu","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=particlegroup('mupm','mu+','mu-')
	process+=particlegroup('mupmnu','mu+','mu-','nu_mu','nu_mubar')
	process+=addProcess(thefactory,"p p mupm mupmnu","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("60","120")
	elif "W-e" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p e+ nu","0","2","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p p e- nu","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=particlegroup('epm','e+','e-')
	process+=addProcess(thefactory,"p p epm nu","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("60","120")

	elif "W-mu" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p mu+ nu","0","2","LeptonPairMassScale",0,0)
	process+=addProcess(thefactory,"p p mu- nu","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=particlegroup('mupm','mu+','mu-')
	process+=addProcess(thefactory,"p p mupm nu","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("60","120")

	elif "Z-e" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("60","120")
	elif "Z-mu" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("60","120")
	elif "Z-jj" in parameterName :
	if(simulation=="Merging"):
	logging.warning("Z-jj not explicitly tested for %s " % simulation)
	sys.exit(0)
	process+=addProcess(thefactory,"p p e+ e- j j","2","2","LeptonPairMassScale",0,0)
	process+=addFirstJet("40")
	process+=addSecondJet("30")
	process+=addLeptonPairCut("60","120")
	elif "Z-LowMass-e" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("20","70")
	elif "Z-MedMass-e" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("40","130")
	elif "Z-LowMass-mu" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",2,2)
	process+=addLeptonPairCut("10","70")
	elif "Z-Mass1" in parameterName :
	process+=addLeptonPairCut("10","35")
	if "-e" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",2,2)
	else :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",2,2)
	elif "Z-Mass2" in parameterName :
	process+=addLeptonPairCut("25","70")
	if "-e" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",2,2)
	else :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",2,2)
	elif "Z-Mass3" in parameterName :
	process+=addLeptonPairCut("60","120")
	if "-e" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",2,2)
	else :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",2,2)
	elif "Z-Mass4" in parameterName :
	process+=addLeptonPairCut("115","8000")
	if "-e" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",2,2)
	else :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",2,2)
	elif "Z-HighMass1" in parameterName :
	process+=addLeptonPairCut("116","400")
	if "-e" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",2,2)
	else :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",2,2)
	elif "Z-HighMass2" in parameterName :
	process+=addLeptonPairCut("400","7000")
	if "-e" in parameterName :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",2,2)
	else :
	if(simulation=="Matchbox"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",0,0)
	elif(simulation=="Merging"):
	process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",2,2)
	elif "W-Jet" in parameterName :
	if(simulation=="Merging"):
	logging.warning("W-Jet not explicitly tested for %s " % simulation)
	sys.exit(0)

	process+=addProcess(thefactory,"p p e+ nu j","1","2","HTScale",0,0)
	process+=addProcess(thefactory,"p p e- nu j","1","2","HTScale",0,0)

	process+=addLeptonPairCut("60","120")
	if "W-Jet-1-e" in parameterName :
	process+=addFirstJet("100")
	parameterName=parameterName.replace("W-Jet-1-e","W-Jet-e")
	elif "W-Jet-2-e" in parameterName :
	process+=addFirstJet("190")
	parameterName=parameterName.replace("W-Jet-2-e","W-Jet-e")
	elif "W-Jet-3-e" in parameterName :
	process+=addFirstJet("270")
	parameterName=parameterName.replace("W-Jet-3-e","W-Jet-e")
	else :
	logging.error("Exit 00005")
	sys.exit(1)
	elif "Z-Jet" in parameterName :
	if(simulation=="Merging"):
	logging.warning("Z-Jet not explicitly tested for %s " % simulation)
	sys.exit(0)


	if "-e" in parameterName :
	process+=addProcess(thefactory,"p p e+ e- j","1","2","HTScale",0,0)
	if "Z-Jet-0-e" in parameterName :
	process+=addFirstJet("35")
	parameterName=parameterName.replace("Z-Jet-0-e","Z-Jet-e")
	elif "Z-Jet-1-e" in parameterName :
	process+=addFirstJet("100")
	parameterName=parameterName.replace("Z-Jet-1-e","Z-Jet-e")
	elif "Z-Jet-2-e" in parameterName :
	process+=addFirstJet("190")
	parameterName=parameterName.replace("Z-Jet-2-e","Z-Jet-e")
	elif "Z-Jet-3-e" in parameterName :
	process+=addFirstJet("270")
	parameterName=parameterName.replace("Z-Jet-3-e","Z-Jet-e")
	else :
	logging.error("Exit 00006")
	sys.exit(1)
	else :
	process+=addProcess(thefactory,"p p mu+ mu- j","1","2","HTScale",0,0)
	process+=addFirstJet("35")
	parameterName=parameterName.replace("Z-Jet-0-mu","Z-Jet-mu")
	process+=addLeptonPairCut("60","120")
	elif "Z-bb" in parameterName :
	if(simulation=="Merging"):
	logging.warning("Z-bb not explicitly tested for %s " % simulation)
	sys.exit(0)
	parameters["bscheme"]=fourFlavour
	process+="set /Herwig/Particles/b:HardProcessMass 4.2GeV\nset /Herwig/Particles/bbar:HardProcessMass 4.2GeV\n"
	process+=addProcess(thefactory,"p p e+ e- b bbar","2","2","FixedScale",0,0)
	process+=addLeptonPairCut("66","116")
	process+=addFirstJet("18")
	process+=addSecondJet("15")
	process+=addLeptonPairCut("60","120")
	elif "Z-b" in parameterName :
	if(simulation=="Merging"):
	logging.warning("Z-b not explicitly tested for %s " % simulation)
	sys.exit(0)
	process+=particlegroup('bjet','b','bbar')
	process+=addProcess(thefactory,"p p e+ e- bjet","1","2","FixedScale",0,0)
	process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 91.2*GeV\n"
	process+=addLeptonPairCut("60","120")
	process+=addFirstJet("15")
	elif "W-b" in parameterName :
	if(simulation=="Merging"):
	logging.warning("W-b not explicitly tested for %s " % simulation)
	sys.exit(0)
	parameters["bscheme"]=fourFlavour
	process += "set /Herwig/Particles/b:HardProcessMass 4.2GeV\nset /Herwig/Particles/bbar:HardProcessMass 4.2GeV\n"
	process+=addProcess(thefactory,"p p e- nu b bbar","2","2","FixedScale",0,0)
	process+=addProcess(thefactory,"p p mu+ nu b bbar","2","2","FixedScale",0,0)
	process += "set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 80.4*GeV\n"
	process+=addFirstJet("30")
	process+=addLeptonPairCut("60","120")
	else :
	logging.error(" Process %s not supported for Matchbox matrix elements" % name)
	sys.exit(1)
	# LHC-GammaGamma
	elif(collider=="LHC-GammaGamma" ) :
	if "-7-" in parameterName : process = StringBuilder(collider_lumi(7000.0))
	elif "-8-" in parameterName : process = StringBuilder(collider_lumi(8000.0))
	else : process = StringBuilder(collider_lumi(7000.0))
	if(simulation=="") :
	if "7" in parameterName : process += insert_ME("MEgg2ff","Muon")
	else :
	logging.error(" Process %s not supported for default matrix elements" % name)
	sys.exit(1)
	else :
	logging.error("LHC-GammaGamma not supported for %s " % simulation)
	sys.exit(1)


	parameters['parameterFile'] = os.path.join(collider,"{c}-{pn}.in".format(c=collider, pn=parameterName))
	parameters['runname'] = 'Rivet-%s' % name
	parameters['process'] = str(process)
	if have_hadronic_collider :
	if collider == "EHS" :
	parameters['collider'] = "PPCollider.in\nread snippets/FixedTarget-PP.in"
	else :
	parameters['collider'] = "PPCollider.in"

	#check if selecteddecaymode and addedBRReweighter is consistent

	if selecteddecaymode and not addedBRReweighter:
	logging.error("Decaymode was selected but no BRReweighter was added.")
	sys.exit(1)

	if addedBRReweighter and not selecteddecaymode:
	logging.error("BRReweighter was added but no Decaymode was selected.")
	sys.exit(1)

	# check that we only add one process if in merging mode:

	if numberOfAddedProcesses > 1 and simulation =="Merging":
	logging.error("In Merging only one process is allowed at the moment. See ticket #403.")
	sys.exit(1)

	# Check if a process was added for Merging or Matchbox:

	if numberOfAddedProcesses == 0 and (simulation =="Merging" or simulation =="Matchbox"):
	logging.error("No process was selected.")
	sys.exit(1)

	# get template and write the file

	with open(os.path.join("Rivet/Templates",templateName), 'r') as f:
	templateText = f.read()

	template = Template( templateText )

	with open(os.path.join("Rivet",name+".in"), 'w') as f:
	f.write( template.substitute(parameters) )
	diff --git a/Tests/python/merge-GammaGamma b/Tests/python/merge-GammaGamma
	new file mode 100755
	--- /dev/null
	+++ b/Tests/python/merge-GammaGamma
	@@ -0,0 +1,124 @@
	+#! /usr/bin/env python
	+import logging,sys
	+import os, yoda, copy
	+
	+if sys.version_info[:3] < (2,4,0):
	+ print "rivet scripts require Python version >= 2.4.0... exiting"
	+ sys.exit(1)
	+
	+
	+#############################################
	+
	+def fillAbove(desthisto, sourcehistosbysqrts):
	+ print "in fill above"
	+ if type(desthisto) is yoda.core.Scatter2D :
	+ for sqrts in sorted(sourcehistosbysqrts.keys()) :
	+ h=sourcehistosbysqrts[sqrts]
	+ for i in range(0,h.numPoints) :
	+ if sqrts==h.points[i].x :
	+ desthisto.addPoint(h.points[i])
	+
	+ elif(type(desthisto)==yoda.core.Profile1D) :
	+ for sqrts, h in sorted(sourcehistosbysqrts.iteritems()) :
	+ for i in range(0,h.numBins) :
	+ if(sqrts>=h.bins[i].xMin and \
	+ sqrts<=h.bins[i].xMax) :
	+ desthisto.bins[i] += h.bins[i]
	+ break
	+ else :
	+ logging.error("Unknown analysis object" + desthisto.path)
	+ sys.exit(1)
	+
	+def merge(hpath):
	+ global inhistos
	+ global outhistos
	+ try:
	+ fillAbove(outhistos[hpath], inhistos[hpath])
	+ except:
	+ pass
	+
	+def useOne(hpath, sqrts):
	+ global inhistos
	+ global outhistos
	+ try:
	+ outhistos[hpath] = inhistos[hpath][float(sqrts)]
	+ except:
	+ pass
	+
	+if __name__ == "__main__":
	+ import logging
	+ from optparse import OptionParser, OptionGroup
	+ parser = OptionParser(usage="%prog name")
	+ verbgroup = OptionGroup(parser, "Verbosity control")
	+ verbgroup.add_option("-v", "--verbose", action="store_const", const=logging.DEBUG, dest="LOGLEVEL",
	+ default=logging.INFO, help="print debug (very verbose) messages")
	+ verbgroup.add_option("-q", "--quiet", action="store_const", const=logging.WARNING, dest="LOGLEVEL",
	+ default=logging.INFO, help="be very quiet")
	+ parser.add_option_group(verbgroup)
	+ (opts, args) = parser.parse_args()
	+ logging.basicConfig(level=opts.LOGLEVEL, format="%(message)s")
	+ ## Check args
	+ if len(args) < 1:
	+ logging.error("Must specify at least the name of the files")
	+ sys.exit(1)
	+
	+#######################################
	+
	+yodafiles=["-mumu-3.5","-mumu-4.5","-mumu-5.5","-mumu-6.5","-mumu-7.5","-mumu-9.0","-mumu-12.5","-mumu-17.5","-mumu-30.0"]
	+
	+# Get histos
	+inhistos = {}
	+outhistos = {}
	+for f in yodafiles:
	+ file = "Rivet-%s%s.yoda" % (args[0], f)
	+ sqrts=float(f.split("-")[-1].replace(".yoda",""))
	+ if not os.access(file, os.R_OK):
	+ logging.error("%s cannot be read" % file)
	+ continue
	+ try:
	+ aos = yoda.read(file)
	+ except:
	+ logging.error("%s cannot be parsed as yoda" % file)
	+ continue
	+ ## Get histos from this YODA file
	+ for aopath, ao in aos.iteritems() :
	+ if("RAW" in aopath or "_XSEC" in aopath or "_EVTCOUNT" in aopath ) :continue
	+ # merge of different energy values
	+ if("L3_2004_I645127" in aopath) :
	+ if not inhistos.has_key(aopath):
	+ inhistos[aopath] = {}
	+ if not inhistos[aopath].has_key(sqrts):
	+ inhistos[aopath][sqrts] = ao
	+ else:
	+ raise Exception("A set with sqrts = %s already exists" % ( sqrts))
	+ else :
	+ outhistos[aopath] = ao
	+## Make empty output histos if needed
	+for hpath,hsets in inhistos.iteritems():
	+ if("L3_2004_I645127" in hpath ) :
	+ if(type(hsets.values()[0])==yoda.core.Scatter2D) :
	+ outhistos[hpath] = yoda.core.Scatter2D(hsets.values()[0].path,
	+ hsets.values()[0].title)
	+ elif(type(hsets.values()[0])==yoda.core.Profile1D) :
	+ outhistos[hpath] = yoda.core.Profile1D(hsets.values()[0].path,
	+ hsets.values()[0].title)
	+ for i in range(0,hsets.values()[0].numBins) :
	+ outhistos[hpath].addBin(hsets.values()[0].bins[i].xMin,
	+ hsets.values()[0].bins[i].xMax)
	+ elif(type(hsets.values()[0])==yoda.core.Histo1D) :
	+ outhistos[hpath] = yoda.core.Histo1D(hsets.values()[0].path,
	+ hsets.values()[0].title)
	+ for i in range(0,hsets.values()[0].numBins) :
	+ outhistos[hpath].addBin(hsets.values()[0].bins[i].xMin,
	+ hsets.values()[0].bins[i].xMax)
	+ else :
	+ logging.error("Histogram %s is of unknown type" % hpath)
	+ sys.exit(1)
	+merge("/L3_2004_I645127:PROCESS=GG/d03-x01-y05")
	+
	+# Choose output file
	+name = args[0]+".yoda"
	+# output the yoda file
	+print "Write yoda to ",name
	+yoda.writeYODA(outhistos,name)
	+sys.exit(0)
	diff --git a/Tests/python/merge-LEP-Gamma b/Tests/python/merge-LEP-Gamma
	new file mode 100755
	--- /dev/null
	+++ b/Tests/python/merge-LEP-Gamma
	@@ -0,0 +1,136 @@
	+#! /usr/bin/env python
	+import logging,sys
	+import os, yoda, copy
	+
	+if sys.version_info[:3] < (2,4,0):
	+ print "rivet scripts require Python version >= 2.4.0... exiting"
	+ sys.exit(1)
	+
	+
	+#############################################
	+
	+def fillAbove(desthisto, sourcehistosbysqrts):
	+ if type(desthisto) is yoda.core.Scatter2D :
	+ for sqrts in sorted(sourcehistosbysqrts.keys()) :
	+ h=sourcehistosbysqrts[sqrts]
	+ for i in range(0,h.numPoints) :
	+ if sqrts==h.points[i].x :
	+ desthisto.addPoint(h.points[i])
	+
	+ elif(type(desthisto)==yoda.core.Profile1D) :
	+ for sqrts, h in sorted(sourcehistosbysqrts.iteritems()) :
	+ for i in range(0,h.numBins) :
	+ if(sqrts>=h.bins[i].xMin and \
	+ sqrts<=h.bins[i].xMax) :
	+ desthisto.bins[i] += h.bins[i]
	+ break
	+ else :
	+ logging.error("Unknown analysis object" + desthisto.path)
	+ sys.exit(1)
	+
	+def merge(hpath):
	+ global inhistos
	+ global outhistos
	+ try:
	+ fillAbove(outhistos[hpath], inhistos[hpath])
	+ except:
	+ pass
	+
	+def useOne(hpath, sqrts):
	+ global inhistos
	+ global outhistos
	+ try:
	+ outhistos[hpath] = inhistos[hpath][float(sqrts)]
	+ except:
	+ pass
	+
	+if __name__ == "__main__":
	+ import logging
	+ from optparse import OptionParser, OptionGroup
	+ parser = OptionParser(usage="%prog name")
	+ verbgroup = OptionGroup(parser, "Verbosity control")
	+ verbgroup.add_option("-v", "--verbose", action="store_const", const=logging.DEBUG, dest="LOGLEVEL",
	+ default=logging.INFO, help="print debug (very verbose) messages")
	+ verbgroup.add_option("-q", "--quiet", action="store_const", const=logging.WARNING, dest="LOGLEVEL",
	+ default=logging.INFO, help="be very quiet")
	+ parser.add_option_group(verbgroup)
	+ (opts, args) = parser.parse_args()
	+ logging.basicConfig(level=opts.LOGLEVEL, format="%(message)s")
	+ ## Check args
	+ if len(args) < 1:
	+ logging.error("Must specify at least the name of the files")
	+ sys.exit(1)
	+
	+#######################################
	+
	+yodafiles=["-Direct-mumu-161","-Direct-mumu-172","-Direct-mumu-183",
	+ "-Direct-mumu-189","-Direct-mumu-196","-Direct-mumu-206",
	+ "-Direct-tautau-189","-Direct-tautau-196","-Direct-tautau-206",
	+ "-Direct-Jets-198","-Single-Resolved-Jets-198","-Double-Resolved-Jets-198",
	+ "-Direct-Jets-206","-Single-Resolved-Jets-206","-Double-Resolved-Jets-206",]
	+
	+# Get histos
	+inhistos = {}
	+outhistos = {}
	+for f in yodafiles:
	+ file = "Rivet-%s%s.yoda" % (args[0], f)
	+ sqrts=float(f.split("-")[-1].replace(".yoda",""))
	+ if not os.access(file, os.R_OK):
	+ logging.error("%s cannot be read" % file)
	+ continue
	+ try:
	+ aos = yoda.read(file)
	+ except:
	+ logging.error("%s cannot be parsed as yoda" % file)
	+ continue
	+ ## Get histos from this YODA file
	+ for aopath, ao in aos.iteritems() :
	+ if("RAW" in aopath or "_XSEC" in aopath or "_EVTCOUNT" in aopath ) :continue
	+ # merge of different energy values
	+ if("L3_2004_I645127" in aopath) :
	+ if(("d01" in aopath and "mu" in file) or
	+ ("d02" in aopath and "tau" in file)) :
	+ if not inhistos.has_key(aopath):
	+ inhistos[aopath] = {}
	+ if not inhistos[aopath].has_key(sqrts):
	+ inhistos[aopath][sqrts] = ao
	+ else:
	+ raise Exception("A set with sqrts = %s already exists" % ( sqrts))
	+ else :
	+ if(aopath in outhistos) :
	+ outhistos[aopath] += ao
	+ else :
	+ outhistos[aopath] = ao
	+## Make empty output histos if needed
	+for hpath,hsets in inhistos.iteritems():
	+ if("L3_2004_I645127" in hpath ) :
	+ if(type(hsets.values()[0])==yoda.core.Scatter2D) :
	+ outhistos[hpath] = yoda.core.Scatter2D(hsets.values()[0].path,
	+ hsets.values()[0].title)
	+ elif(type(hsets.values()[0])==yoda.core.Profile1D) :
	+ outhistos[hpath] = yoda.core.Profile1D(hsets.values()[0].path,
	+ hsets.values()[0].title)
	+ for i in range(0,hsets.values()[0].numBins) :
	+ outhistos[hpath].addBin(hsets.values()[0].bins[i].xMin,
	+ hsets.values()[0].bins[i].xMax)
	+ elif(type(hsets.values()[0])==yoda.core.Histo1D) :
	+ outhistos[hpath] = yoda.core.Histo1D(hsets.values()[0].path,
	+ hsets.values()[0].title)
	+ for i in range(0,hsets.values()[0].numBins) :
	+ outhistos[hpath].addBin(hsets.values()[0].bins[i].xMin,
	+ hsets.values()[0].bins[i].xMax)
	+ else :
	+ logging.error("Histogram %s is of unknown type" % hpath)
	+ sys.exit(1)
	+
	+merge("/L3_2004_I645127/d01-x01-y01")
	+merge("/L3_2004_I645127/d01-x01-y02")
	+merge("/L3_2004_I645127/d02-x01-y01")
	+
	+
	+# Choose output file
	+name = args[0]+".yoda"
	+# output the yoda file
	+print "Write yoda to ",name
	+yoda.writeYODA(outhistos,name)
	+sys.exit(0)

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 11:26 AM (17 h, 34 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5111423
Default Alt Text: (226 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions