#ifndef __JSIRENE__JEVENTSHAPEVARIABLES__
#define __JSIRENE__JEVENTSHAPEVARIABLES__
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"
#include "JMath/JMatrix3S.hh"
#include "JMath/JMathSupportkit.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JPhysics/JPhysicsToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JSirene/JVisibleEnergyToolkit.hh"
/**
* \file
* Event shape variables.
* \author bjjung
*/
namespace JSIRENE {}
namespace JPP { using namespace JSIRENE; }
namespace JSIRENE {
using JMATH::JMatrix3S;
using JGEOMETRY3D::JPosition3D;
/**
* Compute thrust for a given range of tracks.\n
* The definition was taken from the description in section 15.2.2 of arXiv:hep-ph/0603175v2.
*
* \param __begin beginning of track data
* \param __end end of track data
* \return thrust
*/
inline double getThrust(const std::vector::const_iterator __begin,
const std::vector::const_iterator __end)
{
using namespace std;
using namespace JPP;
static const int MAX_STEPS = 20;
static const double epsilon = 1e-9;
Vec axis0(0.0, 0.0, 0.0); //!< Thrust axis
// Set thrust axis start value (= summed final state particle momenta)
for (vector::const_iterator i = __begin; i != __end; ++i) {
if (is_finalstate(*i)) {
axis0 += getKineticEnergy(*i) * i->dir;
}
}
if (axis0.len() < epsilon) { // If total final state momentum is zero, set start value equal to first final state track momentum
vector::const_iterator i = find_if(__begin, __end, is_finalstate);
if (i != __end) {
axis0 = getKineticEnergy(*i) * i->dir;
} else {
THROW(JValueOutOfRange, "getThrust(): No final state tracks given.");
}
}
axis0.normalize();
// Find axis which maximizes the thrust
double res = numeric_limits::max(); //!< residual angle
for (int k = 0; k < MAX_STEPS && res > epsilon; ++k) {
Vec axis1(0.0, 0.0, 0.0);
for (vector::const_iterator i = __begin; i != __end; ++i) {
if (is_finalstate(*i)) {
const Vec p = getKineticEnergy(*i) * i->dir;
axis1 += (p.dot(axis0) > 0.0 ? p : -p);
}
}
axis1.normalize();
res = acos( axis0.dot(axis1) );
axis0 = axis1;
}
// Compute thrust
double pSum = 0.0;
double norm = 0.0;
for (vector::const_iterator i = __begin; i != __end; ++i) {
if (is_finalstate(*i)) {
const Vec p = getKineticEnergy(*i) * i->dir;
pSum += fabs(p.dot(axis0));
norm += p.len();
}
}
return pSum / norm;
}
/**
* Compute thrust for a given event.\n
* The definition was taken from the description in section 15.2.2 of arXiv:hep-ph/0603175v2.
*
* \param event event
* \return thrust
*/
inline double getThrust(const Evt& event)
{
return getThrust(event.mc_trks.begin(), event.mc_trks.end());
}
/**
* Class for computing Fox-Wolfram moments.\n
* The Fox-Wolfram moments are calculated as defined in section 15.2.3 of arXiv:hep-ph/0603175v2.
*/
struct JFoxWolframMoments :
std::vector
{
typedef typename std::vector::reference reference;
/**
* Default constructor.
*/
JFoxWolframMoments() :
std::vector(0)
{}
/**
* Constructor.
*
* \param __begin beginning of track data
* \param __end end of track data
* \param can detector can
* \param N maximum order
*/
JFoxWolframMoments(const std::vector::const_iterator __begin,
const std::vector::const_iterator __end,
const JCylinder3D& can = getMaximumContainmentVolume(),
const int N = 4) :
std::vector(N+1)
{
using namespace std;
using namespace JPP;
const double Evis = getVisibleEnergy(__begin, __end, can);
if (Evis > 0.0) {
for (vector::const_iterator i = __begin; i != __end; ++i) {
if (is_finalstate(*i)) {
const Vec p0 = getKineticEnergy(*i) * i->dir;
for (vector::const_iterator j = __begin; j != __end; ++j) {
if (is_finalstate(*j)) {
const Vec p1 = getKineticEnergy(*j) * j->dir;
for (int k = 0; k < (int) this->size(); ++k) {
(*this)[k] += getContribution(p0, p1, k) / Evis / Evis;
}
}
}
}
}
}
}
/**
* Constructor.
*
* \param event event
* \param can can
* \param N maximum order
*/
JFoxWolframMoments(const Evt& event,
const JCylinder3D& can = getMaximumContainmentVolume(),
const int N = 4) :
JFoxWolframMoments(event.mc_trks.begin(),
event.mc_trks.end(),
can,
N)
{}
private:
/**
* Get contribution to Fox-Wolfram moment from two momenta.
*
* \param p0 first momentum vector
* \param p1 second momentum vector
* \param n Legendre polynomial index
* \return Fox-Wolfram moment contribution
*/
double getContribution(const Vec& p0,
const Vec& p1,
const size_t n)
{
using namespace JPP;
return p0.len() * p1.len() * legendre(n, p0.dot(p1) / p0.len() / p1.len());
}
};
/**
* Class for sphericity tensor calculations.\n
* The tensor is constructed following the definition in section 15.2.1 of arXiv:hep-ph/0603175v2.
*/
struct JSphericityTensor :
public JMatrix3S
{
typedef typename JMatrix3S::eigen_values eigen_values;
/**
* Default constructor.
*/
JSphericityTensor()
{}
/**
* Constructor.
*
* \param __begin beginning of track data
* \param __end end of track data
* \param r the power of the momentum-dependence
*/
JSphericityTensor(std::vector::const_iterator __begin,
std::vector::const_iterator __end,
const double r = 2.0)
{
using namespace std;
using namespace JPP;
double norm = 0.0;
for (vector::const_iterator i = __begin; i != __end; ++i) {
if (is_finalstate(*i)) {
const Vec p = getKineticEnergy(*i) * i->dir;
const double c = pow(p.len(), r-2);
a00 += c * p.x * p.x;
a01 += c * p.x * p.y;
a02 += c * p.x * p.z;
a10 += c * p.y * p.x;
a11 += c * p.y * p.y;
a12 += c * p.y * p.z;
a20 += c * p.z * p.x;
a21 += c * p.z * p.y;
a22 += c * p.z * p.z;
norm += c * p.dot(p);
}
}
this->div(norm);
lambda = JMatrix3S::getEigenValues();
}
/**
* Constructor.
*
* \param event event
* \param r the power of the momentum dependence
*/
JSphericityTensor(const Evt& event,
const double r = 2.0) :
JSphericityTensor(event.mc_trks.begin(), event.mc_trks.end(), r)
{}
/**
* Get eigenvalues.
*
* \return eigenvalues
*/
const eigen_values& getEigenValues() const
{
return lambda;
}
/**
* Get sphericity.
*
* \return sphericity
*/
double getSphericity() const
{
return 3 * (lambda[1] + lambda[2]) / 2.0;
}
/**
* Get aplanarity.
*
* \return aplanarity
*/
double getAplanarity() const
{
return 3 * lambda[2] / 2.0;
}
/**
* Get circularity.
*
* \return circularity
*/
double getCircularity() const
{
return (fabs(lambda[0] + lambda[1]) > 0.0 ?
2 * lambda[1] / (lambda[0] + lambda[1]) : 0.0);
}
/**
* Get planar flow.
*
* \return planar flow
*/
double getPlanarFlow() const
{
return 4 * lambda[0] * lambda[1] / (lambda[0] + lambda[1]) / (lambda[0] + lambda[1]);
}
/**
* Get C-variable.
*
* \return C
*/
double getC() const
{
return 3 * (lambda[0]*lambda[1] + lambda[0]*lambda[2] + lambda[1]*lambda[2]);
}
/**
* Get D-variable.
*
* \return D
*/
double getD() const
{
return 27 * lambda[0] * lambda[1] * lambda[2];
}
private:
eigen_values lambda; //!< sorted eigenvalues
};
/**
* Class for hit inertia tensor calculations.\n
* The given methods are inspired by section 5.2.2 of Stephanie Hickford's thesis.
*/
struct JHitInertiaTensor :
public JMatrix3S
{
typedef typename JMatrix3S::eigen_values eigen_values;
/**
* Default constructor.
*/
JHitInertiaTensor()
{}
/**
* Constructor.
*
* \param __begin beginning of hit data
* \param __end end of hit data
* \param reference reference position (e.g. the location of the interaction vertex)
*/
JHitInertiaTensor(std::vector::const_iterator __begin,
std::vector::const_iterator __end,
const JPosition3D& reference)
{
using namespace std;
using namespace JPP;
for (vector::const_iterator hit = __begin; hit != __end; ++hit) {
const JPosition3D D = sqrt(hit->a) * (getPosition(hit->pos) - reference);
a00 += D.getLengthSquared() - D.getX() * D.getX();
a01 += 0.0 - D.getX() * D.getY();
a02 += 0.0 - D.getX() * D.getZ();
a10 += 0.0 - D.getY() * D.getX();
a11 += D.getLengthSquared() - D.getY() * D.getY();
a12 += 0.0 - D.getY() * D.getZ();
a20 += 0.0 - D.getZ() * D.getX();
a21 += 0.0 - D.getZ() * D.getY();
a22 += D.getLengthSquared() - D.getZ() * D.getZ();
}
lambda = getEigenValues();
}
/**
* Constructor.
*
* \param event event
* \param reference reference position (e.g. the location of the interaction vertex)
*/
JHitInertiaTensor(const Evt& event,
const JPosition3D& reference) :
JHitInertiaTensor(event.hits.begin(), event.hits.end(), reference)
{}
/**
* Get eigenvalue ratio.
*
* \return eigenvalue ratio
*/
double getEigenvalueRatio() const
{
return lambda[2] / (lambda[0] + lambda[1] + lambda[2]);
}
private:
eigen_values lambda; //!< eigenvalues
};
}
#endif