JAAnetToolkit.hh

Go to the documentation of this file.

#ifndef __JAANET__JAANETTOOLKIT__
#define __JAANET__JAANETTOOLKIT__

#include <stdlib.h>
#include <limits>
#include <algorithm>

#include "evt/Hit.hh"
#include "evt/Vec.hh"
#include "evt/Trk.hh"
#include "evt/Evt.hh"

#include "JLang/JException.hh"
#include "JTools/JRange.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 "JTrigger/JHitL0.hh"
#include "JDetector/JTimeRange.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTRouter.hh"
#include "JFit/JFitApplications.hh"
#include "JAAnet/JParticleTypes.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 AAnet data format.
 */
namespace JAANET {

  using JLANG::JIndexOutOfRange;
  using JGEOMETRY3D::JPosition3D;
  using JGEOMETRY3D::JDirection3D;
  using JGEOMETRY3D::JAxis3D;
  using JGEOMETRY3D::JTrack3D;
  using JGEOMETRY3D::JTrack3E;
  using JGEOMETRY3D::JVertex3D;
  using JGEOMETRY3D::JTransformation3D;
  using JTRIGGER::JHitL0;
  using JDETECTOR::JTimeRange;
  using JDETECTOR::JModuleRouter;
  using JDETECTOR::JPMTRouter;
  using JFIT::JPP_RECONSTRUCTION_TYPE;
  using JTOOLS::JRange;

  /**
   * AAshower reconstruction type for AAnet.
   */
  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)
  {
    // if IGAIN = 0 in km3 program, then pure_xx values are not set.
    
    if (hit.pure_a == 0.0)
      return hit.t;
    else
      return hit.pure_t;
  }

  
  /**
   * Get true charge of hit.
   *
   * \param  hit              hit
   * \return                  true charge of hit
   */
  inline double getNPE(const Hit& hit)
  {
    // if IGAIN = 0 in km3 program, then pure_xx values are not set.
    
    if (hit.pure_a == 0.0)
      return hit.a;
    else
      return hit.pure_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. earlist and latest hit time) of Monte Carlo event.
   * Note that the global event time in not taken into account. 
   *
   * \param  event            event
   * \return                  time range
   */
  inline JTimeRange getTimeRange(const Evt& event)
  {
    JTimeRange timeRange(JTimeRange::DEFAULT_RANGE);

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

    return timeRange;
  }
  

  /**
   * Get time range (i.e. earlist and latest hit time) of Monte Carlo event.
   * The given time window is offset so that the lower limit is equal to the earliest hit time.
   * The time window is then applied to all hits.
   * 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)
  {
    using namespace std;

    JTimeRange timeRange(JTimeRange::DEFAULT_RANGE);

    double t0 = numeric_limits<double>::max();

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

    if (t0 != numeric_limits<double>::max()) {

      JTimeRange window(T_ns);
    
      window.fixLowerLimit(t0);
    
      for (std::vector<Hit>::const_iterator hit = event.mc_hits.begin(); hit != event.mc_hits.end(); ++hit) {
        if (!is_noise(*hit) && window(getTime(*hit))) {
          timeRange.include(getTime(*hit));
        }
      }
    }

    return timeRange;
  }


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


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


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

  /**
   * 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 vertex.
   *
   * \param  track          track
   * \return                vertex
   */
  inline JVertex3D getVertex(const Trk& track)
  {
    return JVertex3D(getPosition(track), track.t);
  }


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


  /**
   * Get transformation.
   *
   * \param  hit            hit
   * \return                hit
   */
  inline JHitL0 getHit(const Hit& hit)
  {
    using namespace JPP;

    return JHitL0(JDAQPMTIdentifier(hit.dom_id, hit.channel_id), JAxis3D(getPosition(hit.pos), getDirection(hit.dir)), JHit(hit.t, hit.tot) );
  }


  /**
   * Get transformation.
   *
   * \param  hit            hit
   * \param  router         module router
   * \return                hit
   */
  inline JHitL0 getHit(const Hit& hit, const JModuleRouter& router)
  {
    using namespace JPP;

    const JPMT& pmt = router.getModule(hit.dom_id).getPMT(hit.channel_id);

    return JHitL0(JDAQPMTIdentifier(hit.dom_id, hit.channel_id), pmt.getAxis(), JHit(getTime(hit.tdc, pmt.getCalibration()), getToT(hit.tot, pmt.getCalibration())));
  }


  /**
   * Get transformation.
   *
   * \param  hit            hit
   * \param  router         pmt router
   * \return                hit
   */
  inline JHitL0 getHit(const Hit& hit, const JPMTRouter& router)
  {
    using namespace JPP;

    const JPMT& pmt = router.getPMT(hit.pmt_id);

    return JHitL0(JDAQPMTIdentifier(hit.dom_id, hit.channel_id), pmt.getAxis(), JHit(getTime(hit.tdc, pmt.getCalibration()), getToT(hit.tot, pmt.getCalibration())));
  }


  /**
   * 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 hadron
   *
   * \param track           track
   * \return                true if hadron; else fails
   */
  inline bool is_hadron  (const Trk& track) { return (!is_neutrino(track) &&
                                                      !is_electron(track) &&
                                                      !is_muon    (track) &&
                                                      !is_tau     (track)); }

  /**
   * 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 event has an incoming neutrino.
   *
   * \param  evt            event
   * \return                true if neutrino; else false
   */
  inline bool has_neutrino(const Evt& evt)
  { 
    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 (has_neutrino(evt))
      return evt.mc_trks[0];
    else
      THROW(JIndexOutOfRange, "This event has no neutrino.");
  }

  /**
   * 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 track.id == hit.origin 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 track.id == hit.origin 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));
  }

  /**
   * Auxiliary function to get average direction and total energy of a hadronic cascade
   *
   * \param  evt            event
   * \return                cascade (default if there are no hadrons)
   */
  inline Trk get_hadronic_cascade(const Evt& evt)
  {
    Trk cascade;
    
    for (vector<Trk>::const_iterator track = evt.mc_trks.begin(); track != evt.mc_trks.end(); ++track) {
      if (is_hadron(*track)) {
        cascade.E   += track->E;
        cascade.dir += track->dir * track->E;
        cascade.t    = track->t;
        cascade.pos  = track->pos;
      }
    }

    cascade.dir.normalize();

    return cascade;
  }
  
  /**
   * 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 (vector<Trk>::iterator p = evt.mc_trks.begin(); p != evt.mc_trks.end(); p++) {
      p->t += tOff;
    }
  }

  /**
   * Reconstruction type dependent comparison of track quality.
   *
   * The specialisation of this class should implement the function object operator
   * <pre>
   *    inline bool operator()(const Trk& first, const Trk& second) const
   * </pre>
   * and return true if the first track is better then the second track.
   */
  template<int reconstruction_type>
  struct quality_sorter {
    /**
     * The default comparison is based on:
     *    -# number of reconstruction stages, i.e. <tt>Trk::rec_stages.size()</tt> (the larger the better);
     *    -# likelihood of the final track, i.e. <tt>Trk::lik</tt> (the larger the better).
     *
     * \param  first          first  track
     * \param  second         second track
     * \return                true if first track has better quality than second; else false
     */
    inline bool operator()(const Trk& first, const Trk& second) const
    {
      if (first.rec_stages.size() == second.rec_stages.size())
        return first.lik > second.lik;
      else
        return first.rec_stages.size() > second.rec_stages.size();
    }
  };

  /**
   * Auxiliary class to test whether given track has specified history.\n
   * The history of a track consists of a reconstruction type and a range of application types.
   */
  struct has_history {
    /**
     * Constructor.
     *
     * \param  type           reconstruction type
     * \param  range          range of application types
     */
    has_history(const int type, const JRange<int> range) :
      type (type),
      range(range)
    {}

    /**
     * \param  track          track
     * \return                true if given track has specified history; else false
     */
    inline bool operator()(const Trk& track) const
    {
      if (track.rec_type == type)
        return std::find_if(track.rec_stages.begin(), track.rec_stages.end(), range) != track.rec_stages.end();
      else
        return false;
    }

    const int         type;    //!< reconstruction type
    const JRange<int> range;   //!< range of application types
  };
  
  /**
   * Test whether given track has muon prefit in history.
   *
   * \param  track          track
   * \return                true if muon prefit in history; else false
   */
  inline bool has_muon_prefit(const Trk& track)
  {
    return has_history(JPP_RECONSTRUCTION_TYPE, JFIT::JMUONPREFIT)(track);
  }

  /**
   * Test whether given track has muon simplex fit in history.
   *
   * \param  track          track
   * \return                true if muon simplex fit in history; else false
   */
  inline bool has_muon_simplex(const Trk& track)
  {
    return has_history(JPP_RECONSTRUCTION_TYPE, JFIT::JMUONSIMPLEX)(track);
  }

  /**
   * Test whether given track has muon gandalf fit in history.
   *
   * \param  track          track
   * \return                true if muon gandalf fit in history; else false
   */
  inline bool has_muon_gandalf(const Trk& track)
  {
    return has_history(JPP_RECONSTRUCTION_TYPE, JFIT::JMUONGANDALF)(track);
  }

  /**
   * Test whether given track has muon energy fit in history.
   *
   * \param  track          track
   * \return                true if muon energy fit in history; else false
   */
  inline bool has_muon_energy(const Trk& track)
  {
    return has_history(JPP_RECONSTRUCTION_TYPE, JFIT::JMUONENERGY)(track);
  }

  /**
   * Test whether given track has muon start fit in history.
   *
   * \param  track          track
   * \return                true if muon start fit in history; else false
   */
  inline bool has_muon_start(const Trk& track)
  {
    return has_history(JPP_RECONSTRUCTION_TYPE, JFIT::JMUONSTART)(track);
  }

  /**
   * Test whether given track has default muon fit in history.
   *
   * \param  track          track
   * \return                true if muon fit in history; else false
   */
  inline bool has_muon_fit(const Trk& track)
  {
    return has_history(JPP_RECONSTRUCTION_TYPE, JRange<int>(JFIT::JMUONBEGIN, JFIT::JMUONEND))(track);
  }  

  /**
   * Test whether given track has shower prefit in history.
   *
   * \param  track          track
   * \return                true if shower prefit in history; else false
   */
  inline bool has_shower_prefit(const Trk& track)
  {
    return has_history(JPP_RECONSTRUCTION_TYPE, JFIT::JSHOWERPREFIT)(track);
  }

  /**
   * Test whether given track has shower position fit in history.
   *
   * \param  track          track
   * \return                true if shower position fit in history; else false
   */
  inline bool has_shower_positionfit(const Trk& track)
  {
    return has_history(JPP_RECONSTRUCTION_TYPE, JFIT::JSHOWERPOSITIONFIT)(track);
  }  

  /**
   * Test whether given track has shower complete fit in history.
   *
   * \param  track          track
   * \return                true if shower complete fit in history; else false
   */
  inline bool has_shower_completefit(const Trk& track)
  {
    return has_history(JPP_RECONSTRUCTION_TYPE, JFIT::JSHOWERCOMPLETEFIT)(track);
  }

  /**
   * Test whether given track has default shower fit in history.
   *
   * \param  track          track
   * \return                true if shower fit in history; else false
   */
  inline bool has_shower_fit(const Trk& track)
  {
    return has_history(JPP_RECONSTRUCTION_TYPE, JRange<int>(JFIT::JSHOWERBEGIN, JFIT::JSHOWEREND))(track);
  }

  /**
   * Test whether given track has default shower fit in history.
   *
   * \param  track          track
   * \return                true if shower fit in history; else false
   */
  inline bool has_aashower_fit(const Trk& track)
  {
    return has_history(AASHOWER_RECONSTRUCTION_TYPE, JRange<int>())(track);
  }

  /**
   * Test whether given event has a track according selection.\n
   * The track selector corresponds to the function operator <tt>bool selector(const Trk&);</tt>.
   *
   * \param  evt            event
   * \param  selector       track selector
   * \return                true if at least one corresponding track; else false
   */
  template<class JTrackSelector_t>
  inline bool has_reconstructed_track(const Evt& evt, JTrackSelector_t selector)
  {
    return std::find_if(evt.trks.begin(), evt.trks.end(), selector) != evt.trks.end();
  }

  /**
   * Test whether given event has a track according selection.
   *
   * \param  evt            event
   * \param  range          range of application types
   * \return                true if at least one corresponding track; else false
   */
  template<int reconstruction_type>
  inline bool has_reconstructed_track(const Evt& evt, const JRange<int> range = JRange<int>())
  {
    return has_reconstructed_track(evt, has_history(reconstruction_type, range));
  }
  
  /**
   * Test whether given event has a track with muon reconstruction.
   *
   * \param  evt            event
   * \return                true if at least one reconstructed muon; else false
   */
  inline bool has_reconstructed_muon(const Evt& evt)
  {
    return has_reconstructed_track(evt, has_muon_fit);
  }
  
  /**
   * Test whether given event has a track with shower reconstruction.
   *
   * \param  evt            event
   * \return                true if at least one reconstructed shower; else false
   */
  inline bool has_reconstructed_shower(const Evt& evt)
  {
    return has_reconstructed_track(evt, has_shower_fit);
  }
  
  /**
   * Test whether given event has a track with aashower reconstruction.
   *
   * \param  evt            event
   * \return                true if at least one reconstructed shower; else false
   */
  inline bool has_reconstructed_aashower(const Evt& evt)
  {
    return has_reconstructed_track(evt, has_aashower_fit);
  }

  /**
   * Get best reconstructed track.\n
   * The track selector corresponds to the function operator <tt>bool selector(const Trk&);</tt> and
   * the track comparator to <tt>bool comparator(const Trk&, const Trk&);</tt>.
   *
   * \param  evt            event
   * \param  selector       track selector
   * \param  comparator     track comparator
   * \return                track
   */
  template<class JTrackSelector_t, class JQualitySorter_t>
  inline const Trk& get_best_reconstructed_track(const Evt&       evt,
                                                 JTrackSelector_t selector,
                                                 JQualitySorter_t comparator)
  {
    std::vector<Trk>::const_iterator p = std::find_if(evt.trks.begin(), evt.trks.end(), selector);

    for (std::vector<Trk>::const_iterator i = p; i != evt.trks.end(); ++i) {
      if (selector(*i) && comparator(*i, *p)) {
        p = i;
      }
    }

    if (p != evt.trks.end())
      return *p;
    else
      THROW(JIndexOutOfRange, "This event has no reconstructed track with given selector.");
  }

  /**
   * Get best reconstructed track.
   *
   * \param  evt            event
   * \param  range          range of application types
   * \return                track
   */
  template<int reconstruction_type>
  inline const Trk& get_best_reconstructed_track(const Evt& evt, const JRange<int> range = JRange<int>())
  {
    return get_best_reconstructed_track(evt, has_history(reconstruction_type, range), quality_sorter<reconstruction_type>());
  }


  /**
   * Get best reconstructed muon.
   *
   * \param  evt            event
   * \return                track
   */
  inline const Trk& get_best_reconstructed_muon(const Evt& evt)
  {
    return get_best_reconstructed_track(evt, has_muon_fit, quality_sorter<JPP_RECONSTRUCTION_TYPE>());
  }

  /**
   * Get best reconstructed shower.
   *
   * \param  evt            event
   * \return                track
   */
  inline const Trk& get_best_reconstructed_shower(const Evt& evt)
  {
    return get_best_reconstructed_track(evt, has_shower_fit, quality_sorter<JPP_RECONSTRUCTION_TYPE>());
  }

  /**
   * Get best reconstructed aashower.
   *
   * \param  evt            event
   * \return                track
   */
  inline const Trk& get_best_reconstructed_aashower(const Evt& evt)
  {
    return get_best_reconstructed_track(evt, has_aashower_fit, quality_sorter<AASHOWER_RECONSTRUCTION_TYPE>());
  }
}

#endif