No OneTemporary
Actions

Size

141 KB

Subscribers

None

View Options

	diff --git a/Cuts/Cuts.cc b/Cuts/Cuts.cc
	--- a/Cuts/Cuts.cc
	+++ b/Cuts/Cuts.cc
	@@ -1,578 +1,615 @@
	// -- C++ --
	//
	// Cuts.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the Cuts class.
	//

	#include "Cuts.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/PDT/ParticleData.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/EventRecord/SubProcess.h"
	#include "ThePEG/EventRecord/Collision.h"
	#include "ThePEG/EventRecord/TmpTransform.h"
	#include "ThePEG/Utilities/UtilityBase.h"
	#include "ThePEG/Utilities/HoldFlag.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Config/algorithm.h"

	using namespace ThePEG;

	Cuts::Cuts(Energy MhatMin)
	: theSMax(ZERO), theY(0), theCurrentSHat(-1.0*GeV2),
	theCurrentYHat(0), theMHatMin(MhatMin), theMHatMax(Constants::MaxEnergy),
	theYHatMin(-Constants::MaxRapidity), theYHatMax(Constants::MaxRapidity),
	theX1Min(0.0), theX1Max(1.0), theX2Min(0.0), theX2Max(1.0),
	theScaleMin(ZERO), theScaleMax(Constants::MaxEnergy2),
	- theSubMirror(false) {}
	+ theSubMirror(false), theCutWeight(1.0), theLastCutWeight(1.0) {}

	Cuts::~Cuts() {}

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

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

	void Cuts::doinitrun() {
	Interfaced::doinitrun();
	if ( Debug::level ) {
	describe();
	for_each(theOneCuts, mem_fun(&OneCutBase::describe));
	for_each(theTwoCuts, mem_fun(&TwoCutBase::describe));
	for_each(theMultiCuts, mem_fun(&MultiCutBase::describe));
	}
	}

	void Cuts::describe() const {
	CurrentGenerator::log()
	<< fullName() << ":\n"
	<< "MHat = " << theMHatMin/GeV << " .. " << theMHatMax/GeV << " GeV\n"
	<< "Scale = " << theScaleMin/GeV2 << " .. " << theScaleMax/GeV2 << " GeV2\n"
	<< "YHat = " << theYHatMin << " .. " << theYHatMax << '\n'
	<< "X1 = " << theX1Min << " .. " << theX1Max << '\n'
	<< "X2 = " << theX2Min << " .. " << theX2Max << "\n\n";
	}

	void Cuts::initialize(Energy2 smax, double Y) {
	theSMax = smax;
	theMHatMax = min(theMHatMax, sqrt(smax));
	theY = Y;
	theSubMirror = false;
	}

	void Cuts::initEvent() {
	theCurrentSHat = -1.0*GeV2;
	theCurrentYHat = 0.0;
	theSubMirror = false;
	}

	bool Cuts::initSubProcess(Energy2 shat, double yhat, bool mirror) const {
	+ theCutWeight = 1.0;
	+ theLastCutWeight = 1.0;
	theSubMirror = mirror;
	theCurrentSHat = shat;
	theCurrentYHat = yhat;
	if ( shat <= sHatMin() \|\| shat > sHatMax()(1.0 + 1000.0Constants::epsilon) ) return false;
	if ( yhat <= yHatMin() \|\| yhat >= yHatMax() ) return false;
	double x1 = min(1.0, sqrt(shat/SMax())*exp(yhat));
	if ( x1 <= x1Min() \|\| x1 > x1Max() ) return false;
	double x2 = min(1.0, sqrt(shat/SMax())*exp(-yhat));
	if ( x2 <= x2Min() \|\| x2 > x2Max() ) return false;
	return true;
	}

	bool Cuts::passCuts(const tcPDVector & ptype, const vector<LorentzMomentum> & p,
	tcPDPtr t1, tcPDPtr t2) const {
	if ( subMirror() ) {
	vector<LorentzMomentum> pmir = p;
	for ( int i = 0, N = pmir.size(); i < N; ++i ) pmir[i].setZ(-pmir[i].z());
	swap(t1,t2);
	HoldFlag<> nomir(theSubMirror, false);
	return passCuts(ptype, pmir, t1, t2);
	}

	+ bool pass = true;
	+ theCutWeight = 1.0;
	+ theLastCutWeight = 1.0;
	+
	+ // ATTENTION fuzzy jet finding
	if ( jetFinder() )
	if ( ptype.size() > jetFinder()->minOutgoing() ) {
	vector<LorentzMomentum> jets = p;
	tcPDVector jettype = ptype;
	if ( jetFinder()->cluster(jettype,jets,this,t1,t2) ){
	return passCuts(jettype,jets,t1,t2);
	}
	}

	for ( int i = 0, N = p.size(); i < N; ++i )
	- for ( int j = 0, M = theOneCuts.size(); j < M; ++j )
	- if ( !theOneCuts[j]->passCuts(this, ptype[i], p[i]) ) return false;
	+ for ( int j = 0, M = theOneCuts.size(); j < M; ++j ) {
	+ pass &= theOneCuts[j]->passCuts(this, ptype[i], p[i]);
	+ theCutWeight *= theLastCutWeight;
	+ theLastCutWeight = 1.0;
	+ if ( !pass )
	+ return false;
	+ }

	for ( int i1 = 0, N1 = p.size() - 1; i1 < N1; ++i1 )
	for ( int i2 = i1 + 1, N2 = p.size(); i2 < N2; ++i2 )
	- for ( int j = 0, M = theTwoCuts.size(); j < M; ++j )
	- if ( !theTwoCuts[j]->passCuts(this, ptype[i1], ptype[i2],
	- p[i1], p[i2]) ) return false;
	- for ( int j = 0, M = theMultiCuts.size(); j < M; ++j )
	- if ( !theMultiCuts[j]->passCuts(this, ptype, p) ) return false;
	+ for ( int j = 0, M = theTwoCuts.size(); j < M; ++j ) {
	+ pass &= theTwoCuts[j]->passCuts(this, ptype[i1], ptype[i2],
	+ p[i1], p[i2]);
	+ theCutWeight *= theLastCutWeight;
	+ theLastCutWeight = 1.0;
	+ if ( !pass )
	+ return false;
	+ }
	+
	+ for ( int j = 0, M = theMultiCuts.size(); j < M; ++j ) {
	+ pass &= theMultiCuts[j]->passCuts(this, ptype, p);
	+ theCutWeight *= theLastCutWeight;
	+ theLastCutWeight = 1.0;
	+ if ( !pass )
	+ return false;
	+ }
	+
	if ( t1 ) {
	LorentzMomentum p1(ZERO, ZERO, 0.5*sqrt(currentSHat()),
	0.5*sqrt(currentSHat()));
	for ( int i = 0, N = p.size(); i < N; ++i )
	- for ( int j = 0, M = theTwoCuts.size(); j < M; ++j )
	- if ( !theTwoCuts[j]->passCuts(this, t1, ptype[i], p1, p[i],
	- true, false) )
	+ for ( int j = 0, M = theTwoCuts.size(); j < M; ++j ) {
	+ pass &= theTwoCuts[j]->passCuts(this, t1, ptype[i], p1, p[i],
	+ true, false);
	+ theCutWeight *= theLastCutWeight;
	+ theLastCutWeight = 1.0;
	+ if ( !pass )
	return false;
	+ }
	}
	+
	if ( t2 ) {
	LorentzMomentum p2(ZERO, ZERO,
	-0.5sqrt(currentSHat()), 0.5sqrt(currentSHat()));
	for ( int i = 0, N = p.size(); i < N; ++i )
	- for ( int j = 0, M = theTwoCuts.size(); j < M; ++j )
	- if ( !theTwoCuts[j]->passCuts(this, ptype[i], t2, p[i], p2,
	- false, true) )
	+ for ( int j = 0, M = theTwoCuts.size(); j < M; ++j ) {
	+ pass &= theTwoCuts[j]->passCuts(this, ptype[i], t2, p[i], p2,
	+ false, true);
	+ theCutWeight *= theLastCutWeight;
	+ theLastCutWeight = 1.0;
	+ if ( !pass )
	return false;
	+ }
	}
	- return true;
	+
	+ return pass;
	+
	}

	bool Cuts::passCuts(const tcPVector & p, tcPDPtr t1, tcPDPtr t2) const {
	tcPDVector ptype(p.size());
	vector<LorentzMomentum> mom(p.size());
	for ( int i = 0, N = p.size(); i < N; ++i ) {
	ptype[i] = p[i]->dataPtr();
	mom[i] = p[i]->momentum();
	}
	return passCuts(ptype, mom, t1, t2);
	}

	bool Cuts::passCuts(const SubProcess & sub) const {
	if ( !passCuts(tcPVector(sub.outgoing().begin(), sub.outgoing().end()),
	sub.incoming().first->dataPtr(),
	sub.incoming().second->dataPtr()) ) return false;
	return true;
	}

	bool Cuts::passCuts(const Collision & coll) const {
	tSubProPtr sub = coll.primarySubProcess();
	LorentzMomentum phat = sub->incoming().first->momentum() +
	sub->incoming().second->momentum();
	if ( !initSubProcess(phat.m2(), phat.rapidity()) ) return false;
	TmpTransform<tSubProPtr> tmp(sub, Utilities::getBoostToCM(sub->incoming()));
	if ( !passCuts(*sub) ) return false;
	return true;
	}

	Energy2 Cuts::minS(const tcPDVector & pv) const {
	Energy2 mins = ZERO;
	for ( int i = 0, N = theMultiCuts.size(); i < N; ++i )
	mins = max(mins, theMultiCuts[i]->minS(pv));
	return mins;
	}

	Energy2 Cuts::maxS(const tcPDVector & pv) const {
	Energy2 maxs = SMax();
	for ( int i = 0, N = theMultiCuts.size(); i < N; ++i )
	maxs = min(maxs, theMultiCuts[i]->maxS(pv));
	return maxs;
	}

	Energy2 Cuts::minSij(tcPDPtr pi, tcPDPtr pj) const {
	Energy2 mins = ZERO;
	for ( int i = 0, N = theTwoCuts.size(); i < N; ++i )
	mins = max(mins, theTwoCuts[i]->minSij(pi, pj));
	if ( mins > ZERO ) return mins;
	mins = sqr(pi->massMin() + pj->massMin());
	mins = max(mins, sqr(minKTClus(pi, pj))/4.0);
	mins = max(mins, minDurham(pi, pj)*currentSHat()/2.0);
	mins = max(mins, minKT(pi)minKT(pj)minDeltaR(pi, pj)/4.0);
	return mins;
	}

	Energy2 Cuts::minTij(tcPDPtr pi, tcPDPtr po) const {
	Energy2 mint = ZERO;
	for ( int i = 0, N = theTwoCuts.size(); i < N; ++i )
	mint = max(mint, theTwoCuts[i]->minTij(pi, po));
	if ( mint > ZERO ) return mint;
	mint = max(mint, sqr(minKT(po)));
	return mint;
	}

	double Cuts::minDeltaR(tcPDPtr pi, tcPDPtr pj) const {
	double mindr = 0.0;
	for ( int i = 0, N = theTwoCuts.size(); i < N; ++i )
	mindr = max(mindr, theTwoCuts[i]->minDeltaR(pi, pj));
	return mindr;
	}

	Energy Cuts::minKTClus(tcPDPtr pi, tcPDPtr pj) const {
	Energy minkt = ZERO;
	for ( int i = 0, N = theTwoCuts.size(); i < N; ++i )
	minkt = max(minkt, theTwoCuts[i]->minKTClus(pi, pj));
	return minkt;
	}

	double Cuts::minDurham(tcPDPtr pi, tcPDPtr pj) const {
	double y = 0.0;
	for ( int i = 0, N = theTwoCuts.size(); i < N; ++i )
	y = max(y, theTwoCuts[i]->minDurham(pi, pj));
	return y;
	}

	Energy Cuts::minKT(tcPDPtr p) const {
	Energy minkt = ZERO;
	for ( int i = 0, N = theOneCuts.size(); i < N; ++i )
	minkt = max(minkt, theOneCuts[i]->minKT(p));
	if ( minkt > ZERO ) return minkt;
	minkt = minKTClus(p, tcPDPtr());
	return minkt;
	}

	double Cuts::minEta(tcPDPtr p) const {
	double mineta = -Constants::MaxRapidity;
	for ( int i = 0, N = theOneCuts.size(); i < N; ++i )
	mineta = max(mineta, theOneCuts[i]->minEta(p));
	return mineta;
	}

	double Cuts::maxEta(tcPDPtr p) const {
	double maxeta = Constants::MaxRapidity;
	for ( int i = 0, N = theOneCuts.size(); i < N; ++i )
	maxeta = min(maxeta, theOneCuts[i]->maxEta(p));
	return maxeta;
	}

	double Cuts::minYStar(tcPDPtr p) const {
	if ( currentSHat() < ZERO ) return -Constants::MaxRapidity;
	if ( subMirror() ) {
	HoldFlag<> nomir(theSubMirror, false);
	return -maxYStar(p);
	}
	double etamin = minEta(p);
	double ytot = Y() + currentYHat();
	if ( etamin > 0.0 ) {
	Energy minkt = minKT(p);
	Energy maxm = p->massMax();
	return asinh(minkt*sinh(etamin)/sqrt(sqr(minkt) + sqr(maxm))) - ytot;
	} else {
	return etamin - ytot;
	}
	}

	double Cuts::maxYStar(tcPDPtr p) const {
	if ( currentSHat() < ZERO ) return Constants::MaxRapidity;
	if ( subMirror() ) {
	HoldFlag<> nomir(theSubMirror, false);
	return -minYStar(p);
	}
	double etamax = maxEta(p);
	double ytot = Y() + currentYHat();
	if ( etamax > 0.0 ) {
	return etamax - ytot;
	} else {
	Energy minkt = minKT(p);
	Energy maxm = p->massMax();
	return asinh(minkt*sinh(etamax)/sqrt(sqr(minkt) + sqr(maxm))) - ytot;
	}
	}

	double Cuts::minRapidityMax(tcPDPtr p) const {
	double minRapidityMax = -Constants::MaxRapidity;
	for ( int i = 0, N = theOneCuts.size(); i < N; ++i )
	minRapidityMax = max(minRapidityMax, theOneCuts[i]->minRapidityMax(p));
	return minRapidityMax;
	}

	double Cuts::maxRapidityMin(tcPDPtr p) const {
	double maxRapidityMin = Constants::MaxRapidity;
	for ( int i = 0, N = theOneCuts.size(); i < N; ++i )
	maxRapidityMin = min(maxRapidityMin, theOneCuts[i]->maxRapidityMin(p));
	return maxRapidityMin;
	}

	void Cuts::persistentOutput(PersistentOStream & os) const {
	os << ounit(theSMax, GeV2) << theY << ounit(theCurrentSHat, GeV2)
	<< theCurrentYHat << ounit(theMHatMin, GeV) << ounit(theMHatMax, GeV)
	<< theYHatMin << theYHatMax
	<< theX1Min << theX1Max << theX2Min << theX2Max << ounit(theScaleMin, GeV2)
	<< ounit(theScaleMax, GeV2) << theOneCuts << theTwoCuts << theMultiCuts
	- << theJetFinder << theSubMirror;
	+ << theJetFinder << theSubMirror
	+ << theCutWeight << theLastCutWeight;
	}

	void Cuts::persistentInput(PersistentIStream & is, int) {
	is >> iunit(theSMax, GeV2) >> theY >> iunit(theCurrentSHat, GeV2)
	>> theCurrentYHat >> iunit(theMHatMin, GeV) >> iunit(theMHatMax, GeV)
	>> theYHatMin >> theYHatMax
	>> theX1Min >> theX1Max >> theX2Min >> theX2Max >> iunit(theScaleMin, GeV2)
	>> iunit(theScaleMax, GeV2) >> theOneCuts >> theTwoCuts >> theMultiCuts
	- >> theJetFinder >> theSubMirror;
	+ >> theJetFinder >> theSubMirror
	+ >> theCutWeight >> theLastCutWeight;
	}

	ClassDescription<Cuts> Cuts::initCuts;
	// Definition of the static class description member.

	Energy Cuts::maxMHatMin() const {
	return theMHatMax;
	}

	Energy Cuts::minMHatMax() const {
	return theMHatMin;
	}

	Energy2 Cuts::maxScaleMin() const {
	return theScaleMax;
	}

	Energy2 Cuts::minScaleMax() const {
	return theScaleMin;
	}

	double Cuts::maxYHatMin() const {
	return theYHatMax;
	}

	double Cuts::minYHatMax() const {
	return theYHatMin;
	}

	double Cuts::maxX1Min() const {
	return theX1Max;
	}

	double Cuts::minX1Max() const {
	return theX1Min;
	}

	double Cuts::maxX2Min() const {
	return theX2Max;
	}

	double Cuts::minX2Max() const {
	return theX2Min;
	}

	void Cuts::Init() {

	typedef double (ThePEG::Cuts::*IGFN)() const;
	typedef void (ThePEG::Cuts::*ISFN)(double);

	static ClassDocumentation<Cuts> documentation
	("Cuts is a class for implementing kinematical cuts in ThePEG. The "
	"class itself only implements cuts on the total momentum of the hard "
	"sub-process, implemented as minimum and maximum values of \\f$x_1\\f$ "
	"and \\f$x_2\\f$ (or \\f$\\hat{s}\\f$ and \\f$\\hat{y}\\f$. Further cuts "
	"can be implemented either by inheriting from this base class, in which "
	"the virtual cut() function should be overridden, or by assigning "
	"objects of class OneCutBase, TwoCutBase and MultiCutBase defining "
	"cuts on single particles, pairs of particles and groups of "
	"particles respectively.");

	static Parameter<Cuts,Energy> interfaceMHatMin
	("MHatMin",
	"The minimum allowed value of \\f$\\sqrt{\\hat{s}}\\f$.",
	&Cuts::theMHatMin, GeV, 2.0*GeV, ZERO, Constants::MaxEnergy,
	true, false, Interface::limited,
	0, 0, 0, &Cuts::maxMHatMin, 0);
	interfaceMHatMin.setHasDefault(false);

	static Parameter<Cuts,Energy> interfaceMHatMax
	("MHatMax",
	"The maximum allowed value of \\f$\\sqrt{\\hat{s}}\\f$.",
	&Cuts::theMHatMax, GeV, 100.0*GeV, ZERO, ZERO,
	true, false, Interface::lowerlim,
	0, 0,
	&Cuts::minMHatMax, 0, 0);
	interfaceMHatMax.setHasDefault(false);

	static Parameter<Cuts,Energy2> interfaceScaleMin
	("ScaleMin",
	"The minimum allowed value of the scale to be used in PDFs and "
	"coupling constants.",
	&Cuts::theScaleMin, GeV2, ZERO, ZERO, Constants::MaxEnergy2,
	true, false, Interface::limited,
	0, 0, 0, &Cuts::maxScaleMin, 0);
	interfaceScaleMin.setHasDefault(false);

	static Parameter<Cuts,Energy2> interfaceScaleMax
	("ScaleMax",
	"The maximum allowed value of the scale to be used in PDFs and "
	"coupling constants.",
	&Cuts::theScaleMax, GeV2, 10000.0*GeV2, ZERO, ZERO,
	true, false, Interface::lowerlim,
	0, 0,
	&Cuts::minScaleMax, 0, 0);
	interfaceScaleMax.setHasDefault(false);

	static Parameter<Cuts,double> interfaceYHatMin
	("YHatMin",
	"The minimum value of the rapidity of the hard sub-process "
	"(wrt. the rest system of the colliding particles).",
	&Cuts::theYHatMin, -10.0, 0.0, Constants::MaxRapidity,
	true, false, Interface::upperlim,
	(ISFN)0, (IGFN)0, (IGFN)0, &Cuts::maxYHatMin, (IGFN)0);
	interfaceYHatMin.setHasDefault(false);

	static Parameter<Cuts,double> interfaceYHatMax
	("YHatMax",
	"The maximum value of the rapidity of the hard sub-process "
	"(wrt. the rest system of the colliding particles).",
	&Cuts::theYHatMax, 10.0, -Constants::MaxRapidity, 0.0,
	true, false, Interface::lowerlim,
	(ISFN)0, (IGFN)0, &Cuts::minYHatMax, (IGFN)0, (IGFN)0);
	interfaceYHatMax.setHasDefault(false);

	static Parameter<Cuts,double> interfaceX1Min
	("X1Min",
	"The minimum value of the positive light-cone fraction of the hard "
	"sub-process.",
	&Cuts::theX1Min, 0.0, 0.0, 1.0,
	true, false, Interface::limited,
	(ISFN)0, (IGFN)0, (IGFN)0, &Cuts::maxX1Min, (IGFN)0);
	interfaceX1Min.setHasDefault(false);

	static Parameter<Cuts,double> interfaceX1Max
	("X1Max",
	"The maximum value of the positive light-cone fraction of the hard "
	"sub-process.",
	&Cuts::theX1Max, 0.0, 0.0, 1.0,
	true, false, Interface::limited,
	(ISFN)0, (IGFN)0, &Cuts::minX1Max, (IGFN)0, (IGFN)0);
	interfaceX1Max.setHasDefault(false);

	static Parameter<Cuts,double> interfaceX2Min
	("X2Min",
	"The minimum value of the negative light-cone fraction of the hard "
	"sub-process.",
	&Cuts::theX2Min, 0.0, 0.0, 1.0,
	true, false, Interface::limited,
	(ISFN)0, (IGFN)0, (IGFN)0, &Cuts::maxX2Min, (IGFN)0);
	interfaceX2Min.setHasDefault(false);

	static Parameter<Cuts,double> interfaceX2Max
	("X2Max",
	"The maximum value of the negative light-cone fraction of the hard "
	"sub-process.",
	&Cuts::theX2Max, 0.0, 0.0, 1.0,
	true, false, Interface::limited,
	(ISFN)0, (IGFN)0, &Cuts::minX2Max, (IGFN)0, (IGFN)0);
	interfaceX2Max.setHasDefault(false);

	static RefVector<Cuts,OneCutBase> interfaceOneCuts
	("OneCuts",
	"The objects defining cuts on single outgoing partons from the "
	"hard sub-process.",
	&Cuts::theOneCuts, -1, true, false, true, false, false);

	static RefVector<Cuts,TwoCutBase> interfaceTwoCuts
	("TwoCuts",
	"The objects defining cuts on pairs of particles in the "
	"hard sub-process.",
	&Cuts::theTwoCuts, -1, true, false, true, false, false);

	static RefVector<Cuts,MultiCutBase> interfaceMultiCuts
	("MultiCuts",
	"The objects defining cuts on sets of outgoing particles from the "
	"hard sub-process.",
	&Cuts::theMultiCuts, -1, true, false, true, false, false);


	static Reference<Cuts,JetFinder> interfaceJetFinder
	("JetFinder",
	"Set a JetFinder object used to define cuts on the"
	"level of reconstructed jets as needed for higher order corrections.",
	&Cuts::theJetFinder, false, false, true, true, false);


	interfaceX1Min.rank(10);
	interfaceX1Max.rank(9);
	interfaceX2Min.rank(8);
	interfaceX2Max.rank(7);
	interfaceMHatMin.rank(6);
	interfaceMHatMax.rank(5);
	interfaceYHatMin.rank(4);
	interfaceYHatMax.rank(3);
	interfaceOneCuts.rank(2);
	interfaceTwoCuts.rank(1);

	}

	double Cuts::yHatMin() const {
	return theX1Min > 0.0 && theX2Max > 0.0?
	max(theYHatMin, 0.5*log(theX1Min/theX2Max)): theYHatMin;
	}

	double Cuts::yHatMax() const {
	return theX1Max > 0.0 && theX2Min > 0.0?
	min(theYHatMax, 0.5*log(theX1Max/theX2Min)): theYHatMax;
	}

	bool Cuts::yHat(double y) const {
	return y > yHatMin() && y < yHatMax();
	}

	double Cuts::x1Min() const {
	return max(theX1Min, (theMHatMin/sqrt(SMax()))*exp(theYHatMin));
	}

	double Cuts::x1Max() const {
	return min(theX1Max, (theMHatMax/sqrt(SMax()))*exp(theYHatMax));
	}

	bool Cuts::x1(double x) const {
	return x > x1Min() && x <= x1Max();
	}

	double Cuts::x2Min() const {
	return max(theX2Min, (theMHatMin/sqrt(SMax()))/exp(theYHatMax));
	}

	double Cuts::x2Max() const {
	return min(theX2Max, (theMHatMax/sqrt(SMax()))/exp(theYHatMin));
	}

	bool Cuts::x2(double x) const {
	return x > x2Min() && x <= x2Max();
	}

	template <typename T>
	vector<typename Ptr<T>::transient_const_pointer>
	Cuts::oneCutObjects() const {
	typedef typename Ptr<T>::transient_const_pointer tcPtr;
	vector<tcPtr> ret;
	for ( int i = 0, N = theOneCuts.size(); i < N; ++i )
	if ( dynamic_ptr_cast<tcPtr>(theOneCuts[i]) )
	ret.push_back(dynamic_ptr_cast<tcPtr>(theOneCuts[i]));
	return ret;
	}

	template <typename T>
	vector<typename Ptr<T>::transient_const_pointer>
	Cuts::twoCutObjects() const {
	typedef typename Ptr<T>::transient_const_pointer tcPtr;
	vector<tcPtr> ret;
	for ( int i = 0, N = theTwoCuts.size(); i < N; ++i )
	if ( dynamic_ptr_cast<tcPtr>(theTwoCuts[i]) )
	ret.push_back(dynamic_ptr_cast<tcPtr>(theTwoCuts[i]));
	return ret;
	}

	template <typename T>
	vector<typename Ptr<T>::transient_const_pointer>
	Cuts::multiCutObjects() const {
	typedef typename Ptr<T>::transient_const_pointer tcPtr;
	vector<tcPtr> ret;
	for ( int i = 0, N = theMultiCuts.size(); i < N; ++i )
	if ( dynamic_ptr_cast<tcPtr>(theMultiCuts[i]) )
	ret.push_back(dynamic_ptr_cast<tcPtr>(theMultiCuts[i]));
	return ret;
	}
	diff --git a/Cuts/Cuts.h b/Cuts/Cuts.h
	--- a/Cuts/Cuts.h
	+++ b/Cuts/Cuts.h
	@@ -1,763 +1,785 @@
	// -- C++ --
	//
	// Cuts.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef THEPEG_Cuts_H
	#define THEPEG_Cuts_H
	//
	// This is the declaration of the Cuts class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "Cuts.fh"
	#include "OneCutBase.h"
	#include "TwoCutBase.h"
	#include "MultiCutBase.h"
	#include "JetFinder.h"

	namespace ThePEG {

	/**
	* Cuts is a class for implementing kinematical cuts in ThePEG. The
	* class itself only implements cuts on the total momentum of the hard
	* sub-process, implemented as minimum and maximum values of \f$x_1\f$
	* and \f$x_2\f$ (or \f$\hat{s}=x_1x_2S_{tot}\f$ and
	* \f$\hat{y}=\log(x_1/x_2)/2\f$. Further cuts can be implemented
	* either by inheriting from this base class, in which the virtual
	* cut() function should be overridden, or by assigning objects of
	* class OneCutBase, TwoCutBase and MultiCutBase defining cuts on
	* single particles, pairs of particles and groups of particles in the
	* hard sub-process respectively.
	*
	* The Cuts object must be initialized specifying the overall
	* laboratory frame, giving the total squared invariant mass, \f$S\f$,
	* and the rapidity, \f$Y\f$, of the colliding particles in this
	* frame. The colliding particles are thus assumed to be directed
	* along the \f$z\f$-axis.
	*
	* For each event, the Cuts object must also be initialized giving the
	* squared invarint mass, \f$\hat{s}\f$, and the total rapidity,
	* \f$\hat{y}\f$, of the hard sub-process in the center-of-mass frame
	* of the colliding particles. Note that this means that the
	* transformation between the lab frame and the rest frame of the hard
	* sub-process is assumed to be a simple boost along the z-axis.
	*
	* @see \ref CutsInterfaces "The interfaces"
	* defined for Cuts.
	*/
	class Cuts: public Interfaced {

	public:

	/**
	* A vector of OneCutBase pointers.
	*/
	typedef vector<OneCutPtr> OneCutVector;

	/**
	* A vector of TwoCutBase pointers.
	*/
	typedef vector<TwoCutPtr> TwoCutVector;

	/**
	* A vector of MultiCutBase pointers.
	*/
	typedef vector<MultiCutPtr> MultiCutVector;

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	Cuts(Energy MhatMin=2*GeV);

	/**
	* The destructor.
	*/
	virtual ~Cuts();
	//@}

	public:

	/** @name Initialization functions. */
	//@{
	/**
	* Initialize this object specifying the maximum total invariant
	* mass squared, \a smax, and the total rapidity, \a Y, of the
	* colliding particles (for the maximum invariant mass). A sub-class
	* overriding this function must make sure the base-class function
	* is called. This function should be called once in the beginning
	* of a run.
	*/
	virtual void initialize(Energy2 smax, double Y);

	/**
	* Initialize this object for a new event. A sub-class overriding
	* this function must make sure the base-class function is called.
	* This function is called before the generation of a new
	* sub-process, before the incoming partons have been generated.
	*/
	virtual void initEvent();

	/**
	* Set information about the invariant mass squared, \a shat, and
	* rapidity, \a yhat, of the hard sub-process. The rapidity should
	* be given wrt. the center of mass of the colliding particles. A
	* sub-class overriding this function must make sure the base-class
	* function is called. This function is called before the generation
	* of a new sub-process, after the incoming partons have been
	* generated. If \a mirror is true any questions regarding cuts on
	* the sub-process in the functions minYStar(tcPDPtr),
	* maxYStar(tcPDPtr p), passCuts(const tcPDVector &, const
	* vector<LorentzMomentum> &, tcPDPtr, tcPDPtr) and passCuts(const
	* tcPVector &, tcPDPtr t1, tcPDPtr) will assume that the z-axis is
	* reversed in the sub-process rest frame. Returns false if the
	* given values were outside of the cuts.
	*/
	virtual bool
	initSubProcess(Energy2 shat, double yhat, bool mirror = false) const;
	//@}

	/** @name Check functions to see if a state has passed the cuts or not. */
	//@{
	/**
	* Check if the outgoing particles, with the given types and
	* momenta, from a sub-process passes the cuts. The particles must
	* be given in the rest frame of tha hard sub-process, and the
	* initSubProcess must have been called before. Also the types of
	* the incoming partons, \a t1 and \a t2, may be given if availible.
	*/
	virtual bool passCuts(const tcPDVector & ptype, const vector<LorentzMomentum> & p,
	tcPDPtr t1 = tcPDPtr(), tcPDPtr t2 = tcPDPtr()) const;

	/**
	* Check if the outgoing particles from a sub-process passes the
	* cuts. The particles must be given in the rest frame of tha hard
	* sub-process, and the initSubProcess must have been called
	* before. Also the types of the incoming partons, \a t1 and \a t2,
	* may be given if availible.
	*/
	bool passCuts(const tcPVector & p,
	tcPDPtr t1 = tcPDPtr(), tcPDPtr t2 = tcPDPtr()) const;

	/**
	* Check if the incoming and outgoing particles in the given
	* sub-process passes the cuts. The sub-process must be given in its
	* rest frame, and the initSubProcess must have been called before.
	*/
	bool passCuts(const SubProcess & sub) const;

	/**
	* Check if the given collision passes the cuts. The collision must
	* be given in its rest frame.
	*/
	bool passCuts(const Collision & coll) const;
	//@}

	/** @name Access to cuts of the underlying cut objects. */
	//@{
	/**
	* Return the minimum allowed squared invariant mass of two outgoing
	* partons of type \a pi and \a pj. This function first determines
	* the minimum from the corresponding function from in TwoCutBase
	* objects. If no minimum was found, one is derived from
	* minKTClus(), minDurham(), minKT() and minDeltaR(), if possible.
	*/
	Energy2 minSij(tcPDPtr pi, tcPDPtr pj) const;

	/**
	* Return the minimum allowed value of the negative of the squared
	* invariant mass of an incoming parton of type \a pi and an
	* outgoing parton of type \a po. This function first determines the
	* minimum from the corresponding function from in TwoCutBase
	* objects. If no minimum was found, one is derived from minKT(), if
	* possible.
	*/
	Energy2 minTij(tcPDPtr pi, tcPDPtr po) const;

	/**
	* Return the minimum allowed value of \f$\Delta
	* R_{ij}=\sqrt{\Delta\eta_{ij}^2+\Delta\phi_{ij}^2}\f$ of two
	* outgoing partons of type \a pi and \a pj. Simply returns the
	* maximum of the results from calling the corresponding function in
	* the TwoCutBase objects.
	*/
	double minDeltaR(tcPDPtr pi, tcPDPtr pj) const;

	/**
	* Return the minimum allowed value of the longitudinally invariant
	* \f$k_\perp\f$-algorithms distance measure. This is defined as
	* \f$\min(p_{\perp i}, p_{\perp
	* j})\sqrt{\Delta\eta_{ij}^2+\Delta\phi_{ij}^2}\f$ for two outgoing
	* partons, or simply \f$p_{\perp i}\f$ or \f$p_{\perp j}\f$ for a
	* single outgoing parton. Returns 0 if both partons are incoming. A
	* null pointer indicates an incoming parton, hence the type of the
	* incoming parton is irrelevant. Simply returns the maximum of the
	* results from calling the corresponding function in the TwoCutBase
	* objects.
	*/
	Energy minKTClus(tcPDPtr pi, tcPDPtr pj) const;

	/**
	* Return the minimum allowed value of the Durham
	* \f$k_\perp\f$-algorithms distance measure. This is defined as
	* \f$2\min(E_j^2, E_j^2)(1-\cos\theta_{ij})/\hat{s}\f$ for two
	* outgoing partons. Simply returns the maximum of the results from
	* calling the corresponding function in the TwoCutBase objects.
	*/
	double minDurham(tcPDPtr pi, tcPDPtr pj) const;

	/**
	* Return the minimum allowed value of the transverse momentum of an
	* outgoing parton. This function first determines the minimum from
	* the corresponding function from in OneCutBase objects. If no
	* minimum was found, one is derived from minKTClus(), if possible.
	*/
	Energy minKT(tcPDPtr p) const;

	/**
	* Return the minimum allowed pseudo-rapidity of an outgoing parton
	* of the given type. The pseudo-rapidity is measured in the lab
	* system. Simply returns the maximum of the results from calling
	* the corresponding function in the OneCutBase objects.
	*/
	double minEta(tcPDPtr p) const;

	/**
	* Return the maximum allowed pseudo-rapidity of an outgoing parton
	* of the given type. The pseudo-rapidity is measured in the lab
	* system. Simply returns the minimum of the results from calling
	* the corresponding function in the OneCutBase objects.
	*/
	double maxEta(tcPDPtr p) const;

	/**
	* Return the minimum allowed rapidity of an outgoing parton
	* of the given type. The rapidity is measured in the lab
	* system. Simply returns the maximum of the results from calling
	* the corresponding function in the OneCutBase objects.
	*/
	double minRapidityMax(tcPDPtr p) const;

	/**
	* Return the maximum allowed rapidity of an outgoing parton
	* of the given type. The rapidity is measured in the lab
	* system. Simply returns the minimum of the results from calling
	* the corresponding function in the OneCutBase objects.
	*/
	double maxRapidityMin(tcPDPtr p) const;

	/**
	* Return the minimum allowed rapidity of an outgoing parton of the
	* given type in the center-of-mass system of the hard sub-process.
	* Only available after initSubProcess() has been called.
	*/
	double minYStar(tcPDPtr p) const;

	/**
	* Return the minimum allowed rapidity of an outgoing parton of the
	* given type in the center-of-mass system of the hard sub-process.
	* Only available after initSubProcess() has been called.
	*/
	double maxYStar(tcPDPtr p) const;

	/**
	* Return the minimum allowed value of the squared invariant mass of
	* a set of outgoing partons of the given types. Typically used to
	* cut off the tails of the mass of a resonance for
	* efficiency. Simply returns the maximum of the results from
	* calling the corresponding function in the MultiCutBase objects.
	*/
	Energy2 minS(const tcPDVector & pv) const;

	/**
	* Return the maximum allowed value of the squared invariant mass of
	* a set of outgoing partons of the given types. Typically used to
	* cut off the tails of the mass of a resonance for
	* efficiency. Simply returns the minimum of the results from
	* calling the corresponding function in the MultiCutBase objects.
	*/
	Energy2 maxS(const tcPDVector & pv) const;
	//@}

	/** @name Direct access to underlying cut objects. */
	//@{
	/**
	* Return a vector of pointers to objects of the given class (with
	* base class OneCutBase).
	*/
	template <typename T>
	vector<typename Ptr<T>::transient_const_pointer>
	oneCutObjects() const;

	/**
	* Return a vector of pointers to objects of the given class (with
	* base class TwoCutBase).
	*/
	template <typename T>
	vector<typename Ptr<T>::transient_const_pointer>
	twoCutObjects() const;

	/**
	* Return a vector of pointers to objects of the given class (with
	* base class MultiCutBase).
	*/
	template <typename T>
	vector<typename Ptr<T>::transient_const_pointer>
	multiCutObjects() const;

	/**
	* Return the objects defining cuts on single outgoing partons from the
	* hard sub-process.
	*/
	const OneCutVector& oneCuts() const { return theOneCuts; }

	/**
	* Return the objects defining cuts on pairs of particles in the hard
	* sub-process.
	*/
	const TwoCutVector& twoCuts() const { return theTwoCuts; }

	/**
	* Return the objects defining cuts on sets of outgoing particles from the
	* hard sub-process.
	*/
	const MultiCutVector& multiCuts() const { return theMultiCuts; }

	/**
	* Return the jet finder
	*/
	Ptr<JetFinder>::tptr jetFinder() const { return theJetFinder; }

	/**
	* Add a OneCutBase object.
	*/
	void add(tOneCutPtr c) { theOneCuts.push_back(c); }

	/**
	* Add a TwoCutBase object.
	*/
	void add(tTwoCutPtr c) { theTwoCuts.push_back(c); }

	/**
	* Add a MultiCutBase object.
	*/
	void add(tMultiCutPtr c) { theMultiCuts.push_back(c); }
	//@}

	public:

	/** @name Simple access functions. */
	//@{
	/**
	* The maximum allowed total invariant mass squared allowed for
	* events to be considered.
	*/
	Energy2 SMax() const { return theSMax; }


	/**
	* The total rapidity of the colliding particles corresponding to
	* the maximum invariant mass squared, SMax().
	*/
	double Y() const { return theY; }

	/**
	* The invariant mass squared of the hard sub-process of the event
	* being considered.
	*/
	Energy2 currentSHat() const { return theCurrentSHat; }

	/**
	* The total rapidity of hard sub-process (wrt. the rest system of
	* the colliding particles so that currentYHat() + Y() gives the
	* true rapidity) of the event being considered.
	*/
	double currentYHat() const { return theCurrentYHat; }

	//@}

	/** @name Functions to inquire about specific cuts. */
	//@{
	/**
	* The minimum allowed value of \f$\hat{s}\f$.
	*/
	Energy2 sHatMin() const { return max(sqr(theMHatMin), theX1MintheX2MinSMax()); }

	/**
	* The maximum allowed value of \f$\hat{s}\f$.
	*/
	Energy2 sHatMax() const { return min(sqr(theMHatMax), theX1MaxtheX2MaxSMax()); }

	/**
	* Check if the given \f$\hat{s}\f$ is within the cuts.
	*/
	bool sHat(Energy2 sh) const {
	return sh > sHatMin() && sh <= sHatMax()(1.0 + 1000.0Constants::epsilon);
	}

	/**
	* The minimum allowed value of \f$\sqrt{\hat{s}}\f$.
	*/
	Energy mHatMin() const { return max(theMHatMin, sqrt(theX1MintheX2MinSMax())); }

	/**
	* The maximum allowed value of \f$\sqrt{\hat{s}}\f$.
	*/
	Energy mHatMax() const { return min(theMHatMax, sqrt(theX1MaxtheX2MaxSMax())); }

	/**
	* The minimum value of the rapidity of the hard sub-process
	* (wrt. the rest system of the colliding particles).
	*/
	double yHatMin() const;

	/**
	* The maximum value of the rapidity of the hard sub-process
	* (wrt. the rest system of the colliding particles).
	*/
	double yHatMax() const;

	/**
	* Check if the given \f$\hat{y}\f$ is within the cuts.
	*/
	bool yHat(double y) const;

	/**
	* The minimum value of the positive light-cone fraction of the hard
	* sub-process.
	*/
	double x1Min() const;

	/**
	* The maximum value of the positive light-cone fraction of the hard
	* sub-process.
	*/
	double x1Max() const;

	/**
	* Check if the given \f$x_1\f$ is within the cuts.
	*/
	bool x1(double x) const;

	/**
	* The minimum value of the negative light-cone fraction of the hard
	* sub-process.
	*/
	double x2Min() const;

	/**
	* The maximum value of the negative light-cone fraction of the hard
	* sub-process.
	*/
	double x2Max() const;

	/**
	* Check if the given \f$x_2\f$ is within the cuts.
	*/
	bool x2(double x) const;

	/**
	* The minimum allowed value of the scale to be used in PDF's and
	* coupling constants.
	*/
	Energy2 scaleMin() const { return theScaleMin; }

	/**
	* The maximum allowed value of the scale to be used in PDF's and
	* coupling constants.
	*/
	Energy2 scaleMax() const { return theScaleMax; }

	/**
	* Check if the given scale is within the cuts.
	*/
	bool scale(Energy2 Q2) const { return Q2 > scaleMin() && Q2 < scaleMax(); }

	/**
	* Set true if a matrix element is should be using this cut and is
	* mirrored along the z-axis .
	*/
	bool subMirror() const { return theSubMirror; }
	+
	+ /**
	+ * Return the overall cut weight
	+ */
	+ double cutWeight() const { return theCutWeight; }
	+
	+ /**
	+ * Set the cut weight as appropriate from the call to the last n-cut
	+ * object.
	+ */
	+ void lastCutWeight(double w) const { theLastCutWeight = w; }
	//@}

	public:

	/**
	* Describe the currently active cuts in the log file.
	*/
	virtual void describe() const;

	protected:

	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();
	//@}

	public:

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

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

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

	protected:

	/** @name 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:

	/**
	* Helper function used by the interface.
	*/
	Energy maxMHatMin() const;

	/**
	* Helper function used by the interface.
	*/
	Energy minMHatMax() const;

	/**
	* Helper function used by the interface.
	*/
	double maxYHatMin() const;

	/**
	* Helper function used by the interface.
	*/
	double minYHatMax() const;

	/**
	* Helper function used by the interface.
	*/
	double maxX1Min() const;

	/**
	* Helper function used by the interface.
	*/
	double minX1Max() const;

	/**
	* Helper function used by the interface.
	*/
	double maxX2Min() const;

	/**
	* Helper function used by the interface.
	*/
	double minX2Max() const;

	/**
	* Helper function used by the interface.
	*/
	Energy2 maxScaleMin() const;

	/**
	* Helper function used by the interface.
	*/
	Energy2 minScaleMax() const;

	private:

	/**
	* The maximum allowed total invariant mass squared allowed for
	* events to be considered.
	*/
	Energy2 theSMax;

	/**
	* The total rapidity of the colliding particles corresponding to
	* the maximum invariant mass squared, SMax().
	*/
	double theY;

	/**
	* The invariant mass squared of the hard sub-process of the event
	* being considered.
	*/
	mutable Energy2 theCurrentSHat;

	/**
	* The total rapidity of hard sub-process (wrt. the rest system of
	* the colliding particles so that currentYHat() + Y() gives the
	* true rapidity) of the event being considered.
	*/
	mutable double theCurrentYHat;

	/**
	* The minimum allowed value of \f$\sqrt{\hat{s}}\f$.
	*/
	Energy theMHatMin;

	/**
	* The maximum allowed value of \f$\sqrt{\hat{s}}\f$.
	*/
	Energy theMHatMax;

	/**
	* The minimum value of the rapidity of the hard sub-process
	* (wrt. the rest system of the colliding particles).
	*/
	double theYHatMin;

	/**
	* The maximum value of the rapidity of the hard sub-process
	* (wrt. the rest system of the colliding particles).
	*/
	double theYHatMax;

	/**
	* The minimum value of the positive light-cone fraction of the hard
	* sub-process.
	*/
	double theX1Min;

	/**
	* The maximum value of the positive light-cone fraction of the hard
	* sub-process.
	*/
	double theX1Max;

	/**
	* The minimum value of the negative light-cone fraction of the hard
	* sub-process.
	*/
	double theX2Min;

	/**
	* The maximum value of the negative light-cone fraction of the hard
	* sub-process.
	*/
	double theX2Max;

	/**
	* The minimum allowed value of the scale to be used in PDF's and
	* coupling constants.
	*/
	Energy2 theScaleMin;

	/**
	* The maximum allowed value of the scale to be used in PDF's and
	* coupling constants.
	*/
	Energy2 theScaleMax;

	/**
	* The objects defining cuts on single outgoing partons from the
	* hard sub-process.
	*/
	OneCutVector theOneCuts;

	/**
	* The objects defining cuts on pairs of particles in the hard
	* sub-process.
	*/
	TwoCutVector theTwoCuts;

	/**
	* The objects defining cuts on sets of outgoing particles from the
	* hard sub-process.
	*/
	MultiCutVector theMultiCuts;

	/**
	* An optional jet finder used to define cuts on the level of
	* reconstructed jets.
	*/
	Ptr<JetFinder>::ptr theJetFinder;

	/**
	* Set to true if a matrix element is should be using this cut and is
	* mirrored along the z-axis .
	*/
	mutable bool theSubMirror;

	+ /**
	+ * The overall cut weight
	+ */
	+ mutable double theCutWeight;
	+
	+ /**
	+ * The cut weight as appropriate from the call to the last n-cut
	+ * object.
	+ */
	+ mutable double theLastCutWeight;
	+
	private:

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

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

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the Cuts class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Cuts>
	: public ClassTraitsBase<Cuts> {
	/** Return a platform-independent class name */
	static string className() { return "ThePEG::Cuts"; }
	};

	/** @endcond */

	}

	#endif /* THEPEG_Cuts_H */
	diff --git a/Cuts/JetCuts.cc b/Cuts/JetCuts.cc
	--- a/Cuts/JetCuts.cc
	+++ b/Cuts/JetCuts.cc
	@@ -1,205 +1,232 @@
	// -- C++ --
	//
	// JetCuts.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	// Copyright (C) 2009-2012 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the JetCuts class.
	//

	#include "JetCuts.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/PDT/ParticleData.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Repository/CurrentGenerator.h"

	using namespace ThePEG;

	JetCuts::JetCuts()
	: theOrdering(orderPt) {}

	JetCuts::~JetCuts() {}

	void JetCuts::describe() const {
	CurrentGenerator::log() << name() << " cutting on jets ordered in "
	<< (ordering() == orderPt ? "pt " : "y ")
	<< " with regions:\n";
	for ( vector<Ptr<JetRegion>::ptr>::const_iterator r = jetRegions().begin();
	r != jetRegions().end(); ++r )
	(**r).describe();
	if ( !jetPairRegions().empty() ) {
	CurrentGenerator::log() << "and jet pair cuts:\n";
	for ( vector<Ptr<JetPairRegion>::ptr>::const_iterator r = jetPairRegions().begin();
	r != jetPairRegions().end(); ++r )
	(**r).describe();
	for ( vector<Ptr<MultiJetRegion>::ptr>::const_iterator r = multiJetRegions().begin();
	r != multiJetRegions().end(); ++r )
	(**r).describe();
	}
	if ( !jetVetoRegions().empty() ) {
	CurrentGenerator::log() << "vetoing jets inside:\n";
	for ( vector<Ptr<JetRegion>::ptr>::const_iterator r = jetVetoRegions().begin();
	r != jetVetoRegions().end(); ++r )
	(**r).describe();
	}
	+ CurrentGenerator::log() << "\n";
	}

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

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

	void JetCuts::persistentOutput(PersistentOStream & os) const {
	os << theUnresolvedMatcher
	<< theJetRegions << theJetVetoRegions
	- << theJetPairRegions << theMultiJetRegions << theOrdering;
	+ << theJetPairRegions << theMultiJetRegions
	+ << theOrdering;
	}

	void JetCuts::persistentInput(PersistentIStream & is, int) {
	is >> theUnresolvedMatcher
	>> theJetRegions >> theJetVetoRegions
	- >> theJetPairRegions >> theMultiJetRegions >> theOrdering;
	+ >> theJetPairRegions >> theMultiJetRegions
	+ >> theOrdering;
	}

	struct PtLarger {
	inline bool operator()(const LorentzMomentum& a,
	const LorentzMomentum& b) const {
	return a.perp() > b.perp();
	}
	};

	struct YLess {
	inline bool operator()(const LorentzMomentum& a,
	const LorentzMomentum& b) const {
	return a.rapidity() < b.rapidity();
	}
	};

	bool JetCuts::passCuts(tcCutsPtr parent, const tcPDVector & ptype,
	const vector<LorentzMomentum> & p) const {

	vector<LorentzMomentum> jets;
	tcPDVector::const_iterator ptypeit = ptype.begin();
	vector<LorentzMomentum>::const_iterator pit = p.begin();
	for ( ; ptypeit != ptype.end(); ++ptypeit, ++pit )
	if ( unresolvedMatcher()->check(**ptypeit) )
	jets.push_back(*pit);
	if ( ordering() == orderPt ) {
	sort(jets.begin(),jets.end(),PtLarger());
	} else if ( ordering() == orderY ) {
	sort(jets.begin(),jets.end(),YLess());
	} else assert(false);

	for ( vector<Ptr<JetRegion>::ptr>::const_iterator r = jetRegions().begin();
	r != jetRegions().end(); ++r )
	(**r).reset();

	set<int> matchedJets;

	for ( size_t k = 0; k < jets.size(); ++k ) {
	for ( vector<Ptr<JetRegion>::ptr>::const_iterator r = jetRegions().begin();
	r != jetRegions().end(); ++r ) {
	if ( (**r).matches(parent,k+1,jets[k]) ) {
	matchedJets.insert(k+1);
	break;
	}
	}
	}

	+ double cutWeight = 1.0;
	+ bool pass = true;
	+
	for ( vector<Ptr<JetRegion>::ptr>::const_iterator r = jetRegions().begin();
	- r != jetRegions().end(); ++r )
	- if ( !(**r).didMatch() )
	+ r != jetRegions().end(); ++r ) {
	+ pass &= (**r).didMatch();
	+ cutWeight = (*r).cutWeight();
	+ if ( !pass ) {
	+ parent->lastCutWeight(cutWeight);
	return false;
	+ }
	+ }

	for ( vector<Ptr<JetPairRegion>::ptr>::const_iterator r = jetPairRegions().begin();
	- r != jetPairRegions().end(); ++r )
	- if ( !(**r).matches(parent) )
	+ r != jetPairRegions().end(); ++r ) {
	+ pass &= (**r).matches(parent);
	+ cutWeight = (*r).cutWeight();
	+ if ( !pass ) {
	+ parent->lastCutWeight(cutWeight);
	return false;
	+ }
	+ }

	for ( vector<Ptr<MultiJetRegion>::ptr>::const_iterator r = multiJetRegions().begin();
	- r != multiJetRegions().end(); ++r )
	- if ( !(**r).matches() )
	+ r != multiJetRegions().end(); ++r ) {
	+ pass &= (**r).matches();
	+ cutWeight = (*r).cutWeight();
	+ if ( !pass ) {
	+ parent->lastCutWeight(cutWeight);
	return false;
	+ }
	+ }

	if ( !jetVetoRegions().empty() ) {
	for ( size_t k = 0; k < jets.size(); ++k ) {
	if ( matchedJets.find(k+1) != matchedJets.end() )
	continue;
	for ( vector<Ptr<JetRegion>::ptr>::const_iterator r = jetVetoRegions().begin();
	r != jetVetoRegions().end(); ++r ) {
	- if ( (**r).matches(parent,k+1,jets[k]) )
	+ pass &= !(**r).matches(parent,k+1,jets[k]);
	+ cutWeight = 1. - (*r).cutWeight();
	+ if ( !pass ) {
	+ parent->lastCutWeight(cutWeight);
	return false;
	+ }
	}
	}
	}

	- return true;
	+ parent->lastCutWeight(cutWeight);
	+
	+ return pass;

	}

	ClassDescription<JetCuts> JetCuts::initJetCuts;
	// Definition of the static class description member.

	void JetCuts::Init() {

	static ClassDocumentation<JetCuts> documentation
	("JetCuts combines various JetRegion and JetPairRegion objects into a "
	"cut object.");

	static Reference<JetCuts,MatcherBase> interfaceUnresolvedMatcher
	("UnresolvedMatcher",
	"A matcher identifying unresolved partons",
	&JetCuts::theUnresolvedMatcher, false, false, true, false, false);

	static RefVector<JetCuts,JetRegion> interfaceJetRegions
	("JetRegions",
	"The jet regions to be used.",
	&JetCuts::theJetRegions, -1, false, false, true, false, false);

	static RefVector<JetCuts,JetRegion> interfaceJetVetoRegions
	("JetVetoRegions",
	"The jet veto regions to be used.",
	&JetCuts::theJetVetoRegions, -1, false, false, true, false, false);

	static RefVector<JetCuts,JetPairRegion> interfaceJetPairRegions
	("JetPairRegions",
	"The jet pair regions to be used.",
	&JetCuts::theJetPairRegions, -1, false, false, true, false, false);

	static RefVector<JetCuts,MultiJetRegion> interfaceMultiJetRegions
	("MultiJetRegions",
	"The multi jet regions to be used.",
	&JetCuts::theMultiJetRegions, -1, false, false, true, false, false);


	static Switch<JetCuts,int> interfaceOrdering
	("Ordering",
	"The ordering to apply on jets.",
	&JetCuts::theOrdering, orderPt, false, false);
	static SwitchOption interfaceOrderingOrderPt
	(interfaceOrdering,
	"OrderPt",
	"Order in decreasing transverse momentum.",
	orderPt);
	static SwitchOption interfaceOrderingOrderY
	(interfaceOrdering,
	"OrderY",
	"Order in rapidity.",
	orderY);

	}

	diff --git a/Cuts/JetPairRegion.cc b/Cuts/JetPairRegion.cc
	--- a/Cuts/JetPairRegion.cc
	+++ b/Cuts/JetPairRegion.cc
	@@ -1,190 +1,199 @@
	// -- C++ --
	//
	// JetPairRegion.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 Leif Lonnblad
	// Copyright (C) 2009-2012 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the JetPairRegion class.
	//

	#include "JetPairRegion.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Repository/CurrentGenerator.h"

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

	using namespace ThePEG;

	JetPairRegion::JetPairRegion()
	: theMassMin(0.*GeV), theMassMax(Constants::MaxEnergy),
	theDeltaRMin(0.0), theDeltaRMax(Constants::MaxRapidity),
	theDeltaYMin(0.0), theDeltaYMax(Constants::MaxRapidity),
	- theOppositeHemispheres(false) {}
	+ theOppositeHemispheres(false), theCutWeight(1.0) {}

	JetPairRegion::~JetPairRegion() {}

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

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

	void JetPairRegion::describe() const {

	CurrentGenerator::log()
	<< "JetPairRegion '" << name() << "' matching "
	<< " JetRegions '" << firstRegion()->name()
	<< "' and '" << secondRegion()->name() << "' with\n";

	CurrentGenerator::log()
	<< "m = " << massMin()/GeV << " .. " << massMax()/GeV << " GeV\n"
	<< "dR = " << deltaRMin() << " .. " << deltaRMax() << "\n"
	<< "dy = " << deltaYMin() << " .. " << deltaYMax() << "\n";

	}

	-bool JetPairRegion::matches(tcCutsPtr parent) const {
	+bool JetPairRegion::matches(tcCutsPtr parent) {

	if ( !firstRegion()->didMatch() \|\|
	- !secondRegion()->didMatch() )
	+ !secondRegion()->didMatch() ) {
	+ theCutWeight = 0.0;
	return false;
	+ }
	+
	+ theCutWeight = 1.0;

	const LorentzMomentum& pi = firstRegion()->lastMomentum();
	const LorentzMomentum& pj = secondRegion()->lastMomentum();

	Energy m = (pi+pj).m();
	- if ( m < massMin() \|\| m > massMax() )
	+ if ( !(firstRegion()->lessThanEnergy(massMin(),m,theCutWeight) &&
	+ firstRegion()->lessThanEnergy(m,massMax(),theCutWeight)) )
	return false;

	double dy = abs(pi.rapidity() - pj.rapidity());
	- if ( dy < deltaYMin() \|\| dy > deltaYMax() )
	+ if ( !(firstRegion()->lessThanRapidity(deltaYMin(),dy,theCutWeight) &&
	+ firstRegion()->lessThanRapidity(dy,deltaYMax(),theCutWeight)) )
	return false;

	double dphi = abs(pi.phi() - pj.phi());
	if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
	double dR = sqrt(sqr(dy)+sqr(dphi));
	- if ( dR < deltaRMin() \|\| dR > deltaRMax() )
	+ if ( !(firstRegion()->lessThanRapidity(deltaRMin(),dR,theCutWeight) &&
	+ firstRegion()->lessThanRapidity(dR,deltaRMax(),theCutWeight)) )
	return false;

	double py =
	(pi.rapidity() + parent->currentYHat()) *
	(pj.rapidity() + parent->currentYHat());
	- if ( theOppositeHemispheres && py > 0.0 )
	+ if ( theOppositeHemispheres &&
	+ !firstRegion()->lessThanRapidity(py,0.0,theCutWeight) )
	return false;

	return true;

	}

	// If needed, insert default implementations of virtual function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void JetPairRegion::persistentOutput(PersistentOStream & os) const {
	os << theFirstRegion << theSecondRegion
	<< ounit(theMassMin,GeV) << ounit(theMassMax,GeV)
	<< theDeltaRMin << theDeltaRMax
	<< theDeltaYMin << theDeltaYMax
	- << theOppositeHemispheres;
	+ << theOppositeHemispheres << theCutWeight;
	}

	void JetPairRegion::persistentInput(PersistentIStream & is, int) {
	is >> theFirstRegion >> theSecondRegion
	>> iunit(theMassMin,GeV) >> iunit(theMassMax,GeV)
	>> theDeltaRMin >> theDeltaRMax
	>> theDeltaYMin >> theDeltaYMax
	- >> theOppositeHemispheres;
	+ >> theOppositeHemispheres >> theCutWeight;
	}


	// * Attention * The following static variable is needed for the type
	// description system in ThePEG. Please check that the template arguments
	// are correct (the class and its base class), and that the constructor
	// arguments are correct (the class name and the name of the dynamically
	// loadable library where the class implementation can be found).
	DescribeClass<JetPairRegion,HandlerBase>
	describeThePEGJetPairRegion("ThePEG::JetPairRegion", "JetCuts.so");

	void JetPairRegion::Init() {

	static ClassDocumentation<JetPairRegion> documentation
	("JetPairRegion implements constraints on jets matching two jet regions.");


	static Reference<JetPairRegion,JetRegion> interfaceFirstRegion
	("FirstRegion",
	"The first region to act on.",
	&JetPairRegion::theFirstRegion, false, false, true, false, false);

	static Reference<JetPairRegion,JetRegion> interfaceSecondRegion
	("SecondRegion",
	"The second region to act on.",
	&JetPairRegion::theSecondRegion, false, false, true, false, false);

	static Parameter<JetPairRegion,Energy> interfaceMassMin
	("MassMin",
	"The minimum jet-jet invariant mass.",
	&JetPairRegion::theMassMin, GeV, 0.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<JetPairRegion,Energy> interfaceMassMax
	("MassMax",
	"The maximum jet-jet invariant mass.",
	&JetPairRegion::theMassMax, GeV, Constants::MaxEnergy, 0.0GeV, 0GeV,
	false, false, Interface::lowerlim);

	static Parameter<JetPairRegion,double> interfaceDeltaRMin
	("DeltaRMin",
	"The minimum jet-jet lego-plot separation.",
	&JetPairRegion::theDeltaRMin, 0.0, 0.0, 0,
	false, false, Interface::lowerlim);

	static Parameter<JetPairRegion,double> interfaceDeltaRMax
	("DeltaRMax",
	"The maximum jet-jet lego-plot separation.",
	&JetPairRegion::theDeltaRMax, Constants::MaxRapidity, 0, 0,
	false, false, Interface::lowerlim);

	static Parameter<JetPairRegion,double> interfaceDeltaYMin
	("DeltaYMin",
	"The minimum jet-jet rapidity separation.",
	&JetPairRegion::theDeltaYMin, 0.0, 0.0, 0,
	false, false, Interface::lowerlim);

	static Parameter<JetPairRegion,double> interfaceDeltaYMax
	("DeltaYMax",
	"The maximum jet-jet rapidity separation.",
	&JetPairRegion::theDeltaYMax, Constants::MaxRapidity, 0, 0,
	false, false, Interface::lowerlim);

	static Switch<JetPairRegion,bool> interfaceOppositeHemispheres
	("OppositeHemispheres",
	"Should the jets go into opposite detector hemispheres?",
	&JetPairRegion::theOppositeHemispheres, false, true, false);
	static SwitchOption interfaceOppositeHemispheresTrue
	(interfaceOppositeHemispheres,
	"True",
	"Require jets to be in different hemispheres",
	true);
	static SwitchOption interfaceOppositeHemispheresFalse
	(interfaceOppositeHemispheres,
	"False",
	"Do not require jets to be in different hemispheres",
	false);
	+
	}

	diff --git a/Cuts/JetPairRegion.h b/Cuts/JetPairRegion.h
	--- a/Cuts/JetPairRegion.h
	+++ b/Cuts/JetPairRegion.h
	@@ -1,205 +1,215 @@
	// -- C++ --
	//
	// JetPairRegion.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 Leif Lonnblad
	// Copyright (C) 2009-2012 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef ThePEG_JetPairRegion_H
	#define ThePEG_JetPairRegion_H
	//
	// This is the declaration of the JetPairRegion class.
	//

	#include "ThePEG/Cuts/JetRegion.h"

	namespace ThePEG {

	/**
	* JetPairRegion implements constraints on jets matching two jet regions.
	*
	* @see JetRegion
	* @see JetCuts
	*
	* @see \ref JetPairRegionInterfaces "The interfaces"
	* defined for JetPairRegion.
	*/
	class JetPairRegion: public HandlerBase {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	JetPairRegion();

	/**
	* The destructor.
	*/
	virtual ~JetPairRegion();
	//@}

	public:

	/**
	* Return the first jet region to act on.
	*/
	Ptr<JetRegion>::tptr firstRegion() const { return theFirstRegion; }

	/**
	* Return the second jet region to act on.
	*/
	Ptr<JetRegion>::tptr secondRegion() const { return theSecondRegion; }

	/**
	* Return the minimum jet-jet invariant mass.
	*/
	Energy massMin() const { return theMassMin; }

	/**
	* Return the maximum jet-jet invariant mass.
	*/
	Energy massMax() const { return theMassMax; }

	/**
	* Return the minimum jet-jet lego-plot separation.
	*/
	double deltaRMin() const { return theDeltaRMin; }

	/**
	* Return the maximum jet-jet lego-plot separation.
	*/
	double deltaRMax() const { return theDeltaRMax; }

	/**
	* Return the minimum jet-jet rapidity separation.
	*/
	double deltaYMin() const { return theDeltaYMin; }

	/**
	* Return the maximum jet-jet rapidity separation.
	*/
	double deltaYMax() const { return theDeltaYMax; }

	+ /**
	+ * Return the cut weight encountered from the last call to matches()
	+ */
	+ double cutWeight() const { return theCutWeight; }
	+
	public:

	/**
	* Describe the currently active cuts in the log file.
	*/
	virtual void describe() const;

	/**
	* Return true, if the requirements on the jet regions are fullfilled.
	*/
	- virtual bool matches(tcCutsPtr parent) const;
	+ virtual bool matches(tcCutsPtr parent);

	public:

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

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

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

	protected:

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

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


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


	private:

	/**
	* The first jet region to act on.
	*/
	Ptr<JetRegion>::ptr theFirstRegion;

	/**
	* The second jet region to act on.
	*/
	Ptr<JetRegion>::ptr theSecondRegion;

	/**
	* The minimum jet-jet invariant mass.
	*/
	Energy theMassMin;

	/**
	* The maximum jet-jet invariant mass.
	*/
	Energy theMassMax;

	/**
	* The minimum jet-jet lego-plot separation.
	*/
	double theDeltaRMin;

	/**
	* The maximum jet-jet lego-plot separation.
	*/
	double theDeltaRMax;

	/**
	* The minimum jet-jet rapidity separation.
	*/
	double theDeltaYMin;

	/**
	* The maximum jet-jet rapidity separation.
	*/
	double theDeltaYMax;

	/**
	* Should the jets go into opposite detector hemispheres?
	*/
	bool theOppositeHemispheres;

	/**
	+ * The cut weight encountered from the last call to matches()
	+ */
	+ double theCutWeight;
	+
	+ /**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	JetPairRegion & operator=(const JetPairRegion &);

	};

	}

	#endif /* ThePEG_JetPairRegion_H */
	diff --git a/Cuts/JetRegion.cc b/Cuts/JetRegion.cc
	--- a/Cuts/JetRegion.cc
	+++ b/Cuts/JetRegion.cc
	@@ -1,168 +1,252 @@
	// -- C++ --
	//
	// JetRegion.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 Leif Lonnblad
	// Copyright (C) 2009-2012 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the JetRegion class.
	//

	#include "JetRegion.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Command.h"
	+#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Repository/CurrentGenerator.h"

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

	using namespace ThePEG;

	JetRegion::JetRegion()
	: thePtMin(0.*GeV), thePtMax(Constants::MaxEnergy),
	- theDidMatch(false), theLastNumber(0) {}
	+ theDidMatch(false), theLastNumber(0),
	+ theFuzzy(false), theCutWeight(1.0),
	+ theEnergyCutWidth(1.0*GeV), theRapidityCutWidth(0.1) {}

	JetRegion::~JetRegion() {}

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

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

	string JetRegion::doYRange(string in) {
	istringstream ins(in);
	double first, second;
	ins >> first >> second;
	if ( first > second )
	swap(first,second);
	theYRanges.push_back(make_pair(first,second));
	return "";
	}

	void JetRegion::describe() const {

	CurrentGenerator::log()
	<< "JetRegion '" << name() << "' matching ";
	if ( accepts().empty() )
	CurrentGenerator::log() << "any jets ";
	else {
	CurrentGenerator::log() << "jets ";
	for ( vector<int>::const_iterator k = accepts().begin();
	k != accepts().end(); ++k ) {
	CurrentGenerator::log() << "#" << *k;
	if ( k != --accepts().end() )
	CurrentGenerator::log() << ", ";
	else
	CurrentGenerator::log() << " ";
	}
	}
	CurrentGenerator::log() << " within:\n";

	CurrentGenerator::log()
	<< "pt = " << ptMin()/GeV << " .. " << ptMax()/GeV << " GeV\n";

	for ( vector<pair<double,double> >::const_iterator r = yRanges().begin();
	r != yRanges().end(); ++r ) {
	CurrentGenerator::log() << "y = " << r->first << " .. " << r->second << "\n";
	}

	}

	+double step(double r) {
	+ if ( r < -0.5 )
	+ return 0.0;
	+ if ( r > 0.5 )
	+ return 1.0;
	+ return r+0.5;
	+}
	+
	+bool JetRegion::lessThanEnergy(Energy a, Energy b, double& weight) const {
	+ if ( !fuzzy() ) {
	+ if ( a < b ) {
	+ weight = 1.0;
	+ return true;
	+ }
	+ weight = 0.0;
	+ return false;
	+ }
	+ double w = step((b-a)/theEnergyCutWidth);
	+ if ( w == 0.0 ) {
	+ weight = 0.0;
	+ return false;
	+ }
	+ weight *= w;
	+ return true;
	+}
	+
	+bool JetRegion::lessThanRapidity(double a, double b, double& weight) const {
	+ if ( !fuzzy() ) {
	+ if ( a < b ) {
	+ weight = 1.0;
	+ return true;
	+ }
	+ weight = 0.0;
	+ return false;
	+ }
	+ double w = step((b-a)/theRapidityCutWidth);
	+ if ( w == 0.0 ) {
	+ weight = 0.0;
	+ return false;
	+ }
	+ weight *= w;
	+ return true;
	+}
	+
	bool JetRegion::matches(tcCutsPtr parent, int n, const LorentzMomentum& p) {

	// one jet region can only contain one jet
	if ( didMatch() )
	return false;

	if ( !accepts().empty() && find(accepts().begin(),accepts().end(),n) == accepts().end() )
	return false;

	- if ( p.perp() < ptMin() \|\| p.perp() > ptMax() )
	+ theCutWeight = 1.0;
	+
	+ if ( !(lessThanEnergy(ptMin(),p.perp(),theCutWeight) &&
	+ lessThanEnergy(p.perp(),ptMax(),theCutWeight)) )
	return false;

	bool inRange = false \|\| yRanges().empty();
	for ( vector<pair<double,double> >::const_iterator r = yRanges().begin();
	r != yRanges().end(); ++r ) {
	- if ( p.rapidity() + parent->currentYHat() > r->first && p.rapidity() + parent->currentYHat() < r->second ) {
	+ double rangeWeight = 1.0;
	+ if ( lessThanRapidity(r->first,p.rapidity() + parent->currentYHat(),rangeWeight) &&
	+ lessThanRapidity(p.rapidity() + parent->currentYHat(),r->second,rangeWeight) ) {
	+ theCutWeight *= rangeWeight;
	inRange = true;
	break;
	}
	}
	- if ( !inRange )
	+ if ( !inRange ) {
	+ theCutWeight = 0.0;
	return false;
	+ }

	theDidMatch = true;
	theLastNumber = n;
	theLastMomentum = p;

	return true;

	}

	// If needed, insert default implementations of virtual function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void JetRegion::persistentOutput(PersistentOStream & os) const {
	os << ounit(thePtMin,GeV) << ounit(thePtMax,GeV)
	- << theYRanges << theAccepts;
	+ << theYRanges << theAccepts << theFuzzy << theCutWeight
	+ << ounit(theEnergyCutWidth,GeV) << theRapidityCutWidth;
	}

	void JetRegion::persistentInput(PersistentIStream & is, int) {
	is >> iunit(thePtMin,GeV) >> iunit(thePtMax,GeV)
	- >> theYRanges >> theAccepts;
	+ >> theYRanges >> theAccepts >> theFuzzy >> theCutWeight
	+ >> iunit(theEnergyCutWidth,GeV) >> theRapidityCutWidth;
	}


	// * Attention * The following static variable is needed for the type
	// description system in ThePEG. Please check that the template arguments
	// are correct (the class and its base class), and that the constructor
	// arguments are correct (the class name and the name of the dynamically
	// loadable library where the class implementation can be found).
	DescribeClass<JetRegion,HandlerBase>
	describeThePEGJetRegion("ThePEG::JetRegion", "JetCuts.so");

	void JetRegion::Init() {

	static ClassDocumentation<JetRegion> documentation
	("JetRegion implements the requirement of finding a jet inside a "
	"given range of transverse momenta, and (pseudo-)rapidity.");

	static Parameter<JetRegion,Energy> interfacePtMin
	("PtMin",
	"The minimum pt required.",
	&JetRegion::thePtMin, GeV, 0.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<JetRegion,Energy> interfacePtMax
	("PtMax",
	"The maximum pt allowed.",
	&JetRegion::thePtMax, GeV, Constants::MaxEnergy, 0.0GeV, 0GeV,
	false, false, Interface::lowerlim);

	static Command<JetRegion> interfaceYRange
	("YRange",
	"Insert a rapidity range.",
	&JetRegion::doYRange, false);

	static ParVector<JetRegion,int> interfaceAccepts
	("Accepts",
	"The jet numbers accepted. If empty, any jets are accepted.",
	&JetRegion::theAccepts, -1, 1, 1, 10,
	false, false, Interface::upperlim);

	+ static Switch<JetRegion,bool> interfaceFuzzy
	+ ("Fuzzy",
	+ "Make this jet region a fuzzy cut",
	+ &JetRegion::theFuzzy, false, false, false);
	+ static SwitchOption interfaceFuzzyOn
	+ (interfaceFuzzy,
	+ "Yes",
	+ "",
	+ true);
	+ static SwitchOption interfaceFuzzyOff
	+ (interfaceFuzzy,
	+ "No",
	+ "",
	+ false);
	+
	+ static Parameter<JetRegion,Energy> interfaceEnergyCutWidth
	+ ("EnergyCutWidth",
	+ "The pt cut smearing.",
	+ &JetRegion::theEnergyCutWidth, GeV, 1.0GeV, 0.0GeV, 0*GeV,
	+ false, false, Interface::lowerlim);
	+
	+ static Parameter<JetRegion,double> interfaceRapidityCutWidth
	+ ("RapidityCutWidth",
	+ "The rapidity cut smearing.",
	+ &JetRegion::theRapidityCutWidth, 0.1, 0.0, 0,
	+ false, false, Interface::lowerlim);
	+
	}

	diff --git a/Cuts/JetRegion.h b/Cuts/JetRegion.h
	--- a/Cuts/JetRegion.h
	+++ b/Cuts/JetRegion.h
	@@ -1,207 +1,247 @@
	// -- C++ --
	//
	// JetRegion.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 Leif Lonnblad
	// Copyright (C) 2009-2012 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef ThePEG_JetRegion_H
	#define ThePEG_JetRegion_H
	//
	// This is the declaration of the JetRegion class.
	//

	#include "ThePEG/Handlers/HandlerBase.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Cuts/Cuts.h"

	namespace ThePEG {

	/**
	* JetRegion implements the requirement of finding a jet inside a
	* given range of transverse momenta, and (pseudo-)rapidity.
	*
	* @see JetPairRegion
	* @see JetCuts
	*
	* @see \ref JetRegionInterfaces "The interfaces"
	* defined for JetRegion.
	*/
	class JetRegion: public HandlerBase {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	JetRegion();

	/**
	* The destructor.
	*/
	virtual ~JetRegion();
	//@}

	public:

	/**
	* Return the minimum pt.
	*/
	Energy ptMin() const { return thePtMin; }

	/**
	* Return the maximum pt.
	*/
	Energy ptMax() const { return thePtMax; }

	/**
	* Return the rapidity ranges.
	*/
	const vector<pair<double,double> >& yRanges() const { return theYRanges; }

	/**
	* Return the jets accepted by this region (with respect to the
	* ordering imposed by the JetCuts object). If empty, any jet will
	* be accepted.
	*/
	const vector<int>& accepts() const { return theAccepts; }

	+ /**
	+ * Return true, if this jet region is fuzzy
	+ */
	+ bool fuzzy() const { return theFuzzy; }
	+
	+ /**
	+ * Return the cut weight encountered from the last call to matches()
	+ */
	+ double cutWeight() const { return theCutWeight; }
	+
	+ /**
	+ * Perform a (potentially) fuzzy check on energy-type quantities
	+ */
	+ bool lessThanEnergy(Energy a, Energy b, double& weight) const;
	+
	+ /**
	+ * Perform a (potentially) fuzzy check on angular-type quantities
	+ */
	+ bool lessThanRapidity(double a, double b, double& weight) const;
	+
	public:

	/**
	* Describe the currently active cuts in the log file.
	*/
	virtual void describe() const;

	/**
	* Return true, if the given jet matches this region.
	*/
	virtual bool matches(tcCutsPtr parent, int n, const LorentzMomentum& p);

	/**
	* Return true, if this region matched a jet in the last call to matches().
	*/
	bool didMatch() { return theDidMatch; }

	/**
	* Reset this region to act on a new event.
	*/
	virtual void reset() { theDidMatch = false; }

	/**
	* Return the number of the last jet matching this region.
	*/
	int lastNumber() const { return theLastNumber; }

	/**
	* Return the momentum of the last jet matching this region.
	*/
	const LorentzMomentum& lastMomentum() const { return theLastMomentum; }

	public:

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

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

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

	protected:

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

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


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


	private:

	/**
	* Command to insert a rapidity range
	*/
	string doYRange(string);

	/**
	* The minimum pt.
	*/
	Energy thePtMin;

	/**
	* The maximum pt.
	*/
	Energy thePtMax;

	/**
	* The rapidity ranges.
	*/
	vector<pair<double,double> > theYRanges;

	/**
	* The jets accepted by this region (with respect to the ordering
	* imposed by the JetCuts object). If empty, any jet will be
	* accepted.
	*/
	vector<int> theAccepts;

	/**
	* True, if this region matched a jet in the last call to matches().
	*/
	bool theDidMatch;

	/**
	* The number of the last jet matching this region.
	*/
	int theLastNumber;

	/**
	* Return the momentum of the last jet matching this region.
	*/
	LorentzMomentum theLastMomentum;

	/**
	+ * True if this region is fuzzy
	+ */
	+ bool theFuzzy;
	+
	+ /**
	+ * The cut weight encountered from the last call to matches()
	+ */
	+ double theCutWeight;
	+
	+ /**
	+ * The smearing width for the pt or mass cuts, if fuzzy
	+ */
	+ Energy theEnergyCutWidth;
	+
	+ /**
	+ * The smearing width for the rapidity cut, if fuzzy
	+ */
	+ double theRapidityCutWidth;
	+
	+ /**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	JetRegion & operator=(const JetRegion &);

	};

	}

	#endif /* ThePEG_JetRegion_H */
	diff --git a/Cuts/MultiJetRegion.cc b/Cuts/MultiJetRegion.cc
	--- a/Cuts/MultiJetRegion.cc
	+++ b/Cuts/MultiJetRegion.cc
	@@ -1,182 +1,190 @@
	// -- C++ --
	//
	// MultiJetRegion.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 Leif Lonnblad
	// Copyright (C) 2009-2012 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the MultiJetRegion class.
	//

	#include "MultiJetRegion.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Repository/CurrentGenerator.h"

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

	using namespace ThePEG;

	MultiJetRegion::MultiJetRegion()
	: theMassMin(0.*GeV), theMassMax(Constants::MaxEnergy),
	theDeltaRMin(0.0), theDeltaRMax(Constants::MaxRapidity),
	- theDeltaYMin(0.0), theDeltaYMax(Constants::MaxRapidity) {}
	+ theDeltaYMin(0.0), theDeltaYMax(Constants::MaxRapidity),
	+ theCutWeight(1.0) {}

	MultiJetRegion::~MultiJetRegion() {}

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

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

	void MultiJetRegion::describe() const {

	CurrentGenerator::log()
	<< "MultiJetRegion '" << name() << "' matching JetRegions:\n";

	for ( vector<Ptr<JetRegion>::ptr>::const_iterator r = regions().begin();
	r != regions().end(); ++r ) {
	CurrentGenerator::log()
	<< "'" << (**r).name() << "'\n";
	}

	CurrentGenerator::log() << "with\n";

	CurrentGenerator::log()
	<< "m = " << massMin()/GeV << " .. " << massMax()/GeV << " GeV\n"
	<< "dR = " << deltaRMin() << " .. " << deltaRMax() << "\n"
	<< "dy = " << deltaYMin() << " .. " << deltaYMax() << "\n";

	}

	-bool MultiJetRegion::matches() const {
	+bool MultiJetRegion::matches() {

	int n = regions().size();

	for ( int i = 0; i < n; ++i )
	for ( int j = i+1; j < n; ++j )
	if ( !matches(i,j) )
	return false;

	return true;

	}

	-bool MultiJetRegion::matches(int i, int j) const {
	+bool MultiJetRegion::matches(int i, int j) {

	if ( !regions()[i]->didMatch() \|\|
	- !regions()[j]->didMatch() )
	+ !regions()[j]->didMatch() ) {
	+ theCutWeight = 0.0;
	return false;
	+ }
	+
	+ theCutWeight = 1.0;

	const LorentzMomentum& pi = regions()[i]->lastMomentum();
	const LorentzMomentum& pj = regions()[j]->lastMomentum();

	Energy m = (pi+pj).m();
	- if ( m < massMin() \|\| m > massMax() )
	+ if ( !(regions()[i]->lessThanEnergy(massMin(),m,theCutWeight) &&
	+ regions()[i]->lessThanEnergy(m,massMax(),theCutWeight)) )
	return false;

	double dy = abs(pi.rapidity() - pj.rapidity());
	- if ( dy < deltaYMin() \|\| dy > deltaYMax() )
	+ if ( !(regions()[i]->lessThanRapidity(deltaYMin(),dy,theCutWeight) &&
	+ regions()[i]->lessThanRapidity(dy,deltaYMax(),theCutWeight)) )
	return false;

	double dphi = abs(pi.phi() - pj.phi());
	if ( dphi > Constants::pi ) dphi = 2.0*Constants::pi - dphi;
	double dR = sqrt(sqr(dy)+sqr(dphi));
	- if ( dR < deltaRMin() \|\| dR > deltaRMax() )
	+ if ( !(regions()[i]->lessThanRapidity(deltaRMin(),dR,theCutWeight) &&
	+ regions()[i]->lessThanRapidity(dR,deltaRMax(),theCutWeight)) )
	return false;

	return true;

	}

	// If needed, insert default implementations of virtual function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void MultiJetRegion::persistentOutput(PersistentOStream & os) const {
	os << theRegions
	<< ounit(theMassMin,GeV) << ounit(theMassMax,GeV)
	<< theDeltaRMin << theDeltaRMax
	- << theDeltaYMin << theDeltaYMax;
	+ << theDeltaYMin << theDeltaYMax << theCutWeight;
	}

	void MultiJetRegion::persistentInput(PersistentIStream & is, int) {
	is >> theRegions
	>> iunit(theMassMin,GeV) >> iunit(theMassMax,GeV)
	>> theDeltaRMin >> theDeltaRMax
	- >> theDeltaYMin >> theDeltaYMax;
	+ >> theDeltaYMin >> theDeltaYMax >> theCutWeight;
	}


	// * Attention * The following static variable is needed for the type
	// description system in ThePEG. Please check that the template arguments
	// are correct (the class and its base class), and that the constructor
	// arguments are correct (the class name and the name of the dynamically
	// loadable library where the class implementation can be found).
	DescribeClass<MultiJetRegion,HandlerBase>
	describeThePEGMultiJetRegion("ThePEG::MultiJetRegion", "JetCuts.so");

	void MultiJetRegion::Init() {

	static ClassDocumentation<MultiJetRegion> documentation
	("MultiJetRegion implements pairwise constraints on jets matching several jet regions.");


	static RefVector<MultiJetRegion,JetRegion> interfaceRegions
	("Regions",
	"The jet regions to act on.",
	&MultiJetRegion::theRegions, -1, false, false, true, false, false);

	static Parameter<MultiJetRegion,Energy> interfaceMassMin
	("MassMin",
	"The minimum jet-jet invariant mass.",
	&MultiJetRegion::theMassMin, GeV, 0.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<MultiJetRegion,Energy> interfaceMassMax
	("MassMax",
	"The maximum jet-jet invariant mass.",
	&MultiJetRegion::theMassMax, GeV, Constants::MaxEnergy, 0.0GeV, 0GeV,
	false, false, Interface::lowerlim);

	static Parameter<MultiJetRegion,double> interfaceDeltaRMin
	("DeltaRMin",
	"The minimum jet-jet lego-plot separation.",
	&MultiJetRegion::theDeltaRMin, 0.0, 0.0, 0,
	false, false, Interface::lowerlim);

	static Parameter<MultiJetRegion,double> interfaceDeltaRMax
	("DeltaRMax",
	"The maximum jet-jet lego-plot separation.",
	&MultiJetRegion::theDeltaRMax, Constants::MaxRapidity, 0, 0,
	false, false, Interface::lowerlim);

	static Parameter<MultiJetRegion,double> interfaceDeltaYMin
	("DeltaYMin",
	"The minimum jet-jet rapidity separation.",
	&MultiJetRegion::theDeltaYMin, 0.0, 0.0, 0,
	false, false, Interface::lowerlim);

	static Parameter<MultiJetRegion,double> interfaceDeltaYMax
	("DeltaYMax",
	"The maximum jet-jet rapidity separation.",
	&MultiJetRegion::theDeltaYMax, Constants::MaxRapidity, 0, 0,
	false, false, Interface::lowerlim);

	}

	diff --git a/Cuts/MultiJetRegion.h b/Cuts/MultiJetRegion.h
	--- a/Cuts/MultiJetRegion.h
	+++ b/Cuts/MultiJetRegion.h
	@@ -1,195 +1,205 @@
	// -- C++ --
	//
	// MultiJetRegion.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 Leif Lonnblad
	// Copyright (C) 2009-2012 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef ThePEG_MultiJetRegion_H
	#define ThePEG_MultiJetRegion_H
	//
	// This is the declaration of the MultiJetRegion class.
	//

	#include "ThePEG/Cuts/JetRegion.h"

	namespace ThePEG {

	/**
	* MultiJetRegion implements pairwise constraints on jets matching several jet regions.
	*
	* @see JetRegion
	* @see JetCuts
	*
	* @see \ref MultiJetRegionInterfaces "The interfaces"
	* defined for MultiJetRegion.
	*/
	class MultiJetRegion: public HandlerBase {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	MultiJetRegion();

	/**
	* The destructor.
	*/
	virtual ~MultiJetRegion();
	//@}

	public:

	/**
	* Return the jet regions to act on.
	*/
	const vector<Ptr<JetRegion>::ptr>& regions() const { return theRegions; }

	/**
	* Return the minimum jet-jet invariant mass.
	*/
	Energy massMin() const { return theMassMin; }

	/**
	* Return the maximum jet-jet invariant mass.
	*/
	Energy massMax() const { return theMassMax; }

	/**
	* Return the minimum jet-jet lego-plot separation.
	*/
	double deltaRMin() const { return theDeltaRMin; }

	/**
	* Return the maximum jet-jet lego-plot separation.
	*/
	double deltaRMax() const { return theDeltaRMax; }

	/**
	* Return the minimum jet-jet rapidity separation.
	*/
	double deltaYMin() const { return theDeltaYMin; }

	/**
	* Return the maximum jet-jet rapidity separation.
	*/
	double deltaYMax() const { return theDeltaYMax; }

	+ /**
	+ * Return the cut weight encountered from the last call to matches()
	+ */
	+ double cutWeight() const { return theCutWeight; }
	+
	public:

	/**
	* Describe the currently active cuts in the log file.
	*/
	virtual void describe() const;

	/**
	* Return true, if the requirements on the jet regions are fullfilled.
	*/
	- virtual bool matches() const;
	+ virtual bool matches();

	/**
	* Return true, if the requirements on two jet regions are fullfilled.
	*/
	- virtual bool matches(int i, int j) const;
	+ virtual bool matches(int i, int j);

	public:

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

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

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

	protected:

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

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


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


	private:

	/**
	* The jet regions to act on.
	*/
	vector<Ptr<JetRegion>::ptr> theRegions;

	/**
	* The minimum jet-jet invariant mass.
	*/
	Energy theMassMin;

	/**
	* The maximum jet-jet invariant mass.
	*/
	Energy theMassMax;

	/**
	* The minimum jet-jet lego-plot separation.
	*/
	double theDeltaRMin;

	/**
	* The maximum jet-jet lego-plot separation.
	*/
	double theDeltaRMax;

	/**
	* The minimum jet-jet rapidity separation.
	*/
	double theDeltaYMin;

	/**
	* The maximum jet-jet rapidity separation.
	*/
	double theDeltaYMax;

	/**
	+ * The cut weight encountered from the last call to matches()
	+ */
	+ double theCutWeight;
	+
	+ /**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	MultiJetRegion & operator=(const MultiJetRegion &);

	};

	}

	#endif /* ThePEG_MultiJetRegion_H */
	diff --git a/Handlers/StandardXComb.cc b/Handlers/StandardXComb.cc
	--- a/Handlers/StandardXComb.cc
	+++ b/Handlers/StandardXComb.cc
	@@ -1,677 +1,691 @@
	// -- C++ --
	//
	// StandardXComb.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	// Copyright (C) 2009-2011 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the StandardXComb class.
	//

	#include "StandardXComb.h"
	#include "StdXCombGroup.h"
	#include "ThePEG/Handlers/StandardEventHandler.h"
	#include "ThePEG/Handlers/SubProcessHandler.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "ThePEG/Utilities/Maths.h"
	#include "ThePEG/PDT/ParticleData.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Utilities/SimplePhaseSpace.h"
	#include "ThePEG/Utilities/UtilityBase.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/EventRecord/SubProcessGroup.h"
	#include "ThePEG/Vectors/LorentzRotation.h"
	#include "ThePEG/MatrixElement/MEBase.h"
	#include "ThePEG/MatrixElement/ColourLines.h"
	#include "ThePEG/Handlers/LuminosityFunction.h"
	#include "ThePEG/Handlers/CascadeHandler.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/EventRecord/TmpTransform.h"

	#ifdef ThePEG_TEMPLATES_IN_CC_FILE
	#include "StandardXComb.tcc"
	#endif

	using namespace ThePEG;

	StandardXComb::StandardXComb()
	: XComb(), isMirror(false), theNDim(0),
	partonDims(make_pair(0, 0)), theKinematicsGenerated(false),
	theLastDiagramIndex(0), theLastPDFWeight(0.0),
	theLastCrossSection(ZERO), theLastJacobian(1.0), theLastME2(-1.0),
	theLastMECrossSection(ZERO), theLastMEPDFWeight(1.0),
	- checkedCuts(false), passedCuts(false) {}
	+ checkedCuts(false), passedCuts(false), theCutWeight(1.0) {}

	StandardXComb::
	StandardXComb(Energy newMaxEnergy, const cPDPair & inc,
	tEHPtr newEventHandler, tSubHdlPtr newSubProcessHandler,
	tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
	const PBPair & newPartonBins, tCutsPtr newCuts,
	tMEPtr newME, const DiagramVector & newDiagrams, bool mir,
	tStdXCombPtr newHead)
	: XComb(newMaxEnergy, inc, newEventHandler,
	newExtractor, newCKKW, newPartonBins, newCuts),
	theSubProcessHandler(newSubProcessHandler), theME(newME),
	theDiagrams(newDiagrams), isMirror(mir), theKinematicsGenerated(false),
	theLastDiagramIndex(0), theLastPDFWeight(0.0),
	theLastCrossSection(ZERO), theLastME2(-1.0), theLastMECrossSection(ZERO),
	theLastMEPDFWeight(1.0), theHead(newHead),
	checkedCuts(false), passedCuts(false) {
	partonDims = pExtractor()->nDims(partonBins());
	if ( matrixElement()->haveX1X2() ) {
	partonDims.first = 0;
	partonDims.second = 0;
	}
	theNDim = matrixElement()->nDim() + partonDims.first + partonDims.second;
	mePartonData() = lastDiagram()->partons();
	}

	StandardXComb::StandardXComb(tMEPtr me, const tPVector & parts,
	DiagramIndex indx)
	: theME(me), isMirror(false), theNDim(0), partonDims(make_pair(0, 0)),
	theKinematicsGenerated(false),
	theLastDiagramIndex(0), theLastPDFWeight(0.0), theLastCrossSection(ZERO),
	theLastME2(-1.0), theLastMECrossSection(ZERO), theLastMEPDFWeight(1.0),
	checkedCuts(false), passedCuts(false) {

	subProcess(new_ptr(SubProcess(make_pair(parts[0], parts[1]),
	tCollPtr(), me)));
	for ( int i = 0, N = parts.size(); i < N; ++i ) {
	subProcess()->addOutgoing(parts[i], false);
	theMEPartonData.push_back(parts[i]->dataPtr());
	theMEMomenta.push_back(parts[i]->momentum());
	}
	lastSHat((meMomenta()[0] + meMomenta()[1]).m2());
	string tag = me->diagrams()[indx]->getTag();
	for ( int i = 0, N = me->diagrams().size(); i < N; ++i )
	if ( me->diagrams()[i]->getTag() == tag )
	theDiagrams.push_back(me->diagrams()[i]);
	}

	StandardXComb::StandardXComb(tStdXCombPtr newHead,
	const PBPair & newPartonBins, tMEPtr newME,
	const DiagramVector & newDiagrams)
	: XComb(newHead->maxEnergy(), newHead->particles(),
	newHead->eventHandlerPtr(), newHead->pExtractor(),
	newHead->CKKWHandler(), newPartonBins, newHead->cuts()),
	theSubProcessHandler(const_ptr_cast<tSubHdlPtr>(newHead->subProcessHandler())),
	theME(newME), theDiagrams(newDiagrams), isMirror(newHead->mirror()),
	theKinematicsGenerated(false), theLastDiagramIndex(0), theLastPDFWeight(0.0),
	theLastCrossSection(ZERO), theLastME2(-1.0), theLastMECrossSection(ZERO),
	theLastMEPDFWeight(1.0), theHead(newHead),
	checkedCuts(false), passedCuts(false) {
	partonDims = pExtractor()->nDims(partonBins());
	if ( matrixElement()->haveX1X2() ) {
	partonDims.first = 0;
	partonDims.second = 0;
	}
	theNDim = matrixElement()->nDim() + partonDims.first + partonDims.second;
	mePartonData() = lastDiagram()->partons();
	/* // wait for C++11
	StandardXComb(newHead->maxEnergy(),newHead->particles(),
	newHead->eventHandlerPtr(),
	const_ptr_cast<tSubHdlPtr>(newHead->subProcessHandler()),
	newHead->pExtractor(),newHead->CKKWHandler(),
	newPartonBins,newHead->cuts(),newME,newDiagrams,newHead->mirror(),
	newHead);
	*/
	}

	StandardXComb::~StandardXComb() {}

	void StandardXComb::recreatePartonBinInstances(Energy2 scale) {

	PBIPair newBins;

	Direction<0> dir(true);
	newBins.first =
	new_ptr(PartonBinInstance(lastPartons().first,partonBins().first,scale));

	dir.reverse();
	newBins.second =
	new_ptr(PartonBinInstance(lastPartons().second,partonBins().second,scale));

	resetPartonBinInstances(newBins);
	setPartonBinInfo();

	lastPartons().first->scale(partonBinInstances().first->scale());
	lastPartons().second->scale(partonBinInstances().second->scale());

	}

	void StandardXComb::refillPartonBinInstances(const double* r) {

	pExtractor()->select(this);
	pExtractor()->updatePartonBinInstances(partonBinInstances());
	pExtractor()->generateSHat(lastS(), partonBinInstances(),
	r, r + nDim() - partonDims.second,true);

	}

	void StandardXComb::setIncomingPartons(tStdXCombPtr labHead) {

	if ( lastPartons().first )
	return;

	if ( !labHead )
	labHead = head();

	createPartonBinInstances();
	lastParticles(labHead->lastParticles());
	setPartonBinInfo();

	lastPartons(make_pair(mePartonData()[0]->produceParticle(Lorentz5Momentum()),
	mePartonData()[1]->produceParticle(Lorentz5Momentum())));

	Lorentz5Momentum pFirst = meMomenta()[0];
	Lorentz5Momentum pSecond = meMomenta()[1];

	if ( labHead->matrixElement()->wantCMS() ) {
	Boost toLab = (labHead->lastPartons().first->momentum() +
	labHead->lastPartons().second->momentum()).boostVector();
	if ( toLab.mag2() > Constants::epsilon ) {
	pFirst.boost(toLab);
	pSecond.boost(toLab);
	}
	}

	lastPartons().first->set5Momentum(pFirst);
	lastPartons().second->set5Momentum(pSecond);

	partonBinInstances().first->parton(lastPartons().first);
	partonBinInstances().second->parton(lastPartons().second);

	lastS((lastParticles().first->momentum() +
	lastParticles().second->momentum()).m2());
	lastSHat((lastPartons().first->momentum() +
	lastPartons().second->momentum()).m2());
	lastP1P2(make_pair(labHead->lastP1(),labHead->lastP2()));

	double x1 =
	lastPartons().first->momentum().plus()/
	lastParticles().first->momentum().plus();
	double x2 =
	lastPartons().second->momentum().minus()/
	lastParticles().second->momentum().minus();

	lastX1X2(make_pair(x1,x2));

	lastY((lastPartons().first->momentum()+
	lastPartons().second->momentum()).rapidity());

	}

	void StandardXComb::fill(const PPair& newParticles,
	const PPair& newPartons,
	const vector<Lorentz5Momentum>& newMEMomenta,
	const DVector& newLastRandomNumbers) {

	lastParticles(newParticles);

	lastP1P2(make_pair(0.0, 0.0));
	lastS((lastParticles().first->momentum() +
	lastParticles().second->momentum()).m2());

	lastPartons(newPartons);
	lastSHat((lastPartons().first->momentum() +
	lastPartons().second->momentum()).m2());
	lastX1X2(make_pair(lastPartons().first->momentum().plus()/
	lastParticles().first->momentum().plus(),
	lastPartons().second->momentum().minus()/
	lastParticles().second->momentum().minus()));
	lastY((lastPartons().first->momentum() +
	lastPartons().second->momentum()).rapidity());

	meMomenta().resize(newMEMomenta.size());
	copy(newMEMomenta.begin(),newMEMomenta.end(),
	meMomenta().begin());

	lastRandomNumbers().resize(newLastRandomNumbers.size());
	copy(newLastRandomNumbers.begin(),newLastRandomNumbers.end(),
	lastRandomNumbers().begin());

	}

	bool StandardXComb::checkInit() {
	Energy summin = ZERO;
	for ( int i = 2, N = mePartonData().size(); i < N; ++i ) {
	summin += mePartonData()[i]->massMin();
	}
	return ( summin < min(maxEnergy(), cuts()->mHatMax()) );
	}

	bool StandardXComb::willPassCuts() {

	if ( checkedCuts )
	return passedCuts;

	checkedCuts = true;

	if ( !head() ) {
	if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
	passedCuts = false;
	+ theCutWeight = 0.0;
	return false;
	}
	} else {
	cuts()->initSubProcess(lastSHat(), lastY(), mirror());
	}

	tcPDVector outdata(mePartonData().begin()+2,mePartonData().end());
	vector<LorentzMomentum> outmomenta(meMomenta().begin()+2,meMomenta().end());
	Boost tocm = (meMomenta()[0]+meMomenta()[1]).findBoostToCM();
	if ( tocm.mag2() > Constants::epsilon ) {
	for ( vector<LorentzMomentum>::iterator p = outmomenta.begin();
	p != outmomenta.end(); ++p ) {
	p->boost(tocm);
	}
	}

	if ( !cuts()->passCuts(outdata,outmomenta,mePartonData()[0],mePartonData()[1]) ) {
	+ theCutWeight = cuts()->cutWeight();
	passedCuts = false;
	return false;
	}

	+ theCutWeight = cuts()->cutWeight();
	passedCuts = true;
	return true;

	}

	void StandardXComb::clean() {
	XComb::clean();
	theLastPDFWeight = 0.0;
	theLastCrossSection = ZERO;
	theLastJacobian = 0.0;
	theLastME2 = 0.0;
	theLastMECrossSection = ZERO;
	theLastMEPDFWeight = 0.0;
	theProjectors.clear();
	theProjector = StdXCombPtr();
	theKinematicsGenerated = false;
	checkedCuts = false;
	passedCuts = false;
	+ theCutWeight = 1.0;
	matrixElement()->flushCaches();
	}

	CrossSection StandardXComb::dSigDR(const double * r) {

	matrixElement()->flushCaches();

	matrixElement()->setXComb(this);

	if ( !matrixElement()->apply() ) {
	subProcess(SubProPtr());
	lastCrossSection(ZERO);
	return ZERO;
	}

	meMomenta().resize(mePartonData().size());

	if ( !matrixElement()->generateKinematics(r) ) {
	subProcess(SubProPtr());
	lastCrossSection(ZERO);
	return ZERO;
	}

	setIncomingPartons();

	lastScale(matrixElement()->scale());
	lastAlphaS(matrixElement()->alphaS());
	lastAlphaEM(matrixElement()->alphaEM());

	if ( (!willPassCuts() &&
	!matrixElement()->headCuts() &&
	!matrixElement()->ignoreCuts()) \|\|
	!matrixElement()->apply() ) {
	subProcess(SubProPtr());
	lastCrossSection(ZERO);
	return ZERO;
	}

	lastPDFWeight(head()->lastPDFWeight());

	matrixElement()->setKinematics();
	+
	CrossSection xsec = matrixElement()->dSigHatDR() * lastPDFWeight();

	+ xsec *= cutWeight();
	+
	subProcess(SubProPtr());
	lastCrossSection(xsec);

	return xsec;

	}

	CrossSection StandardXComb::
	dSigDR(const pair<double,double> ll, int nr, const double * r) {

	matrixElement()->flushCaches();

	if ( matrixElement()->keepRandomNumbers() ) {
	lastRandomNumbers().resize(nDim());
	copy(r,r+nDim(),lastRandomNumbers().begin());
	}

	pExtractor()->select(this);
	setPartonBinInfo();
	lastP1P2(ll);
	lastS(sqr(maxEnergy())/exp(lastP1() + lastP2()));

	meMomenta().resize(mePartonData().size());
	matrixElement()->setXComb(this);

	PPair partons;

	if ( !matrixElement()->haveX1X2() ) {

	if ( !pExtractor()->generateL(partonBinInstances(),
	r, r + nr - partonDims.second) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	partons = make_pair(partonBinInstances().first->parton(),
	partonBinInstances().second->parton());
	lastSHat(lastS()/exp(partonBinInstances().first->l() +
	partonBinInstances().second->l()));
	meMomenta()[0] = partons.first->momentum();
	meMomenta()[1] = partons.second->momentum();

	} else {
	if ( !matrixElement()->generateKinematics(r + partonDims.first) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
	matrixElement()->setKinematics();

	lastScale(matrixElement()->scale());

	partons.first = mePartonData()[0]->produceParticle(meMomenta()[0]);
	partons.second = mePartonData()[1]->produceParticle(meMomenta()[1]);

	Direction<0> dir(true);
	partonBinInstances().first =
	new_ptr(PartonBinInstance(lastParticles().first,partons.first,
	partonBins().first,lastScale()));
	dir.reverse();
	partonBinInstances().second =
	new_ptr(PartonBinInstance(lastParticles().second,partons.second,
	partonBins().second,lastScale()));
	}

	lastPartons(partons);

	if ( lastSHat() < cuts()->sHatMin() ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	lastY(0.5*(partonBinInstances().second->l() -
	partonBinInstances().first->l()));
	if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]);
	if ( matrixElement()->wantCMS() &&
	!matrixElement()->haveX1X2() )
	SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat());

	Energy summ = ZERO;
	if ( meMomenta().size() == 3 ) {
	if ( !matrixElement()->haveX1X2() )
	meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat()));
	} else {
	for ( int i = 2, N = meMomenta().size(); i < N; ++i ) {
	if ( !matrixElement()->haveX1X2() )
	meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->mass());
	summ += mePartonData()[i]->massMin();
	}
	if ( sqr(summ) >= lastSHat() ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	}

	matrixElement()->setXComb(this);
	if ( !matrixElement()->haveX1X2() )
	lastScale(max(lastSHat()/4.0, cuts()->scaleMin()));

	lastSHat(pExtractor()->generateSHat(lastS(), partonBinInstances(),
	r, r + nr - partonDims.second,
	matrixElement()->haveX1X2()));

	if ( !cuts()->sHat(lastSHat()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	r += partonDims.first;

	lastX1X2(make_pair(lastPartons().first->momentum().plus()/
	lastParticles().first->momentum().plus(),
	lastPartons().second->momentum().minus()/
	lastParticles().second->momentum().minus()));

	if ( !cuts()->x1(lastX1()) \|\| !cuts()->x2(lastX2()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	lastY((lastPartons().first->momentum() +
	lastPartons().second->momentum()).rapidity());
	if ( !cuts()->yHat(lastY()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	+
	if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	meMomenta()[0] = lastPartons().first->momentum();
	meMomenta()[1] = lastPartons().second->momentum();
	if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]);
	if ( matrixElement()->wantCMS() &&
	!matrixElement()->haveX1X2() )
	SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat());

	if ( meMomenta().size() == 3 ) {
	if ( !matrixElement()->haveX1X2() )
	meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat()));
	} else {
	if ( sqr(summ) >= lastSHat() ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	}

	matrixElement()->setXComb(this);
	if ( !matrixElement()->haveX1X2() ) {
	if ( !matrixElement()->generateKinematics(r) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	}
	+
	lastScale(matrixElement()->scale());
	if ( !cuts()->scale(lastScale()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	+ // get information on cuts; we don't take this into account here for
	+ // reasons of backward compatibility but this will change eventually
	+ willPassCuts();
	+
	pair<bool,bool> evalPDFS =
	make_pair(matrixElement()->havePDFWeight1(),
	matrixElement()->havePDFWeight2());
	if ( mirror() )
	swap(evalPDFS.first,evalPDFS.second);
	lastPDFWeight(pExtractor()->fullFn(partonBinInstances(), lastScale(),
	evalPDFS));
	if ( lastPDFWeight() == 0.0 ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	matrixElement()->setKinematics();
	CrossSection xsec = matrixElement()->dSigHatDR() * lastPDFWeight();
	if ( xsec == ZERO ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	+ xsec *= cutWeight();
	+
	lastAlphaS (matrixElement()->orderInAlphaS () >0 ?
	matrixElement()->alphaS() : -1.);
	lastAlphaEM(matrixElement()->orderInAlphaEW() >0 ?
	matrixElement()->alphaEM() : -1.);

	matrixElement()->fillProjectors();
	if ( !projectors().empty() ) {
	lastProjector(projectors().select(UseRandom::rnd()));
	}

	subProcess(SubProPtr());
	if ( CKKWHandler() && matrixElement()->maxMultCKKW() > 0 &&
	matrixElement()->maxMultCKKW() > matrixElement()->minMultCKKW() ) {
	newSubProcess();
	CKKWHandler()->setXComb(this);
	xsec *= CKKWHandler()->reweightCKKW(matrixElement()->minMultCKKW(),
	matrixElement()->maxMultCKKW());
	}

	if ( matrixElement()->reweighted() ) {
	newSubProcess();
	xsec = matrixElement()->reWeight() matrixElement()->preWeight();
	}

	lastCrossSection(xsec);

	return xsec;

	}

	void StandardXComb::newSubProcess(bool group) {
	if ( subProcess() ) return;
	if ( head() ) {
	// first get the meMomenta in their CMS, as this may
	// not be the case
	Boost cm = (meMomenta()[0] + meMomenta()[1]).findBoostToCM();
	if ( cm.mag2() > Constants::epsilon ) {
	for ( vector<Lorentz5Momentum>::iterator m = meMomenta().begin();
	m != meMomenta().end(); ++m ) {
	*m = m->boost(cm);
	}
	}
	}
	if ( !lastProjector() ) {
	if ( !group )
	subProcess(new_ptr(SubProcess(lastPartons(), tCollPtr(), matrixElement())));
	else
	subProcess(new_ptr(SubProcessGroup(lastPartons(), tCollPtr(), matrixElement())));
	lastDiagramIndex(matrixElement()->diagram(diagrams()));
	const ColourLines & cl = matrixElement()->selectColourGeometry(lastDiagram());
	Lorentz5Momentum p1 = lastPartons().first->momentum();
	Lorentz5Momentum p2 = lastPartons().second->momentum();
	tPPair inc = lastPartons();
	if ( mirror() ) swap(inc.first, inc.second);
	if ( matrixElement()->wantCMS() &&
	!matrixElement()->haveX1X2() ) {
	LorentzRotation r = Utilities::boostToCM(inc);
	lastDiagram()->construct(subProcess(), *this, cl);
	subProcess()->transform(r.inverse());
	lastPartons().first->set5Momentum(p1);
	lastPartons().second->set5Momentum(p2);
	} else {
	lastDiagram()->construct(subProcess(), *this, cl);
	}
	lastPartons().first ->scale(partonBinInstances().first ->scale());
	lastPartons().second->scale(partonBinInstances().second->scale());
	for ( int i = 0, N = subProcess()->outgoing().size(); i < N; ++i )
	subProcess()->outgoing()[i]->scale(lastScale());
	// construct the spin information for the interaction
	matrixElement()->constructVertex(subProcess());
	// set veto scales
	matrixElement()->setVetoScales(subProcess());
	} else {
	lastProjector()->newSubProcess();
	subProcess(lastProjector()->subProcess());
	lastPartons(lastProjector()->lastPartons());
	lastSHat((lastPartons().first->momentum() +
	lastPartons().second->momentum()).m2());
	lastX1X2(make_pair(lastPartons().first->momentum().plus()/
	lastParticles().first->momentum().plus(),
	lastPartons().second->momentum().minus()/
	lastParticles().second->momentum().minus()));
	lastY(log(lastX1()/lastX2())*0.5);
	partonBinInstances().first->parton(lastPartons().first);
	partonBinInstances().second->parton(lastPartons().second);
	if ( !matrixElement()->keepRandomNumbers() )
	throw Exception() << "Matrix element needs to request random number storage "
	<< "for creating projected subprocesses."
	<< Exception::runerror;
	const double * r = &lastRandomNumbers()[0];
	pExtractor()->generateSHat(lastS(), partonBinInstances(),
	r, r + nDim() - partonDims.second,true);
	pExtractor()->updatePartonBinInstances(partonBinInstances());
	}
	}

	tSubProPtr StandardXComb::construct() {

	matrixElement()->setXComb(this);

	if ( !head() ) {
	if ( !cuts()->initSubProcess(lastSHat(), lastY()) ) throw Veto();
	} else {
	cuts()->initSubProcess(lastSHat(), lastY());
	}

	-
	if ( head() )
	if ( !matrixElement()->apply() )
	return tSubProPtr();

	setPartonBinInfo();
	matrixElement()->setKinematics();

	newSubProcess();

	TmpTransform<tSubProPtr>
	tmp(subProcess(), Utilities::getBoostToCM(subProcess()->incoming()));
	if ( !cuts()->passCuts(*subProcess()) ) throw Veto();

	if ( head() ) {
	subProcess()->head(head()->subProcess());
	subProcess()->groupWeight(lastCrossSection()/head()->lastCrossSection());
	}

	return subProcess();

	}

	void StandardXComb::Init() {}

	void StandardXComb::persistentOutput(PersistentOStream & os) const {
	os << theSubProcessHandler << theME << theStats
	<< theDiagrams << isMirror << theNDim << partonDims
	<< theLastDiagramIndex << theMEInfo << theLastRandomNumbers << theMEPartonData
	<< theLastPDFWeight << ounit(theLastCrossSection,nanobarn) << theLastJacobian
	<< theLastME2 << ounit(theLastMECrossSection,nanobarn) << theLastMEPDFWeight
	<< theHead << theProjectors << theProjector << theKinematicsGenerated
	- << checkedCuts << passedCuts;
	+ << checkedCuts << passedCuts << theCutWeight;
	}

	void StandardXComb::persistentInput(PersistentIStream & is, int) {
	is >> theSubProcessHandler >> theME >> theStats
	>> theDiagrams >> isMirror >> theNDim >> partonDims
	>> theLastDiagramIndex >> theMEInfo >> theLastRandomNumbers >> theMEPartonData
	>> theLastPDFWeight >> iunit(theLastCrossSection,nanobarn) >> theLastJacobian
	>> theLastME2 >> iunit(theLastMECrossSection,nanobarn) >> theLastMEPDFWeight
	>> theHead >> theProjectors >> theProjector >> theKinematicsGenerated
	- >> checkedCuts >> passedCuts;
	+ >> checkedCuts >> passedCuts >> theCutWeight;
	}

	ClassDescription<StandardXComb> StandardXComb::initStandardXComb;

	diff --git a/Handlers/StandardXComb.h b/Handlers/StandardXComb.h
	--- a/Handlers/StandardXComb.h
	+++ b/Handlers/StandardXComb.h
	@@ -1,671 +1,681 @@
	// -- C++ --
	//
	// StandardXComb.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	// Copyright (C) 2009-2011 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef ThePEG_StandardXComb_H
	#define ThePEG_StandardXComb_H
	// This is the declaration of the StandardXComb class.

	#include "ThePEG/Config/ThePEG.h"
	#include "SubProcessHandler.fh"
	#include "ThePEG/PDF/PartonExtractor.fh"
	#include "ThePEG/PDF/PartonBin.h"
	#include "ThePEG/PDF/PartonBinInstance.h"
	#include "ThePEG/Utilities/VSelector.h"
	#include "ThePEG/Utilities/ClassDescription.h"
	#include "ThePEG/Utilities/Maths.h"
	#include "ThePEG/Utilities/XSecStat.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/MatrixElement/MEBase.h"
	#include "ThePEG/Handlers/XComb.h"
	#include "ThePEG/Handlers/StandardEventHandler.h"
	#include "ThePEG/Handlers/SubProcessHandler.fh"
	#include "StandardXComb.fh"

	namespace ThePEG {

	/**
	* The StandardXComb class inherits from the more general XComb class
	* which stores all information about the generation of a hard
	* sub-proces for a given pair of incoming particles, a pair of
	* extracted partons, etc. This class stores more information related
	* to thestandard process generation scheme in ThePEG, such as the
	* PartonExtractor and MEBase object used. It also does some of the
	* administration of the process generation.
	*
	* The main function is dSigDR() which returns the differential cross
	* section w.r.t. a given vector of random numbers in the interval
	* ]0,1[. In the initialization this is used to pre-sample the phase
	* space. In the generation phase it is used to give the cross section
	* for a phase space point, and if this StandardXComb is chosen the
	* construct() function is called to generate the actual sub-process.
	*
	* @see ParonExtractor
	* @see MEBase
	* @see Cuts
	* @see StdXCombGroup
	*/
	class StandardXComb: public XComb {

	public:

	/** A vector of DiagramBase objects. */
	typedef MEBase::DiagramVector DiagramVector;

	/** A vector of indices. */
	typedef MEBase::DiagramIndex DiagramIndex;

	/** MEBase needs to be a friend. */
	friend class MEBase;

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Standard constructor.
	*/
	StandardXComb(Energy newMaxEnergy, const cPDPair & inc,
	tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
	tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
	const PBPair & newPartonBins, tCutsPtr newCuts, tMEPtr newME,
	const DiagramVector & newDiagrams, bool mir,
	tStdXCombPtr newHead = tStdXCombPtr());

	/**
	* Constructor given a head xcomb.
	*/
	StandardXComb(tStdXCombPtr newHead,
	const PBPair & newPartonBins, tMEPtr newME,
	const DiagramVector & newDiagrams);

	/**
	* Default constructor.
	*/
	StandardXComb();

	/**
	* Destructor.
	*/
	virtual ~StandardXComb();

	/**
	* Constructor used by MEBase to create a temporary object to store info.
	*/
	StandardXComb(tMEPtr me, const tPVector & parts, DiagramIndex i);

	//@}

	/** @name Utilities for incoming partons. */
	//@{

	/**
	* Properly setup the PartonBinInstance objects provided a sub
	* process has been constructed using this XComb.
	*/
	void recreatePartonBinInstances(Energy2 scale);

	/**
	* Fill the variables needed to generate remnants; momenta will be
	* used from the partons set in this xcomb, but random numbers need
	* to be provided to (re)generate variables not fixed by the
	* incoming partons.
	*/
	void refillPartonBinInstances(const double* r);

	/**
	* Setup information on incoming partons depending
	* on the information previously supplied through the
	* choice of diagram and incoming momenta in the first
	* two entries of meMomenta(). Partons are not actually
	* extracted from the incoming particles, though a subprocess
	* detached from the current Event may be created.
	*/
	void setIncomingPartons(tStdXCombPtr labHead = tStdXCombPtr());

	/**
	* Fill phase space information as far as possible
	*/
	void fill(const PPair& newParticles,
	const PPair& newPartons,
	const vector<Lorentz5Momentum>& newMEMomenta,
	const DVector& newLastRandomNumbers = DVector());

	//@}

	/** @name Access the assigned objects used in the generation. */
	//@{
	/**
	* Return a pointer to the corresponding sub-process handler. May be
	* null if the standard process generation in ThePEG was not used.
	*/
	tcSubHdlPtr subProcessHandler() const { return theSubProcessHandler; }

	/**
	* The matrix element to be used.
	*/
	tMEPtr matrixElement() const { return theME; }

	/**
	* Return a pointer to the head XComb this XComb
	* depends on. May return NULL, if this is not a
	* member of a XComb group.
	*/
	tStdXCombPtr head() const { return theHead; }

	/**
	* Set the head XComb pointer.
	*/
	void head(tStdXCombPtr headXC) { theHead = headXC; }

	/**
	* Return a selector object of xcombs to choose subprocesses
	* different than the one currently integrated.
	*/
	Selector<tStdXCombPtr>& projectors() { return theProjectors; }

	/**
	* Return a selector object of xcombs to choose subprocesses
	* different than the one currently integrated.
	*/
	const Selector<tStdXCombPtr>& projectors() const { return theProjectors; }

	/**
	* Return a pointer to a projector xcomb which will generate a subprocess
	* different from the one just integrated.
	*/
	tStdXCombPtr lastProjector() const { return theProjector; }

	/**
	* Set a pointer to a projector xcomb which will generate a subprocess
	* different from the one just integrated.
	*/
	void lastProjector(tStdXCombPtr pxc) { theProjector = pxc; }

	//@}

	/** @name Main functions used for the generation. */
	//@{
	/**
	* Try to determine if this subprocess is at all possible.
	*/
	virtual bool checkInit();

	/**
	* The number of dimensions of the phase space used to generate this
	* process.
	*/
	int nDim() const { return theNDim; }

	/**
	* Return the parton extraction dimensions
	*/
	const pair<int,int>& partonDimensions() const { return partonDims; }

	/**
	* Return true, if the current configuration will pass the cuts
	*/
	bool willPassCuts();

	/**
	+ * Return the cut weight encountered from fuzzy cuts
	+ */
	+ double cutWeight() const { return theCutWeight; }
	+
	+ /**
	* Reset all saved data about last generated phasespace point;
	*/
	virtual void clean();

	/**
	* Return true, if kinematics have already been generated
	*/
	bool kinematicsGenerated() const { return theKinematicsGenerated; }

	/**
	* Indicate that kinematics have been generated
	*/
	void didGenerateKinematics() { theKinematicsGenerated = true; }

	/**
	* Generate a phase space point from a vector \a r of \a nr numbers
	* in the interval ]0,1[ and return the corresponding differential
	* cross section.
	*/
	virtual CrossSection dSigDR(const pair<double,double> ll, int nr, const double * r);

	/**
	* If this XComb has a head XComb, return the cross section
	* differential in the variables previously supplied. The PDF weight
	* is taken from the lastPDFWeight supplied by the head XComb
	* object.
	*/
	CrossSection dSigDR(const double * r);

	/**
	* Return the PDF weight used in the last call to dSigDR
	*/
	double lastPDFWeight() const { return theLastPDFWeight; }

	/**
	* Return the cross section calculated in the last call to dSigDR
	*/
	CrossSection lastCrossSection() const { return theLastCrossSection; }

	/**
	* Construct a sub-process object from the information available.
	*/
	virtual tSubProPtr construct();
	//@}

	/** @name Functions used for collecting statistics. */
	//@{
	/**
	* The statistics object for this XComb.
	*/
	virtual const XSecStat & stats() const { return theStats; }

	/**
	* Select the current event. It will later be rejected with a
	* probability given by \a weight.
	*/
	virtual void select(double weight) { theStats.select(weight); }

	/**
	* Accept the current event assuming it was previously selcted.
	*/
	virtual void accept() { theStats.accept(); }

	/**
	* Reweight a selected and accepted event.
	*/
	void reweight(double oldWeight, double newWeight) {
	theStats.reweight(oldWeight,newWeight);
	}

	/**
	* Reject the current event assuming it was previously accepted. If
	* weighted events are produced, the \a weight should be the same as
	* the previous call to select(double).
	*/
	virtual void reject(double weight = 1.0) { theStats.reject(weight); }

	/**
	* Reset statistics.
	*/
	virtual void reset() { theStats.reset(); }
	//@}

	/** @name Access information used by the MEBase object. */
	//@{
	/**
	* The diagrams used by the matrix element.
	*/
	const DiagramVector & diagrams() const { return theDiagrams; }

	/**
	* True if the TreeDiagram's for this matrix element should in fact
	* be mirrored before used to create an actual sub-rocess.
	*/
	bool mirror() const { return isMirror; }

	/**
	* Return the momenta of the partons to be used by the matrix
	* element object, in the order specified by the TreeDiagram objects
	* given by the matrix element.
	*/
	const vector<Lorentz5Momentum> & meMomenta() const { return theMEMomenta; }

	/**
	* Return the last selected diagram.
	*/
	tcDiagPtr lastDiagram() const { return diagrams()[lastDiagramIndex()]; }

	/**
	* Return the parton types to be used by the matrix element object,
	* in the order specified by the TreeDiagram objects given by the
	* matrix element.
	*/
	const cPDVector & mePartonData() const { return theMEPartonData; }

	/**
	* Return the index of the last selected diagram.
	*/
	DiagramIndex lastDiagramIndex() const { return theLastDiagramIndex; }

	/**
	* Get information saved by the matrix element in the calculation of
	* the cross section to be used later when selecting diagrams and
	* colour flow.
	*/
	const DVector & meInfo() const { return theMEInfo; }

	/**
	* Set information saved by the matrix element in the calculation of
	* the cross section to be used later when selecting diagrams and
	* colour flow.
	*/
	void meInfo(const DVector & info) { theMEInfo = info; }

	/**
	* Return the random numbers used to generate the
	* last phase space point, if the matrix element
	* requested so.
	*/
	const DVector& lastRandomNumbers() const { return theLastRandomNumbers; }

	/**
	* Get the last jacobian obtained when generating the kinematics
	* for the call to dSigHatDR.
	*/
	double jacobian() const { return theLastJacobian; }

	/**
	* Return the matrix element squared as calculated
	* for the last phase space point. This may optionally
	* be used by a matrix element for caching.
	*/
	double lastME2() const { return theLastME2; }

	/**
	* Return the partonic cross section as calculated
	* for the last phase space point. This may optionally
	* be used by a matrix element for caching.
	*/
	CrossSection lastMECrossSection() const { return theLastMECrossSection; }

	/**
	* Return the PDF weight as calculated
	* for the last phase space point, if the matrix
	* element does supply PDF weights. This may optionally
	* be used by a matrix element for caching.
	*/
	double lastMEPDFWeight() const { return theLastMEPDFWeight; }
	//@}

	protected:

	/**
	* Construct the corresponding SubProcess object if it hasn't been
	* done before.
	*/
	virtual void newSubProcess(bool group = false);

	/**
	* Return the momenta of the partons to be used by the matrix
	* element object, in the order specified by the TreeDiagram objects
	* given by the matrix element.
	*/
	vector<Lorentz5Momentum> & meMomenta() { return theMEMomenta; }

	/**
	* Access the random numbers used to generate the
	* last phase space point, if the matrix element
	* requested so.
	*/
	DVector& lastRandomNumbers() { return theLastRandomNumbers; }

	/**
	* Return the parton types to be used by the matrix element object,
	* in the order specified by the TreeDiagram objects given by the
	* matrix element.
	*/
	cPDVector & mePartonData() { return theMEPartonData; }

	/**
	* Set the last selected diagram.
	*/
	void lastDiagramIndex(DiagramIndex i) { theLastDiagramIndex = i; }

	/**
	* Set the PDF weight used in the last call to dSigDR
	*/
	void lastPDFWeight(double w) { theLastPDFWeight = w; }

	/**
	* Set the cross section calculated in the last call to dSigDR
	*/
	void lastCrossSection(CrossSection s) { theLastCrossSection = s; }

	/**
	* Set the last jacobian obtained when generating the kinematics for
	* the call to dSigHatDR.
	*/
	void jacobian(double j) { theLastJacobian = j; }

	/**
	* Set the matrix element squared as calculated
	* for the last phase space point. This may optionally
	* be used by a matrix element for caching.
	*/
	void lastME2(double v) { theLastME2 = v; }

	/**
	* Set the partonic cross section as calculated
	* for the last phase space point. This may optionally
	* be used by a matrix element for caching.
	*/
	void lastMECrossSection(CrossSection v) { theLastMECrossSection = v; }

	/**
	* Set the PDF weight as calculated
	* for the last phase space point, if the matrix
	* element does supply PDF weights. This may optionally
	* be used by a matrix element for caching.
	*/
	void lastMEPDFWeight(double v) { theLastMEPDFWeight = v; }

	public:

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

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

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

	private:

	/**
	* The corresponding sub-process handler
	*/
	tSubHdlPtr theSubProcessHandler;

	/**
	* The matrix element to be used.
	*/
	tMEPtr theME;

	/**
	* Statistics gathering for this XComb.
	*/
	XSecStat theStats;

	/**
	* The diagrams used by the matrix element.
	*/
	DiagramVector theDiagrams;

	/**
	* True if the TreeDiagram's for this matrix element should in fact
	* be mirrored before used to create an actual sub-rocess.
	*/
	bool isMirror;

	/**
	* The number of dimensions of the phase space used to generate this
	* process.
	*/
	int theNDim;

	protected:

	/**
	* The number of dimensions of the phase space used for each of the
	* incoming partons.
	*/
	pair<int,int> partonDims;

	private:

	/**
	* True, if kinematics have already been generated
	*/
	bool theKinematicsGenerated;

	/**
	* The momenta of the partons to be used by the matrix element
	* object, in the order specified by the TreeDiagram objects given
	* by the matrix element.
	*/
	vector<Lorentz5Momentum> theMEMomenta;

	/**
	* The parton types to be used by the matrix element object, in the
	* order specified by the TreeDiagram objects given by the matrix
	* element.
	*/
	cPDVector theMEPartonData;

	/**
	* The last selected tree diagram.
	*/
	DiagramIndex theLastDiagramIndex;

	/**
	* Information saved by the matrix element in the calculation of the
	* cross section to be used later when selecting diagrams and colour
	* flow.
	*/
	DVector theMEInfo;

	/**
	* The random numbers used to generate the
	* last phase space point, if the matrix element
	* requested so.
	*/
	DVector theLastRandomNumbers;

	/**
	* The PDF weight used in the last call to dSigDR
	*/
	double theLastPDFWeight;

	/**
	* The cross section calculated in the last call to dSigDR
	*/
	CrossSection theLastCrossSection;

	/**
	* Save the last jacobian obtained when generating the kinematics for
	* the call to dSigHatDR.
	*/
	double theLastJacobian;

	/**
	* The matrix element squared as calculated
	* for the last phase space point. This may optionally
	* be used by a matrix element for caching.
	*/
	double theLastME2;

	/**
	* The partonic cross section as calculated
	* for the last phase space point. This may optionally
	* be used by a matrix element for caching.
	*/
	CrossSection theLastMECrossSection;

	/**
	* The PDF weight as calculated
	* for the last phase space point, if the matrix
	* element does supply PDF weights. This may optionally
	* be used by a matrix element for caching.
	*/
	double theLastMEPDFWeight;

	/**
	* A pointer to the head XComb this XComb
	* depends on. May return NULL, if this is not a
	* member of a XComb group.
	*/
	tStdXCombPtr theHead;

	/**
	* A selector object of xcombs to choose subprocesses
	* different than the one currently integrated.
	*/
	Selector<tStdXCombPtr> theProjectors;

	/**
	* A pointer to a projector xcomb which will generate a subprocess
	* different from the one just integrated.
	*/
	tStdXCombPtr theProjector;

	/**
	* True, if cuts have already been checked
	*/
	bool checkedCuts;

	/**
	* The result of the last call to willPassCuts
	*/
	bool passedCuts;

	+ /**
	+ * The cut weight encountered from fuzzy cuts
	+ */
	+ double theCutWeight;
	+
	private:

	/**
	* Describe a concrete class with persistent data.
	*/
	static ClassDescription<StandardXComb> initStandardXComb;

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

	};

	/** @cond TRAITSPECIALIZATIONS */

	/**
	* This template specialization informs ThePEG about the base class of
	* StandardXComb.
	*/
	template <>
	struct BaseClassTrait<StandardXComb,1>: public ClassTraitsType {
	/** Typedef of the base class of StandardXComb. */
	typedef XComb NthBase;
	};

	/**
	* This template specialization informs ThePEG about the name of the
	* StandardXComb class.
	*/
	template <>
	struct ClassTraits<StandardXComb>:
	public ClassTraitsBase<StandardXComb> {
	/** Return the class name. */
	static string className() { return "ThePEG::StandardXComb"; }
	};

	/** @endcond */

	}

	#endif /* ThePEG_StandardXComb_H */
	diff --git a/Handlers/StdXCombGroup.cc b/Handlers/StdXCombGroup.cc
	--- a/Handlers/StdXCombGroup.cc
	+++ b/Handlers/StdXCombGroup.cc
	@@ -1,391 +1,400 @@
	// -- C++ --
	//
	// StdXCombGroup.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 Leif Lonnblad
	// Copyright (C) 2009-2010 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the StdXCombGroup class.
	//

	#include "StdXCombGroup.h"
	#include "ThePEG/MatrixElement/MEGroup.h"
	#include "ThePEG/Handlers/StandardEventHandler.h"
	#include "ThePEG/Handlers/SubProcessHandler.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "ThePEG/Utilities/Maths.h"
	#include "ThePEG/PDT/ParticleData.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Utilities/SimplePhaseSpace.h"
	#include "ThePEG/Utilities/UtilityBase.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/EventRecord/SubProcessGroup.h"
	#include "ThePEG/Vectors/LorentzRotation.h"
	#include "ThePEG/MatrixElement/MEBase.h"
	#include "ThePEG/MatrixElement/ColourLines.h"
	#include "ThePEG/Handlers/LuminosityFunction.h"
	#include "ThePEG/Handlers/CascadeHandler.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/EventRecord/TmpTransform.h"

	using namespace ThePEG;

	StdXCombGroup::StdXCombGroup(Energy newMaxEnergy, const cPDPair & inc,
	tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
	tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
	const PBPair & newPartonBins, tCutsPtr newCuts, tMEGroupPtr newME,
	const DiagramVector & newDiagrams, bool mir, tStdXCombPtr newHead)
	: StandardXComb(newMaxEnergy,inc,newEventHandler,newSubProcessHandler,
	newExtractor, newCKKW, newPartonBins, newCuts,
	newME, newDiagrams, mir, newHead),
	theMEGroup(newME), theDependent(), theLastHeadCrossSection(ZERO) {}

	StdXCombGroup::StdXCombGroup()
	: StandardXComb(), theDependent() {}

	void StdXCombGroup::build(const PartonPairVec& allPBins) {
	for ( MEVector::const_iterator me = theMEGroup->dependent().begin();
	me != theMEGroup->dependent().end(); ++me ) {
	vector<StdXCombPtr> dep =
	theMEGroup->makeDependentXCombs(this,diagrams().front()->partons(),*me,allPBins);
	copy(dep.begin(),dep.end(),back_inserter(theDependent));
	}
	}

	StdXCombGroup::~StdXCombGroup() { }

	void StdXCombGroup::clean() {
	StandardXComb::clean();
	theLastHeadCrossSection = ZERO;
	for ( vector<StdXCombPtr>::const_iterator dep = theDependent.begin();
	dep != theDependent.end(); ++dep )
	(**dep).clean();
	}

	CrossSection StdXCombGroup::dSigDR(const pair<double,double> ll, int nr, const double * r) {

	matrixElement()->flushCaches();

	if ( matrixElement()->keepRandomNumbers() ) {
	lastRandomNumbers().resize(nDim());
	copy(r,r+nDim(),lastRandomNumbers().begin());
	}

	pExtractor()->select(this);
	setPartonBinInfo();
	lastP1P2(ll);
	lastS(sqr(maxEnergy())/exp(lastP1() + lastP2()));

	meMomenta().resize(mePartonData().size());
	matrixElement()->setXComb(this);

	PPair partons;

	if ( !matrixElement()->haveX1X2() ) {

	if ( !pExtractor()->generateL(partonBinInstances(),
	r, r + nr - partonDims.second) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	partons = make_pair(partonBinInstances().first->parton(),
	partonBinInstances().second->parton());
	lastSHat(lastS()/exp(partonBinInstances().first->l() +
	partonBinInstances().second->l()));
	meMomenta()[0] = partons.first->momentum();
	meMomenta()[1] = partons.second->momentum();

	} else {
	if ( !matrixElement()->generateKinematics(r + partonDims.first) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
	matrixElement()->setKinematics();

	lastScale(matrixElement()->scale());

	partons.first = mePartonData()[0]->produceParticle(meMomenta()[0]);
	partons.second = mePartonData()[1]->produceParticle(meMomenta()[1]);

	Direction<0> dir(true);
	partonBinInstances().first =
	new_ptr(PartonBinInstance(lastParticles().first,partons.first,
	partonBins().first,lastScale()));
	dir.reverse();
	partonBinInstances().second =
	new_ptr(PartonBinInstance(lastParticles().second,partons.second,
	partonBins().second,lastScale()));
	}

	lastPartons(partons);

	if ( lastSHat() < cuts()->sHatMin() ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	lastY(0.5*(partonBinInstances().second->l() -
	partonBinInstances().first->l()));
	if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]);
	if ( matrixElement()->wantCMS() &&
	!matrixElement()->haveX1X2() )
	SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat());

	Energy summ = ZERO;
	if ( meMomenta().size() == 3 ) {
	if ( !matrixElement()->haveX1X2() )
	meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat()));
	} else {
	for ( int i = 2, N = meMomenta().size(); i < N; ++i ) {
	if ( !matrixElement()->haveX1X2() )
	meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->mass());
	summ += mePartonData()[i]->massMin();
	}
	if ( sqr(summ) >= lastSHat() ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	}

	matrixElement()->setXComb(this);
	if ( !matrixElement()->haveX1X2() )
	lastScale(max(lastSHat()/4.0, cuts()->scaleMin()));

	lastSHat(pExtractor()->generateSHat(lastS(), partonBinInstances(),
	r, r + nr - partonDims.second,
	matrixElement()->haveX1X2()));

	if ( !cuts()->sHat(lastSHat()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	r += partonDims.first;

	lastX1X2(make_pair(lastPartons().first->momentum().plus()/
	lastParticles().first->momentum().plus(),
	lastPartons().second->momentum().minus()/
	lastParticles().second->momentum().minus()));

	if ( !cuts()->x1(lastX1()) \|\| !cuts()->x2(lastX2()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	lastY((lastPartons().first->momentum() +
	lastPartons().second->momentum()).rapidity());
	if ( !cuts()->yHat(lastY()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	meMomenta()[0] = lastPartons().first->momentum();
	meMomenta()[1] = lastPartons().second->momentum();
	if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]);
	if ( matrixElement()->wantCMS() &&
	!matrixElement()->haveX1X2() )
	SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat());

	if ( meMomenta().size() == 3 ) {
	if ( !matrixElement()->haveX1X2() )
	meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat()));
	} else {
	if ( sqr(summ) >= lastSHat() ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	}

	matrixElement()->setXComb(this);
	if ( !matrixElement()->haveX1X2() ) {
	if ( !matrixElement()->generateKinematics(r) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	}
	lastScale(matrixElement()->scale());
	if ( !cuts()->scale(lastScale()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	pair<bool,bool> evalPDFS =
	make_pair(matrixElement()->havePDFWeight1(),
	matrixElement()->havePDFWeight2());
	if ( mirror() )
	swap(evalPDFS.first,evalPDFS.second);
	lastPDFWeight(pExtractor()->fullFn(partonBinInstances(), lastScale(),
	evalPDFS));
	if ( lastPDFWeight() == 0.0 ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	matrixElement()->setKinematics();
	CrossSection xsec = matrixElement()->dSigHatDR() * lastPDFWeight();

	bool noHeadPass = !willPassCuts() \|\| xsec == ZERO;
	if ( noHeadPass ) {
	xsec = ZERO;
	lastCrossSection(ZERO);
	}

	+ xsec *= cutWeight();
	+
	lastAlphaS(matrixElement()->alphaS());
	lastAlphaEM(matrixElement()->alphaEM());

	lastHeadCrossSection(xsec);

	CrossSection depxsec = ZERO;

	vector<tStdXCombPtr> activeXCombs;

	if ( !theMEGroup->mcSumDependent() ) {
	for ( vector<StdXCombPtr>::const_iterator dep = theDependent.begin();
	dep != theDependent.end(); ++dep ) {
	if ( noHeadPass && (**dep).matrixElement()->headCuts() )
	continue;
	depxsec += (dep).dSigDR(r + theMEGroup->dependentOffset((dep).matrixElement()));
	if ( theMEGroup->groupReweighted() )
	activeXCombs.push_back(*dep);
	}
	} else {
	if ( theMEGroup->lastDependentXComb() ) {
	tStdXCombPtr dxc = theMEGroup->lastDependentXComb();
	if ( !(noHeadPass && dxc->matrixElement()->headCuts()) )
	depxsec = theDependent.size()*
	dxc->dSigDR(r + theMEGroup->dependentOffset(dxc->matrixElement()));
	if ( theMEGroup->groupReweighted() )
	activeXCombs.push_back(dxc);
	}
	}

	- if ( xsec != ZERO ) {
	- double rw = 1.0 + depxsec/xsec;
	- xsec *= rw;
	- } else {
	+ if ( xsec != ZERO &&
	+ depxsec != ZERO ) {
	+ if ( abs(xsec) < abs(depxsec) ) {
	+ volatile double rw = (1.+depxsec/xsec);
	+ xsec = xsec*rw;
	+ } else {
	+ volatile double rw = (1.+xsec/depxsec);
	+ xsec = depxsec*rw;
	+ }
	+ } else if ( xsec == ZERO &&
	+ depxsec != ZERO ) {
	xsec = depxsec;
	}

	lastCrossSection(xsec);

	if ( xsec != ZERO )
	theMEGroup->lastEventStatistics();

	matrixElement()->fillProjectors();
	if ( !projectors().empty() ) {
	lastProjector(projectors().select(UseRandom::rnd()));
	}

	if ( theMEGroup->groupReweighted() && !activeXCombs.empty() ) {
	xsec = theMEGroup->reweightHead(activeXCombs)*lastHeadCrossSection();
	depxsec = ZERO;
	for ( vector<tStdXCombPtr>::const_iterator dep = activeXCombs.begin();
	dep != activeXCombs.end(); ++dep ) {
	depxsec += theMEGroup->reweightDependent(dep,activeXCombs)(**dep).lastCrossSection();
	}
	if ( xsec != ZERO ) {
	double rw = 1.0 + depxsec/xsec;
	xsec *= rw;
	} else {
	xsec = depxsec;
	}
	}

	subProcess(SubProPtr());
	if ( CKKWHandler() && matrixElement()->maxMultCKKW() > 0 &&
	matrixElement()->maxMultCKKW() > matrixElement()->minMultCKKW() ) {
	newSubProcess(theMEGroup->subProcessGroups());
	CKKWHandler()->setXComb(this);
	xsec *= CKKWHandler()->reweightCKKW(matrixElement()->minMultCKKW(),
	matrixElement()->maxMultCKKW());
	}

	if ( matrixElement()->reweighted() ) {
	newSubProcess(theMEGroup->subProcessGroups());
	xsec = matrixElement()->reWeight() matrixElement()->preWeight();
	}

	return xsec;

	}

	void StdXCombGroup::newSubProcess(bool) {

	StandardXComb::newSubProcess(theMEGroup->subProcessGroups());

	if ( !theMEGroup->subProcessGroups() )
	return;

	subProcess()->groupWeight(lastHeadCrossSection()/lastCrossSection());

	Ptr<SubProcessGroup>::tptr group =
	dynamic_ptr_cast<Ptr<SubProcessGroup>::tptr>(subProcess());
	assert(group);

	for ( vector<StdXCombPtr>::iterator dep = theDependent.begin();
	dep != theDependent.end(); ++dep ) {
	if ( (**dep).lastCrossSection() == ZERO )
	continue;
	tSubProPtr ds;
	try {
	ds = (**dep).construct();
	} catch(Veto&) {
	throw Exception() << "A veto was encountered while constructing a dependent sub process.\n"
	<< "This situation should not have happened. Will veto the event."
	<< Exception::eventerror;
	}
	if ( ds )
	group->add(ds);
	}

	}

	tSubProPtr StdXCombGroup::construct() {

	matrixElement()->setXComb(this);

	setPartonBinInfo();
	matrixElement()->setKinematics();

	newSubProcess(theMEGroup->subProcessGroups());

	if ( !theMEGroup->subProcessGroups() ) {
	if ( !cuts()->initSubProcess(lastSHat(), lastY()) ) throw Veto();
	TmpTransform<tSubProPtr>
	tmp(subProcess(), Utilities::getBoostToCM(subProcess()->incoming()));
	if ( !cuts()->passCuts(*subProcess()) ) throw Veto();
	}

	return subProcess();

	}


	void StdXCombGroup::Init() {}

	void StdXCombGroup::persistentOutput(PersistentOStream & os) const {
	os << theDependent << theMEGroup << ounit(theLastHeadCrossSection,nanobarn);
	}

	void StdXCombGroup::persistentInput(PersistentIStream & is, int) {
	is >> theDependent >> theMEGroup >> iunit(theLastHeadCrossSection,nanobarn);
	}

	ClassDescription<StdXCombGroup> StdXCombGroup::initStdXCombGroup;

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 4:42 PM (1 d, 12 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805209
Default Alt Text: (141 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions