|
|
[Rivet-svn] r3068 - in trunk: include/Rivet include/Rivet/Projections src/Projections
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Fri Apr 29 12:33:13 BST 2011
Author: holsch
Date: Fri Apr 29 12:33:13 2011
New Revision: 3068
Log:
Add Spherocity and FParameter Projections fromthe Banfi paper
Added:
trunk/include/Rivet/Projections/FParameter.hh
trunk/include/Rivet/Projections/Spherocity.hh
trunk/src/Projections/FParameter.cc
trunk/src/Projections/Spherocity.cc
Modified:
trunk/include/Rivet/Makefile.am
trunk/src/Projections/Makefile.am
Modified: trunk/include/Rivet/Makefile.am
==============================================================================
--- trunk/include/Rivet/Makefile.am Tue Apr 26 23:13:36 2011 (r3067)
+++ trunk/include/Rivet/Makefile.am Fri Apr 29 12:33:13 2011 (r3068)
@@ -44,6 +44,7 @@
Projections/FinalState.hh \
Projections/FinalStateHCM.hh \
Projections/FoxWolframMoments.hh \
+ Projections/FParameter.hh \
Projections/HadronicFinalState.hh \
Projections/Hemispheres.hh \
Projections/IdentifiedFinalState.hh \
@@ -64,6 +65,7 @@
Projections/ParisiTensor.hh \
Projections/PVertex.hh \
Projections/Sphericity.hh \
+ Projections/Spherocity.hh \
Projections/SVertex.hh \
Projections/Thrust.hh \
Projections/TotalVisibleMomentum.hh \
Added: trunk/include/Rivet/Projections/FParameter.hh
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/include/Rivet/Projections/FParameter.hh Fri Apr 29 12:33:13 2011 (r3068)
@@ -0,0 +1,89 @@
+// -*- C++ -*-
+#ifndef RIVET_Sphericity_HH
+#define RIVET_FParameter_HH
+
+#include "Rivet/Projection.hh"
+#include "Rivet/Projections/FinalState.hh"
+#include "Rivet/Event.hh"
+
+
+namespace Rivet {
+
+ class FParameter : public Projection {
+
+ public:
+
+ /// @name Constructors etc.
+ //@{
+
+ /// Constructor
+ FParameter(const FinalState& fsp);
+
+ /// Clone on the heap.
+ virtual const Projection* clone() const {
+ return new FParameter(*this);
+ }
+
+ //@}
+
+
+ protected:
+
+ /// Perform the projection on the Event
+ void project(const Event& e);
+
+ /// Compare with other projections
+ //int compare(const Projection& p) const;
+ // Taken from Thrust.hh
+ int compare(const Projection& p) const {
+ return mkNamedPCmp(p, "FS");
+ }
+ public:
+
+ /// Reset the projection
+ void clear();
+
+ /// @name Access the event shapes by name
+ /// @{
+ /// F-Parametr
+ double F() const { return lambda1() >= lambda2() ? lambda2()/lambda1() : lambda1()/lambda2(); }
+ /// @}
+
+ /// @name Access the linearised transverse momentum tensor eigenvalues
+ /// @{
+ double lambda1() const { return _lambdas[0]; }
+ double lambda2() const { return _lambdas[1]; }
+ /// @}
+
+
+ /// @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);
+
+ /// Manually calculate the sphericity, without engaging the caching system
+ void calc(const vector<FourMomentum>& fsmomenta);
+
+ /// Manually calculate the sphericity, without engaging the caching system
+ void calc(const vector<Vector3>& fsmomenta);
+
+ //@}
+ private:
+ /// Eigenvalues.
+ vector<double> _lambdas;
+
+ private:
+
+ /// Actually do the calculation
+ void _calcFParameter(const vector<Vector3>& fsmomenta);
+
+ };
+}
+
+
+#endif
Added: trunk/include/Rivet/Projections/Spherocity.hh
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/include/Rivet/Projections/Spherocity.hh Fri Apr 29 12:33:13 2011 (r3068)
@@ -0,0 +1,128 @@
+// -*- C++ -*-
+#ifndef RIVET_Spherocity_HH
+#define RIVET_Spherocity_HH
+
+#include "Rivet/Projection.hh"
+#include "Rivet/Projections/AxesDefinition.hh"
+#include "Rivet/Projections/FinalState.hh"
+#include "Rivet/Event.hh"
+
+namespace Rivet {
+
+
+ /**
+ @brief Get the transverse spherocity scalars for hadron-colliders.
+
+ @author Holger Schulz
+
+ The scalar (maximum) transverse spherocity is defined as
+ \f[
+ T = \mathrm{max}_{\vec{n_\perp}} \frac{\sum_i \left|\vec{p}_{\perp,i} \cdot \vec{n} \right|}{\sum_i |\vec{p}_{\perp,i}|}
+ \f],
+ with the direction of the unit vector \f$ \vec{n_\perp} \f$ which maximises \f$ T \f$
+ being identified as the spherocity axis. The unit vector which maximises the spherocity
+ scalar in the plane perpendicular to \f$ \vec{n} \f$ is the "spherocity major"
+ direction, and the vector perpendicular to both the spherocity and spherocity major directions
+ is the spherocity minor. Both the major and minor directions have associated spherocity
+ scalars.
+
+ Care must be taken in the case of Drell-Yann processes - there we should use the
+ newly proposed observable a_T.
+
+ */
+ class Spherocity : public AxesDefinition {
+ public:
+
+ /// Constructor.
+ Spherocity(const FinalState& fsp)
+ : _calculatedSpherocity(false)
+ {
+ setName("Spherocity");
+ addProjection(fsp, "FS");
+ }
+
+ /// Clone on the heap.
+ virtual const Projection* clone() const {
+ return new Spherocity(*this);
+ }
+
+ protected:
+
+ /// Perform the projection on the Event
+ void project(const Event& e) {
+ const vector<Particle> ps
+ = applyProjection<FinalState>(e, "FS").particles();
+ calc(ps);
+ }
+
+ /// Compare projections
+ int compare(const Projection& p) const {
+ return mkNamedPCmp(p, "FS");
+ }
+
+
+ public:
+
+ ///@{ Spherocity scalar accessors
+ /// The spherocity scalar, \f$ S \f$, (minimum spherocity).
+ double spherocity() const { return _spherocities[0]; }
+ ///@}
+
+ ///@{ Spherocity axis accessors
+ /// The spherocity axis.
+ const Vector3& spherocityAxis() const { return _spherocityAxes[0]; }
+ /// The spherocity major axis (axis of max spherocity perpendicular to spherocity axis).
+ const Vector3& spherocityMajorAxis() const { return _spherocityAxes[1]; }
+ /// The spherocity minor axis (axis perpendicular to spherocity and spherocity major).
+ const Vector3& spherocityMinorAxis() const { return _spherocityAxes[2]; }
+ ///@}
+
+ ///@{ AxesDefinition axis accessors.
+ const Vector3& axis1() const { return spherocityAxis(); }
+ const Vector3& axis2() const { return spherocityMajorAxis(); }
+ const Vector3& axis3() const { return spherocityMinorAxis(); }
+ ///@}
+
+
+ public:
+
+ /// @name Direct methods
+ /// Ways to do the calculation directly, without engaging the caching system
+ //@{
+
+ /// Manually calculate the spherocity, without engaging the caching system
+ void calc(const FinalState& fs);
+
+ /// Manually calculate the spherocity, without engaging the caching system
+ void calc(const vector<Particle>& fsparticles);
+
+ /// Manually calculate the spherocity, without engaging the caching system
+ void calc(const vector<FourMomentum>& fsmomenta);
+
+ /// Manually calculate the spherocity, without engaging the caching system
+ void calc(const vector<Vector3>& threeMomenta);
+
+ //@}
+
+
+ private:
+
+ /// The spherocity scalars.
+ vector<double> _spherocities;
+
+ /// The spherocity axes.
+ vector<Vector3> _spherocityAxes;
+
+ /// Caching flag to avoid costly recalculations.
+ bool _calculatedSpherocity;
+
+ private:
+
+ /// Explicitly calculate the spherocity values.
+ void _calcSpherocity(const vector<Vector3>& fsmomenta);
+
+ };
+
+}
+
+#endif
Added: trunk/src/Projections/FParameter.cc
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Projections/FParameter.cc Fri Apr 29 12:33:13 2011 (r3068)
@@ -0,0 +1,116 @@
+// -*- C++ -*-
+#include "Rivet/Projections/FParameter.hh"
+#include "Rivet/Cmp.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/Tools/Utils.hh"
+
+namespace Rivet {
+
+ FParameter::FParameter(const FinalState& fsp) {
+ setName("FParameter");
+ addProjection(fsp, "FS");
+ clear();
+ }
+
+
+ void FParameter::clear() {
+ _lambdas = vector<double>(2, 0);
+ }
+
+
+ void FParameter::project(const Event& e) {
+ const ParticleVector prts = applyProjection<FinalState>(e, "FS").particles();
+ calc(prts);
+ }
+
+
+ void FParameter::calc(const FinalState& fs) {
+ calc(fs.particles());
+ }
+
+ void FParameter::calc(const vector<Particle>& fsparticles) {
+ vector<Vector3> threeMomenta;
+ threeMomenta.reserve(fsparticles.size());
+ foreach (const Particle& p, fsparticles) {
+ const Vector3 p3 = p.momentum().vector3();
+ threeMomenta.push_back(p3);
+ }
+ _calcFParameter(threeMomenta);
+ }
+
+ void FParameter::calc(const vector<FourMomentum>& fsmomenta) {
+ vector<Vector3> threeMomenta;
+ threeMomenta.reserve(fsmomenta.size());
+ foreach (const FourMomentum& v, fsmomenta) {
+ threeMomenta.push_back(v.vector3());
+ }
+ _calcFParameter(threeMomenta);
+ }
+
+ void FParameter::calc(const vector<Vector3>& fsmomenta) {
+ _calcFParameter(fsmomenta);
+ }
+
+ // Actually do the calculation
+ void FParameter::_calcFParameter(const vector<Vector3>& fsmomenta) {
+
+ // Return (with "safe nonsense" sphericity params) if there are no final state particles.
+ if (fsmomenta.empty()) {
+ getLog() << Log::DEBUG << "No particles in final state..." << endl;
+ clear();
+ return;
+ }
+
+ // A small iteration over full momenta but set z-coord. to 0.0 to get transverse momenta
+ vector <Vector3> fsperpmomenta;
+ foreach (const Vector3& p, fsmomenta) {
+ fsperpmomenta.push_back(Vector3(p.x(), p.y(), 0.0));
+ }
+
+ // Iterate over all the final state particles.
+ Matrix<2> mMom;
+ getLog() << Log::DEBUG << "Number of particles = " << fsperpmomenta.size() << endl;
+ foreach (const Vector3& p3, fsperpmomenta) {
+
+ double prefactor = 1.0/mod(p3);
+
+ Matrix<2> mMomPart;
+ for (size_t i = 0; i < 2; ++i) {
+ for (size_t j = 0; j < 2; ++j) {
+ mMomPart.set(i,j, p3[i]*p3[j]);
+ }
+ }
+ mMom += prefactor * mMomPart;
+ }
+
+ getLog() << Log::DEBUG << "Linearised transverse momentum tensor = " << endl << mMom << endl;
+
+ // Check that the matrix is symmetric.
+ const bool isSymm = mMom.isSymm();
+ if (!isSymm) {
+ getLog() << Log::ERROR << "Error: momentum tensor not symmetric:" << endl;
+ getLog() << Log::ERROR << "[0,1] vs. [1,0]: " << mMom.get(0,1) << ", " << mMom.get(1,0) << endl;
+ getLog() << Log::ERROR << "[0,2] vs. [2,0]: " << mMom.get(0,2) << ", " << mMom.get(2,0) << endl;
+ getLog() << Log::ERROR << "[1,2] vs. [2,1]: " << mMom.get(1,2) << ", " << mMom.get(2,1) << endl;
+ }
+ // If not symmetric, something's wrong (we made sure the error msg appeared first).
+ assert(isSymm);
+
+ // Diagonalize momentum matrix.
+ const EigenSystem<2> eigen2 = diagonalize(mMom);
+ getLog() << Log::DEBUG << "Diag momentum tensor = " << endl << eigen2.getDiagMatrix() << endl;
+
+ // Reset and set eigenvalue parameters.
+ _lambdas.clear();
+ const EigenSystem<2>::EigenPairs epairs = eigen2.getEigenPairs();
+ assert(epairs.size() == 2);
+ for (size_t i = 0; i < 2; ++i) {
+ _lambdas.push_back(epairs[i].first);
+ }
+
+ // Debug output.
+ getLog() << Log::DEBUG << "Lambdas = ("
+ << lambda1() << ", " << lambda2() << ")" << endl;
+ getLog() << Log::DEBUG << "Sum of lambdas = " << lambda1() + lambda2() << endl;
+ }
+}
Modified: trunk/src/Projections/Makefile.am
==============================================================================
--- trunk/src/Projections/Makefile.am Tue Apr 26 23:13:36 2011 (r3067)
+++ trunk/src/Projections/Makefile.am Fri Apr 29 12:33:13 2011 (r3068)
@@ -11,6 +11,7 @@
FinalState.cc \
FinalStateHCM.cc \
FoxWolframMoments.cc \
+ FParameter.cc \
HadronicFinalState.cc \
Hemispheres.cc \
IdentifiedFinalState.cc \
@@ -27,6 +28,7 @@
ParisiTensor.cc \
PVertex.cc \
Sphericity.cc \
+ Spherocity.cc \
SVertex.cc \
Thrust.cc \
TotalVisibleMomentum.cc \
Added: trunk/src/Projections/Spherocity.cc
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Projections/Spherocity.cc Fri Apr 29 12:33:13 2011 (r3068)
@@ -0,0 +1,213 @@
+// -*- C++ -*-
+#include "Rivet/Rivet.hh"
+#include "Rivet/Projections/Spherocity.hh"
+#include "Rivet/Tools/Logging.hh"
+
+
+namespace Rivet {
+
+
+ void Spherocity::calc(const FinalState& fs) {
+ calc(fs.particles());
+ }
+
+ void Spherocity::calc(const vector<Particle>& fsparticles) {
+ vector<Vector3> threeMomenta;
+ threeMomenta.reserve(fsparticles.size());
+ foreach (const Particle& p, fsparticles) {
+ const Vector3 p3 = p.momentum().vector3();
+ threeMomenta.push_back(p3);
+ }
+ _calcSpherocity(threeMomenta);
+ }
+
+ void Spherocity::calc(const vector<FourMomentum>& fsmomenta) {
+ vector<Vector3> threeMomenta;
+ threeMomenta.reserve(fsmomenta.size());
+ foreach (const FourMomentum& v, fsmomenta) {
+ threeMomenta.push_back(v.vector3());
+ }
+ _calcSpherocity(threeMomenta);
+ }
+
+ void Spherocity::calc(const vector<Vector3>& fsmomenta) {
+ _calcSpherocity(fsmomenta);
+ }
+
+
+ /////////////////////////////////////////////////
+
+
+ inline bool mod2Cmp(const Vector<2>& a, const Vector<2>& b) {
+ return a.mod2() > b.mod2();
+ }
+
+ // Inline functions to avoid template hell
+ inline double dot(const Vector<2>& a, const Vector<2>& b) {
+ return a[0]*b[0] + a[1]*b[1];
+ }
+
+ inline Vector<2> unit(const Vector<2>& v) {
+ assert(mod(v) > 0);
+ Vector<2> returnthis;
+ returnthis.set(0, v[0]/mod(v));
+ returnthis.set(1, v[1]/mod(v));
+ return returnthis;
+ }
+
+
+
+ //// Do the general case spherocity calculation
+ void _calcS(const vector<Vector3 >& perpmomenta, double& sphero, Vector3& saxis) {
+ // According to the Salam paper, p5, footnote 4, the
+ // axis n that minimises the Spherocity value ALWAYS coincides with the
+ // direction of one of the transverse momentum vectors of the events particles.
+ // Thus we carry out the calculation of Sphero for all particles and pick the
+ // one that yields the smallerst values
+
+ vector<Vector3> p = perpmomenta;
+ vector<double> sval;
+
+
+ // Prepare vector to store unit vector representation of all particle momenta
+ // and also calculate the transverse momentum sum
+ vector<Vector3> units;
+ double sumperpmomenta = 0.0;
+ foreach (const Vector3& p, perpmomenta) {
+ units.push_back(Vector3(p.x(), p.y(), 0.0).unit());
+ sumperpmomenta += p.mod();
+ }
+
+ // Spherocity calculation
+ //
+ // The outer loop is for iteration over all unit vectors
+ foreach (const Vector3& u, units){
+ double s =0.0;
+ for (unsigned int k=0 ; k<p.size() ; k++)
+ s += fabs(p[k].cross(u).mod() );
+
+ sval.push_back(s);
+ }
+
+
+ // Pick the solution with the smallest spherocity
+ sphero = 999.;
+ for (unsigned int i=0 ; i<units.size() ; i++)
+ if (sval[i] < sphero){
+ sphero = sval[i];
+ saxis = units[i];
+ }
+
+ }
+
+
+
+ // Do the full calculation
+ void Spherocity::_calcSpherocity(const vector<Vector3>& fsmomenta) {
+
+ // Make a vector of the three-momenta in the final state
+ // Explicitly set the z-component (parallel to beam axis) to zero
+ // This creates a 3D-vector representation of the transverse momentum
+ // but take the full information momentum vectors as input
+
+ vector<Vector3> fsperpmomenta;
+ // A small iteration over full momenta but set z-coord. to 0.0 to get transverse momenta
+ foreach (const Vector3& p, fsmomenta) {
+ fsperpmomenta.push_back(Vector3(p.x(), p.y(), 0.0));
+ }
+
+ // This returns the scalar sum of (transverse) momenta
+ double perpmomentumSum(0.0);
+ foreach (const Vector3& p, fsperpmomenta) {
+ perpmomentumSum += mod(p);
+ }
+
+ // Clear the caches
+ _spherocities.clear();
+ _spherocityAxes.clear();
+
+ // If there are fewer than 2 visible particles, we can't do much
+ // This is blindly copied from the Thrust projection
+ if (fsmomenta.size() < 2) {
+ for (int i = 0; i < 3; ++i) {
+ _spherocities.push_back(-1);
+ _spherocityAxes.push_back(Vector3(0,0,0));
+ }
+ return;
+ }
+
+ // Handle special case of spherocity = 1 if there are only 2 particles
+ // This is blindly copied from the Thrust projection
+ if (fsmomenta.size() == 2) {
+ Vector3 axis(0,0,0);
+ _spherocities.push_back(1.0);
+ _spherocities.push_back(0.0);
+ _spherocities.push_back(0.0);
+ axis = fsmomenta[0].unit();
+ if (axis.z() < 0) axis = -axis;
+ _spherocityAxes.push_back(axis);
+ /// @todo Improve this --- special directions bad...
+ /// (a,b,c) _|_ 1/(a^2+b^2) (b,-a,0) etc., but which combination minimises error?
+ if (axis.z() < 0.75)
+ _spherocityAxes.push_back( (axis.cross(Vector3(0,0,1))).unit() );
+ else
+ _spherocityAxes.push_back( (axis.cross(Vector3(0,1,0))).unit() );
+ _spherocityAxes.push_back( _spherocityAxes[0].cross(_spherocityAxes[1]) );
+ return;
+ }
+
+
+ // Temporary variables for calcs
+ Vector3 axis(0,0,0);
+ double val = 0.;
+
+
+ // Get spherocity
+ _calcS(fsperpmomenta, val, axis);
+ getLog() << Log::DEBUG << "Mom sum = " << perpmomentumSum << endl;
+ double spherocity = PI*PI* val*val / (4 * perpmomentumSum*perpmomentumSum);
+ _spherocities.push_back(spherocity);
+
+ // See if calclulated spherocity value makes sense
+ if (spherocity < 0.0 || spherocity > 1.0) {
+ getLog() << Log::WARN << "Spherocity = " << spherocity << endl;
+ }
+
+ getLog() << Log::DEBUG << "Spherocity value = " << spherocity << endl;
+
+ getLog() << Log::DEBUG << "Sperocity axis = " << axis << endl;
+
+ _spherocityAxes.push_back(axis);
+
+
+ //// Get spherocity minor
+ ////
+ ////
+ //// The Beam axis is eZ = (0, 0, 1)
+ ////
+ //double perpMinor = 0.0;
+ //const Vector3 ez = Vector3(0, 0, 1);
+ //foreach (const Vector3& v, fsperpmomenta) {
+ //Vector3 temp = Vector3(v[0], v[1], 0.0);
+ //perpMinor += mod(temp.cross(_spherocityAxes[0]));
+ //}
+ //_spherocitys.push_back(perpMinor / perpmomentumSum);
+
+
+ //// Manual check
+ //double test = 0.;
+ //Vector<2> ex;
+ //ex.set(0, 1.0);
+ //ex.set(1, 0);
+ //foreach (const Vector<2> & v, fsperpmomenta2D) {
+ //test+=fabs(dot(ex, v));
+ //}
+ //if (test/perpmomentumSum < _spherocities[0]) {
+ //std::cout << "Warning: " << test/perpmomentumSum << " > " << _spherocitys[0] << " " << _spherocityAxes[0] << endl;
+ //}
+
+ }
+
+
+
+}
More information about the Rivet-svn
mailing list
|