JReconstruction/JEvtToolkit.hh

Go to the documentation of this file.

#ifndef __JFIT__JEVTTOOLKIT__
#define __JFIT__JEVTTOOLKIT__

#include <istream>
#include <ostream>
#include <algorithm>
#include <cmath>

#include "km3net-dataformat/definitions/reconstruction.hh"

#include "JLang/JException.hh"
#include "JLang/JPredicate.hh"
#include "JLang/JFind_if.hh"
#include "JTools/JRange.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JVersor3D.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JTrack3E.hh"
#include "JMath/JMathToolkit.hh"
#include "JMath/JConstants.hh"
#include "JTools/JRange.hh"
#include "Jeep/JFunctionAdaptor.hh"

#include "JFit/JLine1Z.hh"
#include "JFit/JShower3EZ.hh"

#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JHistory.hh"


/**
 * \author mdejong
 */

using namespace std;

namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }

/**
 * Model fits to data.
 */
namespace JRECONSTRUCTION {

  using JLANG::JIndexOutOfRange;
  using JGEOMETRY3D::JAngle3D;
  using JGEOMETRY3D::JVersor3D;
  using JGEOMETRY3D::JVector3D;
  using JGEOMETRY3D::JPosition3D;
  using JGEOMETRY3D::JDirection3D;
  using JGEOMETRY3D::JAxis3D;
  using JGEOMETRY3D::JTrack3D;
  using JGEOMETRY3D::JTrack3E;
  using JEEP::JFunctionAdaptor;
  using JTOOLS::JRange;
  using JFIT::JLine1Z;
  using JFIT::JShower3Z;


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


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


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


  /**
   * Get track.
   *
   * \param  fit           fit
   * \return               track
   */
  inline JTrack3E getTrack(const JFit& fit)
  {
    return JTrack3E(JTrack3D(getAxis(fit), fit.getT()), fit.getE());
  }

  
  /**
   * Get fit.
   *
   * \param  history       history
   * \param  track         track
   * \param  Q             quality
   * \param  NDF           number of degrees of freedom
   * \param  energy        Energy, which would be set as 0, if the fit does not reconstruct energy
   * \param  status        status of the fit as defined in enum JFitStatus.hh
   * \return               fit
   */
  inline JFit getFit(const JHistory& history,
                     const JTrack3D& track,
                     const double    Q,
                     const int       NDF,
                     const double    energy = 0.0,
                     const int       status = 0)
  {
    return JFit(history,
                track.getX(),  track.getY(),  track.getZ(),
                track.getDX(), track.getDY(), track.getDZ(),
                track.getT(),
                Q, NDF,
                energy, status);
  }


  /**
   * Get fit.
   *
   * \param  history       history
   * \param  track         track
   * \param  angle         angle
   * \param  Q             quality
   * \param  NDF           number of degrees of freedom
   * \param  energy        Energy, which would be set as 0, if the fit does not reconstruct energy
   * \param  status        status of the fit as defined in JFitStatus.hh
   * \return               fit
   */
  inline JFit getFit(const JHistory& history,
                     const JLine1Z&  track,
                     const JAngle3D& angle,
                     const double    Q,
                     const int       NDF,
                     const double    energy = 0.0,
                     const int       status = 0)
  {
    return JFit(history,
                track.getX(),  track.getY(),  track.getZ(),
                angle.getDX(), angle.getDY(), angle.getDZ(),
                track.getT(),
                Q, NDF,
                energy, status);
  }
  

  /**
   * Get dot product.
   *
   * \param  first         first  fit
   * \param  second        second fit
   * \return               dot product
   */
  inline double getDot(const JFit& first, const JFit& second)
  {
    return JMATH::getDot(getDirection(first), getDirection(second));
  }


  /**
   * Get dot product.
   *
   * \param  fit           fit
   * \param  dir           direction
   * \return               dot product
   */
  inline double getDot(const JFit& fit, const JDirection3D& dir)
  {
    return JMATH::getDot(getDirection(fit), dir);
  }


  /**
   * Get quality of fit.\n
   * The larger the quality, the better the fit.
   *
   * \param  chi2          chi2
   * \param  N             number of hits
   * \param  NDF           number of degrees of freedom
   * \return               quality
   */
  inline double getQuality(const double chi2, const int N, const int NDF)
  {
    return N  -  0.25 * chi2 / NDF;
  }


  /**
   * Get quality of fit.\n
   * The larger the quality, the better the fit.
   *
   * \param  chi2          chi2
   * \param  NDF           number of degrees of freedom
   * \return               quality
   */
  inline double getQuality(const double chi2, const int NDF)
  {
    return NDF  -  0.25 * chi2 / NDF;
  }


  /**
   * Get quality of fit.\n
   * The larger the quality, the better the fit.
   *
   * \param  chi2          chi2
   * \return               quality
   */
  inline double getQuality(const double chi2)
  {
    return -chi2;
  }


  /**
   * Comparison of fit results.
   *
   * \param  first         first  fit
   * \param  second        second fit
   * \return               true if first fit has better quality than second; else false
   */
  inline bool qualitySorter(const JFit& first, const JFit& second)
  {
    if (first.getHistory().size() == second.getHistory().size())
      return first.getQ() > second.getQ();
    else
      return first.getHistory().size() > second.getHistory().size();
  }


  /**
   * General purpose sorter of fit results.
   *
   * The default constructor will sort fit results according the method JRECONSTRUCTION::qualitySorter.\n
   * A different method can dynamically be loaded from a (shared) library using class JEEP::JFunctionAdaptor.
   * For the definition of an alternative method, see e.g.\ quality_sorter.cc 
   */
  struct JQualitySorter :
    public JFunctionAdaptor<bool, const JFit&, const JFit&>
  {
    typedef JFunctionAdaptor<bool, const JFit&, const JFit&>  function_adaptor_type;
    typedef function_adaptor_type::pF                         pF;
    
    
    /**
     * Default constructor.
     */
    JQualitySorter() :
      function_adaptor_type(qualitySorter)
    {}
  };


  /**
   * Test whether given fit has specified history.
   *
   * \param  fit            fit
   * \param  type           application type
   * \return                true if type in history; else false
   */
  inline bool has_history(const JFit& fit, const int type)
  {
    return std::find_if(fit.getHistory().begin(),
                        fit.getHistory().end(),
                        JLANG::make_predicate(&JEvent::type, type)) != fit.getHistory().end();
  }


  /**
   * Test whether given fit has specified history.
   *
   * \param  fit            fit
   * \param  range          application type range
   * \return                true if type in history; else false
   */
  inline bool has_history(const JFit& fit, const JRange<int>& range)
  {
    return std::find_if(fit.getHistory().begin(),
                        fit.getHistory().end(),
                        JLANG::make_find_if(&JEvent::type, range)) != fit.getHistory().end();
  }


  /**
   * Test whether given fit has muon prefit in history.
   *
   * \param  fit            fit
   * \return                true if muon prefit in history; else false
   */
  inline bool has_muon_prefit(const JFit& fit)
  {
    return has_history(fit, JMUONPREFIT);
  }


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


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


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


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


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


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


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


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


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


  /**
   * 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 JEvt& evt, JTrackSelector_t selector)
  {
    return std::find_if(evt.begin(), evt.end(), selector) != evt.end();
  }


  /**
   * 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 JEvt& 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 JEvt& evt)
  {
    return has_reconstructed_track(evt, has_shower_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 comprator(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 JFit& get_best_reconstructed_track(const JEvt&      evt,
                                                  JTrackSelector_t selector,
                                                  JQualitySorter_t comparator)
  {
    JEvt::const_iterator p = std::find_if(evt.begin(), evt.begin(), selector);

    for (JEvt::const_iterator i = p; i != evt.end(); ++i) {
      if (selector(*i) && comparator(*i, *p)) {
        p = i;
      }
    }

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


  /**
   * Get best reconstructed muon.
   *
   * \param  evt            event
   * \return                track
   */
  inline const JFit& get_best_reconstructed_muon(const JEvt& evt)
  {
    return get_best_reconstructed_track(evt, has_muon_fit, qualitySorter);
  }


  /**
   * Get best reconstructed shower.
   *
   * \param  evt            event
   * \return                track
   */
  inline const JFit& get_best_reconstructed_shower(const JEvt& evt)
  {
    return get_best_reconstructed_track(evt, has_shower_fit, qualitySorter);
  }


  /**
   * Auxiliary class to compare fit results with respect to a reference direction (e.g.\ true muon).
   * The sort operation results in an ordered set of fit results with increasing angle between
   * the reference direction and that of the fit results.
   */
  class JPointing {
  public:
    /**
     * Default constructor.
     */
    JPointing()
    {}


    /**
     * Constructor.
     *
     * \param  dir           reference direction
     */
    JPointing(const JVersor3D& dir)
    {
      this->dir = dir;
    }


    /**
     * Get direction.
     *
     * \return               direction
     */
    JVersor3D getDirection() const
    {
      return dir;
    }


    /**
     * Get dot product between reference direction and fit result.
     *
     * \param  fit           fit
     * \return               dot product
     */
    inline double getDot(const JFit& fit) const
    {
      return JRECONSTRUCTION::getDot(fit, dir);
    }


    /**
     * Get angle between reference direction and fit result.
     *
     * \param  fit           fit
     * \return               angle [deg]
     */
    inline double getAngle(const JFit& fit) const
    {
      const double dot = getDot(fit);

      if      (dot >= +1.0)
        return   0.0;
      else if (dot <= -1.0)
        return 180.0;
      else
        return acos(dot) * 180.0 / JMATH::PI;
    }


    /**
     * Comparison of fit results.
     *
     * \param  first         first  fit
     * \param  second        second fit
     * \return               true if first fit worse; else false
     */
    inline bool operator()(const JFit& first, const JFit& second) const
    {
      return this->getDot(first) < this->getDot(second);
    }


    /**
     * Select best fit result.
     *
     * \param  __begin       begin of fit results
     * \param  __end         end   of fit results
     * \return               best fit result
     */
    template<class T>
    inline T operator()(T __begin, T __end) const
    {
      return std::max_element(__begin, __end, *this);
    }

  protected:
    JVersor3D dir;
  };

  /**
   * Auxiliary class to compare fit results with respect to a reference position (e.g.\ true muon).
   * The sort operation results in an ordered set of fit results with increasing distance between
   * the reference position and that of the fit results.
   */
  class JPosition {
  public:
    /**
     * Constructor.
     *
     * \param  pos           reference position
     */
    JPosition(JVector3D pos)
    {
      this->pos = pos;
    }

    /**
     * Comparison of fit results.
     *
     * \param  first         first  fit
     * \param  second        second fit
     * \return               true if first fit better; else false
     */
    inline bool operator()(const JFit& first, const JFit& second) const
    {
      return this->pos.getDistance(getPosition(first)) < this->pos.getDistance(getPosition(second));
    }

     /**
     * Select best fit result.
     *
     * \param  __begin       begin of fit results
     * \param  __end         end   of fit results
     * \return               best fit result
     */
    template<class T>
    inline T operator()(T __begin, T __end) const
    {
      return std::min_element(__begin, __end, *this);
    }
    
  protected:
    JVector3D pos;
  };


  class JShowerEnergy {
  public:
    /**                                                                                                                 
     * Constructor.                                                                                                     
     *                                                                                                     
     * \param  E                reference energy
     */
    JShowerEnergy(double E)
    {
      this->energy = E;
    }
    
    /**
     * Comparison of fit results.
     *
     * \param  first         first  fit
     * \param  second        second fit
     * \return               true if first fit better; else false
     */
    inline bool operator()(const JFit& first, const JFit& second) const
    {
      return abs(this->energy - first.getE()) < abs(this->energy - second.getE());
    }

    /**
     * Select best fit result.
     *
     * \param  __begin       begin of fit results
     * \param  __end         end   of fit results
     * \return               best fit result
     */
    template<class T>
    inline T operator()(T __begin, T __end) const
    {
      return std::min_element(__begin, __end, *this);
    }
    
  protected:
    double energy;
  };


  /**
   * Auxiliary class to evaluate atmospheric muon hypothesis.
   * The hypothesis is tested by means of the difference in quality 
   * between the best upward and best downward track.
   */
  class JAtmosphericMuon {
  public:
    /**
     * Default constructor.
     */
    JAtmosphericMuon() :
      dot1(0.0),
      dot2(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  theta1        upper hemisphere angle [deg]
     * \param  theta2        lower hemisphere angle [deg]
     */
    JAtmosphericMuon(const double theta1,
                     const double theta2) :
      dot1(cos(theta1 * JMATH::PI/180.0)),
      dot2(cos(theta2 * JMATH::PI/180.0))
    {}


    /**
     * Test is event is atmospheric muon.
     *
     * \param  __begin       begin of data
     * \param  __end         end   of data
     * \return               negative if preferably down / positive if preferably up
     */
    double operator()(JEvt::const_iterator __begin,
                      JEvt::const_iterator __end) const
    {
      double Qup   = 0.0;
      double Qdown = 0.0;

      for (JEvt::const_iterator i = __begin; i != __end; ++i) {

        if        (i->getDZ() >= dot1) {
            
          if (i->getQ() > Qup) {
            Qup   = i->getQ();
          }
          
        } else if (i->getDZ() <= dot2) {
          
          if (i->getQ() > Qdown) {
            Qdown = i->getQ();
          }
        }
      }

      return Qup - Qdown;
    }


    /**
     * Read atmospheric muon analyser from input.
     *
     * \param  in            input stream
     * \param  object        atmospheric muon analyser
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JAtmosphericMuon& object)
    {
      double theta1, theta2;

      in >> theta1 >> theta2;

      object.dot1 = cos(theta1 * JMATH::PI/180.0);
      object.dot2 = cos(theta2 * JMATH::PI/180.0);

      return in;
    }


    /**
     * Write atmospheric muon analyser to output.
     *
     * \param  out           output stream
     * \param  object        atmospheric muon analyser
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JAtmosphericMuon& object)
    {
      out << acos(object.dot1) * 180.0 / JMATH::PI << ' '
          << acos(object.dot2) * 180.0 / JMATH::PI;

      return out;
    }


    double dot1;
    double dot2;
  };
}

#endif