JEventShapeVariables.hh

Go to the documentation of this file.

#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 "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;
  using JGEOMETRY3D::getMaximumContainmentVolume;


  /**
   * Compute thrust for a given range of tracks.\n
   * The definition was taken from the description in section 15.2.2 of <a href="https://arxiv.org/pdf/hep-ph/0603175.pdf%E2%80%8B">arXiv:hep-ph/0603175v2</a>.
   *
   * \param  __begin           beginning of track data
   * \param  __end             end of track data
   * \return                   thrust
   */
  inline double getThrust(const std::vector<Trk>::const_iterator __begin,
                          const std::vector<Trk>::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<Trk>::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<Trk>::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<double>::max(); //!< residual angle
    
    for (int k = 0; k < MAX_STEPS && res > epsilon; ++k) {

      Vec axis1(0.0, 0.0, 0.0);
      
      for (vector<Trk>::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<Trk>::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 <a href="https://arxiv.org/pdf/hep-ph/0603175.pdf%E2%80%8B">arXiv:hep-ph/0603175v2</a>.
   *
   * \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 <a href="https://arxiv.org/pdf/hep-ph/0603175.pdf%E2%80%8B">arXiv:hep-ph/0603175v2</a>.
   */
  struct JFoxWolframMoments :
    std::vector<double>
  {
    typedef typename std::vector<double>::reference     reference;
    
    /**
     * Default constructor.
     */
    JFoxWolframMoments() :
      std::vector<double>(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<Trk>::const_iterator __begin,
                       const std::vector<Trk>::const_iterator __end,
                       const JCylinder3D&                     can = getMaximumContainmentVolume(),
                       const int                              N   = 4) :
      std::vector<double>(N+1)
    {
      using namespace std;
      using namespace JPP;

      const double Evis = getVisibleEnergy(__begin, __end, can);

      if (Evis > 0.0) {

        for (vector<Trk>::const_iterator i = __begin; i != __end; ++i) {

          if (is_finalstate(*i)) {

            const Vec p0 = getKineticEnergy(*i) * i->dir;
        
            for (vector<Trk>::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 <a href="https://arxiv.org/pdf/hep-ph/0603175.pdf%E2%80%8B">arXiv:hep-ph/0603175v2</a>.
   */
  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<Trk>::const_iterator __begin,
                      std::vector<Trk>::const_iterator __end,
                      const double                     r = 2.0)
    {
      using namespace std;
      using namespace JPP;
      
      double norm = 0.0;

      for (vector<Trk>::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 <a href="https://s3.cern.ch/inspire-prod-files-e/ed2c2d7867bc6d49dd7fd414af45f3d2">Stephanie Hickford's thesis</a>.
   */
  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<Hit>::const_iterator __begin,
                      std::vector<Hit>::const_iterator __end,
                      const JPosition3D&               reference)
    {
      using namespace std;
      using namespace JPP;

      for (vector<Hit>::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