JReconstruction/JEvtToolkit.hh

Go to the documentation of this file.

#ifndef __JRECONSTRUCTION__JEVTTOOLKIT__
#define __JRECONSTRUCTION__JEVTTOOLKIT__
 
#include <string>
#include <istream>
#include <ostream>
#include <map>
#include <algorithm>
#include <cmath>
 
#include "km3net-dataformat/definitions/reconstruction.hh"
#include "km3net-dataformat/definitions/fitparameters.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/JVertex3D.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JTrack3E.hh"
#include "JGeometry3D/JShower3D.hh"
#include "JGeometry3D/JShower3E.hh"
#include "JGeometry3D/JOmega3D.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/JFitStatus.hh"
#include "JReconstruction/JHistory.hh"
 
 
/**
 * \author mdejong
 */
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
 
/**
 * Model fits to data.
 */
namespace JRECONSTRUCTION {
 
  using JLANG::JIndexOutOfRange;
  using JGEOMETRY3D::JAngle3D;
  using JGEOMETRY3D::JVector3D;
  using JGEOMETRY3D::JPosition3D;
  using JGEOMETRY3D::JDirection3D;
  using JGEOMETRY3D::JAxis3D;
  using JGEOMETRY3D::JVertex3D;
  using JGEOMETRY3D::JTrack3D;
  using JGEOMETRY3D::JTrack3E;
  using JGEOMETRY3D::JShower3D;
  using JGEOMETRY3D::JShower3E;
  using JEEP::JFunctionAdaptor;
  using JTOOLS::JRange;
  using JFIT::JLine1Z;
  using JFIT::JShower3Z;
 
 
  static const int JSTART_ZMIN_M                    = 11;        //!< start position of track
  static const int JSTART_ZMAX_M                    = 12;        //!< end position of track
 
 
  /**
   * Auxiliary data structure to get weight of given fit.
   */
  static struct JWeight : 
    public std::map<std::string, int>
  {
    /**
     * Default constructor.
     */
    JWeight() 
    {
#define MAKE_ENTRY(A) std::make_pair(#A, A)
 
      this->insert(MAKE_ENTRY(JGANDALF_BETA0_RAD));
      this->insert(MAKE_ENTRY(JGANDALF_BETA1_RAD));
      this->insert(MAKE_ENTRY(JGANDALF_LIKELIHOOD_RATIO));
      this->insert(MAKE_ENTRY(JGANDALF_NUMBER_OF_HITS));
      this->insert(MAKE_ENTRY(JENERGY_ENERGY));
      this->insert(MAKE_ENTRY(JENERGY_CHI2));
      this->insert(MAKE_ENTRY(JGANDALF_LAMBDA));
      this->insert(MAKE_ENTRY(JGANDALF_NUMBER_OF_ITERATIONS));
      this->insert(MAKE_ENTRY(JSTART_ZMIN_M));
      this->insert(MAKE_ENTRY(JSTART_ZMAX_M));
      this->insert(MAKE_ENTRY(JSTART_NPE_MIP_TOTAL));
      this->insert(MAKE_ENTRY(JSTART_NPE_MIP_MISSED));
      this->insert(MAKE_ENTRY(JSTART_LENGTH_METRES));
      this->insert(MAKE_ENTRY(JSTART_BACKGROUND_LOGP));
      this->insert(MAKE_ENTRY(JVETO_NPE));
      this->insert(MAKE_ENTRY(JVETO_NUMBER_OF_HITS));
      this->insert(MAKE_ENTRY(JENERGY_MUON_RANGE_METRES));
      this->insert(MAKE_ENTRY(JENERGY_NOISE_LIKELIHOOD));
      this->insert(MAKE_ENTRY(JENERGY_NDF));
      this->insert(MAKE_ENTRY(JENERGY_NUMBER_OF_HITS));
      this->insert(MAKE_ENTRY(JCOPY_Z_M));
      this->insert(MAKE_ENTRY(JPP_COVERAGE_ORIENTATION));
      this->insert(MAKE_ENTRY(JPP_COVERAGE_POSITION));
      this->insert(MAKE_ENTRY(JENERGY_MINIMAL_ENERGY));
      this->insert(MAKE_ENTRY(JENERGY_MAXIMAL_ENERGY));
      this->insert(MAKE_ENTRY(JSHOWERFIT_ENERGY));
      this->insert(MAKE_ENTRY(AASHOWERFIT_ENERGY));
      this->insert(MAKE_ENTRY(AASHOWERFIT_NUMBER_OF_HITS));
 
#undef MAKE_ENTRY
    }
 
 
    /**
     * Has weight.
     *
     * \param  key             key
     * \return                 true if valid key; else false
     */
    bool operator()(const std::string& key) const
    {
      return this->count(key) != 0;
    }
 
 
    /**
     * Get weight.
     *
     * \param  fit             fit
     * \param  key             key
     * \param  value           default value
     * \return                 value
     */
    double operator()(const JFit& fit, const std::string& key, const double value) const
    {
      const_iterator p = this->find(key);
 
      if (p != this->end() && fit.hasW(p->second))
        return fit.getW(p->second);
      else
        return value;
    }
 
  } getWeight;
 
 
  /**
   * 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 vertex.
   *
   * \param  fit           fit
   * \return               vertex
   */
  inline JVertex3D getVertex(const JFit& fit)
  {
    return JVertex3D(getPosition(fit), fit.getT());
  }
 
 
  /**
   * 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 shower.
   *
   * \param  fit           fit
   * \return               shower
   */
  inline JShower3E getShower(const JFit& fit)
  {
    return JShower3E(JShower3D(getVertex(fit), getDirection(fit)), 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 = SINGLE_STAGE)
  {
    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 = SINGLE_STAGE)
  {
    return JFit(history,
                track.getX(),  track.getY(),  track.getZ(),
                angle.getDX(), angle.getDY(), angle.getDZ(),
                track.getT(),
                Q, NDF,
                energy, status);
  }
 
  
  /**
   * Get fit.
   *
   * \param  history       history
   * \param  shower        shower
   * \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 JShower3D& shower,
                     const double     Q,
                     const int        NDF,
                     const double     energy = 0.0,
                     const int        status = SINGLE_STAGE)
  {
    return JFit(history,
                shower.getX(),  shower.getY(),  shower.getZ(),
                shower.getDX(), shower.getDY(), shower.getDZ(),
                shower.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 space angle.
   *
   * \param  first         first  fit
   * \param  second        second fit
   * \return               angle [deg]
   */
  inline double getAngle(const JFit& first, const JFit& second)
  {
    return JMATH::getAngle(getDirection(first), getDirection(second));
  }
 
 
  /**
   * Get space angle.
   *
   * \param  fit           fit
   * \param  dir           direction
   * \return               angle [deg]
   */
  inline double getAngle(const JFit& fit, const JDirection3D& dir)
  {
    return JMATH::getAngle(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 JFit&);</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 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 JFit& get_best_reconstructed_track(const JEvt&      evt,
                                                  JTrackSelector_t selector,
                                                  JQualitySorter_t comparator)
  {
    JEvt::const_iterator p = std::find_if(evt.begin(), evt.end(), 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 JDirection3D& dir)
    {
      this->dir = dir;
    }
 
 
    /**
     * Constructor.
     *
     * \param  fit           fit
     */
    JPointing(const JFit& fit)
    {
      this->dir = JRECONSTRUCTION::getDirection(fit);
    }
 
 
    /**
     * Get direction.
     *
     * \return               direction
     */
    JDirection3D getDirection() const
    {
      return 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, dir);
 
      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 better; else false
     */
    inline bool operator()(const JFit& first, const JFit& second) const
    {
      return getDot(first, dir) > getDot(second, dir);
    }
 
 
    /**
     * 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:
    JDirection3D dir;
  };
 
 
  /**
   * Auxiliary class to compare fit results with respect to a reference position.
   * 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;
  };
 
 
  /**
   * Auxiliary class to compare fit results with respect to a reference energy.
   * The sort operation results in an ordered set of fit results with increasing difference between
   * the reference energy and that of the fit results.
   */
  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 fabs(this->energy - first.getE()) < fabs(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;
  };
 
 
  /**
   * Gridify set of fits.
   *
   * \param  __begin         begin of fits
   * \param  __end           end   of fits
   * \param  N               number of directions in grid
   * \return                 end   of fits sorted by quality
   */
  inline JEvt::iterator gridify(JEvt::iterator __begin, JEvt::iterator __end, const int N)
  {
    using namespace std;
    using namespace JPP;
 
    static const double MINIMAL_ANGLE_DEG  =  60.0;
 
    // estimate grid angle
 
    double angle_deg = 10.0 * sqrt((double) 240 / (double) N);
 
    if (angle_deg > MINIMAL_ANGLE_DEG) {
      angle_deg = MINIMAL_ANGLE_DEG;
    }
 
    const JOmega3D ps(angle_deg * PI/180.0);
 
    JEvt::iterator out = __begin;
 
    for (JOmega3D_t::const_iterator dir = ps.begin(); dir != ps.end() && out != __end; ++dir) {
 
      const JPointing pointing(*dir);
 
      sort(out, __end, pointing);
 
      JEvt::iterator p = out;
 
      for (JEvt::iterator i = out; i != __end && pointing.getAngle(*i) <= angle_deg; ++i) {
        if (qualitySorter(*i, *p)) {
          p = i;
        }
      }
 
      swap(*(out++), *p);
    }
 
    sort(__begin, out, qualitySorter);
 
    return out;
  }
 
 
  /**
   * 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