No OneTemporary
Actions

Size

28 KB

Subscribers

None

View Options

	diff --git a/Cuts/Cuts.cc b/Cuts/Cuts.cc
	--- a/Cuts/Cuts.cc
	+++ b/Cuts/Cuts.cc
	@@ -1,615 +1,635 @@
	// -- 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";
	}

	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 ( jetFinder() ) {
	+
	if ( ptype.size() > jetFinder()->minOutgoing() ) {
	- vector<LorentzMomentum> jets = p;
	- tcPDVector jettype = ptype;
	+
	+
	+ 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 )
	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 )
	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 ) {
	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 ) {
	pass &= theTwoCuts[j]->passCuts(this, ptype[i], t2, p[i], p2,
	false, true);
	theCutWeight *= theLastCutWeight;
	theLastCutWeight = 1.0;
	if ( !pass )
	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;
	}

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

	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/JetFinder.cc b/Cuts/JetFinder.cc
	--- a/Cuts/JetFinder.cc
	+++ b/Cuts/JetFinder.cc
	@@ -1,76 +1,113 @@
	// -- 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.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the JetFinder class.
	//

	#include "JetFinder.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	+#include "ThePEG/Interface/Switch.h"
	+#include "ThePEG/Interface/Command.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"

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

	using namespace ThePEG;

	JetFinder::JetFinder()
	- : theMinOutgoing(1) {}
	+ : theMinOutgoing(1), theRestrictConstituents(false),
	+ theConstituentRapidityRange(Constants::MaxRapidity,Constants::MaxRapidity) {}

	JetFinder::~JetFinder() {}

	+string JetFinder::doYRange(string in) {
	+ istringstream ins(in);
	+ double first, second;
	+ ins >> first >> second;
	+ if ( first > second )
	+ swap(first,second);
	+ theConstituentRapidityRange = make_pair(first,second);
	+ return "";
	+}
	+
	+

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

	void JetFinder::persistentOutput(PersistentOStream & os) const {
	- os << theUnresolvedMatcher << theMinOutgoing;
	+ os << theUnresolvedMatcher << theMinOutgoing
	+ << theRestrictConstituents << theConstituentRapidityRange;
	}

	void JetFinder::persistentInput(PersistentIStream & is, int) {
	- is >> theUnresolvedMatcher >> theMinOutgoing;
	+ is >> theUnresolvedMatcher >> theMinOutgoing
	+ >> theRestrictConstituents >> theConstituentRapidityRange;
	}


	// * 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).
	DescribeAbstractClass<JetFinder,Interfaced>
	describeJetFinder("ThePEG::JetFinder", "");

	void JetFinder::Init() {

	static ClassDocumentation<JetFinder> documentation
	("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.");


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


	static Parameter<JetFinder,unsigned int> interfaceMinOutgoing
	("MinOutgoing",
	"The minimum number of outgoing partons to be clustered.",
	&JetFinder::theMinOutgoing, 1, 1, 0,
	false, false, Interface::lowerlim);

	+
	+ static Switch<JetFinder,bool> interfaceRestrictConstituents
	+ ("RestrictConstituents",
	+ "Restrict the constituents for clustering.",
	+ &JetFinder::theRestrictConstituents, false, false, false);
	+ static SwitchOption interfaceRestrictConstituentsYes
	+ (interfaceRestrictConstituents,
	+ "Yes",
	+ "",
	+ true);
	+ static SwitchOption interfaceRestrictConstituentsNo
	+ (interfaceRestrictConstituents,
	+ "No",
	+ "",
	+ false);
	+
	+ static Command<JetFinder> interfaceYRange
	+ ("ConstituentRapidityRange",
	+ "Restrict clustering to a rapidity range.",
	+ &JetFinder::doYRange, false);
	+
	}

	diff --git a/Cuts/JetFinder.h b/Cuts/JetFinder.h
	--- a/Cuts/JetFinder.h
	+++ b/Cuts/JetFinder.h
	@@ -1,129 +1,170 @@
	// -- 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; }
	+
	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 */

File Metadata

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

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions