MeasureDefinition.hh
No OneTemporary
Actions

Size

34 KB

Subscribers

None

MeasureDefinition.hh
View Options

	// Nsubjettiness Package
	// Questions/Comments? jthaler@jthaler.net
	//
	// Copyright (c) 2011-14
	// Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
	//
	// $Id: MeasureDefinition.hh 828 2015-07-20 14:52:06Z jthaler $
	//----------------------------------------------------------------------
	// This file is part of FastJet contrib.
	//
	// It is free software; you can redistribute it and/or modify it under
	// the terms of the GNU General Public License as published by the
	// Free Software Foundation; either version 2 of the License, or (at
	// your option) any later version.
	//
	// It is distributed in the hope that it will be useful, but WITHOUT
	// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
	// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
	// License for more details.
	//
	// You should have received a copy of the GNU General Public License
	// along with this code. If not, see <http://www.gnu.org/licenses/>.
	//----------------------------------------------------------------------

	#ifndef __FASTJET_CONTRIB_MEASUREDEFINITION_HH__
	#define __FASTJET_CONTRIB_MEASUREDEFINITION_HH__

	#include "fastjet/PseudoJet.hh"
	#include <cmath>
	#include <vector>
	#include <list>
	#include <limits>

	#include "TauComponents.hh"

	FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh

	namespace contrib{



	// The following Measures are available (and their relevant arguments):
	// Recommended for usage as jet shapes
	class DefaultMeasure; // Default measure from which next classes derive (should not be called directly)
	class NormalizedMeasure; // (beta,R0)
	class UnnormalizedMeasure; // (beta)
	class NormalizedCutoffMeasure; // (beta,R0,Rcutoff)
	class UnnormalizedCutoffMeasure; // (beta,Rcutoff)

	// New measures as of v2.2
	// Recommended for usage as event shapes (or for jet finding)
	class ConicalMeasure; // (beta,R)
	class OriginalGeometricMeasure; // (R)
	class ModifiedGeometricMeasure; // (R)
	class ConicalGeometricMeasure; // (beta, gamma, R)
	class XConeMeasure; // (beta, R)

	// Formerly GeometricMeasure, now no longer recommended, kept commented out only for cross-check purposes
	//class DeprecatedGeometricMeasure; // (beta)
	//class DeprecatedGeometricCutoffMeasure; // (beta,Rcutoff)


	///////
	//
	// MeasureDefinition
	//
	///////

	//This is a helper class for the Minimum Axes Finders. It is defined later.
	class LightLikeAxis;

	///------------------------------------------------------------------------
	/// \class MeasureDefinition
	/// \brief Base class for measure definitions
	///
	/// This is the base class for measure definitions. Derived classes will calculate
	/// the tau_N of a jet given a specific measure and a set of axes. The measure is
	/// determined by various jet and beam distances (and possible normalization factors).
	///------------------------------------------------------------------------
	class MeasureDefinition {

	public:

	/// Description of measure and parameters
	virtual std::string description() const = 0;

	/// In derived classes, this should return a copy of the corresponding derived class
	virtual MeasureDefinition* create() const = 0;

	//The following five functions define the measure by which tau_N is calculated,
	//and are overloaded by the various measures below

	/// Distanes to jet axis. This is called many times, so needs to be as fast as possible
	/// Unless overloaded, it just calls jet_numerator
	virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
	return jet_numerator(particle,axis);
	}

	/// Distanes to beam. This is called many times, so needs to be as fast as possible
	/// Unless overloaded, it just calls beam_numerator
	virtual double beam_distance_squared(const fastjet::PseudoJet& particle) const {
	return beam_numerator(particle);
	}

	/// The jet measure used in N-(sub)jettiness
	virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const = 0;
	/// The beam measure used in N-(sub)jettiness
	virtual double beam_numerator(const fastjet::PseudoJet& particle) const = 0;

	/// A possible normalization factor
	virtual double denominator(const fastjet::PseudoJet& particle) const = 0;

	/// Run a one-pass minimization routine. There is a generic one-pass minimization that works for a wide variety of measures.
	/// This should be overloaded to create a measure-specific minimization scheme
	virtual std::vector<fastjet::PseudoJet> get_one_pass_axes(int n_jets,
	const std::vector<fastjet::PseudoJet>& inputs,
	const std::vector<fastjet::PseudoJet>& seedAxes,
	int nAttempts = 1000, // cap number of iterations
	double accuracy = 0.0001 // cap distance of closest approach
	) const;

	public:

	/// Returns the tau value for a give set of particles and axes
	double result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const {
	return component_result(particles,axes).tau();
	}

	/// Short-hand for the result() function
	inline double operator() (const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const {
	return result(particles,axes);
	}

	/// Return all of the TauComponents for specific input particles and axes
	TauComponents component_result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const;

	/// Create the partitioning according the jet/beam distances and stores them a TauPartition
	TauPartition get_partition(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const;

	/// Calculate the tau result using an existing partition
	TauComponents component_result_from_partition(const TauPartition& partition, const std::vector<fastjet::PseudoJet>& axes) const;



	/// virtual destructor
	virtual ~MeasureDefinition(){}

	protected:

	/// Flag set by derived classes to choose whether or not to use beam/denominator
	TauMode _tau_mode;

	/// Flag set by derived classes to say whether cheap get_one_pass_axes method can be used (true by default)
	bool _useAxisScaling;

	/// This is the only constructor, which requires _tau_mode and _useAxisScaling to be manually set by derived classes.
	MeasureDefinition() : _tau_mode(UNDEFINED_SHAPE), _useAxisScaling(true) {}


	/// Used by derived classes to set whether or not to use beam/denominator information
	void setTauMode(TauMode tau_mode) {
	_tau_mode = tau_mode;
	}

	/// Used by derived classes to say whether one can use cheap get_one_pass_axes
	void setAxisScaling(bool useAxisScaling) {
	_useAxisScaling = useAxisScaling;
	}

	/// Uses denominator information?
	bool has_denominator() const { return (_tau_mode == NORMALIZED_JET_SHAPE \|\| _tau_mode == NORMALIZED_EVENT_SHAPE);}
	/// Uses beam information?
	bool has_beam() const {return (_tau_mode == UNNORMALIZED_EVENT_SHAPE \|\| _tau_mode == NORMALIZED_EVENT_SHAPE);}

	/// Create light-like axis (used in default one-pass minimization routine)
	fastjet::PseudoJet lightFrom(const fastjet::PseudoJet& input) const {
	double length = sqrt(pow(input.px(),2) + pow(input.py(),2) + pow(input.pz(),2));
	return fastjet::PseudoJet(input.px()/length,input.py()/length,input.pz()/length,1.0);
	}

	/// Shorthand for squaring
	static inline double sq(double x) {return x*x;}

	};


	///////
	//
	// Default Measures
	//
	///////


	///------------------------------------------------------------------------
	/// \enum DefaultMeasureType
	/// \brief Options for default measure
	///
	/// Can be used to switch between pp and ee measure types in the DefaultMeasure
	///------------------------------------------------------------------------
	enum DefaultMeasureType {
	pt_R, /// use transverse momenta and boost-invariant angles,
	E_theta, /// use energies and angles,
	lorentz_dot, /// use dot product inspired measure
	perp_lorentz_dot /// use conical geometric inspired measures
	};

	///------------------------------------------------------------------------
	/// \class DefaultMeasure
	/// \brief Base class for default N-subjettiness measure definitions
	///
	/// This class is the default measure as defined in the original N-subjettiness papers.
	/// Based on the conical measure, but with a normalization factor
	/// This measure is defined as the pT of the particle multiplied by deltaR
	/// to the power of beta. This class includes the normalization factor determined by R0
	///------------------------------------------------------------------------
	class DefaultMeasure : public MeasureDefinition {

	public:

	/// Description
	virtual std::string description() const;
	/// To allow copying around of these objects
	virtual DefaultMeasure* create() const {return new DefaultMeasure(*this);}

	/// fast jet distance
	virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
	return angleSquared(particle, axis);
	}

	/// fast beam distance
	virtual double beam_distance_squared(const fastjet::PseudoJet& /particle/) const {
	return _RcutoffSq;
	}

	/// true jet distance (given by general definitions of energy and angle)
	virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const{
	double jet_dist = angleSquared(particle, axis);
	if (jet_dist > 0.0) {
	return energy(particle) * std::pow(jet_dist,_beta/2.0);
	} else {
	return 0.0;
	}
	}

	/// true beam distance
	virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
	return energy(particle) * std::pow(_Rcutoff,_beta);
	}

	/// possible denominator for normalization
	virtual double denominator(const fastjet::PseudoJet& particle) const {
	return energy(particle) * std::pow(_R0,_beta);
	}

	/// Special minimization routine (from v1.0 of N-subjettiness)
	virtual std::vector<fastjet::PseudoJet> get_one_pass_axes(int n_jets,
	const std::vector<fastjet::PseudoJet>& inputs,
	const std::vector<fastjet::PseudoJet>& seedAxes,
	int nAttempts, // cap number of iterations
	double accuracy // cap distance of closest approach
	) const;

	protected:
	double _beta; ///< Angular exponent
	double _R0; ///< Normalization factor
	double _Rcutoff; ///< Cutoff radius
	double _RcutoffSq; ///< Cutoff radius squared
	DefaultMeasureType _measure_type; ///< Type of measure used (i.e. pp style or ee style)


	/// Constructor is protected so that no one tries to call this directly.
	DefaultMeasure(double beta, double R0, double Rcutoff, DefaultMeasureType measure_type = pt_R)
	: MeasureDefinition(), _beta(beta), _R0(R0), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)), _measure_type(measure_type)
	{
	if (beta <= 0) throw Error("DefaultMeasure: You must choose beta > 0.");
	if (R0 <= 0) throw Error("DefaultMeasure: You must choose R0 > 0.");
	if (Rcutoff <= 0) throw Error("DefaultMeasure: You must choose Rcutoff > 0.");
	}

	/// Added set measure method in case it becomes useful later
	void setDefaultMeasureType(DefaultMeasureType measure_type) {
	_measure_type = measure_type;
	}

	/// Generalized energy value (determined by _measure_type)
	double energy(const PseudoJet& jet) const;
	/// Generalized angle value (determined by _measure_type)
	double angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const;

	/// Name of _measure_type, so description will include the measure type
	std::string measure_type_name() const {
	if (_measure_type == pt_R) return "pt_R";
	else if (_measure_type == E_theta) return "E_theta";
	else if (_measure_type == lorentz_dot) return "lorentz_dot";
	else if (_measure_type == perp_lorentz_dot) return "perp_lorentz_dot";
	else return "Measure Type Undefined";
	}

	/// templated for speed (TODO: probably should remove, since not clear that there is a speed gain)
	template <int N> std::vector<LightLikeAxis> UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
	const std::vector <fastjet::PseudoJet> & inputJets,
	double accuracy) const;

	/// called by get_one_pass_axes to update axes iteratively
	std::vector<LightLikeAxis> UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
	const std::vector <fastjet::PseudoJet> & inputJets,
	double accuracy) const;
	};


	///------------------------------------------------------------------------
	/// \class NormalizedCutoffMeasure
	/// \brief Dimensionless default measure, with radius cutoff
	///
	/// This measure is just a wrapper for DefaultMeasure
	///------------------------------------------------------------------------
	class NormalizedCutoffMeasure : public DefaultMeasure {

	public:

	/// Constructor
	NormalizedCutoffMeasure(double beta, double R0, double Rcutoff, DefaultMeasureType measure_type = pt_R)
	: DefaultMeasure(beta, R0, Rcutoff, measure_type) {
	setTauMode(NORMALIZED_JET_SHAPE);
	}

	/// Description
	virtual std::string description() const;

	/// For copying purposes
	virtual NormalizedCutoffMeasure* create() const {return new NormalizedCutoffMeasure(*this);}

	};

	///------------------------------------------------------------------------
	/// \class NormalizedMeasure
	/// \brief Dimensionless default measure, with no cutoff
	///
	/// This measure is the same as NormalizedCutoffMeasure, with Rcutoff taken to infinity.
	///------------------------------------------------------------------------
	class NormalizedMeasure : public NormalizedCutoffMeasure {

	public:

	/// Constructor
	NormalizedMeasure(double beta, double R0, DefaultMeasureType measure_type = pt_R)
	: NormalizedCutoffMeasure(beta, R0, std::numeric_limits<double>::max(), measure_type) {
	_RcutoffSq = std::numeric_limits<double>::max();
	setTauMode(NORMALIZED_JET_SHAPE);
	}

	/// Description
	virtual std::string description() const;
	/// For copying purposes
	virtual NormalizedMeasure* create() const {return new NormalizedMeasure(*this);}

	};


	///------------------------------------------------------------------------
	/// \class UnnormalizedCutoffMeasure
	/// \brief Dimensionful default measure, with radius cutoff
	///
	/// This class is the unnormalized conical measure. The only difference from NormalizedCutoffMeasure
	/// is that the denominator is defined to be 1.0 by setting _has_denominator to false.
	/// class UnnormalizedCutoffMeasure : public NormalizedCutoffMeasure {
	///------------------------------------------------------------------------
	class UnnormalizedCutoffMeasure : public DefaultMeasure {

	public:

	/// Since all methods are identical, UnnormalizedMeasure inherits directly
	/// from NormalizedMeasure. R0 is a dummy value since the value of R0 is unecessary for this class,
	/// and the "false" flag sets _has_denominator in MeasureDefinition to false so no denominator is used.
	UnnormalizedCutoffMeasure(double beta, double Rcutoff, DefaultMeasureType measure_type = pt_R)
	: DefaultMeasure(beta, std::numeric_limits<double>::quiet_NaN(), Rcutoff, measure_type) {
	setTauMode(UNNORMALIZED_EVENT_SHAPE);
	}

	/// Description
	virtual std::string description() const;
	/// For copying purposes
	virtual UnnormalizedCutoffMeasure* create() const {return new UnnormalizedCutoffMeasure(*this);}

	};


	///------------------------------------------------------------------------
	/// \class UnnormalizedMeasure
	/// \brief Dimensionless default measure, with no cutoff
	///
	/// This measure is the same as UnnormalizedCutoffMeasure, with Rcutoff taken to infinity.
	///------------------------------------------------------------------------
	class UnnormalizedMeasure : public UnnormalizedCutoffMeasure {

	public:
	/// Since all methods are identical, UnnormalizedMeasure inherits directly
	/// from NormalizedMeasure. R0 is a dummy value since the value of R0 is unecessary for this class,
	/// and the "false" flag sets _has_denominator in MeasureDefinition to false so no denominator is used.
	UnnormalizedMeasure(double beta, DefaultMeasureType measure_type = pt_R)
	: UnnormalizedCutoffMeasure(beta, std::numeric_limits<double>::max(), measure_type) {
	_RcutoffSq = std::numeric_limits<double>::max();
	setTauMode(UNNORMALIZED_JET_SHAPE);
	}

	/// Description
	virtual std::string description() const;

	/// For copying purposes
	virtual UnnormalizedMeasure* create() const {return new UnnormalizedMeasure(*this);}

	};


	///------------------------------------------------------------------------
	/// \class ConicalMeasure
	/// \brief Dimensionful event-shape measure, with radius cutoff
	///
	/// Very similar to UnnormalizedCutoffMeasure, but with different normalization convention
	/// and using the new default one-pass minimization algorithm.
	/// Axes are also made to be light-like to ensure sensible behavior
	/// Intended to be used as an event shape.
	///------------------------------------------------------------------------
	class ConicalMeasure : public MeasureDefinition {

	public:

	/// Constructor
	ConicalMeasure(double beta, double Rcutoff)
	: MeasureDefinition(), _beta(beta), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)) {
	if (beta <= 0) throw Error("ConicalMeasure: You must choose beta > 0.");
	if (Rcutoff <= 0) throw Error("ConicalMeasure: You must choose Rcutoff > 0.");
	setTauMode(UNNORMALIZED_EVENT_SHAPE);
	}

	/// Description
	virtual std::string description() const;
	/// For copying purposes
	virtual ConicalMeasure* create() const {return new ConicalMeasure(*this);}

	/// fast jet distance
	virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
	PseudoJet lightAxis = lightFrom(axis);
	return particle.squared_distance(lightAxis);
	}

	/// fast beam distance
	virtual double beam_distance_squared(const fastjet::PseudoJet& /particle/) const {
	return _RcutoffSq;
	}


	/// true jet distance
	virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
	PseudoJet lightAxis = lightFrom(axis);
	double jet_dist = particle.squared_distance(lightAxis)/_RcutoffSq;
	double jet_perp = particle.perp();

	if (_beta == 2.0) {
	return jet_perp * jet_dist;
	} else {
	return jet_perp * pow(jet_dist,_beta/2.0);
	}
	}

	/// true beam distance
	virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
	return particle.perp();
	}

	/// no denominator used for this measure
	virtual double denominator(const fastjet::PseudoJet& /particle/) const {
	return std::numeric_limits<double>::quiet_NaN();
	}

	protected:
	double _beta; ///< angular exponent
	double _Rcutoff; ///< effective jet radius
	double _RcutoffSq;///< effective jet radius squared
	};



	///------------------------------------------------------------------------
	/// \class OriginalGeometricMeasure
	/// \brief Dimensionful event-shape measure, with dot-product distances
	///
	/// This class is the original (and hopefully now correctly coded) geometric measure.
	/// This measure is defined by the Lorentz dot product between
	/// the particle and the axis. This class does not include normalization of tau_N.
	/// New in Nsubjettiness v2.2
	/// NOTE: This is defined differently from the DeprecatedGeometricMeasure which are now commented out.
	///------------------------------------------------------------------------
	class OriginalGeometricMeasure : public MeasureDefinition {

	public:
	/// Constructor
	OriginalGeometricMeasure(double Rcutoff)
	: MeasureDefinition(), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)) {
	if (Rcutoff <= 0) throw Error("OriginalGeometricMeasure: You must choose Rcutoff > 0.");
	setTauMode(UNNORMALIZED_EVENT_SHAPE);
	setAxisScaling(false); // No need to rescale axes (for speed up in one-pass minimization)
	}

	/// Description
	virtual std::string description() const;
	/// For copying purposes
	virtual OriginalGeometricMeasure* create() const {return new OriginalGeometricMeasure(*this);}

	// This class uses the default jet_distance_squared and beam_distance_squared

	/// true jet measure
	virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
	return dot_product(lightFrom(axis), particle)/_RcutoffSq;
	}

	/// true beam measure
	virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
	fastjet::PseudoJet beam_a(0,0,1,1);
	fastjet::PseudoJet beam_b(0,0,-1,1);
	double min_perp = std::min(dot_product(beam_a, particle),dot_product(beam_b, particle));
	return min_perp;
	}

	/// no denominator needed for this measure.
	virtual double denominator(const fastjet::PseudoJet& /particle/) const {
	return std::numeric_limits<double>::quiet_NaN();
	}

	protected:
	double _Rcutoff; ///< Effective jet radius (rho = R^2)
	double _RcutoffSq; ///< Effective jet radius squared

	};


	///------------------------------------------------------------------------
	/// \class ModifiedGeometricMeasure
	/// \brief Dimensionful event-shape measure, with dot-product distances, modified beam measure
	///
	/// This class is the Modified geometric measure. This jet measure is defined by the Lorentz dot product between
	/// the particle and the axis, as in the Original Geometric Measure. The beam measure is defined differently from
	/// the above OriginalGeometric to allow for more conical jets. New in Nsubjettiness v2.2
	///------------------------------------------------------------------------
	class ModifiedGeometricMeasure : public MeasureDefinition {

	public:
	/// Constructor
	ModifiedGeometricMeasure(double Rcutoff)
	: MeasureDefinition(), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)) {
	if (Rcutoff <= 0) throw Error("ModifiedGeometricMeasure: You must choose Rcutoff > 0.");
	setTauMode(UNNORMALIZED_EVENT_SHAPE);
	setAxisScaling(false); // No need to rescale axes (for speed up in one-pass minimization)
	}

	/// Description
	virtual std::string description() const;
	/// For copying purposes
	virtual ModifiedGeometricMeasure* create() const {return new ModifiedGeometricMeasure(*this);}

	// This class uses the default jet_distance_squared and beam_distance_squared

	/// True jet measure
	virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
	return dot_product(lightFrom(axis), particle)/_RcutoffSq;
	}

	/// True beam measure
	virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
	fastjet::PseudoJet lightParticle = lightFrom(particle);
	return 0.5particle.mperp()lightParticle.pt();
	}

	/// This measure does not require a denominator
	virtual double denominator(const fastjet::PseudoJet& /particle/) const {
	return std::numeric_limits<double>::quiet_NaN();
	}

	protected:
	double _Rcutoff; ///< Effective jet radius (rho = R^2)
	double _RcutoffSq; ///< Effective jet radius squared


	};

	///------------------------------------------------------------------------
	/// \class ConicalGeometricMeasure
	/// \brief Dimensionful event-shape measure, basis for XCone jet algorithm
	///
	/// This class is the Conical Geometric measure. This measure is defined by the Lorentz dot product between
	/// the particle and the axis normalized by the axis and particle pT, as well as a factor of cosh(y) to vary
	/// the rapidity depepdence of the beam. New in Nsubjettiness v2.2, and the basis for the XCone jet algorithm
	///------------------------------------------------------------------------
	class ConicalGeometricMeasure : public MeasureDefinition {

	public:

	/// Constructor
	ConicalGeometricMeasure(double jet_beta, double beam_gamma, double Rcutoff)
	: MeasureDefinition(), _jet_beta(jet_beta), _beam_gamma(beam_gamma), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)){
	if (jet_beta <= 0) throw Error("ConicalGeometricMeasure: You must choose beta > 0.");
	if (beam_gamma <= 0) throw Error("ConicalGeometricMeasure: You must choose gamma > 0.");
	if (Rcutoff <= 0) throw Error("ConicalGeometricMeasure: You must choose Rcutoff > 0.");
	setTauMode(UNNORMALIZED_EVENT_SHAPE);
	}

	/// Description
	virtual std::string description() const;
	/// For copying purposes
	virtual ConicalGeometricMeasure* create() const {return new ConicalGeometricMeasure(*this);}

	/// fast jet measure
	virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
	fastjet::PseudoJet lightAxis = lightFrom(axis);
	double pseudoRsquared = 2.0dot_product(lightFrom(axis),particle)/(lightAxis.pt()particle.pt());
	return pseudoRsquared;
	}

	/// fast beam measure
	virtual double beam_distance_squared(const fastjet::PseudoJet& /particle/) const {
	return _RcutoffSq;
	}

	/// true jet measure
	virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
	double jet_dist = jet_distance_squared(particle,axis)/_RcutoffSq;
	if (jet_dist > 0.0) {
	fastjet::PseudoJet lightParticle = lightFrom(particle);
	double weight = (_beam_gamma == 1.0) ? 1.0 : std::pow(0.5*lightParticle.pt(),_beam_gamma - 1.0);
	return particle.pt() * weight * std::pow(jet_dist,_jet_beta/2.0);
	} else {
	return 0.0;
	}
	}

	/// true beam measure
	virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
	fastjet::PseudoJet lightParticle = lightFrom(particle);
	double weight = (_beam_gamma == 1.0) ? 1.0 : std::pow(0.5*lightParticle.pt(),_beam_gamma - 1.0);
	return particle.pt() * weight;
	}

	/// no denominator needed
	virtual double denominator(const fastjet::PseudoJet& /particle/) const {
	return std::numeric_limits<double>::quiet_NaN();
	}

	protected:
	double _jet_beta; ///< jet angular exponent
	double _beam_gamma; ///< beam angular exponent (gamma = 1.0 is recommended)
	double _Rcutoff; ///< effective jet radius
	double _RcutoffSq; ///< effective jet radius squared

	};

	///------------------------------------------------------------------------
	/// \class XConeMeasure
	/// \brief Dimensionful event-shape measure used in XCone jet algorithm
	///
	/// This class is the XCone Measure. This is the default measure for use with the
	/// XCone algorithm. It is identical to the conical geometric measure but with gamma = 1.0.
	///------------------------------------------------------------------------
	class XConeMeasure : public ConicalGeometricMeasure {

	public:
	/// Constructor
	XConeMeasure(double jet_beta, double R)
	: ConicalGeometricMeasure(jet_beta,
	1.0, // beam_gamma, hard coded at gamma = 1.0 default
	R // Rcutoff scale
	) { }

	/// Description
	virtual std::string description() const;
	/// For copying purposes
	virtual XConeMeasure* create() const {return new XConeMeasure(*this);}

	};

	///------------------------------------------------------------------------
	/// \class LightLikeAxis
	/// \brief Helper class to define light-like axes directions
	///
	/// This is a helper class for the minimization routines.
	/// It creates a convenient way of defining axes in order to better facilitate calculations.
	///------------------------------------------------------------------------
	class LightLikeAxis {

	public:
	/// Bare constructor
	LightLikeAxis() : _rap(0.0), _phi(0.0), _weight(0.0), _mom(0.0) {}
	/// Constructor
	LightLikeAxis(double my_rap, double my_phi, double my_weight, double my_mom) :
	_rap(my_rap), _phi(my_phi), _weight(my_weight), _mom(my_mom) {}

	/// Rapidity
	double rap() const {return _rap;}
	/// Azimuth
	double phi() const {return _phi;}
	/// weight factor
	double weight() const {return _weight;}
	/// pt momentum
	double mom() const {return _mom;}

	/// set rapidity
	void set_rap(double my_set_rap) {_rap = my_set_rap;}
	/// set azimuth
	void set_phi(double my_set_phi) {_phi = my_set_phi;}
	/// set weight factor
	void set_weight(double my_set_weight) {_weight = my_set_weight;}
	/// set pt momentum
	void set_mom(double my_set_mom) {_mom = my_set_mom;}
	/// set all kinematics
	void reset(double my_rap, double my_phi, double my_weight, double my_mom) {_rap=my_rap; _phi=my_phi; _weight=my_weight; _mom=my_mom;}

	/// Return PseudoJet version
	fastjet::PseudoJet ConvertToPseudoJet();

	/// Squared distance to PseudoJet
	double DistanceSq(const fastjet::PseudoJet& input) const {
	return DistanceSq(input.rap(),input.phi());
	}

	/// Distance to PseudoJet
	double Distance(const fastjet::PseudoJet& input) const {
	return std::sqrt(DistanceSq(input));
	}

	/// Squared distance to Lightlikeaxis
	double DistanceSq(const LightLikeAxis& input) const {
	return DistanceSq(input.rap(),input.phi());
	}

	/// Distance to Lightlikeaxis
	double Distance(const LightLikeAxis& input) const {
	return std::sqrt(DistanceSq(input));
	}

	private:
	double _rap; ///< rapidity
	double _phi; ///< azimuth
	double _weight; ///< weight factor
	double _mom; ///< pt momentum

	/// Internal squared distance calculation
	double DistanceSq(double rap2, double phi2) const {
	double rap1 = _rap;
	double phi1 = _phi;

	double distRap = rap1-rap2;
	double distPhi = std::fabs(phi1-phi2);
	if (distPhi > M_PI) {distPhi = 2.0*M_PI - distPhi;}
	return distRapdistRap + distPhidistPhi;
	}

	/// Internal distance calculation
	double Distance(double rap2, double phi2) const {
	return std::sqrt(DistanceSq(rap2,phi2));
	}

	};


	////------------------------------------------------------------------------
	///// \class DeprecatedGeometricCutoffMeasure
	//// This class is the old, incorrectly coded geometric measure.
	//// It is kept in case anyone wants to check old code, but should not be used for production purposes.
	//class DeprecatedGeometricCutoffMeasure : public MeasureDefinition {
	//
	//public:
	//
	// // Please, please don't use this.
	// DeprecatedGeometricCutoffMeasure(double jet_beta, double Rcutoff)
	// : MeasureDefinition(),
	// _jet_beta(jet_beta),
	// _beam_beta(1.0), // This is hard coded, since alternative beta_beam values were never checked.
	// _Rcutoff(Rcutoff),
	// _RcutoffSq(sq(Rcutoff)) {
	// setTauMode(UNNORMALIZED_EVENT_SHAPE);
	// setAxisScaling(false);
	// if (jet_beta != 2.0) {
	// throw Error("Geometric minimization is currently only defined for beta = 2.0.");
	// }
	// }
	//
	// virtual std::string description() const;
	//
	// virtual DeprecatedGeometricCutoffMeasure* create() const {return new DeprecatedGeometricCutoffMeasure(*this);}
	//
	// virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
	// fastjet::PseudoJet lightAxis = lightFrom(axis);
	// double pseudoRsquared = 2.0dot_product(lightFrom(axis),particle)/(lightAxis.pt()particle.pt());
	// return pseudoRsquared;
	// }
	//
	// virtual double beam_distance_squared(const fastjet::PseudoJet& /particle/) const {
	// return _RcutoffSq;
	// }
	//
	// virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
	// fastjet::PseudoJet lightAxis = lightFrom(axis);
	// double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(lightAxis.pt(),_beam_beta - 1.0);
	// return particle.pt() * weight * std::pow(jet_distance_squared(particle,axis),_jet_beta/2.0);
	// }
	//
	// virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
	// double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(particle.pt()/particle.e(),_beam_beta - 1.0);
	// return particle.pt() * weight * std::pow(_Rcutoff,_jet_beta);
	// }
	//
	// virtual double denominator(const fastjet::PseudoJet& /particle/) const {
	// return std::numeric_limits<double>::quiet_NaN();
	// }
	//
	// virtual std::vector<fastjet::PseudoJet> get_one_pass_axes(int n_jets,
	// const std::vector<fastjet::PseudoJet>& inputs,
	// const std::vector<fastjet::PseudoJet>& seedAxes,
	// int nAttempts, // cap number of iterations
	// double accuracy // cap distance of closest approach
	// ) const;
	//
	//protected:
	// double _jet_beta;
	// double _beam_beta;
	// double _Rcutoff;
	// double _RcutoffSq;
	//
	//};
	//
	//// ------------------------------------------------------------------------
	//// / \class DeprecatedGeometricMeasure
	//// Same as DeprecatedGeometricMeasureCutoffMeasure, but with Rcutoff taken to infinity.
	//// NOTE: This class should not be used for production purposes.
	//class DeprecatedGeometricMeasure : public DeprecatedGeometricCutoffMeasure {
	//
	//public:
	// DeprecatedGeometricMeasure(double beta)
	// : DeprecatedGeometricCutoffMeasure(beta,std::numeric_limits<double>::max()) {
	// _RcutoffSq = std::numeric_limits<double>::max();
	// setTauMode(UNNORMALIZED_JET_SHAPE);
	// }
	//
	// virtual std::string description() const;
	//
	// virtual DeprecatedGeometricMeasure* create() const {return new DeprecatedGeometricMeasure(*this);}
	//};


	} //namespace contrib

	FASTJET_END_NAMESPACE

	#endif // __FASTJET_CONTRIB_MEASUREDEFINITION_HH__

File Metadata

Mime Type: text/x-c++
Expires: Thu, Apr 3, 8:02 PM (2 h, 57 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4711896
Default Alt Text: MeasureDefinition.hh (34 KB)

MeasureDefinition.hhNo OneTemporaryActions

MeasureDefinition.hhView Options

File Metadata

Event Timeline

MeasureDefinition.hh
No OneTemporary
Actions

MeasureDefinition.hh
View Options