JAAnetToolkit.hh

Go to the documentation of this file.

#ifndef __JAANET__JAANETTOOLKIT__
#define __JAANET__JAANETTOOLKIT__

#include <stdlib.h>
#include <limits>
#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/fitparameters.hh"
#include "km3net-dataformat/definitions/reconstruction.hh"
#include "km3net-dataformat/tools/reconstruction.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 "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 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.\ 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.getLength() < T_ns.getLength()) {
      time_range.setLength(T_ns.getLength());
    }

    return time_range;
  }


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

#endif