No OneTemporary
Actions

Size

50 KB

Referenced Files

None

Subscribers

None

View Options

	diff --git a/include/Rivet/Math/MathUtils.hh b/include/Rivet/Math/MathUtils.hh
	--- a/include/Rivet/Math/MathUtils.hh
	+++ b/include/Rivet/Math/MathUtils.hh
	@@ -1,506 +1,510 @@
	// -- C++ --
	#ifndef RIVET_MathUtils_HH
	#define RIVET_MathUtils_HH

	#include "Rivet/Math/MathHeader.hh"
	#include "Rivet/Tools/RivetBoost.hh"
	#include <cassert>

	namespace Rivet {


	/// @name Comparison functions for safe floating point equality tests
	//@{

	/// Compare a floating point number to zero with a degree
	/// of fuzziness expressed by the absolute @a tolerance parameter.
	inline bool isZero(double val, double tolerance=1E-8) {
	return (fabs(val) < tolerance);
	}

	/// Compare an integral-type number to zero.
	///
	/// Since there is no risk of floating point error, this function just exists
	/// in case @c isZero is accidentally used on an integer type, to avoid
	/// implicit type conversion. The @a tolerance parameter is ignored.
	inline bool isZero(long val, double UNUSED(tolerance)=1E-8) {
	return val == 0;
	}


	/// @brief Compare two floating point numbers for equality with a degree of fuzziness
	///
	/// The @a tolerance parameter is fractional, based on absolute values of the args.
	inline bool fuzzyEquals(double a, double b, double tolerance=1E-5) {
	const double absavg = (fabs(a) + fabs(b))/2.0;
	const double absdiff = fabs(a - b);
	const bool rtn = (isZero(a) && isZero(b)) \|\| absdiff < tolerance*absavg;
	// cout << a << " == " << b << "? " << rtn << endl;
	return rtn;
	}

	/// @brief Compare two integral-type numbers for equality with a degree of fuzziness.
	///
	/// Since there is no risk of floating point error with integral types,
	/// this function just exists in case @c fuzzyEquals is accidentally
	/// used on an integer type, to avoid implicit type conversion. The @a
	/// tolerance parameter is ignored, even if it would have an
	/// absolute magnitude greater than 1.
	inline bool fuzzyEquals(long a, long b, double UNUSED(tolerance)=1E-5) {
	return a == b;
	}


	/// @brief Compare two floating point numbers for >= with a degree of fuzziness
	///
	/// The @a tolerance parameter on the equality test is as for @c fuzzyEquals.
	inline bool fuzzyGtrEquals(double a, double b, double tolerance=1E-5) {
	return a > b \|\| fuzzyEquals(a, b, tolerance);
	}

	/// @brief Compare two integral-type numbers for >= with a degree of fuzziness.
	///
	/// Since there is no risk of floating point error with integral types,
	/// this function just exists in case @c fuzzyGtrEquals is accidentally
	/// used on an integer type, to avoid implicit type conversion. The @a
	/// tolerance parameter is ignored, even if it would have an
	/// absolute magnitude greater than 1.
	inline bool fuzzyGtrEquals(long a, long b, double UNUSED(tolerance)=1E-5) {
	return a >= b;
	}


	/// @brief Compare two floating point numbers for <= with a degree of fuzziness
	///
	/// The @a tolerance parameter on the equality test is as for @c fuzzyEquals.
	inline bool fuzzyLessEquals(double a, double b, double tolerance=1E-5) {
	return a < b \|\| fuzzyEquals(a, b, tolerance);
	}

	/// @brief Compare two integral-type numbers for <= with a degree of fuzziness.
	///
	/// Since there is no risk of floating point error with integral types,
	/// this function just exists in case @c fuzzyLessEquals is accidentally
	/// used on an integer type, to avoid implicit type conversion. The @a
	/// tolerance parameter is ignored, even if it would have an
	/// absolute magnitude greater than 1.
	inline bool fuzzyLessEquals(long a, long b, double UNUSED(tolerance)=1E-5) {
	return a <= b;
	}

	//@}


	/// @name Ranges and intervals
	//@{

	/// Represents whether an interval is open (non-inclusive) or closed (inclusive).
	///
	/// For example, the interval \f$ [0, \pi) \f$ is closed (an inclusive
	/// boundary) at 0, and open (a non-inclusive boundary) at \f$ \pi \f$.
	enum RangeBoundary { OPEN=0, SOFT=0, CLOSED=1, HARD=1 };


	/// @brief Determine if @a value is in the range @a low to @a high, for floating point numbers
	///
	/// Interval boundary types are defined by @a lowbound and @a highbound.
	/// @todo Optimise to one-line at compile time?
	template<typename NUM>
	inline bool inRange(NUM value, NUM low, NUM high,
	RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
	if (lowbound == OPEN && highbound == OPEN) {
	return (value > low && value < high);
	} else if (lowbound == OPEN && highbound == CLOSED) {
	return (value > low && fuzzyLessEquals(value, high));
	} else if (lowbound == CLOSED && highbound == OPEN) {
	return (fuzzyGtrEquals(value, low) && value < high);
	} else { // if (lowbound == CLOSED && highbound == CLOSED) {
	return (fuzzyGtrEquals(value, low) && fuzzyLessEquals(value, high));
	}
	}

	/// Alternative version of inRange for doubles, which accepts a pair for the range arguments.
	template<typename NUM>
	inline bool inRange(NUM value, pair<NUM, NUM> lowhigh,
	RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
	return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
	}


	/// @brief Determine if @a value is in the range @a low to @a high, for integer types
	///
	/// Interval boundary types are defined by @a lowbound and @a highbound.
	/// @todo Optimise to one-line at compile time?
	inline bool inRange(int value, int low, int high,
	RangeBoundary lowbound=CLOSED, RangeBoundary highbound=CLOSED) {
	if (lowbound == OPEN && highbound == OPEN) {
	return (value > low && value < high);
	} else if (lowbound == OPEN && highbound == CLOSED) {
	return (value > low && value <= high);
	} else if (lowbound == CLOSED && highbound == OPEN) {
	return (value >= low && value < high);
	} else { // if (lowbound == CLOSED && highbound == CLOSED) {
	return (value >= low && value <= high);
	}
	}

	/// Alternative version of @c inRange for ints, which accepts a pair for the range arguments.
	inline bool inRange(int value, pair<int, int> lowhigh,
	RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
	return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
	}

	//@}


	/// @name Miscellaneous numerical helpers
	//@{

	/// Named number-type squaring operation.
	template <typename NUM>
	inline NUM sqr(NUM a) {
	return a*a;
	}

	- /// Named number-type addition in quadrature operation.
	+ /// @brief Named number-type addition in quadrature operation.
	+ ///
	+ /// @note Result has the sqrt operation applied.
	template <typename Num>
	inline Num add_quad(Num a, Num b) {
	return sqrt(aa + bb);
	}

	/// Named number-type addition in quadrature operation.
	+ ///
	+ /// @note Result has the sqrt operation applied.
	template <typename Num>
	inline Num add_quad(Num a, Num b, Num c) {
	return sqrt(aa + bb + c*c);
	}

	/// Return a/b, or @a fail if b = 0
	inline double divide(double num, double den, double fail=0.0) {
	return (!isZero(den)) ? num/den : fail;
	}

	/// A more efficient version of pow for raising numbers to integer powers.
	template <typename Num>
	inline Num intpow(Num val, unsigned int exp) {
	assert(exp >= 0);
	if (exp == 0) return (Num) 1;
	else if (exp == 1) return val;
	return val * intpow(val, exp-1);
	}

	/// Find the sign of a number
	inline int sign(double val) {
	if (isZero(val)) return ZERO;
	const int valsign = (val > 0) ? PLUS : MINUS;
	return valsign;
	}

	/// Find the sign of a number
	inline int sign(int val) {
	if (val == 0) return ZERO;
	return (val > 0) ? PLUS : MINUS;
	}

	/// Find the sign of a number
	inline int sign(long val) {
	if (val == 0) return ZERO;
	return (val > 0) ? PLUS : MINUS;
	}

	//@}


	/// @name Binning helper functions
	//@{

	/// @brief Make a list of @a nbins + 1 values equally spaced between @a start and @a end inclusive.
	///
	/// NB. The arg ordering and the meaning of the nbins variable is "histogram-like",
	/// as opposed to the Numpy/Matlab version.
	inline vector<double> linspace(size_t nbins, double start, double end, bool include_end=true) {
	assert(end >= start);
	assert(nbins > 0);
	vector<double> rtn;
	const double interval = (end-start)/static_cast<double>(nbins);
	double edge = start;
	while (inRange(edge, start, end, CLOSED, CLOSED)) {
	rtn.push_back(edge);
	edge += interval;
	}
	assert(rtn.size() == nbins+1);
	if (!include_end) rtn.pop_back();
	return rtn;
	}


	/// @brief Make a list of @a nbins + 1 values exponentially spaced between @a start and @a end inclusive.
	///
	/// NB. The arg ordering and the meaning of the nbins variable is "histogram-like",
	/// as opposed to the Numpy/Matlab version, and the start and end arguments are expressed
	/// in "normal" space, rather than as the logarithms of the start/end values as in Numpy/Matlab.
	inline vector<double> logspace(size_t nbins, double start, double end, bool include_end=true) {
	assert(end >= start);
	assert(start > 0);
	assert(nbins > 0);
	const double logstart = std::log(start);
	const double logend = std::log(end);
	const vector<double> logvals = linspace(nbins, logstart, logend);
	assert(logvals.size() == nbins+1);
	vector<double> rtn;
	foreach (double logval, logvals) {
	rtn.push_back(std::exp(logval));
	}
	assert(rtn.size() == nbins+1);
	if (!include_end) rtn.pop_back();
	return rtn;
	}

	namespace BWHelpers {
	/// @brief CDF for the Breit-Wigner distribution
	inline double CDF(double x, double mu, double gamma) {
	// normalize to (0;1) distribution
	const double xn = (x - mu)/gamma;
	return std::atan(xn)/M_PI + 0.5;
	}

	/// @brief inverse CDF for the Breit-Wigner distribution
	inline double antiCDF(double p, double mu, double gamma) {
	const double xn = std::tan(M_PI*(p-0.5));
	return gamma*xn + mu;
	}
	}

	/// @brief Make a list of @a nbins + 1 values spaced for equal area
	/// Breit-Wigner binning between @a start and @a end inclusive. @a
	/// mu and @a gamma are the Breit-Wigner parameters.
	///
	/// NB. The arg ordering and the meaning of the nbins variable is "histogram-like",
	/// as opposed to the Numpy/Matlab version, and the start and end arguments are expressed
	/// in "normal" space.
	inline vector<double> BWspace(size_t nbins, double start, double end,
	double mu, double gamma) {
	assert(end >= start);
	assert(nbins > 0);
	const double pmin = BWHelpers::CDF(start,mu,gamma);
	const double pmax = BWHelpers::CDF(end, mu,gamma);
	const vector<double> edges = linspace(nbins, pmin, pmax);
	assert(edges.size() == nbins+1);
	vector<double> rtn;
	foreach (double edge, edges) {
	rtn.push_back(BWHelpers::antiCDF(edge,mu,gamma));
	}
	assert(rtn.size() == nbins+1);
	return rtn;
	}


	/// @brief Return the bin index of the given value, @a val, given a vector of bin edges
	///
	/// NB. The @a binedges vector must be sorted
	template <typename NUM>
	inline int index_between(const NUM& val, const vector<NUM>& binedges) {
	if (!inRange(val, binedges.front(), binedges.back())) return -1; //< Out of histo range
	int index = -1;
	for (size_t i = 1; i < binedges.size(); ++i) {
	if (val < binedges[i]) {
	index = i-1;
	break;
	}
	}
	assert(inRange(index, -1, binedges.size()-1));
	return index;
	}

	//@}


	/// @name Statistics functions
	//@{

	/// Calculate the mean of a sample
	inline double mean(const vector<int>& sample) {
	double mean = 0.0;
	for (size_t i=0; i<sample.size(); ++i) {
	mean += sample[i];
	}
	return mean/sample.size();
	}

	// Calculate the error on the mean, assuming poissonian errors
	inline double mean_err(const vector<int>& sample) {
	double mean_e = 0.0;
	for (size_t i=0; i<sample.size(); ++i) {
	mean_e += sqrt(sample[i]);
	}
	return mean_e/sample.size();
	}

	/// Calculate the covariance (variance) between two samples
	inline double covariance(const vector<int>& sample1, const vector<int>& sample2) {
	const double mean1 = mean(sample1);
	const double mean2 = mean(sample2);
	const size_t N = sample1.size();
	double cov = 0.0;
	for (size_t i = 0; i < N; i++) {
	const double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
	cov += cov_i;
	}
	if (N > 1) return cov/(N-1);
	else return 0.0;
	}

	/// Calculate the error on the covariance (variance) of two samples, assuming poissonian errors
	inline double covariance_err(const vector<int>& sample1, const vector<int>& sample2) {
	const double mean1 = mean(sample1);
	const double mean2 = mean(sample2);
	const double mean1_e = mean_err(sample1);
	const double mean2_e = mean_err(sample2);
	const size_t N = sample1.size();
	double cov_e = 0.0;
	for (size_t i = 0; i < N; i++) {
	const double cov_i = (sqrt(sample1[i]) - mean1_e)*(sample2[i] - mean2) +
	(sample1[i] - mean1)*(sqrt(sample2[i]) - mean2_e);
	cov_e += cov_i;
	}
	if (N > 1) return cov_e/(N-1);
	else return 0.0;
	}


	/// Calculate the correlation strength between two samples
	inline double correlation(const vector<int>& sample1, const vector<int>& sample2) {
	const double cov = covariance(sample1, sample2);
	const double var1 = covariance(sample1, sample1);
	const double var2 = covariance(sample2, sample2);
	const double correlation = cov/sqrt(var1*var2);
	const double corr_strength = correlation*sqrt(var2/var1);
	return corr_strength;
	}

	/// Calculate the error of the correlation strength between two samples assuming Poissonian errors
	inline double correlation_err(const vector<int>& sample1, const vector<int>& sample2) {
	const double cov = covariance(sample1, sample2);
	const double var1 = covariance(sample1, sample1);
	const double var2 = covariance(sample2, sample2);
	const double cov_e = covariance_err(sample1, sample2);
	const double var1_e = covariance_err(sample1, sample1);
	const double var2_e = covariance_err(sample2, sample2);

	// Calculate the correlation
	const double correlation = cov/sqrt(var1*var2);
	// Calculate the error on the correlation
	const double correlation_err = cov_e/sqrt(var1*var2) -
	cov/(2pow(3./2., var1var2)) * (var1_e * var2 + var1 * var2_e);


	// Calculate the error on the correlation strength
	const double corr_strength_err = correlation_err*sqrt(var2/var1) +
	correlation/(2sqrt(var2/var1)) (var2_e/var1 - var2*var1_e/pow(2, var2));

	return corr_strength_err;
	}
	//@}


	/// @name Angle range mappings
	//@{

	/// @brief Reduce any number to the range [-2PI, 2PI]
	///
	/// Achieved by repeated addition or subtraction of 2PI as required. Used to
	/// normalise angular measures.
	inline double _mapAngleM2PITo2Pi(double angle) {
	double rtn = fmod(angle, TWOPI);
	if (isZero(rtn)) return 0;
	assert(rtn >= -TWOPI && rtn <= TWOPI);
	return rtn;
	}

	/// Map an angle into the range (-PI, PI].
	inline double mapAngleMPiToPi(double angle) {
	double rtn = _mapAngleM2PITo2Pi(angle);
	if (isZero(rtn)) return 0;
	rtn = (rtn > PI ? rtn-TWOPI :
	rtn <= -PI ? rtn+TWOPI : rtn);
	assert(rtn > -PI && rtn <= PI);
	return rtn;
	}

	/// Map an angle into the range [0, 2PI).
	inline double mapAngle0To2Pi(double angle) {
	double rtn = _mapAngleM2PITo2Pi(angle);
	if (isZero(rtn)) return 0;
	if (rtn < 0) rtn += TWOPI;
	if (rtn == TWOPI) rtn = 0;
	assert(rtn >= 0 && rtn < TWOPI);
	return rtn;
	}

	/// Map an angle into the range [0, PI].
	inline double mapAngle0ToPi(double angle) {
	double rtn = fabs(mapAngleMPiToPi(angle));
	if (isZero(rtn)) return 0;
	assert(rtn > 0 && rtn <= PI);
	return rtn;
	}

	/// Map an angle into the enum-specified range.
	inline double mapAngle(double angle, PhiMapping mapping) {
	switch (mapping) {
	case MINUSPI_PLUSPI:
	return mapAngleMPiToPi(angle);
	case ZERO_2PI:
	return mapAngle0To2Pi(angle);
	case ZERO_PI:
	return mapAngle0To2Pi(angle);
	default:
	throw Rivet::UserError("The specified phi mapping scheme is not implemented");
	}
	}

	//@}


	/// @name Phase space measure helpers
	//@{

	/// @brief Calculate the difference between two angles in radians
	///
	/// Returns in the range [0, PI].
	inline double deltaPhi(double phi1, double phi2) {
	return mapAngle0ToPi(phi1 - phi2);
	}

	/// Calculate the difference between two pseudorapidities,
	/// returning the unsigned value.
	inline double deltaEta(double eta1, double eta2) {
	return fabs(eta1 - eta2);
	}

	/// Calculate the distance between two points in 2D rapidity-azimuthal
	/// ("\f$ \eta-\phi \f$") space. The phi values are given in radians.
	inline double deltaR(double rap1, double phi1, double rap2, double phi2) {
	const double dphi = deltaPhi(phi1, phi2);
	return sqrt( sqr(rap1-rap2) + sqr(dphi) );
	}

	/// Calculate a rapidity value from the supplied energy @a E and longitudinal momentum @a pz.
	inline double rapidity(double E, double pz) {
	if (isZero(E - pz)) {
	throw std::runtime_error("Divergent positive rapidity");
	return MAXDOUBLE;
	}
	if (isZero(E + pz)) {
	throw std::runtime_error("Divergent negative rapidity");
	return -MAXDOUBLE;
	}
	return 0.5*log((E+pz)/(E-pz));
	}

	//@}


	}


	#endif
	diff --git a/include/Rivet/Math/Vector4.hh b/include/Rivet/Math/Vector4.hh
	--- a/include/Rivet/Math/Vector4.hh
	+++ b/include/Rivet/Math/Vector4.hh
	@@ -1,936 +1,936 @@
	#ifndef RIVET_MATH_VECTOR4
	#define RIVET_MATH_VECTOR4

	#include "Rivet/Math/MathHeader.hh"
	#include "Rivet/Math/MathUtils.hh"
	#include "Rivet/Math/VectorN.hh"
	#include "Rivet/Math/Vector3.hh"

	namespace Rivet {


	class FourVector;
	class FourMomentum;
	class LorentzTransform;
	typedef FourVector Vector4;
	FourVector transform(const LorentzTransform& lt, const FourVector& v4);


	/// @brief Specialisation of VectorN to a general (non-momentum) Lorentz 4-vector.
	class FourVector : public Vector<4> {
	friend FourVector multiply(const double a, const FourVector& v);
	friend FourVector multiply(const FourVector& v, const double a);
	friend FourVector add(const FourVector& a, const FourVector& b);
	friend FourVector transform(const LorentzTransform& lt, const FourVector& v4);

	public:

	FourVector() : Vector<4>() { }

	template<typename V4>
	FourVector(const V4& other) {
	this->setT(other.t());
	this->setX(other.x());
	this->setY(other.y());
	this->setZ(other.z());
	}

	FourVector(const Vector<4>& other)
	: Vector<4>(other) { }

	FourVector(const double t, const double x, const double y, const double z) {
	this->setT(t);
	this->setX(x);
	this->setY(y);
	this->setZ(z);
	}

	virtual ~FourVector() { }

	public:

	double t() const { return get(0); }
	double x() const { return get(1); }
	double y() const { return get(2); }
	double z() const { return get(3); }
	FourVector& setT(const double t) { set(0, t); return *this; }
	FourVector& setX(const double x) { set(1, x); return *this; }
	FourVector& setY(const double y) { set(2, y); return *this; }
	FourVector& setZ(const double z) { set(3, z); return *this; }

	double invariant() const {
	// Done this way for numerical precision
	return (t() + z())(t() - z()) - x()x() - y()*y();
	}

	/// Angle between this vector and another
	double angle(const FourVector& v) const {
	return vector3().angle( v.vector3() );
	}

	/// Angle between this vector and another (3-vector)
	double angle(const Vector3& v3) const {
	return vector3().angle(v3);
	}

	/// @brief Square of the projection of the 3-vector on to the \f$ x-y \f$ plane
	/// This is a more efficient function than @c polarRadius, as it avoids the square root.
	/// Use it if you only need the squared value, or e.g. an ordering by magnitude.
	double polarRadius2() const {
	return vector3().polarRadius2();
	}

	/// Synonym for polarRadius2
	double perp2() const {
	return vector3().perp2();
	}

	/// Synonym for polarRadius2
	double rho2() const {
	return vector3().rho2();
	}

	/// Projection of 3-vector on to the \f$ x-y \f$ plane
	double polarRadius() const {
	return vector3().polarRadius();
	}

	/// Synonym for polarRadius
	double perp() const {
	return vector3().perp();
	}

	/// Synonym for polarRadius
	double rho() const {
	return vector3().rho();
	}

	/// Angle subtended by the 3-vector's projection in x-y and the x-axis.
	double azimuthalAngle(const PhiMapping mapping = ZERO_2PI) const {
	return vector3().azimuthalAngle(mapping);
	}

	/// Synonym for azimuthalAngle.
	double phi(const PhiMapping mapping = ZERO_2PI) const {
	return vector3().phi(mapping);
	}

	/// Angle subtended by the 3-vector and the z-axis.
	double polarAngle() const {
	return vector3().polarAngle();
	}

	/// Synonym for polarAngle.
	double theta() const {
	return vector3().theta();
	}

	/// Pseudorapidity (defined purely by the 3-vector components)
	double pseudorapidity() const {
	return vector3().pseudorapidity();
	}

	/// Synonym for pseudorapidity.
	double eta() const {
	return vector3().eta();
	}

	/// Get the spatial part of the 4-vector as a 3-vector.
	Vector3 vector3() const {
	return Vector3(get(1), get(2), get(3));
	}


	public:

	/// Contract two 4-vectors, with metric signature (+ - - -).
	double contract(const FourVector& v) const {
	const double result = t()v.t() - x()v.x() - y()v.y() - z()v.z();
	return result;
	}

	/// Contract two 4-vectors, with metric signature (+ - - -).
	double dot(const FourVector& v) const {
	return contract(v);
	}

	/// Contract two 4-vectors, with metric signature (+ - - -).
	double operator*(const FourVector& v) const {
	return contract(v);
	}

	/// Multiply by a scalar.
	FourVector& operator*=(double a) {
	_vec = multiply(a, *this)._vec;
	return *this;
	}

	/// Divide by a scalar.
	FourVector& operator/=(double a) {
	_vec = multiply(1.0/a, *this)._vec;
	return *this;
	}

	/// Add to this 4-vector.
	FourVector& operator+=(const FourVector& v) {
	_vec = add(*this, v)._vec;
	return *this;
	}

	/// Subtract from this 4-vector. NB time as well as space components are subtracted.
	FourVector& operator-=(const FourVector& v) {
	_vec = add(*this, -v)._vec;
	return *this;
	}

	/// Multiply all components (space and time) by -1.
	FourVector operator-() const {
	FourVector result;
	result._vec = -_vec;
	return result;
	}

	};


	/// Contract two 4-vectors, with metric signature (+ - - -).
	inline double contract(const FourVector& a, const FourVector& b) {
	return a.contract(b);
	}

	/// Contract two 4-vectors, with metric signature (+ - - -).
	inline double dot(const FourVector& a, const FourVector& b) {
	return contract(a, b);
	}

	inline FourVector multiply(const double a, const FourVector& v) {
	FourVector result;
	result._vec = a * v._vec;
	return result;
	}

	inline FourVector multiply(const FourVector& v, const double a) {
	return multiply(a, v);
	}

	inline FourVector operator*(const double a, const FourVector& v) {
	return multiply(a, v);
	}

	inline FourVector operator*(const FourVector& v, const double a) {
	return multiply(a, v);
	}

	inline FourVector operator/(const FourVector& v, const double a) {
	return multiply(1.0/a, v);
	}

	inline FourVector add(const FourVector& a, const FourVector& b) {
	FourVector result;
	result._vec = a._vec + b._vec;
	return result;
	}

	inline FourVector operator+(const FourVector& a, const FourVector& b) {
	return add(a, b);
	}

	inline FourVector operator-(const FourVector& a, const FourVector& b) {
	return add(a, -b);
	}

	/// Calculate the Lorentz self-invariant of a 4-vector.
	/// \f$ v_\mu v^\mu = g_{\mu\nu} x^\mu x^\nu \f$.
	inline double invariant(const FourVector& lv) {
	return lv.invariant();
	}

	/// Angle (in radians) between spatial parts of two Lorentz vectors.
	inline double angle(const FourVector& a, const FourVector& b) {
	return a.angle(b);
	}

	/// Angle (in radians) between spatial parts of two Lorentz vectors.
	inline double angle(const Vector3& a, const FourVector& b) {
	return angle( a, b.vector3() );
	}

	/// Angle (in radians) between spatial parts of two Lorentz vectors.
	inline double angle(const FourVector& a, const Vector3& b) {
	return a.angle(b);
	}

	/// Calculate transverse length sq. \f$ \rho^2 \f$ of a Lorentz vector.
	inline double polarRadius2(const FourVector& v) {
	return v.polarRadius2();
	}
	/// Synonym for polarRadius2.
	inline double perp2(const FourVector& v) {
	return v.perp2();
	}
	/// Synonym for polarRadius2.
	inline double rho2(const FourVector& v) {
	return v.rho2();
	}

	/// Calculate transverse length \f$ \rho \f$ of a Lorentz vector.
	inline double polarRadius(const FourVector& v) {
	return v.polarRadius();
	}
	/// Synonym for polarRadius.
	inline double perp(const FourVector& v) {
	return v.perp();
	}
	/// Synonym for polarRadius.
	inline double rho(const FourVector& v) {
	return v.rho();
	}

	/// Calculate azimuthal angle of a Lorentz vector.
	inline double azimuthalAngle(const FourVector& v, const PhiMapping mapping = ZERO_2PI) {
	return v.azimuthalAngle(mapping);
	}
	/// Synonym for azimuthalAngle.
	inline double phi(const FourVector& v, const PhiMapping mapping = ZERO_2PI) {
	return v.phi(mapping);
	}


	/// Calculate polar angle of a Lorentz vector.
	inline double polarAngle(const FourVector& v) {
	return v.polarAngle();
	}
	/// Synonym for polarAngle.
	inline double theta(const FourVector& v) {
	return v.theta();
	}

	/// Calculate pseudorapidity of a Lorentz vector.
	inline double pseudorapidity(const FourVector& v) {
	return v.pseudorapidity();
	}
	/// Synonym for pseudorapidity.
	inline double eta(const FourVector& v) {
	return v.eta();
	}



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



	/// Specialized version of the FourVector with momentum/energy functionality.
	class FourMomentum : public FourVector {
	friend FourMomentum multiply(const double a, const FourMomentum& v);
	friend FourMomentum multiply(const FourMomentum& v, const double a);
	friend FourMomentum add(const FourMomentum& a, const FourMomentum& b);
	friend FourMomentum transform(const LorentzTransform& lt, const FourMomentum& v4);

	public:
	FourMomentum() { }

	template<typename V4>
	FourMomentum(const V4& other) {
	this->setE(other.t());
	this->setPx(other.x());
	this->setPy(other.y());
	this->setPz(other.z());
	}

	FourMomentum(const Vector<4>& other)
	: FourVector(other) { }

	FourMomentum(const double E, const double px, const double py, const double pz) {
	this->setE(E);
	this->setPx(px);
	this->setPy(py);
	this->setPz(pz);
	}

	~FourMomentum() {}

	public:
	/// Get energy \f$ E \f$ (time component of momentum).
	double E() const { return t(); }

	/// Get 3-momentum part, \f$ p \f$.
	Vector3 p() const { return vector3(); }

	/// Get x-component of momentum \f$ p_x \f$.
	double px() const { return x(); }

	/// Get y-component of momentum \f$ p_y \f$.
	double py() const { return y(); }

	/// Get z-component of momentum \f$ p_z \f$.
	double pz() const { return z(); }

	/// Set energy \f$ E \f$ (time component of momentum).
	FourMomentum& setE(double E) { setT(E); return *this; }

	/// Set x-component of momentum \f$ p_x \f$.
	FourMomentum& setPx(double px) { setX(px); return *this; }

	/// Set y-component of momentum \f$ p_y \f$.
	FourMomentum& setPy(double py) { setY(py); return *this; }

	/// Set z-component of momentum \f$ p_z \f$.
	FourMomentum& setPz(double pz) { setZ(pz); return *this; }

	/// @brief Get the mass \f$ m = \sqrt{E^2 - p^2} \f$ (the Lorentz self-invariant).
	///
	/// For spacelike momenta, the mass will be -sqrt(\|mass2\|).
	double mass() const {
	// assert(Rivet::isZero(mass2()) \|\| mass2() > 0);
	// if (Rivet::isZero(mass2())) {
	// return 0.0;
	// } else {
	// return sqrt(mass2());
	// }
	return sign(mass2()) * sqrt(fabs(mass2()));
	}

	/// Get the squared mass \f$ m^2 = E^2 - p^2 \f$ (the Lorentz self-invariant).
	double mass2() const {
	return invariant();
	}

	/// Calculate the rapidity.
	double rapidity() const {
	return 0.5 * std::log( (E() + pz()) / (E() - pz()) );
	}

	/// Calculate the squared transverse momentum \f$ p_T^2 \f$.
	double pT2() const {
	return vector3().polarRadius2();
	}

	/// Calculate the transverse momentum \f$ p_T \f$.
	double pT() const {
	return sqrt(pT2());
	}

	/// Calculate the transverse energy \f$ E_T^2 = E^2 \sin^2{\theta} \f$.
	double Et2() const {
	return Et() * Et();
	}

	/// Calculate the transverse energy \f$ E_T = E \sin{\theta} \f$.
	double Et() const {
	return E() * sin(polarAngle());
	}

	/// Calculate the boost vector (in units of \f$ \beta \f$).
	Vector3 boostVector() const {
	// const Vector3 p3 = vector3();
	// const double m2 = mass2();
	// if (Rivet::isZero(m2)) return p3.unit();
	// else {
	// // Could also do this via beta = tanh(rapidity), but that's
	// // probably more messy from a numerical hygiene point of view.
	// const double p2 = p3.mod2();
	// const double beta = sqrt( p2 / (m2 + p2) );
	// return beta * p3.unit();
	// }
	/// @todo Be careful about c=1 convention...
	return Vector3(px()/E(), py()/E(), pz()/E());
	}

	/// Struct for sorting by increasing energy
	struct byEAscending {
	bool operator()(const FourMomentum& left, const FourMomentum& right) const{
	double pt2left = left.E();
	double pt2right = right.E();
	return pt2left < pt2right;
	}

	bool operator()(const FourMomentum* left, const FourMomentum* right) const{
	return (*this)(left, right);
	}
	};

	/// Struct for sorting by decreasing energy
	struct byEDescending {
	bool operator()(const FourMomentum& left, const FourMomentum& right) const{
	return byEAscending()(right, left);
	}

	bool operator()(const FourMomentum* left, const FourVector* right) const{
	return (*this)(left, right);
	}
	};


	/// Multiply by a scalar
	FourMomentum& operator*=(double a) {
	_vec = multiply(a, *this)._vec;
	return *this;
	}

	/// Divide by a scalar
	FourMomentum& operator/=(double a) {
	_vec = multiply(1.0/a, *this)._vec;
	return *this;
	}

	- /// Subtract from this 4-vector. NB time as well as space components are subtracted.
	+ /// Add to this 4-vector. NB time as well as space components are added.
	FourMomentum& operator+=(const FourMomentum& v) {
	_vec = add(*this, v)._vec;
	return *this;
	}

	/// Subtract from this 4-vector. NB time as well as space components are subtracted.
	FourMomentum& operator-=(const FourMomentum& v) {
	_vec = add(*this, -v)._vec;
	return *this;
	}

	/// Multiply all components (time and space) by -1.
	FourMomentum operator-() const {
	FourMomentum result;
	result._vec = -_vec;
	return result;
	}


	};


	inline FourMomentum multiply(const double a, const FourMomentum& v) {
	FourMomentum result;
	result._vec = a * v._vec;
	return result;
	}

	inline FourMomentum multiply(const FourMomentum& v, const double a) {
	return multiply(a, v);
	}

	inline FourMomentum operator*(const double a, const FourMomentum& v) {
	return multiply(a, v);
	}

	inline FourMomentum operator*(const FourMomentum& v, const double a) {
	return multiply(a, v);
	}

	inline FourMomentum operator/(const FourMomentum& v, const double a) {
	return multiply(1.0/a, v);
	}

	inline FourMomentum add(const FourMomentum& a, const FourMomentum& b) {
	FourMomentum result;
	result._vec = a._vec + b._vec;
	return result;
	}

	inline FourMomentum operator+(const FourMomentum& a, const FourMomentum& b) {
	return add(a, b);
	}

	inline FourMomentum operator-(const FourMomentum& a, const FourMomentum& b) {
	return add(a, -b);
	}



	/// Get the mass \f$ m = \sqrt{E^2 - p^2} \f$ (the Lorentz self-invariant) of a momentum 4-vector.
	inline double mass(const FourMomentum& v) {
	return v.mass();
	}

	/// Get the squared mass \f$ m^2 = E^2 - p^2 \f$ (the Lorentz self-invariant) of a momentum 4-vector.
	inline double mass2(const FourMomentum& v) {
	return v.mass2();
	}

	/// Calculate the rapidity of a momentum 4-vector.
	inline double rapidity(const FourMomentum& v) {
	return v.rapidity();
	}

	/// Calculate the squared transverse momentum \f$ p_T^2 \f$ of a momentum 4-vector.
	inline double pT2(const FourMomentum& v) {
	return v.pT2();
	}

	/// Calculate the transverse momentum \f$ p_T \f$ of a momentum 4-vector.
	inline double pT(const FourMomentum& v) {
	return v.pT();
	}

	/// Calculate the transverse energy squared, \f$ E_T^2 = E^2 \sin^2{\theta} \f$ of a momentum 4-vector.
	inline double Et2(const FourMomentum& v) {
	return v.Et2();}

	/// Calculate the transverse energy \f$ E_T = E \sin{\theta} \f$ of a momentum 4-vector.
	inline double Et(const FourMomentum& v) {
	return v.Et();
	}

	/// Calculate the velocity boost vector of a momentum 4-vector.
	inline Vector3 boostVector(const FourMomentum& v) {
	return v.boostVector();
	}


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


	/// @name \f$ \Delta R \f$ calculations from 4-vectors
	//@{

	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors.
	/// There is a scheme ambiguity for momentum-type four vectors as to whether
	/// the pseudorapidity (a purely geometric concept) or the rapidity (a
	/// relativistic energy-momentum quantity) is to be used: this can be chosen
	/// via the optional scheme parameter. Use of this scheme option is
	/// discouraged in this case since @c RAPIDITY is only a valid option for
	/// vectors whose type is really the FourMomentum derived class.
	inline double deltaR(const FourVector& a, const FourVector& b,
	RapScheme scheme = PSEUDORAPIDITY) {
	switch (scheme) {
	case PSEUDORAPIDITY :
	return deltaR(a.vector3(), b.vector3());
	case RAPIDITY:
	{
	const FourMomentum* ma = dynamic_cast<const FourMomentum*>(&a);
	const FourMomentum* mb = dynamic_cast<const FourMomentum*>(&b);
	if (!ma \|\| !mb) {
	string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
	throw std::runtime_error(err);
	}
	return deltaR(ma, mb, scheme);
	}
	default:
	throw std::runtime_error("The specified deltaR scheme is not yet implemented");
	}
	}


	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors.
	/// There is a scheme ambiguity for momentum-type four vectors
	/// as to whether the pseudorapidity (a purely geometric concept) or the
	/// rapidity (a relativistic energy-momentum quantity) is to be used: this can
	/// be chosen via the optional scheme parameter.
	inline double deltaR(const FourVector& v,
	double eta2, double phi2,
	RapScheme scheme = PSEUDORAPIDITY) {
	switch (scheme) {
	case PSEUDORAPIDITY :
	return deltaR(v.vector3(), eta2, phi2);
	case RAPIDITY:
	{
	const FourMomentum* mv = dynamic_cast<const FourMomentum*>(&v);
	if (!mv) {
	string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
	throw std::runtime_error(err);
	}
	return deltaR(*mv, eta2, phi2, scheme);
	}
	default:
	throw std::runtime_error("The specified deltaR scheme is not yet implemented");
	}
	}


	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors.
	/// There is a scheme ambiguity for momentum-type four vectors
	/// as to whether the pseudorapidity (a purely geometric concept) or the
	/// rapidity (a relativistic energy-momentum quantity) is to be used: this can
	/// be chosen via the optional scheme parameter.
	inline double deltaR(double eta1, double phi1,
	const FourVector& v,
	RapScheme scheme = PSEUDORAPIDITY) {
	switch (scheme) {
	case PSEUDORAPIDITY :
	return deltaR(eta1, phi1, v.vector3());
	case RAPIDITY:
	{
	const FourMomentum* mv = dynamic_cast<const FourMomentum*>(&v);
	if (!mv) {
	string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
	throw std::runtime_error(err);
	}
	return deltaR(eta1, phi1, *mv, scheme);
	}
	default:
	throw std::runtime_error("The specified deltaR scheme is not yet implemented");
	}
	}


	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors.
	/// There is a scheme ambiguity for momentum-type four vectors
	/// as to whether the pseudorapidity (a purely geometric concept) or the
	/// rapidity (a relativistic energy-momentum quantity) is to be used: this can
	/// be chosen via the optional scheme parameter.
	inline double deltaR(const FourMomentum& a, const FourMomentum& b,
	RapScheme scheme = PSEUDORAPIDITY) {
	switch (scheme) {
	case PSEUDORAPIDITY:
	return deltaR(a.vector3(), b.vector3());
	case RAPIDITY:
	return deltaR(a.rapidity(), a.azimuthalAngle(), b.rapidity(), b.azimuthalAngle());
	default:
	throw std::runtime_error("The specified deltaR scheme is not yet implemented");
	}
	}

	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors.
	/// There is a scheme ambiguity for momentum-type four vectors
	/// as to whether the pseudorapidity (a purely geometric concept) or the
	/// rapidity (a relativistic energy-momentum quantity) is to be used: this can
	/// be chosen via the optional scheme parameter.
	inline double deltaR(const FourMomentum& v,
	double eta2, double phi2,
	RapScheme scheme = PSEUDORAPIDITY) {
	switch (scheme) {
	case PSEUDORAPIDITY:
	return deltaR(v.vector3(), eta2, phi2);
	case RAPIDITY:
	return deltaR(v.rapidity(), v.azimuthalAngle(), eta2, phi2);
	default:
	throw std::runtime_error("The specified deltaR scheme is not yet implemented");
	}
	}


	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors.
	/// There is a scheme ambiguity for momentum-type four vectors
	/// as to whether the pseudorapidity (a purely geometric concept) or the
	/// rapidity (a relativistic energy-momentum quantity) is to be used: this can
	/// be chosen via the optional scheme parameter.
	inline double deltaR(double eta1, double phi1,
	const FourMomentum& v,
	RapScheme scheme = PSEUDORAPIDITY) {
	switch (scheme) {
	case PSEUDORAPIDITY:
	return deltaR(eta1, phi1, v.vector3());
	case RAPIDITY:
	return deltaR(eta1, phi1, v.rapidity(), v.azimuthalAngle());
	default:
	throw std::runtime_error("The specified deltaR scheme is not yet implemented");
	}
	}

	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors.
	/// There is a scheme ambiguity for momentum-type four vectors
	/// as to whether the pseudorapidity (a purely geometric concept) or the
	/// rapidity (a relativistic energy-momentum quantity) is to be used: this can
	/// be chosen via the optional scheme parameter.
	inline double deltaR(const FourMomentum& a, const FourVector& b,
	RapScheme scheme = PSEUDORAPIDITY) {
	switch (scheme) {
	case PSEUDORAPIDITY:
	return deltaR(a.vector3(), b.vector3());
	case RAPIDITY:
	return deltaR(a.rapidity(), a.azimuthalAngle(), FourMomentum(b).rapidity(), b.azimuthalAngle());
	default:
	throw std::runtime_error("The specified deltaR scheme is not yet implemented");
	}
	}

	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two four-vectors.
	/// There is a scheme ambiguity for momentum-type four vectors
	/// as to whether the pseudorapidity (a purely geometric concept) or the
	/// rapidity (a relativistic energy-momentum quantity) is to be used: this can
	/// be chosen via the optional scheme parameter.
	inline double deltaR(const FourVector& a, const FourMomentum& b,
	RapScheme scheme = PSEUDORAPIDITY) {
	switch (scheme) {
	case PSEUDORAPIDITY:
	return deltaR(a.vector3(), b.vector3());
	case RAPIDITY:
	return deltaR(FourMomentum(a).rapidity(), a.azimuthalAngle(), b.rapidity(), b.azimuthalAngle());
	default:
	throw std::runtime_error("The specified deltaR scheme is not yet implemented");
	}
	}

	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a
	/// three-vector and a four-vector.
	inline double deltaR(const FourMomentum& a, const Vector3& b) {
	return deltaR(a.vector3(), b);
	}

	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a
	/// three-vector and a four-vector.
	inline double deltaR(const Vector3& a, const FourMomentum& b) {
	return deltaR(a, b.vector3());
	}

	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a
	/// three-vector and a four-vector.
	inline double deltaR(const FourVector& a, const Vector3& b) {
	return deltaR(a.vector3(), b);
	}

	/// @brief Calculate the 2D rapidity-azimuthal ("eta-phi") distance between a
	/// three-vector and a four-vector.
	inline double deltaR(const Vector3& a, const FourVector& b) {
	return deltaR(a, b.vector3());
	}

	//@}

	/// @name \f$ \Delta phi \f$ calculations from 4-vectors
	//@{

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(const FourMomentum& a, const FourMomentum& b) {
	return deltaPhi(a.vector3(), b.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(const FourMomentum& v, double phi2) {
	return deltaPhi(v.vector3(), phi2);
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(double phi1, const FourMomentum& v) {
	return deltaPhi(phi1, v.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(const FourVector& a, const FourVector& b) {
	return deltaPhi(a.vector3(), b.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(const FourVector& v, double phi2) {
	return deltaPhi(v.vector3(), phi2);
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(double phi1, const FourVector& v) {
	return deltaPhi(phi1, v.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(const FourVector& a, const FourMomentum& b) {
	return deltaPhi(a.vector3(), b.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(const FourMomentum& a, const FourVector& b) {
	return deltaPhi(a.vector3(), b.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(const FourVector& a, const Vector3& b) {
	return deltaPhi(a.vector3(), b);
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(const Vector3& a, const FourVector& b) {
	return deltaPhi(a, b.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(const FourMomentum& a, const Vector3& b) {
	return deltaPhi(a.vector3(), b);
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaPhi(const Vector3& a, const FourMomentum& b) {
	return deltaPhi(a, b.vector3());
	}

	//@}

	/// @name \f$ \|\Delta eta\| \f$ calculations from 4-vectors
	//@{

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(const FourMomentum& a, const FourMomentum& b) {
	return deltaEta(a.vector3(), b.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(const FourMomentum& v, double eta2) {
	return deltaEta(v.vector3(), eta2);
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(double eta1, const FourMomentum& v) {
	return deltaEta(eta1, v.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(const FourVector& a, const FourVector& b) {
	return deltaEta(a.vector3(), b.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(const FourVector& v, double eta2) {
	return deltaEta(v.vector3(), eta2);
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(double eta1, const FourVector& v) {
	return deltaEta(eta1, v.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(const FourVector& a, const FourMomentum& b) {
	return deltaEta(a.vector3(), b.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(const FourMomentum& a, const FourVector& b) {
	return deltaEta(a.vector3(), b.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(const FourVector& a, const Vector3& b) {
	return deltaEta(a.vector3(), b);
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(const Vector3& a, const FourVector& b) {
	return deltaEta(a, b.vector3());
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(const FourMomentum& a, const Vector3& b) {
	return deltaEta(a.vector3(), b);
	}

	/// Calculate the difference in azimuthal angle between two spatial vectors.
	inline double deltaEta(const Vector3& a, const FourMomentum& b) {
	return deltaEta(a, b.vector3());
	}

	//@}

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


	/// @name 4-vector string representations
	//@{

	/// Render a 4-vector as a string.
	inline std::string toString(const FourVector& lv) {
	ostringstream out;
	out << "(" << (fabs(lv.t()) < 1E-30 ? 0.0 : lv.t())
	<< "; " << (fabs(lv.x()) < 1E-30 ? 0.0 : lv.x())
	<< ", " << (fabs(lv.y()) < 1E-30 ? 0.0 : lv.y())
	<< ", " << (fabs(lv.z()) < 1E-30 ? 0.0 : lv.z())
	<< ")";
	return out.str();
	}

	/// Write a 4-vector to an ostream.
	inline std::ostream& operator<<(std::ostream& out, const FourVector& lv) {
	out << toString(lv);
	return out;
	}

	//@}


	}

	#endif
	diff --git a/include/Rivet/Projections/MissingMomentum.hh b/include/Rivet/Projections/MissingMomentum.hh
	--- a/include/Rivet/Projections/MissingMomentum.hh
	+++ b/include/Rivet/Projections/MissingMomentum.hh
	@@ -1,90 +1,92 @@
	// -- C++ --
	#ifndef RIVET_MissingMomentum_HH
	#define RIVET_MissingMomentum_HH

	#include "Rivet/Rivet.hh"
	#include "Rivet/Projection.hh"
	#include "Rivet/Projections/VisibleFinalState.hh"
	#include "Rivet/Particle.hh"
	#include "Rivet/Event.hh"

	namespace Rivet {


	/// @brief Calculate missing \f$ E \f$, \f$ E_\perp \f$ etc.
	///
	/// Project out the total visible energy vector, allowing missing
	/// \f$ E \f$, \f$ E_\perp \f$ etc. to be calculated. Final state
	/// visibility restrictions are automatic.
	class MissingMomentum : public Projection {
	public:

	/// Default constructor with uncritical FS.
	MissingMomentum()
	{
	setName("MissingMomentum");
	FinalState fs;
	addProjection(fs, "FS");
	addProjection(VisibleFinalState(fs), "VisibleFS");
	}


	/// Constructor.
	MissingMomentum(const FinalState& fs)
	{
	setName("MissingMomentum");
	addProjection(fs, "FS");
	addProjection(VisibleFinalState(fs), "VisibleFS");
	}


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


	public:

	/// The vector-summed visible four-momentum in the event.
	+ /// @note Reverse this vector with operator- to get the missing momentum vector.
	const FourMomentum& visibleMomentum() const { return _momentum; }

	- /// The vector-summed (in)visible transverse energy in the event
	+ /// The vector-summed visible transverse energy in the event
	+ /// @note Reverse this vector with operator- to get the missing ET vector.
	const Vector3& vectorEt() const { return _vet; }

	/// The scalar-summed visible transverse energy in the event.
	double scalarEt() const { return _set; }


	protected:

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

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


	public:

	/// Clear the projection results.
	void clear();


	private:

	/// The total visible momentum
	FourMomentum _momentum;

	/// Scalar transverse energy
	double _set;

	/// Vector transverse energy
	Vector3 _vet;

	};


	}

	#endif

File Metadata

Mime Type: text/x-diff
Expires: Tue, Sep 30, 5:41 AM (1 h, 9 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6566238
Default Alt Text: (50 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions