No OneTemporary
Actions

Size

36 KB

Subscribers

None

View Options

	diff --git a/include/Rivet/Makefile.am b/include/Rivet/Makefile.am
	--- a/include/Rivet/Makefile.am
	+++ b/include/Rivet/Makefile.am
	@@ -1,132 +1,133 @@
	## Internal headers - not to be installed
	nobase_dist_noinst_HEADERS =
	## Public headers - to be installed
	nobase_pkginclude_HEADERS =


	## Rivet interface
	nobase_pkginclude_HEADERS += \
	Rivet.hh Exceptions.hh \
	Event.hh Run.hh \
	ParticleBase.hh ParticleName.hh \
	Particle.fhh Particle.hh \
	Jet.fhh Jet.hh \
	ProjectionApplier.hh ProjectionHandler.hh \
	Projection.fhh Projection.hh \
	Cmp.fhh Cmp.hh \
	BeamConstraint.hh AnalysisHandler.hh \
	Analysis.hh AnalysisInfo.hh \
	AnalysisBuilder.hh AnalysisLoader.hh \
	Cuts.hh

	## Build config stuff
	nobase_pkginclude_HEADERS += \
	Config/RivetCommon.hh \
	Config/RivetConfig.hh \
	Config/BuildOptions.hh


	## Projections
	nobase_pkginclude_HEADERS += \
	Projections/AxesDefinition.hh \
	Projections/Beam.hh \
	Projections/BeamThrust.hh \
	Projections/CentralEtHCM.hh \
	Projections/ChargedFinalState.hh \
	Projections/ChargedLeptons.hh \
	Projections/ConstLossyFinalState.hh \
	Projections/DISFinalState.hh \
	Projections/DISKinematics.hh \
	Projections/DISLepton.hh \
	Projections/DressedLeptons.hh \
	Projections/FastJets.hh \
	Projections/FinalPartons.hh \
	Projections/FinalState.hh \
	Projections/FoxWolframMoments.hh \
	Projections/FParameter.hh \
	Projections/HadronicFinalState.hh \
	Projections/HeavyHadrons.hh \
	Projections/Hemispheres.hh \
	Projections/IdentifiedFinalState.hh \
	Projections/InitialQuarks.hh \
	Projections/InvMassFinalState.hh \
	Projections/JetAlg.hh \
	Projections/JetShape.hh \
	Projections/LeadingParticlesFinalState.hh \
	Projections/LossyFinalState.hh \
	Projections/MergedFinalState.hh \
	Projections/MissingMomentum.hh \
	Projections/NeutralFinalState.hh \
	Projections/NonHadronicFinalState.hh \
	Projections/ParticleFinder.hh \
	Projections/ParisiTensor.hh \
	Projections/PrimaryHadrons.hh \
	Projections/PromptFinalState.hh \
	+ Projections/SmearedJets.hh \
	Projections/Sphericity.hh \
	Projections/Spherocity.hh \
	Projections/TauFinder.hh \
	Projections/Thrust.hh \
	Projections/TriggerCDFRun0Run1.hh \
	Projections/TriggerCDFRun2.hh \
	Projections/TriggerUA5.hh \
	Projections/UnstableFinalState.hh \
	Projections/VetoedFinalState.hh \
	Projections/VisibleFinalState.hh \
	Projections/WFinder.hh \
	Projections/ZFinder.hh


	## Analysis base class headers
	nobase_pkginclude_HEADERS += \
	Analyses/MC_ParticleAnalysis.hh \
	Analyses/MC_JetAnalysis.hh \
	Analyses/MC_JetSplittings.hh


	## Tools
	nobase_pkginclude_HEADERS += \
	Tools/Logging.hh \
	Tools/Utils.hh \
	Tools/RivetPaths.hh \
	Tools/RivetSTL.hh \
	Tools/RivetBoost.hh \
	Tools/RivetFastJet.hh \
	Tools/RivetHepMC.hh \
	Tools/RivetYODA.hh \
	Tools/RivetMT2.hh \
	Tools/BinnedHistogram.hh \
	Tools/ParticleUtils.hh \
	Tools/ParticleIdUtils.hh \
	Tools/TypeTraits.hh
	nobase_dist_noinst_HEADERS += \
	Tools/osdir.hh


	## Maths
	nobase_pkginclude_HEADERS += \
	Math/Matrices.hh \
	Math/Vector3.hh \
	Math/VectorN.hh \
	Math/MatrixN.hh \
	Math/MatrixDiag.hh \
	Math/MathHeader.hh \
	Math/Vectors.hh \
	Math/LorentzTrans.hh \
	Math/Matrix3.hh \
	Math/MathUtils.hh \
	Math/Vector4.hh \
	Math/Math.hh \
	Math/Units.hh \
	Math/Constants.hh \
	Math/eigen/util.h \
	Math/eigen/regressioninternal.h \
	Math/eigen/regression.h \
	Math/eigen/vector.h \
	Math/eigen/ludecompositionbase.h \
	Math/eigen/ludecomposition.h \
	Math/eigen/linearsolver.h \
	Math/eigen/linearsolverbase.h \
	Math/eigen/matrix.h \
	Math/eigen/vectorbase.h \
	Math/eigen/projective.h \
	Math/eigen/matrixbase.h
	diff --git a/include/Rivet/Projections/FastJets.hh b/include/Rivet/Projections/FastJets.hh
	--- a/include/Rivet/Projections/FastJets.hh
	+++ b/include/Rivet/Projections/FastJets.hh
	@@ -1,256 +1,249 @@
	// -- C++ --
	#ifndef RIVET_FastJets_HH
	#define RIVET_FastJets_HH

	#include "Rivet/Jet.hh"
	#include "Rivet/Particle.hh"
	#include "Rivet/Projection.hh"
	#include "Rivet/Projections/JetAlg.hh"
	#include "Rivet/Projections/FinalState.hh"
	#include "Rivet/Tools/RivetFastJet.hh"

	#include "fastjet/SISConePlugin.hh"
	#include "fastjet/ATLASConePlugin.hh"
	#include "fastjet/CMSIterativeConePlugin.hh"
	#include "fastjet/CDFJetCluPlugin.hh"
	#include "fastjet/CDFMidPointPlugin.hh"
	#include "fastjet/D0RunIIConePlugin.hh"
	#include "fastjet/TrackJetPlugin.hh"
	#include "fastjet/JadePlugin.hh"
	//#include "fastjet/PxConePlugin.hh"

	namespace Rivet {


	/// Project out jets found using the FastJet package jet algorithms.
	class FastJets : public JetAlg {
	public:

	- /// Wrapper enum for selected Fastjet jet algorithms.
	+ /// Wrapper enum for selected FastJet jet algorithms.
	+ /// @todo Move to JetAlg and alias here?
	enum JetAlgName { KT, CAM, SISCONE, ANTIKT,
	// PXCONE,
	ATLASCONE, CMSCONE,
	CDFJETCLU, CDFMIDPOINT, D0ILCONE,
	JADE, DURHAM, TRACKJET };


	/// @name Constructors etc.
	//@{

	/// "Wrapped" argument constructor using Rivet enums for most common
	/// jet alg choices (including some plugins). For the built-in algs,
	/// E-scheme recombination is used. For full control of
	/// FastJet built-in jet algs, use the native arg constructor.
	FastJets(const FinalState& fsp, JetAlgName alg, double rparameter,
	JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
	JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES,
	double seed_threshold=1.0)
	: JetAlg(fsp, usemuons, useinvis), _adef(0) { _init1(alg, rparameter, seed_threshold); }

	/// Native argument constructor, using FastJet alg/scheme enums.
	FastJets(const FinalState& fsp, fastjet::JetAlgorithm type,
	fastjet::RecombinationScheme recom, double rparameter,
	JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
	JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES)
	: JetAlg(fsp, usemuons, useinvis), _adef(0) { _init2(type, recom, rparameter); }

	/// Explicitly pass in a JetDefinition
	FastJets(const FinalState& fsp, const fastjet::JetDefinition& jdef,
	JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
	JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES)
	: JetAlg(fsp, usemuons, useinvis), _adef(0) { _init3(jdef); }

	/// Explicitly pass in an externally-constructed plugin (must be heap-allocated, Rivet will delete)
	FastJets(const FinalState& fsp, fastjet::JetDefinition::Plugin* plugin,
	JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
	JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES)
	: JetAlg(fsp, usemuons, useinvis), _adef(0) { _init4(plugin); }

	/// Explicitly pass in an externally-constructed plugin (must be heap-allocated, Rivet will delete)
	/// @deprecated Use the pointer version -- it makes the lifetime & ownership more obvious
	FastJets(const FinalState& fsp, fastjet::JetDefinition::Plugin& plugin,
	JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
	JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES)
	: JetAlg(fsp, usemuons, useinvis), _adef(0) { _init4(&plugin); }


	/// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method)
	FastJets(JetAlgName alg, double rparameter, double seed_threshold=1.0)
	: _adef(0) { _init1(alg, rparameter, seed_threshold); }
	/// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method)
	FastJets(fastjet::JetAlgorithm type, fastjet::RecombinationScheme recom, double rparameter)
	: _adef(0) { _init2(type, recom, rparameter); }
	/// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method)
	FastJets(fastjet::JetDefinition::Plugin* plugin)
	: _adef(0) { _init3(plugin); }
	/// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method)
	FastJets(fastjet::JetDefinition::Plugin& plugin)
	: _adef(0) { _init3(&plugin); }


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

	//@}


	public:

	/// Reset the projection. Jet def, etc. are unchanged.
	void reset();

	/// @brief Use provided jet area definition
	///
	/// @warning @a adef is NOT copied, the user must ensure that it remains valid!
	///
	/// Provide an adef null pointer to re-disable jet area calculation
	void useJetArea(fastjet::AreaDefinition* adef) {
	_adef = adef;
	}

	/// @name Access to the jets
	//@{

	- /// Number of jets above the \f$ p_\perp \f$ cut.
	- size_t numJets(double ptmin=0.0) const;
	-
	- /// Number of jets.
	- size_t size() const {
	- return numJets();
	- }
	-
	- /// Get the jets (unordered).
	- Jets _jets(double ptmin=0.0) const;
	+ /// Get the jets (unordered) with pT > ptmin.
	+ Jets _jets() const;

	/// Get the pseudo jets (unordered).
	PseudoJets pseudoJets(double ptmin=0.0) const;
	/// Alias
	PseudoJets pseudojets(double ptmin=0.0) const { return pseudoJets(ptmin); }

	/// Get the pseudo jets, ordered by \f$ p_T \f$.
	PseudoJets pseudoJetsByPt(double ptmin=0.0) const {
	return sorted_by_pt(pseudoJets(ptmin));
	}
	/// Alias
	PseudoJets pseudojetsByPt(double ptmin=0.0) const { return pseudoJetsByPt(ptmin); }

	/// Get the pseudo jets, ordered by \f$ E \f$.
	PseudoJets pseudoJetsByE(double ptmin=0.0) const {
	return sorted_by_E(pseudoJets(ptmin));
	}
	/// Alias
	PseudoJets pseudojetsByE(double ptmin=0.0) const { return pseudoJetsByE(ptmin); }

	/// Get the pseudo jets, ordered by rapidity.
	PseudoJets pseudoJetsByRapidity(double ptmin=0.0) const {
	return sorted_by_rapidity(pseudoJets(ptmin));
	}
	/// Alias
	PseudoJets pseudojetsByRapidity(double ptmin=0.0) const { return pseudoJetsByRapidity(ptmin); }

	//@}

	/// Trim (filter) a jet, keeping tag and constituent info in the resulting jet
	Jet trimJet(const Jet& input, const fastjet::Filter& trimmer) const;

	/// @name Access to the FastJet clustering objects such as jet def, area def, and cluster
	//@{

	/// Return the cluster sequence.
	const fastjet::ClusterSequence* clusterSeq() const {
	return _cseq.get();
	}

	/// Return the area-enabled cluster sequence (if an area defn exists, otherwise returns a null ptr).
	const fastjet::ClusterSequenceArea* clusterSeqArea() const {
	if (areaDef() == NULL) return (fastjet::ClusterSequenceArea*) 0;
	return dynamic_cast<fastjet::ClusterSequenceArea*>(_cseq.get());
	}

	/// Return the jet definition.
	const fastjet::JetDefinition& jetDef() const {
	return _jdef;
	}

	/// Return the area definition.. May be null.
	const fastjet::AreaDefinition* areaDef() const {
	return _adef;
	}

	//@}


	/// @name Access to jet splitting variables: DISABLED FROM 2.3.0, USE FASTJET OBJECTS DIRECTLY INSTEAD
	//@{

	// /// Get the subjet splitting variables for the given jet.
	// vector<double> ySubJet(const fastjet::PseudoJet& jet) const;

	// /// @brief Split a jet a la PRL100,242001(2008).
	// ///
	// /// Based on code from G.Salam, A.Davison.
	// fastjet::PseudoJet splitJet(fastjet::PseudoJet jet, double& last_R) const;

	// /// @brief Filter a jet a la PRL100,242001(2008).
	// ///
	// /// Based on code from G.Salam, A.Davison.
	// fastjet::PseudoJet filterJet(fastjet::PseudoJet jet, double& stingy_R, const double def_R) const;

	//@}


	private:

	/// Shared utility functions to implement constructor behaviour
	/// @todo Replace with calls between constructors when C++11 available?
	void _initBase();
	void _init1(JetAlgName alg, double rparameter, double seed_threshold);
	void _init2(fastjet::JetAlgorithm type, fastjet::RecombinationScheme recom, double rparameter);
	void _init3(const fastjet::JetDefinition& plugin);
	void _init4(fastjet::JetDefinition::Plugin* plugin);

	/// Function to make Rivet::Jet from fastjet::PseudoJet, including constituent and tag info
	Jet _makeJetFromPseudoJet(const PseudoJet& pj) const;
	-
	+
	protected:

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

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

	public:

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


	private:

	/// Jet definition
	fastjet::JetDefinition _jdef;

	/// Pointer to user-handled area definition
	fastjet::AreaDefinition* _adef;

	/// Cluster sequence
	shared_ptr<fastjet::ClusterSequence> _cseq;

	/// FastJet external plugin
	shared_ptr<fastjet::JetDefinition::Plugin> _plugin;

	/// Map of vectors of y scales. This is mutable so we can use caching/lazy evaluation.
	mutable map<int, vector<double> > _yscales;

	/// set of particles sorted by their PT2
	//set<Particle, ParticleBase::byPTAscending> _particles;
	map<int, Particle> _particles;

	};

	}

	#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,246 +1,245 @@
	// -- 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:

	/// Enum for the treatment of muons: whether to include all, some, or none in jet-finding
	enum MuonsStrategy { NO_MUONS, DECAY_MUONS, ALL_MUONS };

	/// Enum for the treatment of invisible particles: whether to include all, some, or none in jet-finding
	enum InvisiblesStrategy { NO_INVISIBLES, DECAY_INVISIBLES, ALL_INVISIBLES };



	/// Constructor
	JetAlg(const FinalState& fs, MuonsStrategy usemuons=JetAlg::ALL_MUONS, InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES);

	/// Default constructor
	JetAlg() {};

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

	/// Destructor
	virtual ~JetAlg() { }


	/// @name Control the treatment of muons and invisible particles
	///
	/// Since MC-based jet calibration (and/or particle flow) can add back in
	/// particles that weren't seen in calorimeters/trackers.
	//@{

	/// @brief Include (some) muons in jet construction.
	///
	/// The default behaviour is that jets are only constructed from visible
	/// particles. Some jet studies, including those from ATLAS, use a definition
	/// in which neutrinos from hadron decays are included via MC-based calibrations.
	/// Setting this flag to true avoids the automatic restriction to a VisibleFinalState.
	void useMuons(MuonsStrategy usemuons=ALL_MUONS) {
	_useMuons = usemuons;
	}

	/// @brief Include (some) invisible particles in jet construction.
	///
	/// The default behaviour is that jets are only constructed from visible
	/// particles. Some jet studies, including those from ATLAS, use a definition
	/// in which neutrinos from hadron decays are included via MC-based calibrations.
	/// Setting this flag to true avoids the automatic restriction to a VisibleFinalState.
	void useInvisibles(InvisiblesStrategy useinvis=DECAY_INVISIBLES) {
	_useInvisibles = useinvis;
	}

	/// @brief Include (some) invisible particles in jet construction.
	///
	/// The default behaviour is that jets are only constructed from visible
	/// particles. Some jet studies, including those from ATLAS, use a definition
	/// in which neutrinos from hadron decays are included via MC-based calibrations.
	/// Setting this flag to true avoids the automatic restriction to a VisibleFinalState.
	///
	/// @deprecated Use the enum-arg version instead
	void useInvisibles(bool useinvis) {
	_useInvisibles = useinvis ? DECAY_INVISIBLES : NO_INVISIBLES;
	}

	//@}


	/// @name Access to jet objects
	//@{

	/// 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(const Cut & c = Cuts::open()) const {
	- const Jets rawjets = _jets(0.0); // arg means no pT cut
	+ virtual Jets jets(const Cut& c=Cuts::open()) const {
	+ const Jets rawjets = _jets();
	// Just return a copy of rawjets if the cut is open
	if (c == Cuts::open()) return rawjets;
	// If there is a non-trivial cut...
	+ /// @todo Use an STL erase(remove_if) and lambda function for this
	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, const 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(const 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(const 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.
	DEPRECATED("Use the version with a sorter function argument.")
	- Jets jetsByP(const 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.
	DEPRECATED("Use the version with a sorter function argument.")
	- Jets jetsByE(const 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.
	DEPRECATED("Use the version with a sorter function argument.")
	- Jets jetsByEt(const 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
	DEPRECATED("Use the version with a Cut argument.")
	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;
	+ virtual Jets _jets() const = 0;


	public:

	- /// Number of jets.
	- virtual size_t size() const = 0;
	- /// Determine if the jet collection is empty.
	+ /// Number of jets passing the provided Cut.
	+ size_t numJets(const Cut& c=Cuts::open()) const { return jets(c).size(); }
	+
	+ /// Number of jets (without cuts).
	+ size_t size() const { return jets().size(); }
	+ /// Whether the inclusive 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 to exclude (some) muons from the would-be constituents.
	MuonsStrategy _useMuons;

	/// Flag to determine whether or not to exclude (some) invisible particles from the would-be constituents.
	InvisiblesStrategy _useInvisibles;


	};


	}

	#endif
	diff --git a/include/Rivet/Projections/SmearedJets.hh b/include/Rivet/Projections/SmearedJets.hh
	new file mode 100644
	--- /dev/null
	+++ b/include/Rivet/Projections/SmearedJets.hh
	@@ -0,0 +1,136 @@
	+// -- C++ --
	+#ifndef RIVET_SmearedJets_HH
	+#define RIVET_SmearedJets_HH
	+
	+#include "Rivet/Jet.hh"
	+#include "Rivet/Particle.hh"
	+#include "Rivet/Projection.hh"
	+#include "Rivet/Projections/JetAlg.hh"
	+#include "boost/function.hpp"
	+
	+namespace Rivet {
	+
	+
	+ void JET_SMEAR_IDENTITY(Jet&) { cout << "JSMEAR!" << endl; }
	+ void PARTICLE_SMEAR_IDENTITY(Particle&) { }
	+ double JET_FN0(const Jet& p) { return 0; }
	+ double JET_FN1(const Jet& p) { cout << "JEFF1" << endl; return 1; }
	+ double PARTICLE_FN0(const Particle& p) { return 0; }
	+ double PARTICLE_FN1(const Particle& p) { return 1; }
	+ double P4_FN0(const FourMomentum& p) { return 0; }
	+ double P4_FN1(const FourMomentum& p) { return 1; }
	+
	+
	+ /// Wrapper projection for smearing {@link Jet}s with detector resolutions and efficiencies
	+ class SmearedJets : public JetAlg {
	+ public:
	+
	+ /// @name Constructors etc.
	+ //@{
	+
	+
	+ // SmearedJets(const JetAlg& ja,
	+ // boost::function<double(const Jet&)>& jetEffFn, boost::function<void(Jet&)>& jetSmearFn=jet_smear_identity)
	+ // boost::function<double(const Jet&)>& bTagEffFn=fn1, boost::function<double(const Jet&)>& bTagMistagFn=fn0,
	+ // boost::function<double(const Jet&)>& cTagEffFn=fn1, boost::function<double(const Jet&)>& cTagMistagFn=fn0)
	+ //: _jetEffFn(jetEffFn), _jetSmearFn(jetSmearFn)
	+ // _bTagEffFn(bTagEffFn), _bTagMistagFn(bTagMistagFn),
	+ // _cTagEffFn(cTagEffFn), _cTagMistagFn(cTagMistagFn)
	+
	+
	+ // typedef double (*J2DFN)(Jet&);
	+ // typedef void (*J2VFN)(Jet&);
	+
	+
	+ /// Constructor with jet efficiency function arg
	+ template <typename J2DFN>
	+ SmearedJets(const JetAlg& ja, J2DFN jetEffFn)
	+ : _jetEffFn(jetEffFn), _jetSmearFn(JET_SMEAR_IDENTITY)
	+ {
	+ setName("SmearedJets");
	+ addProjection(ja, "TruthJets");
	+ }
	+
	+
	+ /// Constructor with jet efficiency and smearing function args
	+ template <typename J2DFN, typename J2VFN>
	+ SmearedJets(const JetAlg& ja, J2DFN jetEffFn, J2VFN jetSmearFn)
	+ : _jetEffFn(jetEffFn), _jetSmearFn(jetSmearFn)
	+ {
	+ setName("SmearedJets");
	+ addProjection(ja, "TruthJets");
	+ }
	+
	+
	+ /// Clone on the heap.
	+ virtual const Projection* clone() const {
	+ return new SmearedJets(*this);
	+ }
	+
	+ //@}
	+
	+
	+ int compare(const Projection& p) const {
	+ /// @todo Maybe do the comparison via void* pointer comparisons, extracted before the function<> wrapping
	+ //const SmearedJets& other = dynamic_cast<const SmearedJets&>(p);
	+ //return cmp(_jetEffFn, other._jetEffFn) \|\| cmp(_jetSmearFn, other._jetSmearFn);
	+ return UNDEFINED;
	+ }
	+
	+
	+ /// @todo Move!!! Improve random numbers...
	+ double rand01() const {
	+ return rand() / (double)RAND_MAX;
	+ }
	+
	+
	+ /// @todo Remove from JetAlg API? I think calc() doesn't work well on projections which chain others
	+ void calc(const Particles& constituents, const Particles& tagparticles=Particles()) {
	+ ///
	+ }
	+
	+
	+ void project(const Event& e) {
	+ _recojets = applyProjection<JetAlg>(e, "TruthJets").jetsByPt();
	+ // Filtering
	+ /// @todo Use erase-remove_if?
	+ if (_jetEffFn) {
	+ foreach (Jet& j, _recojets) {
	+ const double jeff = _jetEffFn(j);
	+ MSG_INFO(j.pT()/GeV << " " << jeff);
	+ if (jeff == 1) continue; //< don't remove; no need to roll expensive dice
	+ if (jeff == 0 \|\| jeff > rand01()) MSG_INFO("REMOVE JET!!!"); ///< @todo Actually remove it...
	+ }
	+ }
	+ // Smearing
	+ if (_jetSmearFn) {
	+ foreach (Jet& j, _recojets) {
	+ _jetSmearFn(j); //< smear the jet momentum
	+ }
	+ }
	+ }
	+
	+
	+ Jets _jets() const { return _recojets; }
	+
	+
	+ /// Reset the projection. Smearing functions will be unchanged.
	+ void reset() { _recojets.clear(); }
	+
	+
	+ /// @todo Modified tagger?
	+
	+
	+ private:
	+
	+ Jets _recojets;
	+ std::function<double(const Jet&)> _jetEffFn;
	+ std::function<void(Jet&)> _jetSmearFn; //, _bTagEffFn, _bTagMistagFn, _cTagEffFn, _cTagMistagFn;
	+
	+ };
	+
	+
	+}
	+
	+
	+#endif
	diff --git a/src/Projections/FastJets.cc b/src/Projections/FastJets.cc
	--- a/src/Projections/FastJets.cc
	+++ b/src/Projections/FastJets.cc
	@@ -1,223 +1,222 @@
	// -- C++ --
	#include "Rivet/Config/RivetCommon.hh"
	#include "Rivet/Tools/Logging.hh"
	#include "Rivet/Projections/FastJets.hh"
	#include "Rivet/Projections/HeavyHadrons.hh"
	#include "Rivet/Projections/TauFinder.hh"

	namespace Rivet {


	void FastJets::_initBase() {
	setName("FastJets");
	addProjection(HeavyHadrons(), "HFHadrons");
	addProjection(TauFinder(TauFinder::HADRONIC), "Taus");
	}

	+
	void FastJets::_init1(JetAlgName alg, double rparameter, double seed_threshold) {
	_initBase();
	MSG_DEBUG("JetAlg = " << alg);
	MSG_DEBUG("R parameter = " << rparameter);
	MSG_DEBUG("Seed threshold = " << seed_threshold);
	if (alg == KT) {
	_jdef = fastjet::JetDefinition(fastjet::kt_algorithm, rparameter, fastjet::E_scheme);
	} else if (alg == CAM) {
	_jdef = fastjet::JetDefinition(fastjet::cambridge_algorithm, rparameter, fastjet::E_scheme);
	} else if (alg == ANTIKT) {
	_jdef = fastjet::JetDefinition(fastjet::antikt_algorithm, rparameter, fastjet::E_scheme);
	} else if (alg == DURHAM) {
	_jdef = fastjet::JetDefinition(fastjet::ee_kt_algorithm, fastjet::E_scheme);
	} else {
	// Plugins:
	if (alg == SISCONE) {
	const double OVERLAP_THRESHOLD = 0.75;
	_plugin.reset(new fastjet::SISConePlugin(rparameter, OVERLAP_THRESHOLD));
	// } else if (alg == PXCONE) {
	// string msg = "PxCone currently not supported, since FastJet doesn't install it by default. ";
	// msg += "Please notify the Rivet authors if this behaviour should be changed.";
	// throw Error(msg);
	//_plugin.reset(new fastjet::PxConePlugin(rparameter));
	} else if (alg == ATLASCONE) {
	const double OVERLAP_THRESHOLD = 0.5;
	_plugin.reset(new fastjet::ATLASConePlugin(rparameter, seed_threshold, OVERLAP_THRESHOLD));
	} else if (alg == CMSCONE) {
	_plugin.reset(new fastjet::CMSIterativeConePlugin(rparameter, seed_threshold));
	} else if (alg == CDFJETCLU) {
	const double OVERLAP_THRESHOLD = 0.75;
	_plugin.reset(new fastjet::CDFJetCluPlugin(rparameter, OVERLAP_THRESHOLD, seed_threshold));
	} else if (alg == CDFMIDPOINT) {
	const double OVERLAP_THRESHOLD = 0.5;
	_plugin.reset(new fastjet::CDFMidPointPlugin(rparameter, OVERLAP_THRESHOLD, seed_threshold));
	} else if (alg == D0ILCONE) {
	const double min_jet_Et = 6.0;
	_plugin.reset(new fastjet::D0RunIIConePlugin(rparameter, min_jet_Et));
	} else if (alg == JADE) {
	_plugin.reset(new fastjet::JadePlugin());
	} else if (alg == TRACKJET) {
	_plugin.reset(new fastjet::TrackJetPlugin(rparameter));
	}
	_jdef = fastjet::JetDefinition(_plugin.get());
	}
	}

	+
	void FastJets::_init2(fastjet::JetAlgorithm type,
	fastjet::RecombinationScheme recom, double rparameter) {
	_initBase();
	_jdef = fastjet::JetDefinition(type, rparameter, recom);
	}

	+
	void FastJets::_init3(const fastjet::JetDefinition& jdef) {
	_initBase();
	_jdef = jdef;
	}

	+
	void FastJets::_init4(fastjet::JetDefinition::Plugin* plugin) {
	_initBase();
	_plugin.reset(plugin);
	_jdef = fastjet::JetDefinition(_plugin.get());
	}



	int FastJets::compare(const Projection& p) const {
	const FastJets& other = dynamic_cast<const FastJets&>(p);
	return \
	cmp(_useMuons, other._useMuons) \|\|
	cmp(_useInvisibles, other._useInvisibles) \|\|
	mkNamedPCmp(other, "FS") \|\|
	cmp(_jdef.jet_algorithm(), other._jdef.jet_algorithm()) \|\|
	cmp(_jdef.recombination_scheme(), other._jdef.recombination_scheme()) \|\|
	cmp(_jdef.plugin(), other._jdef.plugin()) \|\|
	cmp(_jdef.R(), other._jdef.R()) \|\|
	cmp(_adef, other._adef);
	}


	namespace {
	+ /// @todo Replace with C++11 lambdas
	bool isPromptInvisible(const Particle& p) { return !(p.isVisible() \|\| p.fromDecay()); }
	// bool isMuon(const Particle& p) { return p.abspid() == PID::MUON; }
	bool isPromptMuon(const Particle& p) { return isMuon(p) && !p.fromDecay(); }
	}


	void FastJets::project(const Event& e) {
	// Assemble final state particles
	const string fskey = (_useInvisibles == JetAlg::NO_INVISIBLES) ? "VFS" : "FS";
	Particles fsparticles = applyProjection<FinalState>(e, fskey).particles();
	// Remove prompt invisibles if needed (already done by VFS if using NO_INVISIBLES)
	if (_useInvisibles == JetAlg::DECAY_INVISIBLES)
	fsparticles.erase( std::remove_if(fsparticles.begin(), fsparticles.end(), isPromptInvisible), fsparticles.end() );
	// Remove prompt/all muons if needed
	if (_useMuons == JetAlg::DECAY_MUONS)
	fsparticles.erase( std::remove_if(fsparticles.begin(), fsparticles.end(), isPromptMuon), fsparticles.end() );
	else if (_useMuons == JetAlg::NO_MUONS)
	fsparticles.erase( std::remove_if(fsparticles.begin(), fsparticles.end(), isMuon), fsparticles.end() );

	// Tagging particles
	/// @todo Allow the user to specify tag particle kinematic thresholds
	const Particles chadrons = applyProjection<HeavyHadrons>(e, "HFHadrons").cHadrons();
	const Particles bhadrons = applyProjection<HeavyHadrons>(e, "HFHadrons").bHadrons();
	const Particles taus = applyProjection<FinalState>(e, "Taus").particles();
	calc(fsparticles, chadrons+bhadrons+taus);
	}


	void FastJets::calc(const Particles& fsparticles, const Particles& tagparticles) {
	_particles.clear();
	vector<fastjet::PseudoJet> pjs;

	MSG_DEBUG("Finding jets from " << fsparticles.size() << " input particles");

	/// @todo Use FastJet3's UserInfo system

	// Store 4 vector data about each particle into FastJet's PseudoJets
	int counter = 1;
	foreach (const Particle& p, fsparticles) {
	const FourMomentum fv = p.momentum();
	fastjet::PseudoJet pj(fv.px(), fv.py(), fv.pz(), fv.E()); ///< @todo Eliminate with implicit cast?
	pj.set_user_index(counter);
	pjs.push_back(pj);
	_particles[counter] = p;
	counter += 1;
	}
	// And the same for ghost tagging particles (with negative user indices)
	counter = 1;
	foreach (const Particle& p, tagparticles) {
	const FourMomentum fv = 1e-20 * p.momentum(); ///< Ghostify the momentum
	fastjet::PseudoJet pj(fv.px(), fv.py(), fv.pz(), fv.E()); ///< @todo Eliminate with implicit cast?
	pj.set_user_index(-counter);
	pjs.push_back(pj);
	_particles[-counter] = p;
	counter += 1;
	}

	MSG_DEBUG("Running FastJet ClusterSequence construction");
	// Choose cseq as basic or area-calculating
	if (_adef == NULL) {
	_cseq.reset(new fastjet::ClusterSequence(pjs, _jdef));
	} else {
	_cseq.reset(new fastjet::ClusterSequenceArea(pjs, _jdef, *_adef));
	}
	}


	void FastJets::reset() {
	_yscales.clear();
	_particles.clear();
	/// @todo _cseq = fastjet::ClusterSequence();
	}


	- size_t FastJets::numJets(double ptmin) const {
	- if (_cseq.get() == NULL) return 0;
	- return _cseq->inclusive_jets(ptmin).size();
	- }
	-
	-
	- Jets FastJets::_jets(double ptmin) const {
	- Jets rtn;
	+ Jets FastJets::_jets() const {
	+ Jets rtn; rtn.reserve(pseudojets().size());
	foreach (const fastjet::PseudoJet& pj, pseudojets()) {
	rtn.push_back(_makeJetFromPseudoJet(pj));
	}
	/// @todo Cache?
	return rtn;
	}


	Jet FastJets::trimJet(const Jet &input, const fastjet::Filter &trimmer)const{
	assert(input.pseudojet().associated_cluster_sequence() == clusterSeq());
	PseudoJet pj = trimmer(input);
	return _makeJetFromPseudoJet(pj);
	}


	PseudoJets FastJets::pseudoJets(double ptmin) const {
	return (clusterSeq() != NULL) ? clusterSeq()->inclusive_jets(ptmin) : PseudoJets();
	}


	Jet FastJets::_makeJetFromPseudoJet(const PseudoJet &pj)const{
	assert(clusterSeq() != NULL);
	-
	+
	// take the constituents from the cluster sequence, unless the jet was not
	// associated with the cluster sequence (could be the case for trimmed jets)
	const PseudoJets parts = (pj.associated_cluster_sequence() == clusterSeq())?
	clusterSeq()->constituents(pj): pj.constituents();
	-
	+
	vector<Particle> constituents, tags;
	constituents.reserve(parts.size());
	foreach (const fastjet::PseudoJet& p, parts) {
	map<int, Particle>::const_iterator found = _particles.find(p.user_index());
	// assert(found != _particles.end());
	if (found == _particles.end() && p.is_pure_ghost()) continue; //< Pure FJ ghosts are ok
	assert(found != _particles.end()); //< Anything else must be known
	assert(found->first != 0); //< All mapping IDs are pos-def (particles) or neg-def (tags)
	if (found->first > 0) constituents.push_back(found->second);
	else if (found->first < 0) tags.push_back(found->second);
	}
	-
	+
	return Jet(pj, constituents, tags);
	}


	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 4:09 PM (1 d, 13 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3799375
Default Alt Text: (36 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions