Evt.hh

Go to the documentation of this file.

#ifndef EVT_HH_INCLUDED
#define EVT_HH_INCLUDED

#include "AAObject.hh"
#include "Hit.hh"
#include "Trk.hh"
#include "Hit.hh"
#include "Exception.hh"

#include "TTimeStamp.h"

#include <vector>

/**
 * The Evt class respresent a Monte Carlo (MC) event as well as an offline event.\n
 * Some data from the online (DAQ) event are copied.
 */

struct Evt: public AAObject
{
  int id;                       ///< offline event identifier 
  int det_id;                   ///< detector identifier from DAQ
  int mc_id;                    ///< identifier of the MC event (as found in ascii or antcc file).
  
  int run_id;                   ///< DAQ run identifier
  int mc_run_id;                ///< MC  run identifier
  
  int frame_index;              ///< from the raw data
  ULong64_t trigger_mask;       ///< trigger mask from raw data (i.e. the trigger bits)
  ULong64_t trigger_counter;    ///< trigger counter
  unsigned int overlays;        ///< number of overlaying triggered events
  TTimeStamp t;                 ///< UTC time of the start of the timeslice the event came from
  
  //hits and tracks
  std::vector<Hit> hits;        ///< list of hits
  std::vector<Trk> trks;        ///< list of reconstructed tracks (can be several because of prefits,showers, etc).   

  //Monte carlo 
  std::vector<double> w;         ///< MC: Weights w[0]=w1, w[1]=w2, w[2]]=w3 (see e.g. <a href="https://simulation.pages.km3net.de/taglist/taglist.pdf">Tag list</a>)
  std::vector<double> w2list;    ///< MC: factors that make up w[1]=w2       (see e.g. <a href="https://simulation.pages.km3net.de/taglist/taglist.pdf">Tag list</a>)
  std::vector<double> w3list;    ///< MC: atmospheric flux information

  double mc_t;                   ///< MC: time of the mc event
  std::vector<Hit> mc_hits;      ///< MC: list of MC truth hits
  std::vector<Trk> mc_trks;      ///< MC: list of MC truth tracks  

  // --- place to store user info ---
  TString  comment;              ///< user can use this as he/she likes
  int      index;                ///< user can use this as he/she likes 
  int      flags;                ///< user can use this as he/she likes 

  
  // basic ctor and trival io
  Evt() :
    id(0), det_id(0), mc_id(0), run_id(0), mc_run_id(0), frame_index(0),
    trigger_mask(0), trigger_counter(0),
    overlays(0), t(0), mc_t(0), index(0), flags(0) {}


  /**
   * Print event.
   *
   * \param  out           output stream
   */
  void print(std::ostream& out = std::cout) const
  {
    out << "Evt: id=" << id << 
      " run_id=" << run_id << 
      " #hits=" << hits.size() <<
      " #mc_hits=" << mc_hits.size() <<
      " #trks=" << trks.size() <<
      " #mc_trks=" << mc_trks.size() << std::endl;
  }


  /**
   * Get UTC time.
   *
   * \return               time [ns]
   */
  double get_ns()
  {
    return this->t.GetSec()*1e9 + this->t.GetNanoSec(); 
  }

  /**
   * Reset event.
   */
  void clear()
  {
    *this = Evt();
  }

  /**
   * Clear the hit vectors and all the references to hits in the tracks
   */
  void clear_hits() 
  {
    hits.clear();
    mc_hits.clear();
    for (auto& t: trks )    t.hit_ids.clear();
    for (auto& t: mc_trks ) t.hit_ids.clear(); 
  }


  /**
   * Get the original Monte-Carlo-time of the hit.
   *
   * Apply a correction to the hit time from the DAQ, so that the 
   * time returned by this fucntion matches the original time used 
   * in the Monte Carlo event.
   *
   * \param  h             hit
   * \return               time [ns]
   */
  double getMCtime(const Hit& h) 
  {
    return h.t+this->get_ns() - this->mc_t;
  }

  /**
   * Get pointer to hit with given identifier.\n
   * This method throws a run-time exception if no hit is available.
   *
   * \param  req_id        hit identifier
   * \param  use_mc_hits   option to use MC hit
   * \return               pointer to hit.
   */
  Hit* get_hit_by_id( int req_id , bool use_mc_hits = false )
  {
    std::vector<Hit>& col = use_mc_hits? mc_hits : hits;

    if ( req_id < int(col.size()) ) // try to see if we are lucky
      {
        if ( col[req_id].id == req_id ) return &(col[id]);
      }

    for( auto& h: col ) 
      {
        if ( h.id == req_id ) return &(col[id]);
      }

    THROW(Exception, "Could not find hit with id=" << req_id);
  }


  /**
   * Get all the tracks that are 'primary'.\n
   * Here, primary means the tracks have no parents.
   *
   * Note that this method only works if MC parent-child relations are availabe.
   *
   * \return               list of pointers to primary tracks
   */
  std::vector<Trk*> primary_trks()
  {
    std::vector<Trk*> r;
    for(auto& t: mc_trks )
      {
        if ( t.is_primary() ) r.push_back(&t);
      }
    return r;
  }

  /**
   * Get a pointer to the (first) neutrino from the MC track list. 
   *
   * \return               pointers to neutrino (may be NULL)
   */
  Trk* neutrino()
  {
    for(auto& t: mc_trks )
      {
        if ( t.is_neutrino() ) return &t;
      }
    return nullptr;
  }

  /**
   * Get a pointer to primary neutrino from the MC track list. 
   *
   * Note that this method only works if MC parent-child relations are availabe.
   *
   * \return               pointer to primary neutrino (may be NULL)
   */
  Trk* primary_neutrino() 
  {
    for( auto& t : mc_trks ) 
      {
        if ( t.is_neutrino() and t.is_primary() ) return &t;
      }
    return nullptr;
  }

  /**
   * Get a pointer to leading lepton from the MC track list. 
   * Here, leading means the lepton that has a neutrino as mother.
   *
   * Note that this method only works if MC parent-child relations are availabe.
   *
   * \return               pointer to leadig lepton (may be NULL)
   */
  Trk* leading_lepton() 
  {
    int nu_id = primary_neutrino()->id;

    for(auto& t: mc_trks )
      {
        if ( t.is_lepton()        &&
             t.mother_id == nu_id &&
             !t.is_orphan() ) return &t;
      }
    return nullptr;
  }

  /**
   * Get the (first) parent of the track 'child'.\n
   * This method throws a run-time exception if no parent is available.
   *
   * \param  child         child particle
   * \return               pointer to parent
   */
  Trk* get_parent_of( Trk& child ) 
  {  
    for(auto& t: mc_trks )
      {
        if (child.mother_id == t.id ) return &t;
      }
    THROW(Exception, "Failed to find parent of track" << child.id << " mother_id="<< child.mother_id);
  }  

  ClassDef(Evt,12)
};

#endif