No OneTemporary
Actions

Size

58 KB

Subscribers

None

View Options

	diff --git a/include/Rivet/Cuts.hh b/include/Rivet/Cuts.hh
	--- a/include/Rivet/Cuts.hh
	+++ b/include/Rivet/Cuts.hh
	@@ -1,83 +1,83 @@
	#ifndef RIVET_Cuts_HH
	#define RIVET_Cuts_HH
	#include <boost/smart_ptr.hpp>

	namespace Rivet {

	/// @internal Forward declaration of helper class. Not for end users.
	class CuttableBase;

	/// @internal Base class for cut objects.
	/// @note End users should always use the @ref Cut typedef instead.
	class CutBase;
	/// Main cut object
	typedef boost::shared_ptr<CutBase> Cut;

	class CutBase {
	public:
	/// Main work method.
	/// @internal Forwards the received object to @ref accept_,
	/// wrapped in the Cuttable converter
	template <typename ClassToCheck> bool accept(const ClassToCheck &) const;
	/// Comparison to another Cut
	virtual bool operator==(const Cut & c) const = 0;
	/// Default destructor
	virtual ~CutBase() {}
	protected:
	/// @internal Actual accept implementation, overloadable by various cut combiners
	virtual bool _accept(const CuttableBase &) const = 0;
	};

	// compare two cuts for equality, forwards to the cut-specific implementation
	inline bool operator==(const Cut & a, const Cut & b) { return *a == b; }

	/// Namespace used for ambiguous identifiers.
	namespace Cuts {
	/// Available categories of cut objects
	enum Quantity { pT=0, pt=0, Et=1, et=1, mass, rap, absrap, eta, abseta, phi };
	/// Fully open cut singleton, accepts everything
	const Cut & open();

	/// @name Shortcuts for common cuts
	//@{
	Cut range(Quantity, double m, double n);
	inline Cut etaIn(double m, double n) { return range(eta,m,n); }
	inline Cut rapIn(double m, double n) { return range(rap,m,n); }
	inline Cut ptIn(double m, double n) { return range(pT,m,n); }
	inline Cut etIn(double m, double n) { return range(Et,m,n); }
	inline Cut massIn(double m, double n) { return range(mass,m,n); }
	/// @todo Can do more here: absetaLess, ptGtr, ...?
	//@}
	}

	/// @name Cut constructors
	//@{
	Cut operator < (Cuts::Quantity, double);
	Cut operator > (Cuts::Quantity, double);
	Cut operator <= (Cuts::Quantity, double);
	Cut operator >= (Cuts::Quantity, double);

	/// @internal Overload helpers for integer arguments
	//@{
	inline Cut operator < (Cuts::Quantity qty, int i) { return qty < double(i); }
	inline Cut operator > (Cuts::Quantity qty, int i) { return qty > double(i); }
	inline Cut operator <= (Cuts::Quantity qty, int i) { return qty <= double(i); }
	inline Cut operator >= (Cuts::Quantity qty, int i) { return qty >= double(i); }
	//@}

	//@}

	/// @name Cut combiners
	//@{
	/// Logical AND operation on two cuts
	- Cut operator & (const Cut aptr, const Cut bptr);
	+ Cut operator & (const Cut & aptr, const Cut & bptr);
	/// Logical OR operation on two cuts
	- Cut operator \| (const Cut aptr, const Cut bptr);
	+ Cut operator \| (const Cut & aptr, const Cut & bptr);
	/// Logical NOT operation on a cut
	- Cut operator ~ (const Cut cptr);
	+ Cut operator ~ (const Cut & cptr);
	/// Logical XOR operation on two cuts
	- Cut operator ^ (const Cut aptr, const Cut bptr);
	+ Cut operator ^ (const Cut & aptr, const Cut & bptr);
	//@}
	}

	#endif
	diff --git a/include/Rivet/Projections/FinalState.hh b/include/Rivet/Projections/FinalState.hh
	--- a/include/Rivet/Projections/FinalState.hh
	+++ b/include/Rivet/Projections/FinalState.hh
	@@ -1,51 +1,51 @@
	// -- C++ --
	#ifndef RIVET_FinalState_HH
	#define RIVET_FinalState_HH

	#include "Rivet/Projections/ParticleFinder.hh"
	#include "Rivet/Cuts.hh"

	namespace Rivet {


	/// @brief Project out all final-state particles in an event.
	/// Probably the most important projection in Rivet!
	class FinalState : public ParticleFinder {
	public:

	/// @name Standard constructors and destructors.
	//@{
	/// @deprecated Keep for backwards compatibility for now
	/// The default constructor. May specify the minimum and maximum
	/// pseudorapidity \f$ \eta \f$ and the min \f$ p_T \f$ (in GeV).
	FinalState(double mineta,
	double maxeta,
	double minpt = 0.0);

	/// Construction using Cuts object
	- FinalState(Cut c = Cuts::open());
	+ FinalState(const Cut & c = Cuts::open());

	/// Clone on the heap.
	virtual const Projection* clone() const {
	return new FinalState(*this);
	}

	//@}


	/// Apply the projection to the event.
	virtual void project(const Event& e);

	/// Compare projections.
	virtual int compare(const Projection& p) const;

	/// Decide if a particle is to be accepted or not.
	/// @todo Rename to _accept or acceptFinal?
	virtual bool accept(const Particle& p) const;

	};


	}

	#endif
	diff --git a/include/Rivet/Projections/JetAlg.hh b/include/Rivet/Projections/JetAlg.hh
	--- a/include/Rivet/Projections/JetAlg.hh
	+++ b/include/Rivet/Projections/JetAlg.hh
	@@ -1,193 +1,193 @@
	// -- C++ --
	#ifndef RIVET_JetAlg_HH
	#define RIVET_JetAlg_HH

	#include "Rivet/Projection.hh"
	#include "Rivet/Projections/FinalState.hh"
	#include "Rivet/Projections/VisibleFinalState.hh"
	#include "Rivet/Particle.hh"
	#include "Rivet/Jet.hh"

	namespace Rivet {


	/// Abstract base class for projections which can return a set of {@link Jet}s.
	class JetAlg : public Projection {
	public:

	/// Constructor
	JetAlg(const FinalState& fs);

	JetAlg() {};

	/// Clone on the heap.
	virtual const Projection* clone() const = 0;

	/// Destructor
	virtual ~JetAlg() { }

	/// @brief Include invisible particles in jet construction.
	/// The default behaviour is that jets are only constructed from visible
	/// (i.e. charged under an SM gauge group) particles. Some jet studies,
	/// including those from ATLAS, use a definition in which neutrinos from hadron
	/// decays are included (via MC correction) in the experimental jet definition.
	/// Setting this flag to true avoids the automatic restriction to a VisibleFinalState.
	void useInvisibles(bool useinvis=true) {
	_useInvisibles = useinvis;
	}


	/// Get jets in no guaranteed order, with optional cuts on \f$ p_\perp \f$ and rapidity.
	/// @note Returns a copy rather than a reference, due to cuts
	/// @todo Can't this be a const Cut& arg?
	- virtual Jets jets(Cut c=Cuts::open()) const {
	+ virtual Jets jets(const Cut & c = Cuts::open()) const {
	const Jets rawjets = _jets(0.0); // arg means no pT cut
	// Just return a copy of rawjets if the cut is open
	if (c == Cuts::open()) return rawjets;
	// If there is a non-trivial cut...
	Jets rtn;
	rtn.reserve(size());
	foreach (const Jet& j, rawjets)
	if (c->accept(j)) rtn.push_back(j);
	return rtn;
	}

	/// Get the jets, ordered by supplied sorting function object, with optional cuts on \f$ p_\perp \f$ and rapidity.
	/// @note Returns a copy rather than a reference, due to cuts and sorting
	template <typename F>
	- Jets jets(F sorter, Cut c=Cuts::open()) const {
	+ Jets jets(F sorter, const Cut & c = Cuts::open()) const {
	/// @todo Will the vector be efficiently std::move'd by value through this function chain?
	return sortBy(jets(c), sorter);
	}

	/// Get the jets, ordered by supplied sorting function object, with optional cuts on \f$ p_\perp \f$ and rapidity.
	/// @note Returns a copy rather than a reference, due to cuts and sorting
	template <typename F>
	- Jets jets(Cut c, F sorter) const {
	+ Jets jets(const Cut & c , F sorter) const {
	/// @todo Will the vector be efficiently std::move'd by value through this function chain?
	return sortBy(jets(c), sorter);
	}


	/// Get the jets, ordered by \f$ p_T \f$, with optional cuts.
	///
	/// @note Returns a copy rather than a reference, due to cuts and sorting
	///
	/// This is a very common use-case, so is available as syntatic sugar for jets(c, cmpMomByPt).
	/// @todo The other sorted accessors should be removed in a cleanup.
	- Jets jetsByPt(Cut c=Cuts::open()) const {
	+ Jets jetsByPt(const Cut & c = Cuts::open()) const {
	return jets(c, cmpMomByPt);
	}


	/// @name Old sorted jet accessors
	/// @deprecated Use the versions with sorter function arguments
	//@{

	/// Get the jets, ordered by \f$ \|p\| \f$, with optional cuts on \f$ p_\perp \f$ and rapidity.
	/// @note Returns a copy rather than a reference, due to cuts and sorting
	/// @deprecated Use the version with a sorter function argument.
	- Jets jetsByP(Cut c=Cuts::open()) const {
	+ Jets jetsByP(const Cut & c = Cuts::open()) const {
	return jets(c, cmpMomByP);
	}

	/// Get the jets, ordered by \f$ E \f$, with optional cuts on \f$ p_\perp \f$ and rapidity.
	/// @note Returns a copy rather than a reference, due to cuts and sorting
	/// @deprecated Use the version with a sorter function argument.
	- Jets jetsByE(Cut c=Cuts::open()) const {
	+ Jets jetsByE(const Cut & c = Cuts::open()) const {
	return jets(c, cmpMomByE);
	}

	/// Get the jets, ordered by \f$ E_T \f$, with optional cuts on \f$ p_\perp \f$ and rapidity.
	/// @note Returns a copy rather than a reference, due to cuts and sorting
	/// @deprecated Use the version with a sorter function argument.
	- Jets jetsByEt(Cut c=Cuts::open()) const {
	+ Jets jetsByEt(const Cut & c = Cuts::open()) const {
	return jets(c, cmpMomByEt);
	}

	//@}


	/// @name Old jet accessors
	/// @deprecated Use the versions with Cut arguments
	//@{

	/// Get jets in no guaranteed order, with optional cuts on \f$ p_\perp \f$ and rapidity.
	///
	/// @deprecated Use the version with a Cut argument
	/// @note Returns a copy rather than a reference, due to cuts
	Jets jets(double ptmin, double ptmax=MAXDOUBLE,
	double rapmin=-MAXDOUBLE, double rapmax=MAXDOUBLE,
	RapScheme rapscheme=PSEUDORAPIDITY) const {
	if (rapscheme == PSEUDORAPIDITY) {
	return jets((Cuts::pT >= ptmin) & (Cuts::pT < ptmax) & (Cuts::rapIn(rapmin, rapmax)));
	} else if (rapscheme == RAPIDITY) {
	return jets((Cuts::pT >= ptmin) & (Cuts::pT < ptmax) & (Cuts::etaIn(rapmin, rapmax)));
	}
	throw LogicError("Unknown rapidity scheme. This shouldn't be possible!");
	}

	/// Get the jets, ordered by \f$ p_T \f$, with a cut on \f$ p_\perp \f$.
	///
	/// @deprecated Use the version with a Cut argument
	/// @note Returns a copy rather than a reference, due to cuts and sorting
	///
	/// This is a very common use-case, so is available as syntatic sugar for jets(Cuts::pT >= ptmin, cmpMomByPt).
	/// @todo The other sorted accessors should be removed in a cleanup.
	Jets jetsByPt(double ptmin) const {
	return jets(Cuts::pT >= ptmin, cmpMomByPt);
	}

	//@}


	protected:

	/// @brief Internal pure virtual method for getting jets in no guaranteed order.
	/// An optional cut on min \f$ p_\perp \f$ is applied in this function, since that is
	/// directly supported by FastJet and it seems a shame to not make use of that. But
	/// all other jet cuts are applied at the @c ::jets() function level.
	/// @todo Remove the ptmin cut
	virtual Jets _jets(double ptmin=0) const = 0;


	public:

	/// Number of jets.
	virtual size_t size() const = 0;
	/// Determine if the jet collection is empty.
	bool empty() const { return size() != 0; }

	/// Clear the projection.
	virtual void reset() = 0;

	typedef Jet entity_type;
	typedef Jets collection_type;

	/// Template-usable interface common to FinalState.
	collection_type entities() const { return jets(); }

	/// Do the calculation locally (no caching).
	virtual void calc(const Particles& constituents, const Particles& tagparticles=Particles()) = 0;


	protected:

	/// Perform the projection on the Event.
	virtual void project(const Event& e) = 0;

	/// Compare projections.
	virtual int compare(const Projection& p) const = 0;


	protected:

	/// Flag to determine whether or not the VFS wrapper is to be used.
	bool _useInvisibles;

	};


	}

	#endif
	diff --git a/include/Rivet/Projections/ParticleFinder.hh b/include/Rivet/Projections/ParticleFinder.hh
	--- a/include/Rivet/Projections/ParticleFinder.hh
	+++ b/include/Rivet/Projections/ParticleFinder.hh
	@@ -1,190 +1,190 @@
	// -- C++ --
	#ifndef RIVET_ParticleFinder_HH
	#define RIVET_ParticleFinder_HH

	#include "Rivet/Projection.hh"
	#include "Rivet/Cuts.hh"

	namespace Rivet {


	/// @brief Base class for projections which return subsets of an event's particles
	class ParticleFinder : public Projection {
	public:

	/// @name Object lifetime management
	//@{

	/// Construction using Cuts object
	- ParticleFinder(Cut c = Cuts::open()) : _cuts(c), _theParticles() {}
	+ ParticleFinder(const Cut & c = Cuts::open()) : _cuts(c), _theParticles() {}

	/// Virtual destructor for inheritance
	virtual ~ParticleFinder() {}

	/// Clone on the heap.
	virtual const Projection* clone() const = 0;

	//@}


	/// @name Particle accessors
	//@{

	/// Get the final-state particles in no particular order, with no cuts.
	virtual const Particles& particles() const { return _theParticles; }

	/// Access the projected final-state particles.
	size_t size() const { return particles().size(); }

	/// Is this final state empty?
	bool empty() const { return particles().empty(); }
	/// @deprecated Is this final state empty?
	bool isEmpty() const { return particles().empty(); }


	/// @brief Get the final-state particles, with optional cuts.
	/// @note Returns a copy rather than a reference, due to cuts
	/// @todo Can't this be a const Cut& arg?
	- Particles particles(Cut c) const {
	+ Particles particles(const Cut & c) const {
	// Just return a copy of particles() if the cut is open
	if (c == Cuts::open()) return particles();
	// If there is a non-trivial cut...
	Particles rtn;
	rtn.reserve(size());
	foreach (const Particle& p, particles())
	if (c->accept(p)) rtn.push_back(p);
	return rtn;
	}

	/// Get the final-state particles, ordered by supplied sorting function object.
	/// @note Returns a copy rather than a reference, due to cuts and sorting
	/// @todo Can't this be a const Cut& arg?
	/// @todo Use a std::function instead of typename F?
	template <typename F>
	- Particles particles(F sorter, Cut c=Cuts::open()) const {
	+ Particles particles(F sorter, const Cut & c=Cuts::open()) const {
	/// @todo Will the vector be efficiently std::move'd by value through this function chain?
	return sortBy(particles(c), sorter);
	}

	/// Get the final-state particles, ordered by supplied sorting function object.
	/// @note Returns a copy rather than a reference, due to cuts and sorting
	/// @todo Can't this be a const Cut& arg?
	/// @todo Use a std::function instead of typename F?
	template <typename F>
	- Particles particles(Cut c, F sorter) const {
	+ Particles particles(const Cut & c, F sorter) const {
	/// @todo Will the vector be efficiently std::move'd by value through this function chain?
	return sortBy(particles(c), sorter);
	}

	/// Get the final-state particles, ordered by decreasing \f$ p_T \f$ and with optional cuts.
	///
	/// This is a very common use-case, so is available as syntatic sugar for particles(c, cmpMomByPt).
	- Particles particlesByPt(Cut c=Cuts::open()) const {
	+ Particles particlesByPt(const Cut & c=Cuts::open()) const {
	return particles(c, cmpMomByPt);
	}

	/// Get the final-state particles, ordered by decreasing \f$ p_T \f$ and with a cut on minimum \f$ p_T \f$.
	///
	/// This is a very common use-case, so is available as syntatic sugar for particles(Cuts::pT >= ptmin, cmpMomByPt).
	Particles particlesByPt(double ptmin) const {
	return particles(Cuts::pT >= ptmin, cmpMomByPt);
	}


	/// @name Little-used sorted accessors
	/// @deprecated Use the versions with a sorter function argument
	//@{

	/// Get the final-state particles, ordered by decreasing \f$ p \f$.
	/// @todo Remove, since there is the templated method or sortByX methods available for these unusual cases?
	/// @deprecated Use the version with a sorter function argument
	- Particles particlesByP(Cut c=Cuts::open()) const {
	+ Particles particlesByP(const Cut & c=Cuts::open()) const {
	return particles(c, cmpMomByP);
	}

	/// Get the final-state particles, ordered by decreasing \f$ E \f$.
	/// @todo Remove, since there is the templated method or sortByX methods available for these unusual cases?
	/// @deprecated Use the version with a sorter function argument
	- Particles particlesByE(Cut c=Cuts::open()) const {
	+ Particles particlesByE(const Cut & c=Cuts::open()) const {
	return particles(c, cmpMomByE);
	}

	/// Get the final-state particles, ordered by decreasing \f$ E_T \f$.
	/// @todo Remove, since there is the templated method or sortByX methods available for these unusual cases?
	/// @deprecated Use the version with a sorter function argument
	- Particles particlesByEt(Cut c=Cuts::open()) const {
	+ Particles particlesByEt(const Cut & c=Cuts::open()) const {
	return particles(c, cmpMomByEt);
	}

	/// Get the final-state particles, ordered by increasing \f$ \eta \f$.
	/// @todo Remove, since there is the templated method or sortByX methods available for these unusual cases?
	/// @deprecated Use the version with a sorter function argument
	- Particles particlesByEta(Cut c=Cuts::open()) const {
	+ Particles particlesByEta(const Cut & c=Cuts::open()) const {
	return particles(c, cmpMomByAscPseudorapidity);
	}

	/// Get the final-state particles, ordered by increasing \f$ \|\eta\| \f$.
	/// @todo Remove, since there is the templated method or sortByX methods available for these unusual cases?
	/// @deprecated Use the version with a sorter function argument
	- Particles particlesByModEta(Cut c=Cuts::open()) const {
	+ Particles particlesByModEta(const Cut & c=Cuts::open()) const {
	return particles(c, cmpMomByAscAbsPseudorapidity);
	}

	/// Get the final-state particles, ordered by increasing \f$ y \f$.
	/// @todo Remove, since there is the templated method or sortByX methods available for these unusual cases?
	/// @deprecated Use the version with a sorter function argument
	- Particles particlesByRapidity(Cut c=Cuts::open()) const {
	+ Particles particlesByRapidity(const Cut & c=Cuts::open()) const {
	return particles(c, cmpMomByAscRapidity);
	}

	/// Get the final-state particles, ordered by increasing \f$ \|y\| \f$.
	/// @todo Remove, since there is the templated method or sortByX methods available for these unusual cases?
	/// @deprecated Use the version with a sorter function argument
	- Particles particlesByModRapidity(Cut c=Cuts::open()) const {
	+ Particles particlesByModRapidity(const Cut & c=Cuts::open()) const {
	return particles(c, cmpMomByAscAbsRapidity);
	}

	//@}

	//@}


	/// @todo Replace with cuts() accessor
	///virtual Cut cuts() const { return _cuts; }


	/// @name For JetAlg compatibility
	//@{

	typedef Particle entity_type;
	typedef Particles collection_type;

	/// Template-usable interface common to JetAlg.
	const collection_type& entities() const {
	return particles();
	}

	//@}


	protected:

	/// Apply the projection to the event.
	virtual void project(const Event& e) = 0;

	/// Compare projections.
	virtual int compare(const Projection& p) const;

	protected:
	/// The applicable cuts
	Cut _cuts;

	/// The final-state particles.
	Particles _theParticles;

	};


	}

	#endif
	diff --git a/include/Rivet/Projections/PromptFinalState.hh b/include/Rivet/Projections/PromptFinalState.hh
	--- a/include/Rivet/Projections/PromptFinalState.hh
	+++ b/include/Rivet/Projections/PromptFinalState.hh
	@@ -1,66 +1,66 @@
	// -- C++ --
	#ifndef RIVET_PromptFinalState_HH
	#define RIVET_PromptFinalState_HH

	#include "Rivet/Projections/FinalState.hh"

	namespace Rivet {


	/// @brief Find final state particles directly connected to the hard process.
	///
	/// The definition of "prompt" used in Rivet is that from high-scale physics, i.e.
	/// particles directly connected to the hard process in an interaction, regardless
	/// of realistic reconstructibility of displaced vertices, etc. By construction
	/// hadrons cannot be considered prompt as they will be colour connected to other
	/// parts of the event through non-perturbative effects: this projection can
	/// return electrons, muons, photons, and exotic particles which do not have a
	/// hadron in their post-hadronization ancestor chain. Flags exist to choose
	/// whether intermediate tau or muon decays invalidate a particle's promptness.
	///
	/// @todo Decide how to treat brem photons off prompt leptons -- are they also prompt? "Decay" does not change the lepton PID...
	class PromptFinalState : public FinalState {
	public:

	/// @name Constructors
	//@{
	// Final State
	PromptFinalState(const FinalState& fsp);

	/// Cut constructor.
	- PromptFinalState(Cut c);
	+ PromptFinalState(const Cut & c);

	/// Clone on the heap.
	virtual const Projection* clone() const {
	return new PromptFinalState(*this);
	}
	//@}

	/// Accept particles from decays of prompt muons as themselves being prompt?
	void acceptMuonDecays(bool acc=true) { _acceptMuDecays = acc; }
	/// Accept particles from decays of prompt taus as themselves being prompt?
	void acceptTauDecays(bool acc=true) { _acceptTauDecays = acc; }

	/// Decide if a given particle is prompt based on set definition flags
	/// @todo Move into ParticleUtils / MCUtils
	/// @note This one doesn't make any judgements about final-stateness
	bool isPrompt(const Particle& p) const;

	protected:

	/// Apply the projection on the supplied event.
	void project(const Event& e);

	/// Compare projections.
	int compare(const Projection& p) const;

	private:

	bool _acceptMuDecays, _acceptTauDecays;

	};

	}


	#endif
	diff --git a/include/Rivet/Projections/WFinder.hh b/include/Rivet/Projections/WFinder.hh
	--- a/include/Rivet/Projections/WFinder.hh
	+++ b/include/Rivet/Projections/WFinder.hh
	@@ -1,157 +1,157 @@
	// -- C++ --
	#ifndef RIVET_WFinder_HH
	#define RIVET_WFinder_HH

	#include "Rivet/Projections/FinalState.hh"

	namespace Rivet {


	/// @brief Convenience finder of leptonically decaying Ws
	///
	/// Chain together different projections as convenience for finding W's
	/// from two leptons in the final state, including photon clustering.
	///
	/// @todo Inherit directly from ParticleFinder, not FinalState
	class WFinder : public FinalState {
	public:

	enum ClusterPhotons { NOCLUSTER=0, CLUSTERNODECAY=1, CLUSTERALL };
	enum PhotonTracking { NOTRACK=0, TRACK=1 };
	enum MassWindow { MASS=0, TRANSMASS=1 };

	/// @name Constructors
	//@{

	/// Constructor taking cuts object
	/// @param inputfs Input final state
	/// @param cuts charged lepton cuts
	/// @param pid type of the charged lepton
	/// @param minmass,maxmass (transverse) mass window
	/// @param missingET minimal amount of missing ET (neutrinos) required
	/// @param dRmax maximum dR of photons around charged lepton to take into account
	/// for W reconstruction (only relevant if one of the following are true)
	/// @param clusterPhotons whether such photons are supposed to be
	/// clustered to the lepton object and thus W mom
	/// @param trackPhotons whether such photons should be added to _theParticles
	/// @param masstype whether mass window should be applied using m or mT
	WFinder(const FinalState& inputfs,
	- Cut cuts,
	+ const Cut & cuts,
	PdgId pid,
	double minmass, double maxmass,
	double missingET,
	double dRmax=0.1,
	ClusterPhotons clusterPhotons=CLUSTERNODECAY,
	PhotonTracking trackPhotons=NOTRACK,
	MassWindow masstype=MASS,
	double masstarget=80.4*GeV);

	/// Clone on the heap.
	virtual const Projection* clone() const {
	return new WFinder(*this);
	}

	//@}


	/// Access to the found bosons
	///
	/// @note Currently either 0 or 1 boson can be found.
	const Particles& bosons() const { return _bosons; }

	/// Access to the found boson (assuming it exists)
	const Particle& boson() const { return _bosons[0]; }


	/// Access to the W constituent clustered leptons
	///
	/// @note Either size 0 if no boson was found or 1 if one boson was found
	const Particles& constituentLeptons() const { return _constituentLeptons; }

	/// Access to the W constituent clustered lepton (assuming it exists)
	const Particle& constituentLepton() const { return _constituentLeptons[0]; }


	/// Access to the W constituent neutrinos
	///
	/// @note Either size 0 if no boson was found or 1 if one boson was found
	const Particles& constituentNeutrinos() const { return _constituentNeutrinos; }

	/// Access to the W constituent neutrinos
	const Particle& constituentNeutrino() const { return _constituentNeutrinos[0]; }


	/// Access to the particles other than the W leptons and clustered photons
	///
	/// Useful for e.g. input to a jet finder
	const FinalState& remainingFinalState() const;

	/// @brief Calculate the transverse mass of the W, from the charged lepton and neutrino
	///
	/// Defined as sqrt(2 pT_l pT_nu (1.0 - cos(dphi_lnu))). Return -1 if no boson found.
	double mT() const {
	if (bosons().empty()) return -1;
	const FourMomentum& l = constituentLeptons()[0].momentum();
	const FourMomentum& nu = constituentNeutrinos()[0].momentum();
	return sqrt( 2 * l.pT() * nu.pT() * (1 - cos(deltaPhi(l, nu))) );
	}


	protected:

	/// Apply the projection on the supplied event.
	void project(const Event& e);

	/// Compare projections.
	int compare(const Projection& p) const;


	public:

	/// Clear the projection
	void clear() {
	_theParticles.clear();
	_bosons.clear();
	_constituentLeptons.clear();
	_constituentNeutrinos.clear();
	}

	private:

	/// Transverse mass cuts
	double _minmass, _maxmass, _masstarget;
	bool _useTransverseMass;

	/// Missing ET cut
	double _etMiss;

	/// Switch for tracking of photons (whether to add them to _theParticles)
	/// This is relevant when the ZFinder::_theParticles are to be excluded
	/// from e.g. the input to a jet finder, to specify whether the clustered
	/// photons are to be excluded as well.
	/// (Yes, some experiments make a difference between clusterPhotons and
	/// trackPhotons!)
	PhotonTracking _trackPhotons;

	/// Lepton flavour
	PdgId _pid;

	/// Neutrino flavour
	PdgId _nu_pid;

	/// list of found bosons (currently either 0 or 1)
	Particles _bosons;

	/// Constituent leptons (currently either 0 or 1)
	Particles _constituentLeptons;

	/// Constituent neutrinos (currently either 0 or 1)
	Particles _constituentNeutrinos;

	};


	}


	#endif
	diff --git a/include/Rivet/Projections/ZFinder.hh b/include/Rivet/Projections/ZFinder.hh
	--- a/include/Rivet/Projections/ZFinder.hh
	+++ b/include/Rivet/Projections/ZFinder.hh
	@@ -1,131 +1,131 @@
	// -- C++ --
	#ifndef RIVET_ZFinder_HH
	#define RIVET_ZFinder_HH

	#include "Rivet/Projections/FinalState.hh"

	namespace Rivet {


	/// @brief Convenience finder of leptonically decaying Zs
	///
	/// Chain together different projections as convenience for finding Z's
	/// from two leptons in the final state, including photon clustering.
	///
	/// @todo Inherit directly from ParticleFinder, not FinalState
	class ZFinder : public FinalState {
	public:

	enum ClusterPhotons { NOCLUSTER=0, CLUSTERNODECAY=1, CLUSTERALL };
	enum PhotonTracking { NOTRACK=0, TRACK=1 };

	/// @name Constructors
	//@{

	/// Constructor taking cuts object
	/// @param inputfs Input final state
	/// @param cuts lepton cuts
	/// @param pid type of the leptons
	/// @param minmass,maxmass mass window
	/// @param dRmax maximum dR of photons around leptons to take into account
	/// for Z reconstruction (only relevant if one of the following are true)
	/// @param clusterPhotons whether such photons are supposed to be
	/// clustered to the lepton objects and thus Z mom
	/// @param trackPhotons whether such photons should be added to _theParticles
	/// (cf. _trackPhotons)
	ZFinder(const FinalState& inputfs,
	- Cut cuts,
	+ const Cut & cuts,
	PdgId pid,
	double minmass, double maxmass,
	double dRmax=0.1,
	ClusterPhotons clusterPhotons=CLUSTERNODECAY,
	PhotonTracking trackPhotons=NOTRACK,
	double masstarget=91.2*GeV);

	/// Clone on the heap.
	virtual const Projection* clone() const {
	return new ZFinder(*this);
	}

	//@}


	/// Access to the found bosons
	///
	/// @note Currently either 0 or 1 boson can be found.
	const Particles& bosons() const { return _bosons; }

	/// Access to the found boson (assuming it exists).
	const Particle boson() const { return _bosons[0]; }

	/// Access to the Z constituent clustered leptons
	///
	/// For example, to make more fine-grained cuts on the clustered leptons.
	/// The positive charge constituent is first in the list (if not empty), and
	/// the negative one second.
	const Particles& constituents() const { return _constituents; }

	/// Access to the Z constituent clustered leptons, sorted by a comparison functor
	///
	/// Unlike the no-arg version, this returns by value (i.e. is less efficient)
	template <typename CMP>
	Particles constituents(const CMP& cmp) const {
	Particles rtn = constituents();
	std::sort(rtn.begin(), rtn.end(), cmp);
	return rtn;
	}

	/// Access to the particles other than the Z leptons and clustered photons
	///
	/// Useful for e.g. input to a jet finder
	const FinalState& remainingFinalState() const;


	protected:

	/// Apply the projection on the supplied event.
	void project(const Event& e);

	/// Compare projections.
	int compare(const Projection& p) const;


	public:

	/// Clear the projection
	void clear() {
	_theParticles.clear();
	_bosons.clear();
	_constituents.clear();
	}


	private:
	/// Mass cuts to apply to clustered leptons (cf. InvMassFinalState)
	double _minmass, _maxmass, _masstarget;

	/// Switch for tracking of photons (whether to add them to _theParticles)
	/// This is relevant when the ZFinder::_theParticles are to be excluded
	/// from e.g. the input to a jet finder, to specify whether the clustered
	/// photons are to be excluded as well.
	/// (Yes, some experiments make a difference between clusterPhotons and
	/// trackPhotons!)
	PhotonTracking _trackPhotons;

	/// Lepton flavour
	PdgId _pid;

	/// list of found bosons (currently either 0 or 1)
	Particles _bosons;

	/// Clustered leptons
	vector<Particle> _constituents;

	};


	}



	#endif
	diff --git a/src/Core/Cuts.cc b/src/Core/Cuts.cc
	--- a/src/Core/Cuts.cc
	+++ b/src/Core/Cuts.cc
	@@ -1,375 +1,375 @@
	#include <Rivet/Cuts.hh>
	#include <boost/make_shared.hpp>

	// headers for converters
	#include <Rivet/Particle.hh>
	#include <Rivet/Jet.hh>
	#include <Rivet/Math/Vectors.hh>
	#include <fastjet/PseudoJet.hh>
	#include <HepMC/SimpleVector.h>

	// todo Sort out what can go into anonymous namespace{}

	namespace Rivet {

	// Base for all wrapper classes that translate ClassToCheck to Cuttable
	class CuttableBase {
	public:
	virtual double getValue(Cuts::Quantity) const = 0;
	virtual ~CuttableBase() {}
	};


	// Cuttables can be directly passed to @ref _accept
	template <>
	bool CutBase::accept<CuttableBase>(const CuttableBase & t) const {
	return _accept(t);
	}

	// bool operator==(const Cut & a, const Cut & b) {
	// return *a == b;
	// }


	// Open Cut singleton
	class Open_Cut : public CutBase {
	public:
	bool operator==(const Cut & c) const {
	shared_ptr<Open_Cut> cc = dynamic_pointer_cast<Open_Cut>(c);
	return bool(cc);
	}
	protected:
	// open cut accepts everything
	bool _accept(const CuttableBase &) const { return true; }
	};

	const Cut & Cuts::open() {
	// only ever need one static open cut object
	static const Cut open = boost::make_shared<Open_Cut>();
	return open;
	}





	// Cut constructor for >=
	class Cut_GtrEq : public CutBase {
	public:
	Cut_GtrEq(const Cuts::Quantity qty, const double low) : qty_(qty), low_(low) {}
	bool operator==(const Cut & c) const {
	shared_ptr<Cut_GtrEq> cc = dynamic_pointer_cast<Cut_GtrEq>(c);
	return cc && qty_ == cc->qty_ && low_ == cc->low_;
	}
	protected:
	bool _accept(const CuttableBase & o) const { return o.getValue(qty_) >= low_; }
	private:
	Cuts::Quantity qty_;
	double low_;
	};

	// Cut constructor for <
	class Cut_Less : public CutBase {
	public:
	Cut_Less(const Cuts::Quantity qty, const double high) : qty_(qty), high_(high) {}
	bool operator==(const Cut & c) const {
	shared_ptr<Cut_Less> cc = dynamic_pointer_cast<Cut_Less>(c);
	return cc && qty_ == cc->qty_ && high_ == cc->high_;
	}
	protected:
	bool _accept(const CuttableBase & o) const { return o.getValue(qty_) < high_; }
	private:
	Cuts::Quantity qty_;
	double high_;
	};

	// Cut constructor for >=
	class Cut_Gtr : public CutBase {
	public:
	Cut_Gtr(const Cuts::Quantity qty, const double low) : qty_(qty), low_(low) {}
	bool operator==(const Cut & c) const {
	shared_ptr<Cut_Gtr> cc = dynamic_pointer_cast<Cut_Gtr>(c);
	return cc && qty_ == cc->qty_ && low_ == cc->low_;
	}
	protected:
	bool _accept(const CuttableBase & o) const { return o.getValue(qty_) > low_; }
	private:
	Cuts::Quantity qty_;
	double low_;
	};

	// Cut constructor for <
	class Cut_LessEq : public CutBase {
	public:
	Cut_LessEq(const Cuts::Quantity qty, const double high) : qty_(qty), high_(high) {}
	bool operator==(const Cut & c) const {
	shared_ptr<Cut_LessEq> cc = dynamic_pointer_cast<Cut_LessEq>(c);
	return cc && qty_ == cc->qty_ && high_ == cc->high_;
	}
	protected:
	bool _accept(const CuttableBase & o) const { return o.getValue(qty_) <= high_; }
	private:
	Cuts::Quantity qty_;
	double high_;
	};


	template <typename T>
	inline Cut make_cut(T t) {
	return boost::make_shared<T>(t);
	}

	Cut operator < (Cuts::Quantity qty, double n) {
	return make_cut(Cut_Less(qty, n));
	}

	Cut operator >= (Cuts::Quantity qty, double n) {
	return make_cut(Cut_GtrEq(qty, n));
	}

	Cut operator <= (Cuts::Quantity qty, double n) {
	return make_cut(Cut_LessEq(qty, n));
	}

	Cut operator > (Cuts::Quantity qty, double n) {
	return make_cut(Cut_Gtr(qty, n));
	}

	Cut Cuts::range(Cuts::Quantity qty, double m, double n) {
	if (m > n) swap(m,n);
	return (qty >= m) & (qty < n);
	}


	//////////////
	/// Combiners

	/// AND, OR, NOT, and XOR objects for combining cuts

	class CutsOr : public CutBase {
	public:
	- CutsOr(const Cut c1, const Cut c2) : cut1(c1), cut2(c2) {}
	+ CutsOr(const Cut & c1, const Cut & c2) : cut1(c1), cut2(c2) {}
	bool operator==(const Cut & c) const {
	shared_ptr<CutsOr> cc = dynamic_pointer_cast<CutsOr>(c);
	return cc && ( ( cut1 == cc->cut1 && cut2 == cc->cut2 )
	\|\| ( cut1 == cc->cut2 && cut2 == cc->cut1 ));
	}
	protected:
	bool _accept(const CuttableBase & o) const {
	return cut1->accept(o) \|\| cut2->accept(o);
	}
	private:
	const Cut cut1;
	const Cut cut2;
	};

	class CutsAnd : public CutBase {
	public:
	- CutsAnd(const Cut c1, const Cut c2) : cut1(c1), cut2(c2) {}
	+ CutsAnd(const Cut & c1, const Cut & c2) : cut1(c1), cut2(c2) {}
	bool operator==(const Cut & c) const {
	shared_ptr<CutsAnd> cc = dynamic_pointer_cast<CutsAnd>(c);
	return cc && ( ( cut1 == cc->cut1 && cut2 == cc->cut2 )
	\|\| ( cut1 == cc->cut2 && cut2 == cc->cut1 ));
	}
	protected:
	bool _accept(const CuttableBase & o) const {
	return cut1->accept(o) && cut2->accept(o);
	}
	private:
	const Cut cut1;
	const Cut cut2;
	};

	class CutInvert : public CutBase {
	public:
	- CutInvert(const Cut c1) : cut(c1) {}
	+ CutInvert(const Cut & c1) : cut(c1) {}
	bool operator==(const Cut & c) const {
	shared_ptr<CutInvert> cc = dynamic_pointer_cast<CutInvert>(c);
	return cc && cut == cc->cut;
	}
	protected:
	bool _accept(const CuttableBase & o) const {
	return !cut->accept(o);
	}
	private:
	const Cut cut;
	};

	class CutsXor : public CutBase {
	public:
	- CutsXor(const Cut c1, const Cut c2) : cut1(c1), cut2(c2) {}
	+ CutsXor(const Cut & c1, const Cut & c2) : cut1(c1), cut2(c2) {}
	bool operator==(const Cut & c) const {
	shared_ptr<CutsXor> cc = dynamic_pointer_cast<CutsXor>(c);
	return cc && ( ( cut1 == cc->cut1 && cut2 == cc->cut2 )
	\|\| ( cut1 == cc->cut2 && cut2 == cc->cut1 ));
	}
	protected:
	bool _accept(const CuttableBase & o) const {
	bool A_and_B = cut1->accept(o) && cut2->accept(o);
	bool A_or_B = cut1->accept(o) \|\| cut2->accept(o);
	return A_or_B && (! A_and_B);
	}
	private:
	const Cut cut1;
	const Cut cut2;
	};

	////////////
	///Operators

	- Cut operator & (const Cut aptr, const Cut bptr) {
	+ Cut operator & (const Cut & aptr, const Cut & bptr) {
	return make_cut(CutsAnd(aptr,bptr));
	}

	- Cut operator \| (const Cut aptr, const Cut bptr) {
	+ Cut operator \| (const Cut & aptr, const Cut & bptr) {
	return make_cut(CutsOr(aptr,bptr));
	}

	- Cut operator ~ (const Cut cptr) {
	+ Cut operator ~ (const Cut & cptr) {
	return make_cut(CutInvert(cptr));
	}

	- Cut operator ^ (const Cut aptr, const Cut bptr) {
	+ Cut operator ^ (const Cut & aptr, const Cut & bptr) {
	return make_cut(CutsXor(aptr,bptr));
	}

	///////////////////////
	/// Cuts

	// Non-functional Cuttable template class. Must be specialized below.
	template <typename T> class Cuttable;

	// Non-cuttables need to be wrapped into a Cuttable first
	#define SPECIALISE_ACCEPT(TYPENAME) \
	template <> \
	bool CutBase::accept<TYPENAME>(const TYPENAME & t) const { \
	return _accept(Cuttable<TYPENAME>(t)); \
	} \


	void qty_not_found() {
	throw Exception("Missing implementation for a Cuts::Quantity.");
	}


	template<>
	class Cuttable<Particle> : public CuttableBase {
	public:
	Cuttable(const Particle& p) : p_(p) {}
	double getValue(Cuts::Quantity qty) const {
	switch ( qty ) {
	case Cuts::pT: return p_.momentum().pT();
	case Cuts::Et: return p_.momentum().Et();
	case Cuts::mass: return p_.momentum().mass();
	case Cuts::rap: return p_.momentum().rapidity();
	case Cuts::absrap: return std::abs(p_.momentum().rapidity());
	case Cuts::eta: return p_.momentum().pseudorapidity();
	case Cuts::abseta: return std::abs(p_.momentum().pseudorapidity());
	case Cuts::phi: return p_.momentum().phi();
	default: qty_not_found();
	}
	return -999.;
	}
	private:
	const Particle & p_;
	};
	SPECIALISE_ACCEPT(Particle)


	template<>
	class Cuttable<FourMomentum> : public CuttableBase {
	public:
	Cuttable(const FourMomentum& fm) : fm_(fm) {}
	double getValue(Cuts::Quantity qty) const {
	switch ( qty ) {
	case Cuts::pT: return fm_.pT();
	case Cuts::Et: return fm_.Et();
	case Cuts::mass: return fm_.mass();
	case Cuts::rap: return fm_.rapidity();
	case Cuts::absrap: return std::abs(fm_.rapidity());
	case Cuts::eta: return fm_.pseudorapidity();
	case Cuts::abseta: return std::abs(fm_.pseudorapidity());
	case Cuts::phi: return fm_.phi();
	default: qty_not_found();
	}
	return -999.;
	}
	private:
	const FourMomentum & fm_;
	};
	SPECIALISE_ACCEPT(FourMomentum)


	template<>
	class Cuttable<Jet> : public CuttableBase {
	public:
	Cuttable(const Jet& jet) : jet_(jet) {}
	double getValue(Cuts::Quantity qty) const {
	switch ( qty ) {
	case Cuts::pT: return jet_.momentum().pT();
	case Cuts::Et: return jet_.momentum().Et();
	case Cuts::mass: return jet_.momentum().mass();
	case Cuts::rap: return jet_.momentum().rapidity();
	case Cuts::absrap: return std::abs(jet_.momentum().rapidity());
	case Cuts::eta: return jet_.momentum().pseudorapidity();
	case Cuts::abseta: return std::abs(jet_.momentum().pseudorapidity());
	case Cuts::phi: return jet_.momentum().phi();
	default: qty_not_found();
	}
	return -999.;
	}
	private:
	const Jet & jet_;
	};
	SPECIALISE_ACCEPT(Jet)


	template<>
	class Cuttable<fastjet::PseudoJet> : public CuttableBase {
	public:
	Cuttable(const fastjet::PseudoJet& pjet) : pjet_(pjet) {}
	double getValue(Cuts::Quantity qty) const {
	switch ( qty ) {
	case Cuts::pT: return pjet_.perp();
	case Cuts::Et: return pjet_.Et();
	case Cuts::mass: return pjet_.m();
	case Cuts::rap: return pjet_.rap();
	case Cuts::absrap: return std::abs(pjet_.rap());
	case Cuts::eta: return pjet_.eta();
	case Cuts::abseta: return std::abs(pjet_.eta());
	case Cuts::phi: return pjet_.phi();
	default: qty_not_found();
	}
	return -999.;
	}
	private:
	const fastjet::PseudoJet & pjet_;
	};
	SPECIALISE_ACCEPT(fastjet::PseudoJet)


	template<>
	class Cuttable<HepMC::FourVector> : public CuttableBase {
	public:
	Cuttable(const HepMC::FourVector & vec) : vec_(vec) {}
	double getValue(Cuts::Quantity qty) const {
	switch ( qty ) {
	case Cuts::pT: return vec_.perp();
	/// @todo case Cuts::Et: return vec_.perp();
	case Cuts::mass: return vec_.m();
	case Cuts::rap: return 0.5*std::log((vec_.t()+vec_.z())/(vec_.t()-vec_.z()));
	case Cuts::absrap: return std::abs(getValue(Cuts::rap));
	case Cuts::eta: return vec_.pseudoRapidity();
	case Cuts::abseta: return std::abs(vec_.pseudoRapidity());
	case Cuts::phi: return vec_.phi();
	default: qty_not_found();
	}
	return -999.;
	}
	private:
	const HepMC::FourVector & vec_;
	};
	SPECIALISE_ACCEPT(HepMC::FourVector)


	}
	diff --git a/src/Projections/FinalState.cc b/src/Projections/FinalState.cc
	--- a/src/Projections/FinalState.cc
	+++ b/src/Projections/FinalState.cc
	@@ -1,107 +1,107 @@
	// -- C++ --
	#include "Rivet/Projections/FinalState.hh"

	namespace Rivet {


	// @deprecated, keep for backwards compatibility for now.
	FinalState::FinalState(double mineta, double maxeta, double minpt)
	{
	using namespace Cuts;
	setName("FinalState");
	const bool openpt = isZero(minpt);
	const bool openeta = (mineta <= -MAXDOUBLE && maxeta >= MAXDOUBLE);
	MSG_TRACE("Check for open FS conditions:" << std::boolalpha
	<< " eta=" << openeta
	<< ", pt=" << openpt);
	if ( openpt && openeta ) {
	_cuts = open();
	}
	else {
	addProjection(FinalState(), "OpenFS");
	if ( openeta )
	_cuts = pT >= minpt;
	else if ( openpt )
	_cuts = etaIn(mineta,maxeta);
	else
	_cuts = etaIn(mineta,maxeta) & (pT >= minpt);
	}
	}


	- FinalState::FinalState(Cut c)
	+ FinalState::FinalState(const Cut & c)
	: ParticleFinder(c)
	{
	setName("FinalState");
	const bool open = ( c == Cuts::open() );
	MSG_TRACE("Check for open FS conditions: " << std::boolalpha << open);
	if ( !open ) {
	addProjection(FinalState(), "OpenFS");
	}
	}




	/// @todo HOW DO WE COMPARE CUTS OBJECTS?
	int FinalState::compare(const Projection& p) const {
	const FinalState& other = dynamic_cast<const FinalState&>(p);

	//MSG_TRACE("FS::compare: " << 1 << " " << this << " " << &p);
	// std::vector<std::pair<double, double> > eta1(_etaRanges);
	//std::vector<std::pair<double, double> > eta2(other._etaRanges);
	//std::sort(eta1.begin(), eta1.end());
	//std::sort(eta2.begin(), eta2.end());

	//MSG_TRACE("FS::compare: " << 2 << " " << this << " " << &p);
	//if (eta1 < eta2) return ORDERED;
	//else if (eta2 < eta1) return UNORDERED;

	//MSG_TRACE("FS::compare: " << 3 << " " << this << " " << &p);
	//return cmp(_ptmin, other._ptmin);
	return _cuts == other._cuts ? EQUIVALENT : UNDEFINED;
	}



	void FinalState::project(const Event& e) {
	_theParticles.clear();

	// Handle "open FS" special case
	if ( _cuts == Cuts::open() ) {
	//MSG_TRACE("Open FS processing: should only see this once per event ("
	// << e.genEvent().event_number() << ")");
	foreach (const GenParticle* p, Rivet::particles(e.genEvent())) {
	if (p->status() == 1) {
	//MSG_TRACE("FS GV = " << p->production_vertex());
	_theParticles.push_back(Particle(*p));
	}
	}
	return;
	}

	// If this is not itself the "open" FS, base the calculations on the open FS' results
	/// @todo In general, we'd like to calculate a restrictive FS based on the most restricted superset FS.
	const Particles allstable = applyProjection<FinalState>(e, "OpenFS").particles();
	foreach (const Particle& p, allstable) {
	const bool passed = accept(p);
	MSG_TRACE("Choosing: ID = " << p.pdgId()
	<< ", pT = " << p.pT()
	<< ", eta = " << p.eta()
	<< ": result = " << std::boolalpha << passed);
	if (passed) _theParticles.push_back(p);
	}
	//MSG_DEBUG("Number of final-state particles = " << _theParticles.size());
	}


	/// Decide if a particle is to be accepted or not.
	bool FinalState::accept(const Particle& p) const {
	// Not having s.c. == 1 should never happen!
	assert(p.genParticle() == NULL \|\| p.genParticle()->status() == 1);

	return _cuts->accept(p);
	}


	}
	diff --git a/src/Projections/PromptFinalState.cc b/src/Projections/PromptFinalState.cc
	--- a/src/Projections/PromptFinalState.cc
	+++ b/src/Projections/PromptFinalState.cc
	@@ -1,66 +1,66 @@
	// -- C++ --
	#include "Rivet/Projections/PromptFinalState.hh"

	namespace Rivet {


	PromptFinalState::PromptFinalState(const FinalState& fsp)
	: _acceptMuDecays(false), _acceptTauDecays(false)
	{
	setName("PromptFinalState");
	addProjection(fsp, "FS");
	}


	- PromptFinalState::PromptFinalState(Cut c)
	+ PromptFinalState::PromptFinalState(const Cut & c)
	: _acceptMuDecays(false), _acceptTauDecays(false)
	{
	setName("PromptFinalState");
	addProjection(FinalState(c), "FS");
	}

	int PromptFinalState::compare(const Projection& p) const {
	const PCmp fscmp = mkNamedPCmp(p, "FS");
	if (fscmp != EQUIVALENT) return fscmp;
	const PromptFinalState& other = dynamic_cast<const PromptFinalState&>(p);
	return cmp(_acceptMuDecays, other._acceptMuDecays) \|\| cmp(_acceptTauDecays, other._acceptTauDecays);
	}


	bool PromptFinalState::isPrompt(const Particle& p) const {
	if (p.genParticle() == NULL) return false; // no HepMC connection, give up! Throw UserError exception?
	/// @todo Shouldn't a const vertex be being returned? Ah, HepMC...
	GenVertex* prodVtx = p.genParticle()->production_vertex();
	if (prodVtx == NULL) return false; // orphaned particle, has to be assume false
	const pair<GenParticle, GenParticle> beams = prodVtx->parent_event()->beam_particles();

	/// @todo Would be nicer to be able to write this recursively up the chain, exiting as soon as a parton or string/cluster is seen
	foreach (const GenParticle* ancestor, Rivet::particles(prodVtx, HepMC::ancestors)) {
	const PdgId pid = ancestor->pdg_id();
	if (ancestor->status() != 2) continue; // no non-standard statuses or beams to be used in decision making
	if (ancestor == beams.first \|\| ancestor == beams.second) continue; // PYTHIA6 uses status 2 for beams, I think... (sigh)
	if (PID::isParton(pid)) continue; // PYTHIA6 also uses status 2 for some partons, I think... (sigh)
	if (PID::isHadron(pid)) return false; // prompt particles can't be from hadron decays
	if (abs(pid) == PID::TAU && p.abspid() != PID::TAU && !_acceptTauDecays) return false; // allow or ban particles from tau decays (permitting tau copies)
	if (abs(pid) == PID::MUON && p.abspid() != PID::MUON && !_acceptMuDecays) return false; // allow or ban particles from muon decays (permitting muon copies)
	}
	return true;
	}


	void PromptFinalState::project(const Event& e) {
	_theParticles.clear();

	const Particles& particles = applyProjection<FinalState>(e, "FS").particles();
	foreach (const Particle& p, particles)
	if (isPrompt(p)) _theParticles.push_back(p);
	MSG_DEBUG("Number of final state particles not from hadron decays = " << _theParticles.size());

	if (getLog().isActive(Log::TRACE)) {
	foreach (const Particle& p, _theParticles)
	MSG_TRACE("Selected: " << p.pid() << ", charge = " << p.charge());
	}
	}


	}
	diff --git a/src/Projections/WFinder.cc b/src/Projections/WFinder.cc
	--- a/src/Projections/WFinder.cc
	+++ b/src/Projections/WFinder.cc
	@@ -1,155 +1,155 @@
	// -- C++ --
	#include "Rivet/Projections/WFinder.hh"
	#include "Rivet/Projections/ChargedFinalState.hh"
	#include "Rivet/Projections/InvMassFinalState.hh"
	#include "Rivet/Projections/MissingMomentum.hh"
	#include "Rivet/Projections/MergedFinalState.hh"
	#include "Rivet/Projections/DressedLeptons.hh"
	#include "Rivet/Projections/VetoedFinalState.hh"

	namespace Rivet {

	WFinder::WFinder(const FinalState& inputfs,
	- Cut fsCut,
	+ const Cut & fsCut,
	PdgId pid,
	double minmass, double maxmass,
	double missingET,
	double dRmax,
	ClusterPhotons clusterPhotons,
	PhotonTracking trackPhotons,
	MassWindow masstype,
	double masstarget) {
	setName("WFinder");

	_minmass = minmass;
	_maxmass = maxmass;
	_masstarget = masstarget;
	_pid = pid;
	_trackPhotons = trackPhotons;
	_useTransverseMass = (masstype == MASS);

	// Check that the arguments are legal
	assert(abs(_pid) == PID::ELECTRON \|\| abs(_pid) == PID::MUON);
	_nu_pid = abs(_pid) + 1;
	assert(abs(_nu_pid) == PID::NU_E \|\| abs(_nu_pid) == PID::NU_MU);

	// Don't make pT or eta cuts on the neutrino
	IdentifiedFinalState neutrinos(inputfs);
	neutrinos.acceptNeutrinos();
	addProjection(neutrinos, "Neutrinos");

	// Lepton clusters
	IdentifiedFinalState bareleptons(inputfs);
	bareleptons.acceptIdPair(pid);
	const bool doClustering = (clusterPhotons != NOCLUSTER);
	const bool useDecayPhotons = (clusterPhotons == CLUSTERALL);
	DressedLeptons leptons(inputfs, bareleptons, dRmax, doClustering, fsCut, useDecayPhotons);
	addProjection(leptons, "DressedLeptons");

	// Add MissingMomentum proj to calc MET
	MissingMomentum vismom(inputfs);
	addProjection(vismom, "MissingET");
	// Set ETmiss
	_etMiss = missingET;

	VetoedFinalState remainingFS;
	remainingFS.addVetoOnThisFinalState(*this);
	addProjection(remainingFS, "RFS");
	}


	/////////////////////////////////////////////////////


	const FinalState& WFinder::remainingFinalState() const {
	return getProjection<FinalState>("RFS");
	}


	int WFinder::compare(const Projection& p) const {
	PCmp dlcmp = mkNamedPCmp(p, "DressedLeptons");
	if (dlcmp != EQUIVALENT) return dlcmp;

	const WFinder& other = dynamic_cast<const WFinder&>(p);
	return (cmp(_minmass, other._minmass) \|\| cmp(_maxmass, other._maxmass) \|\|
	cmp(_useTransverseMass, other._useTransverseMass) \|\|
	cmp(_etMiss, other._etMiss) \|\|
	cmp(_pid, other._pid) \|\| cmp(_trackPhotons, other._trackPhotons));
	}


	void WFinder::project(const Event& e) {
	clear();

	const DressedLeptons& leptons = applyProjection<DressedLeptons>(e, "DressedLeptons");
	const FinalState& neutrinos = applyProjection<FinalState>(e, "Neutrinos");

	// Make and register an invariant mass final state for the W decay leptons
	vector<pair<PdgId, PdgId> > l_nu_ids;
	l_nu_ids += make_pair(abs(_pid), -abs(_nu_pid));
	l_nu_ids += make_pair(-abs(_pid), abs(_nu_pid));
	InvMassFinalState imfs(l_nu_ids, _minmass, _maxmass, _masstarget);
	imfs.useTransverseMass(_useTransverseMass);
	Particles tmp;
	tmp.insert(tmp.end(), leptons.clusteredLeptons().begin(), leptons.clusteredLeptons().end());
	tmp.insert(tmp.end(), neutrinos.particles().begin(), neutrinos.particles().end());
	imfs.calc(tmp);

	if (imfs.particlePairs().size() < 1) return;

	ParticlePair Wconstituents(imfs.particlePairs()[0]);
	Particle p1(Wconstituents.first), p2(Wconstituents.second);

	if (threeCharge(p1) == 0) {
	_constituentLeptons += p2;
	_constituentNeutrinos += p1;
	} else {
	_constituentLeptons += p1;
	_constituentNeutrinos += p2;
	}

	FourMomentum pW = p1.momentum() + p2.momentum();
	const int w3charge = threeCharge(p1) + threeCharge(p2);
	assert(abs(w3charge) == 3);
	const int wcharge = w3charge/3;

	stringstream msg;
	string wsign = (wcharge == 1) ? "+" : "-";
	string wstr = "W" + wsign;
	msg << wstr << " reconstructed from: " << "\n"
	<< " " << p1.momentum() << " " << p1.pid() << "\n"
	<< " + " << p2.momentum() << " " << p2.pid();
	MSG_DEBUG(msg.str());

	// Check missing ET
	const MissingMomentum& vismom = applyProjection<MissingMomentum>(e, "MissingET");
	/// @todo Restrict missing momentum eta range? Use vectorET()?
	if (vismom.scalarEt() < _etMiss) {
	MSG_DEBUG("Not enough missing ET: " << vismom.scalarEt()/GeV
	<< " GeV vs. " << _etMiss/GeV << " GeV");
	return;
	}

	// Make W Particle and insert into particles list
	const PdgId wpid = (wcharge == 1) ? PID::WPLUSBOSON : PID::WMINUSBOSON;
	_bosons.push_back(Particle(wpid, pW));

	// Find the DressedLeptons and neutrinos which survived the IMFS cut such that we can
	// extract their original particles
	foreach (const Particle& p, _constituentNeutrinos) {
	_theParticles.push_back(p);
	}
	foreach (const Particle& p, _constituentLeptons) {
	foreach (const DressedLepton& l, leptons.clusteredLeptons()) {
	if (p.pid() == l.pid() && p.momentum() == l.momentum()) {
	_theParticles.push_back(l.constituentLepton());
	if (_trackPhotons) {
	_theParticles.insert(_theParticles.end(), l.constituentPhotons().begin(), l.constituentPhotons().end());
	}
	}
	}
	}
	}


	}
	diff --git a/src/Projections/ZFinder.cc b/src/Projections/ZFinder.cc
	--- a/src/Projections/ZFinder.cc
	+++ b/src/Projections/ZFinder.cc
	@@ -1,102 +1,102 @@
	// -- C++ --
	#include "Rivet/Projections/ZFinder.hh"
	#include "Rivet/Projections/InvMassFinalState.hh"
	#include "Rivet/Projections/DressedLeptons.hh"
	#include "Rivet/Projections/VetoedFinalState.hh"

	namespace Rivet {


	ZFinder::ZFinder(const FinalState& inputfs,
	- Cut fsCut,
	+ const Cut & fsCut,
	PdgId pid,
	double minmass, double maxmass,
	double dRmax,
	ClusterPhotons clusterPhotons,
	PhotonTracking trackPhotons,
	double masstarget)
	{
	setName("ZFinder");

	_minmass = minmass;
	_maxmass = maxmass;
	_masstarget = masstarget;
	_pid = pid;
	_trackPhotons = trackPhotons;

	IdentifiedFinalState bareleptons(inputfs);
	bareleptons.acceptIdPair(pid);
	const bool doClustering = (clusterPhotons != NOCLUSTER);
	const bool useDecayPhotons = (clusterPhotons == CLUSTERALL);
	DressedLeptons leptons(inputfs, bareleptons, dRmax, doClustering, fsCut, useDecayPhotons);
	addProjection(leptons, "DressedLeptons");

	VetoedFinalState remainingFS;
	remainingFS.addVetoOnThisFinalState(*this);
	addProjection(remainingFS, "RFS");
	}


	/////////////////////////////////////////////////////


	const FinalState& ZFinder::remainingFinalState() const {
	return getProjection<FinalState>("RFS");
	}


	int ZFinder::compare(const Projection& p) const {
	PCmp LCcmp = mkNamedPCmp(p, "DressedLeptons");
	if (LCcmp != EQUIVALENT) return LCcmp;

	const ZFinder& other = dynamic_cast<const ZFinder&>(p);
	return (cmp(_minmass, other._minmass) \|\| cmp(_maxmass, other._maxmass) \|\|
	cmp(_pid, other._pid) \|\| cmp(_trackPhotons, other._trackPhotons));
	}


	void ZFinder::project(const Event& e) {
	clear();

	const DressedLeptons& leptons = applyProjection<DressedLeptons>(e, "DressedLeptons");

	InvMassFinalState imfs(std::make_pair(_pid, -_pid), _minmass, _maxmass, _masstarget);
	Particles tmp;
	tmp.insert(tmp.end(), leptons.clusteredLeptons().begin(), leptons.clusteredLeptons().end());
	imfs.calc(tmp);

	if (imfs.particlePairs().size() < 1) return;
	ParticlePair Zconstituents(imfs.particlePairs()[0]);
	Particle l1(Zconstituents.first), l2(Zconstituents.second);
	if (threeCharge(l1) > 0.0) {
	_constituents += l1, l2;
	} else {
	_constituents += l2, l1;
	}
	FourMomentum pZ = l1.momentum() + l2.momentum();
	assert(threeCharge(l1) + threeCharge(l2) == 0);

	stringstream msg;
	msg << "Z reconstructed from: \n"
	<< " " << l1.momentum() << " " << l1.pid() << "\n"
	<< " + " << l2.momentum() << " " << l2.pid();
	MSG_DEBUG(msg.str());
	_bosons.push_back(Particle(PID::ZBOSON, pZ));

	// Find the DressedLeptons which survived the IMFS cut such that we can
	// extract their original particles
	foreach (const Particle& p, _constituents) {
	foreach (const DressedLepton& l, leptons.clusteredLeptons()) {
	if (p.pid()==l.pid() && p.momentum()==l.momentum()) {
	_theParticles.push_back(l.constituentLepton());
	if (_trackPhotons) {
	_theParticles.insert(_theParticles.end(),
	l.constituentPhotons().begin(), l.constituentPhotons().end());
	}
	}
	}
	}
	}


	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 6:16 AM (1 d, 13 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4981189
Default Alt Text: (58 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions