No OneTemporary
Actions

Size

47 KB

Subscribers

None

View Options

	diff --git a/Cuts/Cuts.cc b/Cuts/Cuts.cc
	--- a/Cuts/Cuts.cc
	+++ b/Cuts/Cuts.cc
	@@ -1,661 +1,663 @@
	// -- 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), 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";
	+ if ( theJetFinder )
	+ theJetFinder->describe();
	}

	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) ) {
	theCutWeight = 0.0;
	return false;
	}
	if ( yhat <= yHatMin() \|\| yhat >= yHatMax() ) {
	theCutWeight = 0.0;
	return false;
	}
	double x1 = min(1.0, sqrt(shat/SMax())*exp(yhat));
	if ( x1 <= x1Min() \|\| x1 > x1Max() ) {
	theCutWeight = 0.0;
	return false;
	}
	double x2 = min(1.0, sqrt(shat/SMax())*exp(-yhat));
	if ( x2 <= x2Min() \|\| x2 > x2Max() ) {
	theCutWeight = 0.0;
	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;

	if ( jetFinder() ) {

	if ( ptype.size() > jetFinder()->minOutgoing() ) {


	vector<LorentzMomentum> jets;
	tcPDVector jettype;

	if ( !jetFinder()->restrictConsitutents() ) {
	jets = p;
	jettype = ptype;
	} else {
	tcPDVector::const_iterator pd = ptype.begin();
	vector<LorentzMomentum>::const_iterator pm = p.begin();
	for ( ; pd != ptype.end(); ++pd, ++pm ) {
	if ( pm->rapidity() > jetFinder()->constituentRapidityRange().first &&
	pm->rapidity() < jetFinder()->constituentRapidityRange().second ) {
	jets.push_back(*pm);
	jettype.push_back(*pd);
	}
	}
	}

	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 ) {
	pass &= theOneCuts[j]->passCuts(this, ptype[i], p[i]);
	theCutWeight *= theLastCutWeight;
	theLastCutWeight = 1.0;
	if ( !pass ) {
	theCutWeight = 0.0;
	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 ) {
	pass &= theTwoCuts[j]->passCuts(this, ptype[i1], ptype[i2],
	p[i1], p[i2]);
	theCutWeight *= theLastCutWeight;
	theLastCutWeight = 1.0;
	if ( !pass ) {
	theCutWeight = 0.0;
	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 ) {
	theCutWeight = 0.0;
	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 ) {
	pass &= theTwoCuts[j]->passCuts(this, t1, ptype[i], p1, p[i],
	true, false);
	theCutWeight *= theLastCutWeight;
	theLastCutWeight = 1.0;
	if ( !pass ) {
	theCutWeight = 0.0;
	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 ) {
	pass &= theTwoCuts[j]->passCuts(this, ptype[i], t2, p[i], p2,
	false, true);
	theCutWeight *= theLastCutWeight;
	theLastCutWeight = 1.0;
	if ( !pass ) {
	theCutWeight = 0.0;
	return false;
	}
	}
	}

	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
	<< theCutWeight << theLastCutWeight << theFuzzyTheta;
	}

	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
	>> theCutWeight >> theLastCutWeight >> theFuzzyTheta;
	}

	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);

	static Reference<Cuts,FuzzyTheta> interfaceFuzzy
	("Fuzzy",
	"The fuzziness to be applied to cuts (may not be supported by all cut objects).",
	&Cuts::theFuzzyTheta, 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/FastJetFinder.cc b/Cuts/FastJetFinder.cc
	--- a/Cuts/FastJetFinder.cc
	+++ b/Cuts/FastJetFinder.cc
	@@ -1,269 +1,297 @@
	// -- C++ --
	//
	// FastJetFinder.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.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the FastJetFinder class.
	//

	#include "FastJetFinder.h"
	#include "ThePEG/Interface/ClassDocumentation.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/Repository/CurrentGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "fastjet/ClusterSequence.hh"


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

	using namespace ThePEG;

	FastJetFinder::FastJetFinder()
	: theDCut(ZERO), theConeRadius(0.7),
	theVariant(kt), theMode(inclusive),
	theRecombination(recoE) {}

	FastJetFinder::~FastJetFinder() {}

	void FastJetFinder::doinit() {
	fastjet::ClusterSequence::set_fastjet_banner_stream(0);
	}

	void FastJetFinder::doinitrun() {
	fastjet::ClusterSequence::set_fastjet_banner_stream(& CurrentGenerator::log());
	// Force banner to be printed right now; suppresses later output.
	fastjet::ClusterSequence::print_banner();
	}

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

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

	+void FastJetFinder::describe() const {
	+ generator()->log()
	+ << "'" << name() << "' clustering jets from constituents matched by '"
	+ << unresolvedMatcher()->name() << "'\n"
	+ << "into " << (theMode == inclusive ? "inclusive" : "exclusive")
	+ << " jets recombining with the "
	+ << (theRecombination == recoPt ? "pt" : "E")
	+ << " scheme\n";
	+ generator()->log() << "The measure used is ";
	+ switch(theVariant) {
	+ case 1: generator()->log() << "kt"; break;
	+ case 2: generator()->log() << "CA"; break;
	+ case 3: generator()->log() << "antiKt"; break;
	+ case 4: generator()->log() << "sphericalKt"; break;
	+ case 5: generator()->log() << "sphericalCA"; break;
	+ case 6: generator()->log() << "sphericalAntiKt"; break;
	+ default: assert(false);
	+ }
	+ generator()->log() << "\n";
	+ generator()->log() << "The cone radius is R = "
	+ << theConeRadius << "\n";
	+ if ( theMode == exclusive ) {
	+ generator()->log() << "The exclusive resolution scale in GeV is D = "
	+ << sqrt(theDCut/GeV2) << "\n";
	+ }
	+ generator()->log() << flush;
	+}
	+
	bool FastJetFinder::cluster(tcPDVector & ptype, vector<LorentzMomentum> & p,
	tcCutsPtr, tcPDPtr, tcPDPtr) const {
	if ( ptype.size() <= minOutgoing() ){
	return false;
	}

	tcPDVector::iterator di = ptype.begin();
	vector<LorentzMomentum>::iterator pi = p.begin();
	size_t index = 0;

	vector<fastjet::PseudoJet> recombinables;
	tcPDVector ptypeBuffer;
	vector<LorentzMomentum> pBuffer;
	for ( ; di != ptype.end(); ++di, ++pi, index++ ) {
	if ( !unresolvedMatcher()->check(**di) ) {
	ptypeBuffer.push_back(*di);
	pBuffer.push_back(*pi);
	continue;
	}
	recombinables.push_back(fastjet::PseudoJet( (pi).x()/GeV,(pi).y()/GeV,(pi).z()/GeV, (pi).t()/GeV ));
	recombinables.back().set_user_index(index);
	}

	fastjet::Strategy strategy = fastjet::Best;

	fastjet::RecombinationScheme recomb_scheme;
	if ( theRecombination == recoE )
	recomb_scheme = fastjet::E_scheme;
	else if ( theRecombination == recoPt)
	recomb_scheme = fastjet::pt_scheme;
	else assert(false);

	fastjet::JetAlgorithm jet_algorithm = fastjet::kt_algorithm;
	if ( theVariant == CA ) {
	jet_algorithm = fastjet::cambridge_algorithm;
	} else if ( theVariant == antiKt ) {
	jet_algorithm = fastjet::antikt_algorithm;
	} else if ( theVariant > 3 ) {
	jet_algorithm = fastjet::ee_genkt_algorithm;
	}

	fastjet::JetDefinition jet_def;
	if ( theVariant < 4 ) {
	jet_def = fastjet::JetDefinition(jet_algorithm, theConeRadius, recomb_scheme, strategy);
	} else {
	int power = 1;
	if ( theVariant == sphericalCA ) {
	power = 0;
	}
	if ( theVariant == sphericalAntiKt ) {
	power = -1;
	}
	jet_def = fastjet::JetDefinition(jet_algorithm, theConeRadius, power, recomb_scheme, strategy);
	}

	fastjet::ClusterSequence clust_seq(recombinables, jet_def);

	double dcut = 0.0;

	if ( theVariant != antiKt &&
	theVariant != sphericalAntiKt ) {
	dcut = theDCut/GeV2;
	} else {
	dcut = theDCut != ZERO ? GeV2/theDCut : ZERO;
	}

	vector<fastjet::PseudoJet> recoJets;
	if ( theMode == inclusive )
	recoJets = clust_seq.inclusive_jets();
	else if ( theMode == exclusive )
	recoJets = clust_seq.exclusive_jets(dcut);


	if ( recoJets.size() + pBuffer.size() == p.size() ){
	return false;
	}
	else {
	tcPDVector ptypeNew;
	vector<LorentzMomentum> pNew;
	for (vector<fastjet::PseudoJet>::const_iterator iter = recoJets.begin();
	iter != recoJets.end(); iter++){
	ptypeNew.push_back(ptype[iter->constituents().begin()->user_index()]);
	pNew.push_back(LorentzMomentum(iter->px()GeV, iter->py()GeV, iter->pz()GeV, iter->E()GeV));
	}

	di = ptypeBuffer.begin();
	pi = pBuffer.begin();
	for (; di != ptypeBuffer.end(); di++, pi++){
	ptypeNew.push_back(*di);
	pNew.push_back(*pi);
	}
	ptype = ptypeNew;
	p = pNew;
	return true;
	}
	}


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


	void FastJetFinder::persistentOutput(PersistentOStream & os) const {
	os << ounit(theDCut,GeV2) << theConeRadius << theVariant << theMode
	<< theRecombination;
	}

	void FastJetFinder::persistentInput(PersistentIStream & is, int) {
	is >> iunit(theDCut,GeV2) >> theConeRadius >> theVariant >> theMode
	>> theRecombination;
	}


	// * 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<FastJetFinder,JetFinder>
	describeFastJetFinder("ThePEG::FastJetFinder", "FastJetFinder.so");

	void FastJetFinder::Init() {

	static ClassDocumentation<FastJetFinder> documentation
	("FastJetFinder implements the class of longitudinally invariant kt "
	"jet clustering algorithms, as relevant for cuts on the real "
	"emission contribution to a NLO calculation. Recombination is "
	"exclusively performed using the pt scheme.");


	static Parameter<FastJetFinder,Energy2> interfaceDCut
	("DCut",
	"The distance cut, when acting exclusively. "
	"The inverse is taken for the anti-kt algorithm, "
	"while for the Cambridge/Aachen variant dCut/GeV2 is used.",
	&FastJetFinder::theDCut, GeV2, 0.0GeV2, 0.0GeV2, 0*GeV2,
	false, false, Interface::lowerlim);


	static Parameter<FastJetFinder,double> interfaceConeRadius
	("ConeRadius",
	"The cone radius R used in inclusive mode.",
	&FastJetFinder::theConeRadius, 0.7, 0.0, 10.0,
	false, false, Interface::limited);

	static Switch<FastJetFinder,int> interfaceVariant
	("Variant",
	"The variant to use.",
	&FastJetFinder::theVariant, kt, false, false);
	static SwitchOption interfaceVariantKt
	(interfaceVariant,
	"Kt",
	"Kt algorithm.",
	kt);
	static SwitchOption interfaceVariantCA
	(interfaceVariant,
	"CA",
	"Cambridge/Aachen algorithm.",
	CA);
	static SwitchOption interfaceVariantAntiKt
	(interfaceVariant,
	"AntiKt",
	"Anti kt algorithm.",
	antiKt);
	static SwitchOption interfaceVariantSphericalKt
	(interfaceVariant,
	"SphericalKt",
	"Spherical kt algorithm.",
	sphericalKt);
	static SwitchOption interfaceVariantSphericalCA
	(interfaceVariant,
	"SphericalCA",
	"Spherical Cambridge/Aachen algorithm.",
	sphericalCA);
	static SwitchOption interfaceVariantSphericalAntiKt
	(interfaceVariant,
	"SphericalAntiKt",
	"Spherical anti kt algorithm.",
	sphericalAntiKt);

	static Switch<FastJetFinder,int> interfaceMode
	("Mode",
	"The mode to use.",
	&FastJetFinder::theMode, inclusive, false, false);
	static SwitchOption interfaceModeInclusive
	(interfaceMode,
	"Inclusive",
	"Find inclusive jets.",
	inclusive);
	static SwitchOption interfaceModeExclusive
	(interfaceMode,
	"Exclusive",
	"Find exclusive jets.",
	exclusive);

	static Switch<FastJetFinder,int> interfaceRecombination
	("RecombinationScheme",
	"The recombination scheme to use.",
	&FastJetFinder::theRecombination, recoE, false, false);
	static SwitchOption interfaceRecombinationPt
	(interfaceRecombination,
	"Pt",
	"Add transverse momenta",
	recoPt);
	static SwitchOption interfaceRecombinationE
	(interfaceRecombination,
	"E",
	"Add the four-momenta",
	recoE);

	}

	diff --git a/Cuts/FastJetFinder.h b/Cuts/FastJetFinder.h
	--- a/Cuts/FastJetFinder.h
	+++ b/Cuts/FastJetFinder.h
	@@ -1,182 +1,187 @@
	// -- C++ --
	//
	// FastJetFinder.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_FastJetFinder_H
	#define THEPEG_FastJetFinder_H
	//
	// This is the declaration of the FastJetFinder class.
	//

	#include "ThePEG/Cuts/JetFinder.h"

	namespace ThePEG {

	/**
	* FastJetFinder implements the class of longitudinally invariant kt
	* jet clustering algorithms.
	*
	* @see \ref FastJetFinderInterfaces "The interfaces"
	* defined for FastJetFinder.
	*/
	class FastJetFinder: public JetFinder {

	public:

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

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

	public:

	/**
	* Perform jet clustering on the given outgoing particles.
	* Optionally, information on the incoming particles is provided.
	* Return true, if a clustering has been performed.
	*/
	virtual bool cluster(tcPDVector & ptype, vector<LorentzMomentum> & p,
	tcCutsPtr parent, tcPDPtr t1 = tcPDPtr(),
	tcPDPtr t2 = tcPDPtr()) const;

	+ /**
	+ * Describe this jet fined.
	+ */
	+ virtual void describe() const;
	+
	public:
	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object. Called in the read phase.
	*/
	virtual void doinit();

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

	/** @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 resolution cut.
	*/
	Energy2 theDCut;

	/**
	* The `cone radius' R.
	*/
	double theConeRadius;

	/**
	* The possible variants.
	*/
	enum variants {
	kt = 1,
	CA = 2,
	antiKt = 3,
	sphericalKt = 4,
	sphericalCA = 5,
	sphericalAntiKt = 6
	};

	/**
	* The variant.
	*/
	int theVariant;

	/**
	* The possible modes.
	*/
	enum modes {
	inclusive = 1,
	exclusive = 2
	};

	/**
	* The mode.
	*/
	int theMode;

	/**
	* The possible recombination schemes.
	*/
	enum recombinations {
	recoPt = 1,
	recoE = 2
	};

	/**
	* The recombination scheme
	*/
	int theRecombination;

	private:

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

	};

	}

	#endif /* THEPEG_FastJetFinder_H */
	diff --git a/Cuts/JetFinder.h b/Cuts/JetFinder.h
	--- a/Cuts/JetFinder.h
	+++ b/Cuts/JetFinder.h
	@@ -1,170 +1,175 @@
	// -- C++ --
	//
	// JetFinder.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 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_JetFinder_H
	#define THEPEG_JetFinder_H
	//
	// This is the declaration of the JetFinder class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "ThePEG/PDT/MatcherBase.h"
	#include "Cuts.fh"

	namespace ThePEG {

	/**
	* JetFinder defines an interface to jet finders to be used when cuts
	* should actually be defined on the level of reconstructed jets such
	* as typically encountered in higher order corrections.
	*
	* @see \ref JetFinderInterfaces "The interfaces"
	* defined for JetFinder.
	*/
	class JetFinder: public Interfaced {

	public:

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

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

	public:

	/**
	* Perform jet clustering on the given outgoing particles.
	* Optionally, information on the incoming particles is provided.
	* Return true, if a clustering has been performed.
	*/
	virtual bool cluster(tcPDVector & ptype, vector<LorentzMomentum> & p,
	tcCutsPtr parent, tcPDPtr t1 = tcPDPtr(),
	tcPDPtr t2 = tcPDPtr()) const = 0;

	/**
	* Return the matcher for unresolved partons.
	*/
	Ptr<MatcherBase>::tptr unresolvedMatcher() const { return theUnresolvedMatcher; }

	/**
	* Set the minimum number of outgoing partons on which clustering
	* should be performed.
	*/
	void minOutgoing(unsigned int n) { theMinOutgoing = n; }

	/**
	* Return the minimum number of outgoing partons on which clustering
	* should be performed.
	*/
	unsigned int minOutgoing() const { return theMinOutgoing; }

	/**
	* Return true, if jets should only be constructed from matching
	* objects inside a given rapidity interval.
	*/
	bool restrictConsitutents() const { return theRestrictConstituents; }

	/**
	* Jets should only be constructed from matching
	* objects inside a given rapidity interval.
	*/
	void restrictConsitutents(bool on) { theRestrictConstituents = on; }

	/**
	* Return the rapidity interval within objects should be considered
	* for clustering, if appropriate
	*/
	const pair<double,double>& constituentRapidityRange() const { return theConstituentRapidityRange; }

	/**
	* Set the rapidity interval within objects should be considered
	* for clustering, if appropriate
	*/
	void constituentRapidityRange(const pair<double,double>& r) { theConstituentRapidityRange = r; }

	+ /**
	+ * Describe this jet finder
	+ */
	+ virtual void describe() const = 0;
	+
	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();


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

	/**
	* A matcher for unresolved partons.
	*/
	Ptr<MatcherBase>::ptr theUnresolvedMatcher;

	/**
	* The minimum number of outgoing partons on which clustering
	* should be performed.
	*/
	unsigned int theMinOutgoing;

	/**
	* True, if jets should only be constructed from matching
	* objects inside a given rapidity interval.
	*/
	bool theRestrictConstituents;

	/**
	* The rapidity interval within objects should be considered
	* for clustering, if appropriate
	*/
	pair<double,double> theConstituentRapidityRange;

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

	};

	}

	#endif /* THEPEG_JetFinder_H */
	diff --git a/LesHouches/LesHouchesFileReader.h b/LesHouches/LesHouchesFileReader.h
	--- a/LesHouches/LesHouchesFileReader.h
	+++ b/LesHouches/LesHouchesFileReader.h
	@@ -1,333 +1,335 @@
	// -- C++ --
	//
	// LesHouchesFileReader.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_LesHouchesFileReader_H
	#define THEPEG_LesHouchesFileReader_H
	// This is the declaration of the LesHouchesFileReader class.

	#include "LesHouchesReader.h"
	#include "LesHouchesFileReader.fh"
	#include "ThePEG/PDT/Decayer.h"
	#include "ThePEG/Utilities/CFileLineReader.h"
	#include <string>
	#include <stdio.h>

	namespace ThePEG {


	/**
	* LesHouchesFileReader is an base class to be used for objects which
	* reads event files from matrix element generators. It inherits from
	* LesHouchesReader and extends it by defining a file handle to be
	* read from, which is opened and closed by the open() and close()
	* functions. Note that the file handle is a standard C filehandle and
	* not a C++ stream. This is because there is no standard way in C++
	* to connect a pipe to a stream for reading eg. gzipped files. This
	* class is able to read plain event files conforming to the Les
	* Houches Event File accord.
	*
	* @see \ref LesHouchesFileReaderInterfaces "Th1e interfaces"
	* defined for LesHouchesFileReader.
	* @see Event
	* @see LesHouchesReader
	*/
	class LesHouchesFileReader: public LesHouchesReader {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Default constructor.
	*/
	- LesHouchesFileReader() : neve(0), ieve(0), theQNumbers(false) {}
	+ LesHouchesFileReader() : neve(0), ieve(0), theQNumbers(false),
	+ theIncludeFxFxTags(false),
	+ theIncludeCentral(false) {}

	/**
	* Copy-constructor. Note that a file which is opened in the object
	* copied from will have to be reopened in this.
	*/
	LesHouchesFileReader(const LesHouchesFileReader &);

	/**
	* Destructor.
	*/
	virtual ~LesHouchesFileReader();
	//@}

	public:

	/** @name Virtual functions specified by the LesHouchesReader base class. */
	//@{
	/**
	* Initialize. This function is called by the LesHouchesEventHandler
	* to which this object is assigned.
	*/
	virtual void initialize(LesHouchesEventHandler & eh);

	/**
	* Open a file with events. Derived classes should overwrite it and
	* first calling it before reading in the run information into the
	* corresponding protected variables.
	*/
	virtual void open();

	/**
	* Close the file from which events have been read.
	*/
	virtual void close();


	/**
	* Read the next event from the file or stream into the
	* corresponding protected variables. Return false if there is no
	* more events or if this was not a LHF event file.
	*/
	virtual bool doReadEvent();
	//@}

	/**
	* Return the name of the file from where to read events.
	*/
	string filename() const { return theFileName; }

	/**
	* Return the optional weights information string ("Names")
	*/

	virtual vector<string> optWeightsNamesFunc();


	public:

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

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

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


	/**
	* Erases all occurences of a substring from a string
	*/

	void erase_substr(std::string& subject, const std::string& search);


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

	/** @name Standard (and non-standard) Interfaced functions. */
	//@{
	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit();

	/**
	* Return true if this object needs to be initialized before all
	* other objects because it needs to extract PDFs from the event file.
	*/
	virtual bool preInitialize() const;
	//@

	protected:

	/**
	* The wrapper around the C FILE stream from which to read
	*/
	CFileLineReader cfile;

	protected:

	/**
	* The number of events in this file.
	*/
	long neve;

	/**
	* The current event number.
	*/
	long ieve;

	/**
	* If the file is a standard Les Houches formatted file (LHF) this
	* is its version number. If empty, this is not a Les Houches
	* formatted file
	*/
	string LHFVersion;

	/**
	* If LHF. All lines (since the last open() or readEvent()) outside
	* the header, init and event tags.
	*/
	string outsideBlock;

	/**
	* If LHF. All lines from the header block.
	*/
	string headerBlock;

	/**
	* If LHF. Additional comments found in the init block.
	*/
	string initComments;

	/**
	* If LHF. Map of attributes (name-value pairs) found in the init
	* tag.
	*/
	map<string,string> initAttributes;

	/**
	* If LHF. Additional comments found with the last read event.
	*/
	string eventComments;

	/**
	* If LHF. Map of attributes (name-value pairs) found in the last
	* event tag.
	*/
	map<string,string> eventAttributes;

	private:

	/**
	* The name of the file from where to read events.
	*/
	string theFileName;

	/**
	* Whether or not to search for QNUMBERS stuff
	*/
	bool theQNumbers;

	/**
	* Include/Read FxFx tags
	*/
	bool theIncludeFxFxTags;

	/**
	* Include central weight (for backup use)
	*/
	bool theIncludeCentral;

	/**
	* Decayer for any decay modes read from the file
	*/
	DecayerPtr theDecayer;

	/**
	* Further information on the weights
	*/
	map<string,string> scalemap;

	/**
	* Temporary holder for optional weights
	*/

	map<string,double> optionalWeightsTemp;


	private:

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

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

	public:

	/** @cond EXCEPTIONCLASSES */
	/** Exception class used by LesHouchesFileReader if reading the file
	* fails. */
	class LesHouchesFileError: public Exception {};
	/** @endcond */

	};

	}


	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/**
	* This template specialization informs ThePEG about the name of the
	* LesHouchesFileReader class and the shared object where it is
	* defined.
	*/
	template <>
	struct ClassTraits<LesHouchesFileReader>
	: public ClassTraitsBase<LesHouchesFileReader> {
	/**
	* Return the class name.
	*/
	static string className() { return "ThePEG::LesHouchesFileReader"; }
	/**
	* Return the name of the shared library to be loaded to get access
	* to the LesHouchesFileReader class and every other class it uses
	* (except the base class).
	*/
	static string library() { return "LesHouches.so"; }

	};

	/** @endcond */

	}

	#endif /* THEPEG_LesHouchesFileReader_H */

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 10:34 AM (1 d, 11 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5111202
Default Alt Text: (47 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions