JAAnetToolkit.hh

Go to the documentation of this file.

#ifndef __JAANET__JAANETTOOLKIT__
#define __JAANET__JAANETTOOLKIT__

#include <cstdlib>
#include <algorithm>

#include "km3net-dataformat/offline/Hit.hh"
#include "km3net-dataformat/offline/Vec.hh"
#include "km3net-dataformat/offline/Trk.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/definitions/trkmembers.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
#include "km3net-dataformat/definitions/reconstruction.hh"
#include "km3net-dataformat/tools/reconstruction.hh"

#include "JLang/JException.hh"
#include "JLang/JPredicate.hh"

#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JTrack3E.hh"
#include "JGeometry3D/JVertex3D.hh"
#include "JGeometry3D/JTransformation3D.hh"

#include "JPhysics/JConstants.hh"
#include "JPhysics/JTimeRange.hh"
#include "JPhysics/JGeane.hh"

#include "JAAnet/JParticleTypes.hh"
#include "JAAnet/JPDB.hh"


/**
 * \file
 *
 * Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
 * \author mdejong
 */
namespace JAANET {}
namespace JPP { using namespace JAANET; }

/**
 * Extensions to Evt data format.
 */
namespace JAANET {

  using JLANG::make_predicate;
  using JLANG::JIndexOutOfRange;
  using JLANG::JValueOutOfRange;
  
  using JGEOMETRY3D::JPosition3D;
  using JGEOMETRY3D::JDirection3D;
  using JGEOMETRY3D::JAxis3D;
  using JGEOMETRY3D::JTrack3D;
  using JGEOMETRY3D::JTrack3E;
  using JGEOMETRY3D::JVertex3D;
  using JGEOMETRY3D::JTransformation3D;
  
  using JPHYSICS::JTimeRange;
  using JPHYSICS::MASS_PROTON;
  using JPHYSICS::MASS_NEUTRON;
  
  /**
   * AAnet shower fit reconstruction type.
   */
  static const int AASHOWER_RECONSTRUCTION_TYPE  = 101;

  /**
   * Enumeration of hit types based on km3 codes.
   */
  enum JHitType_t {
    HIT_TYPE_MUON_DIRECT              =   +5,    //!< Direct    light from muon
    HIT_TYPE_MUON_SCATTERED           =   -5,    //!< Scattered light from muon
    HIT_TYPE_DELTARAYS_DIRECT         =   +4,    //!< Direct    light from delta-rays
    HIT_TYPE_DELTARAYS_SCATTERED      =   -4,    //!< Scattered light from delta-rays
    HIT_TYPE_BREMSSTRAHLUNG_DIRECT    =   +3,    //!< Direct    light from Bremsstrahlung
    HIT_TYPE_BREMSSTRAHLUNG_SCATTERED =   -3,    //!< Scattered light from Bremsstrahlung
    HIT_TYPE_SHOWER_DIRECT            =  +99,    //!< Direct    light from primary shower
    HIT_TYPE_SHOWER_SCATTERED         =  -99,    //!< Scattered light from primary shower
    HIT_TYPE_NOISE                    =   -1,    //!< Random noise
    HIT_TYPE_UNKNOWN                  =    0     //!< Unknown source
  };


  /**
   * Get true time of hit.
   *
   * \param  hit              hit
   * \return                  true time of hit
   */
  inline double getTime(const Hit& hit)
  {
    return hit.t;
  }


  /**
   * Get true charge of hit.
   *
   * \param  hit              hit
   * \return                  true charge of hit
   */
  inline double getNPE(const Hit& hit)
  {
    return hit.a;
  }


  /**
   * Verify hit origin.
   *
   * \param  hit              hit
   * \return                  true if noise; else false
   */
  inline bool is_noise(const Hit& hit)
  {
    return hit.type == HIT_TYPE_NOISE;
  }
  

  /**
   * Get time range (i.e.\ time between earliest and latest hit) of Monte Carlo event.\n
   * Note that the global event time in not taken into account. 
   *
   * \param  event            event
   * \return                  time range
   */
  inline JTimeRange getTimeRange(const Evt& event)
  {
    JTimeRange time_range(JTimeRange::DEFAULT_RANGE());

    for (std::vector<Hit>::const_iterator hit = event.mc_hits.begin(); hit != event.mc_hits.end(); ++hit) {
      if (!is_noise(*hit)) {
        time_range.include(getTime(*hit));
      }
    }

    return time_range;
  }
  

  /**
   * Get time range (i.e.\ time between earliest and latest hit) of Monte Carlo event.\n
   * The resulting time range is larger than or equal to the given time window.\n
   * Note that the global event time in not taken into account. 
   *
   * \param  event            event
   * \param  T_ns             time window [ns]
   * \return                  time range
   */
  inline JTimeRange getTimeRange(const Evt& event, const JTimeRange& T_ns)
  {
    JTimeRange time_range = getTimeRange(event);

    if (!time_range.is_valid()) {
      time_range = T_ns;
    } 

    if (time_range.getLength() < T_ns.getLength()) {
      time_range.setLength(T_ns.getLength());
    }

    return time_range;
  }


  /**
   * Get position.
   *
   * \param  pos            position
   * \return                position
   */
  inline JPosition3D getPosition(const Vec& pos)
  {
    return JPosition3D(pos.x, pos.y, pos.z);
  }


  /**
   * Get position.
   *
   * \param  pos            position
   * \return                position
   */
  inline Vec getPosition(const JPosition3D& pos)
  {
    return Vec(pos.getX(), pos.getY(), pos.getZ());
  }


  /**
   * Get position.
   *
   * \param  track          track
   * \return                position
   */
  inline JPosition3D getPosition(const Trk& track)
  {
    return getPosition(track.pos);
  }


  /**
   * Get direction.
   *
   * \param  dir            direction
   * \return                direction
   */
  inline JDirection3D getDirection(const Vec& dir)
  {
    return JDirection3D(dir.x, dir.y, dir.z);
  }


  /**
   * Get direction.
   *
   * \param  dir            direction
   * \return                direction
   */
  inline Vec getDirection(const JDirection3D& dir)
  {
    return Vec(dir.getDX(), dir.getDY(), dir.getDZ());
  }


  /**
   * Get direction.
   *
   * \param  track          track
   * \return                direction
   */
  inline JDirection3D getDirection(const Trk& track)
  {
    return getDirection(track.dir);
  }


  /**
   * Get axis.
   *
   * \param  track          track
   * \return                axis
   */
  inline JAxis3D getAxis(const Trk& track)
  {
    return JAxis3D(getPosition(track), getDirection(track));
  }


  /**
   * Get track.
   *
   * \param  track          track
   * \return                track
   */
  inline JTrack3E getTrack(const Trk& track)
  {
    return JTrack3E(JTrack3D(getAxis(track), track.t), track.E);
  }


  /**
   * Get transformation.
   *
   * \param  track          track
   * \return                transformation
   */
  inline JTransformation3D getTransformation(const Trk& track)
  {
    return JTransformation3D(getPosition(track), getDirection(track));
  }


  /**
   * Get track information.
   *
   * \param  track          track
   * \param  index          index
   * \param  value          default value
   * \return                actual  value
   */
  inline double getW(const Trk& track, const size_t index, const double value)
  {
    if (index < track.fitinf.size())
      return track.fitinf[index];
    else
      return value;
  }

  
  /**
   * Test whether given track is a photon or neutral pion.
   *
   * \param  track          track
   * \return                true if photon or neutral pion; else false
   */
  inline bool is_photon  (const Trk& track) { return (    track.type  ==  TRACK_TYPE_PHOTON || 
                                                      abs(track.type) ==  TRACK_TYPE_NEUTRAL_PION); }

  /**
   * Test whether given track is a neutrino.
   *
   * \param  track          track
   * \return                true if neutrino; else false
   */
  inline bool is_neutrino(const Trk& track) { return (abs(track.type) == TRACK_TYPE_NUE  ||
                                                      abs(track.type) == TRACK_TYPE_NUMU ||
                                                      abs(track.type) == TRACK_TYPE_NUTAU); }

  /**
   * Test whether given track is a (anti-)electron.
   *
   * \param  track          track
   * \return                true if (anti-)electron; else false
   */
  inline bool is_electron(const Trk& track) { return (abs(track.type) == TRACK_TYPE_ELECTRON); }

  /**
   * Test whether given track is a (anti-)muon.
   *
   * \param  track          track
   * \return                true if (anti-)muon; else false
   */
  inline bool is_muon    (const Trk& track) { return (abs(track.type) == TRACK_TYPE_MUON); }

  /**
   * Test whether given track is a (anti-)tau.
   *
   * \param  track          track
   * \return                true if (anti-)tau; else false
   */
  inline bool is_tau     (const Trk& track) { return (abs(track.type) == TRACK_TYPE_TAU); }

  /**
   * Test whether given track is a charged pion.
   *
   * \param  track          track
   * \return                true if charged pion; else false
   */
  inline bool is_pion    (const Trk& track) { return (abs(track.type) == TRACK_TYPE_CHARGED_PION_PLUS); }

  /**
   * Test whether given track is a (anti-)proton.
   *
   * \param  track          track
   * \return                true if (anti-)proton; else false
   */
  inline bool is_proton  (const Trk& track) { return (abs(track.type) == TRACK_TYPE_PROTON); }

  /**
   * Test whether given track is a (anti-)neutron.
   *
   * \param  track          track
   * \return                true if (anti-)neutron; else false
   */
  inline bool is_neutron (const Trk& track) { return (abs(track.type) == TRACK_TYPE_NEUTRON); }

  /**
   * Test whether given track is a lepton
   *
   * \param track           track
   * \return                true if lepton; else fails
   */
  inline bool is_lepton  (const Trk& track) { return (is_neutrino(track) ||
                                                      is_electron(track) ||
                                                      is_muon    (track) ||
                                                      is_tau     (track)); }

  /**
   * Test whether given track is a charged lepton
   *
   * \param track           track
   * \return                true if lepton; else fails
   */
  inline bool is_charged_lepton  (const Trk& track) { return (is_electron(track) ||
                                                              is_muon    (track) ||
                                                              is_tau     (track)); }
  
  /**
   * Test whether given track is a hadron
   *
   * \param track           track
   * \return                true if hadron; else fails
   */
  inline bool is_hadron  (const Trk& track) { return (abs(track.type) != TRACK_TYPE_PHOTON &&
                                                      !is_lepton(track)); }
  
  /**
   * Test whether given track is a meson (is hadron and third digit of PDG code is zero)
   *
   * \param track           track
   * \return                true if hadron; else fails
   */
  inline bool is_meson   (const Trk& track) { return (is_hadron(track) &&
                                                      ((int)(track.type / 1000)) % 10 == 0); }
  
  /**
   * Test whether given track is a baryon (is hadron and third digit of PDG code is not zero)
   *
   * \param track           track
   * \return                true if hadron; else fails
   */
  inline bool is_baryon  (const Trk& track) { return (is_hadron(track) &&
                                                      ((int)(track.type / 1000)) % 10 != 0); }
  
  /**
   * Test whether given track corresponds to given particle type.
   *
   * \param  track          track
   * \param  type           particle type
   * \return                true if matches the corresponding particle; else false
   */
  inline bool has_particleID (const Trk& track, const int type) { return (track.type == type); }
  
  /** 
   * Test whether given track corresponds to an initial state particle.
   *
   * \param  track          track
   * \return                true if track corresponds to an initial state particle; else false
   */
  inline bool is_initialstate(const Trk& track)
  {
    if        (Evt::ROOT_IO_VERSION >= 14) {
      
      return (track.status == TRK_ST_PRIMARYNEUTRINO ||
              track.status == TRK_ST_PRIMARYCOSMIC   ||
              track.status == TRK_ST_MUONBUNDLE      ||
              track.status == TRK_ST_ININUCLEI);
      
    } else if (Evt::ROOT_IO_VERSION >= 11) {

      return (track.mother_id == TRK_MOTHER_UNDEFINED ||
              track.mother_id == TRK_MOTHER_NONE);

    } else {
      
      return is_neutrino(track) && track.id == 0;
    }
  }
  
  /** 
   * Test whether given track corresponds to a final state particle.
   *
   * \param  track          track
   * \return                true if track corresponds to a final state particle; else false
   */
  inline bool is_finalstate(const Trk& track)
  {
    if (Evt::ROOT_IO_VERSION >= 16) {
      return  track.status == TRK_ST_FINALSTATE;
    } else {
      return !is_initialstate(track);
    }
  }
  
  /**
   * Test whether given event has an incoming neutrino.
   *
   * \param  evt            event
   * \return                true if neutrino; else false
   */
  inline bool has_neutrino(const Evt& evt)
  {
    if (Evt::ROOT_IO_VERSION >= 14) {

      std::vector<Trk>::const_iterator i = find_if(evt.mc_trks.cbegin(), evt.mc_trks.cend(),
                                                   make_predicate(&Trk::status, TRK_ST_PRIMARYNEUTRINO));
      return i != evt.mc_trks.cend();
      
    } else {
      
      return !evt.mc_trks.empty() && is_neutrino(evt.mc_trks[0]);
    }
  }
  
  /**
   * Get incoming neutrino.
   *
   * \param  evt            event
   * \return                neutrino
   */
  inline const Trk& get_neutrino(const Evt& evt)
  {
    if (Evt::ROOT_IO_VERSION >= 14) {
      
      std::vector<Trk>::const_iterator i = find_if(evt.mc_trks.cbegin(), evt.mc_trks.cend(),
                                                   make_predicate(&Trk::status, TRK_ST_PRIMARYNEUTRINO));

      if (i != evt.mc_trks.cend()) { return *i; }
      
    } else if (!evt.mc_trks.empty() && is_neutrino(evt.mc_trks[0])) {
      
      return evt.mc_trks[0];
    }
      
    THROW(JIndexOutOfRange, "JAANET::get_neutrino(): Event " << evt.id << " has no neutrino."); 
  }

  /**
   * Get primary.
   *
   * \param  evt            event
   * \return                primary track
   */
  inline const Trk& get_primary(const Evt& evt)
  {
    using namespace std;
    using namespace JPP;

    for (vector<Trk>::const_iterator i = evt.mc_trks.cbegin(); i != evt.mc_trks.cend(); ++i) {
      if (is_initialstate(*i) && i->status != TRK_ST_ININUCLEI) {
        return *i;
      }
    }

    THROW(JIndexOutOfRange, "JAANET::get_primary(): Cannot retrieve primary track for event " << evt.id << ".");
  }


  /**
   * Get vertex.
   *
   * \param  track          track
   * \return                vertex
   */
  inline JVertex3D getVertex(const Trk& track)
  {
    return JVertex3D(getPosition(track), track.t);
  }

  /**
   * Get event vertex.
   *
   * \param  event          event
   * \return                event vertex
   */
  inline JVertex3D getVertex(const Evt& event)
  {
    using namespace std;
    using namespace JPP;
    
    if (has_neutrino(event)) {

      const Trk& neutrino = get_neutrino(event);

      return getVertex(neutrino);
      
    } else if (!event.mc_trks.empty()) {

      vector<Trk>::const_iterator i = find_if(event.mc_trks.begin(), event.mc_trks.end(), &is_initialstate);

      if (i != event.mc_trks.end()) {
        return getVertex(*i);
      } else {
        THROW(JValueOutOfRange, "getVertex(): No initial state track found.");
      }
      
    } else if (!event.trks.empty()) { // For reconstructed DAQ events
      
      const Trk& track = get_best_reconstructed_track<JPP_RECONSTRUCTION_TYPE>(event);

      return getVertex(track);
      
    } else {

      THROW(JValueOutOfRange, "getVertex(): Could not find a valid event vertex.");
    }
  }

  /**
   * Count the number of electrons in a given event
   *
   * \param  evt            event
   * \return                number of electrons
   */
  inline int count_electrons(const Evt& evt)
  {
    return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_electron);
  }

  /**
   * Test whether given event has an electron.
   *
   * \param  evt            event
   * \return                true if event has electron; else false
   */
  inline bool has_electron(const Evt& evt)
  {
    return count_electrons(evt) != 0;
  }

  /**
   * Get first electron from the event tracklist
   *
   * \param  evt            event
   * \return                electron
   */
  inline const Trk& get_electron(const Evt& evt)
  {
    if (count_electrons(evt) > 0)
      return *std::find_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_electron);
    else
      THROW(JIndexOutOfRange, "This event has no electron.");
  }

  /**
   * Test whether given hit was produced by an electron
   *
   * Warning: This method will only work with the output of KM3Sim.
   * For JSirene or KM3, you should check if <tt>track.id == hit.origin</tt> instead.
   *
   * \param  hit            hit
   * \return                true if hit was produced by an electron; else false
   */
  inline bool from_electron(const Hit& hit)
  {
    return hit.type == GEANT4_TYPE_ELECTRON || hit.type == GEANT4_TYPE_ANTIELECTRON;
  }

  /**
   * Count the number of muons in a given event
   *
   * \param  evt            event
   * \return                number of muons
   */
  inline int count_muons(const Evt& evt)
  {
    return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_muon);
  }

  /**
   * Test whether given event has a muon.
   *
   * \param  evt            event
   * \return                true if event has muon; else false
   */
  inline bool has_muon(const Evt& evt)
  {
    return count_muons(evt) != 0;
  }

  /**
   * Get first muon from the event tracklist
   *
   * \param  evt            event
   * \return                muon
   */
  inline const Trk& get_muon(const Evt& evt)
  {
    if (count_muons(evt) > 0)
      return *std::find_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_muon);
    else
      THROW(JIndexOutOfRange, "This event has no muon.");
  }

  /**
   * Test whether given hit was produced by a muon
   *
   * Warning: This method will only work with the output of KM3Sim.
   * For JSirene or KM3, you should check if <tt>track.id == hit.origin</tt> instead.
   *
   * \param  hit            hit
   * \return                true if hit was produced by a muon; else false
   */
  inline bool from_muon(const Hit& hit)
  {
    return hit.type == GEANT4_TYPE_MUON || hit.type == GEANT4_TYPE_ANTIMUON;
  }

  /**
   * Count the number of taus in a given event
   *
   * \param  evt            event
   * \return                number of taus
   */
  inline int count_taus(const Evt& evt)
  {
    return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_tau);
  }

  /**
   * Test whether given event has a tau.
   *
   * \param  evt            event
   * \return                true if event has tau; else false
   */
  inline bool has_tau(const Evt& evt)
  {
    return count_taus(evt) != 0;
  }

  /**
   * Get first tau from the event tracklist
   *
   * \param  evt            event
   * \return                tau
   */
  inline const Trk& get_tau(const Evt& evt)
  {
    if (count_taus(evt) > 0)
      return *std::find_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_tau);
    else
      THROW(JIndexOutOfRange, "This event has no tau.");
  }

  /**
   * Test whether given hit was produced by a tau
   *
   * Warning: This method will only work with the output of KM3Sim.
   * For JSirene or KM3, you should check if <tt>track.id == hit.origin</tt> instead.
   *
   * \param  hit            hit
   * \return                true if hit was produced by a tau; else false
   */
  inline bool from_tau(const Hit& hit)
  {
    return hit.type == GEANT4_TYPE_TAU || hit.type == GEANT4_TYPE_ANTITAU;
  }

  /**
   * Count the number of hadrons in a given event (not including neutral pions)
   *
   * \param  evt            event
   * \return                number of hadrons
   */
  inline int count_hadrons(const Evt& evt)
  {
    return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_hadron);
  }
  
  /**
   * Test whether given hit was produced by a hadronic shower
   *
   * Warning: This method will only work with the output of KM3Sim.
   * For JSirene or KM3, you should check if <tt>track.id == hit.origin</tt> instead.
   *
   * \param  hit            hit
   * \return                true if hit was produced by a hadron; else false
   */
  inline bool from_hadron(const Hit& hit)
  {
    return (!from_electron(hit) && !from_muon(hit) && !from_tau(hit));
  }

  
  /**
   * Add time offset to mc event;
   * according to current definition, the absolute time of the event is defined by the track "t" attribute;
   * this could change in the future if the global attribute mc_t is assigned to this purpose.
   *
   * \param  evt            event
   * \param  tOff           time offset [ns]
   */
  inline void add_time_offset(Evt& evt, const double tOff)
  {
    for (std::vector<Trk>::iterator p = evt.mc_trks.begin(); p != evt.mc_trks.end(); p++) {
      p->t += tOff;
    }
  }


  /**
   * Get kinetic energy of particle with given mass.
   *
   * \param  E         energy [GeV]
   * \param  m         mass   [GeV]
   * \return           energy [GeV]
   */
  inline double getKineticEnergy(const double E, const double m)
  {
    if (E > m) {
      return sqrt((E - m) * (E + m));
    } else {
      return 0.0;
    }
  }
  
  
  /**
   * Get track kinetic energy.
   *
   * \param  trk          track
   * \return              kinetic energy [GeV]
   */
  inline double getKineticEnergy(const Trk& trk)
  {
    const double energy = trk.E;
    const double mass   = JPDB::getInstance().getPDG(trk.type).mass;

    return getKineticEnergy(energy, mass);
  }
  
  
  /**
   * Get initial state energy of a neutrino interaction.\n
   * This method includes baryon number conservation.
   *
   * \param  evt         event
   * \return             energy [GeV]
   */
  inline double getE0(const Evt& evt)
  {
    using namespace std;
    
    const Trk& neutrino = get_neutrino(evt);
    
    double E0 = neutrino.E;
    
    for (vector<Trk>::const_iterator track = evt.mc_trks.cbegin(); track != evt.mc_trks.end(); ++track) {
      
      if (!is_finalstate(*track)) { continue; }
      
      if        (track->type == TRACK_TYPE_NEUTRON     ||
                 track->type == TRACK_TYPE_SIGMA_MINUS ||
                 track->type == TRACK_TYPE_NEUTRAL_SIGMA) {
          
        E0 += MASS_NEUTRON;
          
      } else if (track->type == TRACK_TYPE_PROTON     ||
                 track->type == TRACK_TYPE_LAMBDA     ||
                 track->type == TRACK_TYPE_SIGMA_PLUS) {
        
        E0 += MASS_PROTON;
        
      } else if (track->type == TRACK_TYPE_ANTINEUTRON) {
        
        E0 -= MASS_NEUTRON;
        
      } else if (track->type == TRACK_TYPE_ANTIPROTON ||
                 track->type == TRACK_TYPE_ANTILAMBDA) {
          
        E0 -= MASS_PROTON;
      }
    }

    return E0;
  }

  
  /**
   * Get final state energy of a neutrino interaction.\n
   * This method includes muon energy loss.
   *
   * \param  evt         event
   * \return             energy [GeV]
   */
  inline double getE1(const Evt& evt)
  {
    using namespace std;
    using namespace JPP;
    
    double E1 = 0.0;
    
    const Trk& neutrino = get_neutrino(evt);

    for (vector<Trk>::const_iterator track = evt.mc_trks.cbegin(); track != evt.mc_trks.cend(); ++track) {

      if (!is_finalstate(*track)) { continue; }      
        
      if (is_muon(*track)) {

        const Trk& mother = ( track->mother_id > TRK_MOTHER_UNDEFINED ?
                              evt.mc_trks[track->mother_id] : neutrino );

        const double distance = (track->pos - mother.pos).len();

        E1 += gWater(track->E, -distance);
          
      } else {

        E1 += track->E;
      }
    }

    return E1;
  }

  
  /**
   * Get momentum vector of the initial state of a neutrino interaction.\n
   * This method assumes neutrino DIS on a stationary nucleus
   *
   * \param  evt         event
   * \return             energy [GeV]
   */
  inline Vec getP0(const Evt& evt)
  {
    const Trk& neutrino = get_neutrino(evt);

    return neutrino.E * neutrino.dir;
  }


  /**
   * Get momentum vector of the final state of a neutrino interaction.\n
   * This method includes muon energy losses.
   *
   * \param  evt         event
   * \return             final state momentum vector [GeV]
   */
  inline Vec getP1(const Evt& evt)
  {
    using namespace std;
    using namespace JPP;
    
    Vec P1(0,0,0);

    const Trk& neutrino = get_neutrino(evt);

    for (vector<Trk>::const_iterator track = evt.mc_trks.cbegin(); track != evt.mc_trks.cend(); ++track) {

      if (!is_finalstate(*track)) { continue; }
      
      double kineticEnergy = 0.0;

      if (is_muon(*track)) {

        const Trk& mother = ( track->mother_id > TRK_MOTHER_UNDEFINED ?
                              evt.mc_trks[track->mother_id] : neutrino );

        const double distance = (track->pos - mother.pos).len();
        const double energy   = gWater(track->E, -distance);

        kineticEnergy = getKineticEnergy(energy, MASS_MUON);
          
      } else {

        kineticEnergy = getKineticEnergy(*track);
      }

      P1 += kineticEnergy * track->dir;
    }
    
    return P1;
  }


  /**
   * Get final state invariant mass.
   *
   * \param  event       event
   * \return             invariant mass [GeV]
   */
  inline double getInvariantMass(const Evt& event)
  {
    using namespace std;
    using namespace JPP;

    double Minv = 0.0;

    for (vector<Trk>::const_iterator track = event.mc_trks.cbegin(); track != event.mc_trks.cend(); ++track) {

      if (!is_finalstate(*track)) { continue; }

      const double mass = JPDB::getInstance().getPDG(track->type).mass;

      Minv += mass;
    }

    return Minv;
  }
}

#endif