Trk.hh

Go to the documentation of this file.

#ifndef TRK_HH_INCLUDED
#define TRK_HH_INCLUDED
 
#include "AAObject.hh"
#include "Vec.hh"
#include <vector>
#include "TDatabasePDG.h"
 
/** 
 * The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower. 
 */
struct Trk: public AAObject
{
  int    id;                          ///< track identifier
  Vec    pos;                         ///< postion of the track at time t
  Vec    dir;                         ///< track direction
  double t;                           ///< track time (when the particle is at pos )
  double E;                           ///< Energy (either MC truth or reconstructed)
 
  double len;                         ///< length, if applicable
  double lik;                         ///< likelihood or lambda value (for aafit, lambda)
  int    type;                        ///< MC: particle type in PDG encoding.
  int               rec_type;         ///< identifyer for the overall fitting algorithm/chain/strategy
  std::vector<int>  rec_stages;       ///< list of identifyers of succesfull fitting stages resulting in this track
 
  int status;                         ///< MC status code
  int mother_id;                      ///< MC id of the parent particle 
 
  std::vector<double> fitinf;         ///< place to store additional fit info, for jgandalf, see JFitParameters.hh
  std::vector<int>    hit_ids;        ///< list of associated hit-ids (corresponds to Hit::id).
  std::vector<double> error_matrix;   ///< (5x5) error covariance matrix (stored as linear vector)
  std::string    comment;             ///< use as you like
 
  /**
   * Default constructor.
   */
  Trk(): id(0),t(0),E(0),len(0),lik(0), type(0), rec_type(0), status(0), mother_id(-1) {}
 
  /**
   * Read track (useful in python).
   *
   * \param  t             track
   */
  void read(const Trk& t) { *this = t;}
 
  /**
   * Write track (useful in python).
   *
   * \param  t             track
   */
  void write(Trk& t)      const { t = *this; }
 
 
  /**
   * Get distance of closes approach to another track.
   *
   * \param  t             track
   * \return               distance
   */
  double distance_to (const Trk& t) const
  {
    Vec v = dir.cross( t.dir );
    v /= v.len();
    const double d = (pos - t.pos ).dot( v );
    return d>0?d:-d;
  }
 
  /**
   * Get distance to given position.
   *
   * \param  p             position
   * \return               distance
   */
  double distance_to( const Vec& p = Vec() ) const
  {
    Vec v = p - pos;
    double dz = v.dot(dir);
    return v.dot(v) - dz*dz;
  }
 
  /**
   * Get the name of the MC particle type.
   *
   * \return               name
   */
  std::string name() const
  {
    TParticlePDG* p = TDatabasePDG::Instance()->GetParticle( type );
    if (!p) return "unnamed state ("+ std::to_string(type)+")";
    return p->GetName();
  }
 
  /**
   * Check if this is a primary particle (i.e. has not mother).
   *
   * \return               true if primary; else false
   */  
  bool is_primary() const
  {
    return mother_id == -1; // and/or status==0 ?
  }
  
  /**
   * Check if this is a netrino.
   *
   * Note that its is checked if the PDG type is a nu-e, nu-mu or nu-tau.
   *
   * \return               true if neutrino; else false
   */  
  bool is_neutrino() const
  {
    const unsigned t = abs(type);
    return t == 12 || t == 14 || t == 16;
  }
 
  /**
   * Check if this is a charged lepton.
   *
   * Note that its is checked if the PDG type is a (anti-)electron, (anti-)muon or (anti-)tua.
   *
   * \return               true if charged lepton; else false
   */
  bool is_lepton() const 
  {
    const unsigned t = abs(type);
    return t == 11 || t == 13 || t == 15;
  }
  
  /**
   * Check if this is a muon.
   *
   * Note that its is checked if the PDG type is a (anti-)muon.
   *
   * \return               true if muon; else false
   */
  bool is_muon() const { return abs(type)==13; }
 
  /**
   * Check if this is an orphan.
   *
   * If a particle's (intermediate) mother has been deleted from the track-list (the 3rd digit of status==1).
   *
   * \return               true if orphan; else false
   */  
  bool is_orphan() const
  {
    return (status / 100) % 10 == 1;
  }
  
  /**
   * Get list of of pointers to tracks, all of which have their mother identifier set to identifier of this track.
   *
   * \param  mctrks        list of input tracks
   * \return               list of pointers to tracks
   */
  std::vector< Trk* > get_daughters ( std::vector<Trk>& mctrks )
  {
    std::vector<Trk*> r;
 
    for( auto& t : mctrks )
      {
        if ( t.mother_id == id ) r.push_back( &t );
      }
    return r;
  }
 
  /**
   * Print track.
   *
   * \param  out           output stream
   */
  void print(std::ostream& out) const
  {
    out << "Trk: id=" << id << " pos="; pos.print(out);
    out << " dir="; dir.print(out);
    out << " t=" << t << " E=" << E << " pdg-type=" << type;
  }
 
  ClassDef(Trk,11)
};
 
#endif