Page MenuHomeHEPForge
Files F7879050

No OneTemporary
Actions

View Options

diff --git a/include/Rivet/Projections/Sphericity.hh b/include/Rivet/Projections/Sphericity.hh
--- a/include/Rivet/Projections/Sphericity.hh
+++ b/include/Rivet/Projections/Sphericity.hh
@@ -1,157 +1,163 @@
// -*- C++ -*-
#ifndef RIVET_Sphericity_HH
#define RIVET_Sphericity_HH
#include "Rivet/Projection.hh"
#include "Rivet/Projections/AxesDefinition.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Event.hh"
+#include "Rivet/Jet.fhh"
namespace Rivet {
+
/// @brief Calculate the sphericity event shape.
///
/// The sphericity tensor (or quadratic momentum tensor) is defined as
/// \f[
/// S^{\alpha \beta} = \frac{\sum_i p_i^\alpha p_i^\beta}{\sum_i |\mathbf{p}_i|^2}
/// \f],
/// where the Greek indices are spatial components and the Latin indices are used
/// for sums over particles. From this, the sphericity, aplanarity and planarity can be
/// calculated by combinations of eigenvalues.
///
/// Defining the three eigenvalues
/// \f$ \lambda_1 \ge \lambda_2 \ge \lambda_3 \f$, with \f$ \lambda_1 + \lambda_2 + \lambda_3 = 1 \f$,
/// the sphericity is
/// \f[
/// S = \frac{3}{2} (\lambda_2 + \lambda_3)
/// \f]
///
/// The aplanarity is \f$ A = \frac{3}{2}\lambda_3 \f$ and the planarity
/// is \f$ P = \frac{2}{3}(S-2A) = \lambda_2 - \lambda_3 \f$. The eigenvectors define a
/// set of spatial axes comparable with the thrust axes, but more sensitive to
/// high momentum particles due to the quadratic sensitivity of the tensor to
/// the particle momenta.
///
/// Since the sphericity is quadratic in the particle momenta, it is not an
/// infrared safe observable in perturbative QCD. This can be fixed by adding
/// a regularizing power of \f$r\f$ to the definition:
/// \f[
/// S^{\alpha \beta} =
/// \frac{\sum_i |\mathbf{p}_i|^{r-2} p_i^\alpha p_i^\beta}
/// {\sum_i |\mathbf{p}_i|^r}
/// \f]
///
/// \f$r\f$ is available as a constructor argument on this class and will be
/// taken into account by the Cmp<Projection> operation, so a single analysis
/// can use several sphericity projections with different \f$r\f$ values without
/// fear of a clash.
///
class Sphericity : public AxesDefinition {
public:
/// @name Constructors etc.
//@{
/// Constructor
Sphericity(double rparam=2.0): _regparam(rparam){}
Sphericity(const FinalState& fsp, double rparam=2.0);
/// Clone on the heap.
DEFAULT_RIVET_PROJ_CLONE(Sphericity);
//@}
protected:
/// Perform the projection on the Event
void project(const Event& e);
/// Compare with other projections
int compare(const Projection& p) const;
+
public:
/// Reset the projection
void clear();
/// @name Access the event shapes by name
- /// @{
+ //@{
/// Sphericity
double sphericity() const { return 3.0 / 2.0 * (lambda2() + lambda3()); }
- /// Transverse Sphericity
+ /// Transverse sphericity
double transSphericity() const { return 2.0 * lambda2() / ( lambda1() + lambda2() ); }
/// Planarity
double planarity() const { return 2 * (sphericity() - 2 * aplanarity()) / 3.0; }
/// Aplanarity
double aplanarity() const { return 3 / 2.0 * lambda3(); }
- /// @}
+ //@}
+
/// @name Access the sphericity basis vectors
- /// @{
+ //@{
/// Sphericity axis
const Vector3& sphericityAxis() const { return _sphAxes[0]; }
/// Sphericity major axis
const Vector3& sphericityMajorAxis() const { return _sphAxes[1]; }
/// Sphericity minor axis
const Vector3& sphericityMinorAxis() const { return _sphAxes[2]; }
- /// @}
+ //@}
- ///@{ AxesDefinition axis accessors.
+
+ /// @name AxesDefinition axis accessors
+ //@{
const Vector3& axis1() const { return sphericityAxis(); }
const Vector3& axis2() const { return sphericityMajorAxis(); }
const Vector3& axis3() const { return sphericityMinorAxis(); }
- ///@}
+ //@}
/// @name Access the momentum tensor eigenvalues
- /// @{
+ //@{
double lambda1() const { return _lambdas[0]; }
double lambda2() const { return _lambdas[1]; }
double lambda3() const { return _lambdas[2]; }
- /// @}
+ //@}
/// @name Direct methods
/// Ways to do the calculation directly, without engaging the caching system
//@{
/// Manually calculate the sphericity, without engaging the caching system
void calc(const FinalState& fs);
/// Manually calculate the sphericity, without engaging the caching system
- void calc(const vector<Particle>& fsparticles);
+ void calc(const Particles& particles);
/// Manually calculate the sphericity, without engaging the caching system
- void calc(const vector<FourMomentum>& fsmomenta);
+ void calc(const Jets& jets);
/// Manually calculate the sphericity, without engaging the caching system
- void calc(const vector<Vector3>& fsmomenta);
+ void calc(const vector<FourMomentum>& momenta);
+
+ /// @brief Manually calculate the sphericity, without engaging the caching system
+ ///
+ /// This one actually does the calculation
+ void calc(const vector<Vector3>& momenta);
//@}
+ private:
- private:
/// Eigenvalues.
vector<double> _lambdas;
/// Sphericity axes.
vector<Vector3> _sphAxes;
/// Regularizing parameter, used to force infra-red safety.
const double _regparam;
- private:
+ };
- /// Actually do the calculation
- void _calcSphericity(const vector<Vector3>& fsmomenta);
-
- };
}
-
#endif
diff --git a/src/Projections/Sphericity.cc b/src/Projections/Sphericity.cc
--- a/src/Projections/Sphericity.cc
+++ b/src/Projections/Sphericity.cc
@@ -1,137 +1,137 @@
// -*- C++ -*-
#include "Rivet/Projections/Sphericity.hh"
+#include "Rivet/Jet.hh"
namespace Rivet {
Sphericity::Sphericity(const FinalState& fsp, double rparam)
: _regparam(rparam)
{
setName("Sphericity");
addProjection(fsp, "FS");
clear();
}
void Sphericity::clear() {
_lambdas = vector<double>(3, 0);
_sphAxes = vector<Vector3>(3, Vector3());
}
int Sphericity::compare(const Projection& p) const {
PCmp fscmp = mkNamedPCmp(p, "FS");
if (fscmp != EQUIVALENT) return fscmp;
const Sphericity& other = dynamic_cast<const Sphericity&>(p);
if (fuzzyEquals(_regparam, other._regparam)) return 0;
return cmp(_regparam, other._regparam);
}
void Sphericity::project(const Event& e) {
const Particles prts = applyProjection<FinalState>(e, "FS").particles();
calc(prts);
}
void Sphericity::calc(const FinalState& fs) {
calc(fs.particles());
}
- void Sphericity::calc(const vector<Particle>& fsparticles) {
+
+ void Sphericity::calc(const Particles& particles) {
vector<Vector3> threeMomenta;
- threeMomenta.reserve(fsparticles.size());
- foreach (const Particle& p, fsparticles) {
- threeMomenta.push_back( p.momentum().vector3() );
- }
- _calcSphericity(threeMomenta);
+ transform(particles, threeMomenta, p3);
+ calc(threeMomenta);
}
- void Sphericity::calc(const vector<FourMomentum>& fsmomenta) {
+
+ void Sphericity::calc(const Jets& jets) {
vector<Vector3> threeMomenta;
- threeMomenta.reserve(fsmomenta.size());
- foreach (const FourMomentum& v, fsmomenta) {
- threeMomenta.push_back(v.vector3());
- }
- _calcSphericity(threeMomenta);
+ transform(jets, threeMomenta, p3);
+ calc(threeMomenta);
}
- void Sphericity::calc(const vector<Vector3>& fsmomenta) {
- _calcSphericity(fsmomenta);
+
+ void Sphericity::calc(const vector<FourMomentum>& momenta) {
+ vector<Vector3> threeMomenta;
+ transform(momenta, threeMomenta, [](const FourMomentum& p4){return p4.vector3();});
+ calc(threeMomenta);
}
- // Actually do the calculation
- void Sphericity::_calcSphericity(const vector<Vector3>& fsmomenta) {
+ void Sphericity::calc(const vector<Vector3>& momenta) {
MSG_DEBUG("Calculating sphericity with r = " << _regparam);
- // Return (with "safe nonsense" sphericity params) if there are no final state particles.
- if (fsmomenta.empty()) {
- MSG_DEBUG("No particles in final state...");
+ // Return (with "safe nonsense" sphericity params) if there are no final state particles
+ if (momenta.empty()) {
+ MSG_DEBUG("No momenta given...");
clear();
return;
}
// Iterate over all the final state particles.
Matrix3 mMom;
double totalMomentum = 0.0;
- MSG_DEBUG("Number of particles = " << fsmomenta.size());
- foreach (const Vector3& p3, fsmomenta) {
+ MSG_DEBUG("Number of particles = " << momenta.size());
+ for (const Vector3& p3 : momenta) {
// Build the (regulated) normalising factor.
totalMomentum += pow(p3.mod(), _regparam);
// Build (regulated) quadratic momentum components.
const double regfactor = pow(p3.mod(), _regparam-2);
if (!fuzzyEquals(regfactor, 1.0)) {
MSG_TRACE("Regfactor (r=" << _regparam << ") = " << regfactor);
}
Matrix3 mMomPart;
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; ++j) {
mMomPart.set(i,j, p3[i]*p3[j]);
}
}
mMom += regfactor * mMomPart;
}
// Normalise to total (regulated) momentum.
mMom /= totalMomentum;
MSG_DEBUG("Momentum tensor = " << "\n" << mMom);
// Check that the matrix is symmetric.
const bool isSymm = mMom.isSymm();
if (!isSymm) {
MSG_ERROR("Error: momentum tensor not symmetric (r=" << _regparam << ")");
MSG_ERROR("[0,1] vs. [1,0]: " << mMom.get(0,1) << ", " << mMom.get(1,0));
MSG_ERROR("[0,2] vs. [2,0]: " << mMom.get(0,2) << ", " << mMom.get(2,0));
MSG_ERROR("[1,2] vs. [2,1]: " << mMom.get(1,2) << ", " << mMom.get(2,1));
}
// If not symmetric, something's wrong (we made sure the error msg appeared first).
assert(isSymm);
// Diagonalize momentum matrix.
const EigenSystem<3> eigen3 = diagonalize(mMom);
MSG_DEBUG("Diag momentum tensor = " << "\n" << eigen3.getDiagMatrix());
// Reset and set eigenvalue/vector parameters.
_lambdas.clear();
_sphAxes.clear();
const EigenSystem<3>::EigenPairs epairs = eigen3.getEigenPairs();
assert(epairs.size() == 3);
for (size_t i = 0; i < 3; ++i) {
_lambdas.push_back(epairs[i].first);
_sphAxes.push_back(Vector3(epairs[i].second));
}
// Debug output.
MSG_DEBUG("Lambdas = ("
<< lambda1() << ", " << lambda2() << ", " << lambda3() << ")");
MSG_DEBUG("Sum of lambdas = " << lambda1() + lambda2() + lambda3());
MSG_DEBUG("Vectors = "
<< sphericityAxis() << ", "
<< sphericityMajorAxis() << ", "
<< sphericityMinorAxis() << ")");
}
+
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:23 PM (1 d, 10 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3798643
Default Alt Text
(10 KB)

Event Timeline