JMuonBundleCategory.hh

Go to the documentation of this file.

#ifndef __JAANET__JMUONBUNDLECATEGORY__
#define __JAANET__JMUONBUNDLECATEGORY__
 
#include <iostream>
 
#include "km3net-dataformat/definitions/trkmembers.hh"
 
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"
 
#include "JLang/JClonable.hh"
 
#include "JTools/JRange.hh"
 
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JTrack3E.hh"
 
#include "JMath/JMathToolkit.hh"
 
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JEvtCategory.hh"
#include "JAAnet/JAAnetToolkit.hh"
 
 
/**
 * \file
 *
 * Classes and methods for defining muon bundle categories.
 * \author bjung
 */
namespace JAANET {}
namespace JPP { using namespace JAANET; }
 
namespace JAANET {
 
  using JLANG::JClonable;
  
  using JTOOLS::JRange;
  
 
  static const JRange<double> DEFAULT_BUNDLE_ENERGY_RANGE              = JRange<double>(0, 1e11);  //!< Default bundle energy range              [GeV]
  static const JRange<double> DEFAULT_BUNDLE_COSINE_ZENITH_ANGLE_RANGE = JRange<double>(0.0, 1.0); //!< Default bundle cosine zenith angle range [-]
  static const JRange<int>    DEFAULT_BUNDLE_MULTIPLICITY_RANGE        = JRange<int>   (0, 1000);  //!< Default bundle multiplicity range        [-]
  static const double         DEFAULT_BUNDLE_RADIAL_EXTENT             = 1e5;                      //!< Default bundle radial extent             [m]
 
 
  /**
   * Auxiliary function to check if given PDG code corresponds to a valid muon bundle primary type.
   *
   * \param  type        PDG type
   * \return             true if given PDG type corresponds to a valid muon bundle primary type; else false
   */
  inline bool is_muon_bundle_primary(const int type)
  {
    using namespace std;
    
    const int    absType = abs(type);
    const string typeStr = to_string(absType);
      
    return ( (absType == PDG_MUONBUNDLE         ||                    // MUPAGE
              absType == TRACK_TYPE_MUON)                          ||
             (typeStr.size() == 10 && typeStr.substr(0,2) == "10") || // CORSIKA nuclei
             (absType != TRACK_TYPE_PHOTON      &&                    // CORSIKA single baryon primaries (e.g. proton or neutron)
              absType != TRACK_TYPE_ELECTRON    &&
              absType != TRACK_TYPE_MUON        &&
              absType != TRACK_TYPE_TAU         &&
              absType != TRACK_TYPE_NUE         &&
              absType != TRACK_TYPE_NUMU        &&
              absType != TRACK_TYPE_NUTAU       &&
              ((int)(absType / 1000)) % 10 != 0) );
  }
  
  
  /**
   * Class for muon bundle categories.
   */
  struct JMuonBundleCategory final :
    public JClonable<JEvtCategory, JMuonBundleCategory>
  {
    /**
     * Default constructor.
     */
    JMuonBundleCategory() :
      Erange     (DEFAULT_BUNDLE_ENERGY_RANGE),
      coszRange  (DEFAULT_BUNDLE_COSINE_ZENITH_ANGLE_RANGE),
      Mrange     (DEFAULT_BUNDLE_MULTIPLICITY_RANGE),
      Rmax       (DEFAULT_BUNDLE_RADIAL_EXTENT)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  primaryType     cosmic-ray primary PDG identifier
     */
    JMuonBundleCategory(const int primaryType) :
      Erange     (DEFAULT_BUNDLE_ENERGY_RANGE),
      coszRange  (DEFAULT_BUNDLE_COSINE_ZENITH_ANGLE_RANGE),
      Mrange     (DEFAULT_BUNDLE_MULTIPLICITY_RANGE),
      Rmax       (DEFAULT_BUNDLE_RADIAL_EXTENT)
    {
      this->primaryType = primaryType;
      
      check_validity();
    }
    
 
    /**
     * Constructor.
     *
     * \param  primaryType     cosmic-ray primary PDG identifier
     * \param  Erange          energy              range
     * \param  coszRange       cosine zenith-angle range
     * \param  Mrange          multiplicity        range
     * \param  Rmax            maximum radial extent
     */
    JMuonBundleCategory(const int             primaryType,
                        const JRange<double>& Erange,
                        const JRange<double>& coszRange,
                        const JRange<int>&    Mrange,
                        const double          Rmax) :
      Erange     (Erange),
      coszRange  (coszRange),
      Mrange     (Mrange),
      Rmax       (Rmax)
    {
      this->primaryType = primaryType;
      
      check_validity();
    }
    
    
    /**
     * Constructor.
     *
     * \param  header             MC header
     */
    JMuonBundleCategory(const JHead& header) :
      JMuonBundleCategory()
    {
      using namespace std;
      using namespace JPP;
 
      if        (is_mupage(header)) {
        this->primaryType = PDG_MUONBUNDLE;
        this->Erange      = header.cut_in.E;
        this->coszRange   = header.cut_in.cosT;
      } else if (is_corsika(header)) {
        this->primaryType = header.primary.type;
        this->Erange      = header.cut_primary.E;
        this->coszRange   = header.cut_primary.cosT;
      }
 
      check_validity();
    }
 
 
    /**
     * Check if muon bundle category is valid.
     *
     * \return                    true if muon bundle category is valid; else false
     */
    inline bool is_valid() const override final
    {
      return (is_muon_bundle_primary(this->primaryType) &&
              Erange   .is_valid()                      &&
              coszRange.is_valid()                      &&
              Mrange.getLowerLimit() >= 0               &&
              Mrange   .is_valid()                      &&
              Rmax > 0);
    }
 
 
    /**
     * Check whether this event category matches with the given MC header.
     *
     * \param  header            header
     * \return                   true if matching; else false
     */
    bool match(const JHead& header) const override final
    {
      if (is_mupage(header) || is_corsika(header)) {
        
        const JMuonBundleCategory cat(header);
 
        return (this->primaryType == cat.getPrimaryType() &&
                this->Erange   .overlap (cat.Erange)      &&
                this->coszRange.overlap (cat.coszRange));
        
      } else {
        
        return false;
      }
    }
    
 
    /**
     * Check whether given event matches with this event category.
     *
     * \param  event             event
     * \return                   true if matching; else false
     */
    bool match(const Evt& event) const override final
    {
      using namespace std;
      using namespace JPP;
 
      const Trk& primary = get_primary(event);
 
      if (!(is_muonbundle(primary) || is_muon(primary)) ||
          !Erange   ( primary.E)                        ||
          !coszRange(-primary.dir.z)) { return false; }
 
      if (Mrange == DEFAULT_BUNDLE_MULTIPLICITY_RANGE &&
          Rmax   == DEFAULT_BUNDLE_RADIAL_EXTENT) { return true; }
 
      const JTrack3E&     tp = getTrack(primary);
      const JDirection3D& dp = tp.getDirection();
      const JPosition3D&  pp = tp.getPosition();
        
      int    multiplicity = 0;
      double maxR         = 0.0;
        
      for (vector<Trk>::const_iterator i = event.mc_trks.cbegin(); i != event.mc_trks.cend(); ++i) {
        if (is_finalstate(*i) && is_muon(*i)) {
 
          const JTrack3E&    ti = getTrack(*i);
          const JPosition3D& pi = ti.getPosition();
            
          const double R = getCross<JVector3D>(dp, pi - pp).getLength();
            
          if (R > maxR) { maxR = R; }
            
          ++multiplicity;
        }
      }
 
      return (Mrange(multiplicity) && maxR < Rmax);
    }
 
 
    /**
     * Less-than method.
     *
     * Note: The energy, cosine zenith-angle and multiplicity ranges are compared based on the lower limit first.\n
     *       If the lower limits are equal, the upper limits are compared.
     *
     * \param  category           event category
     * \return                    true if first this event category less than given event category; else false
     */
    inline bool less(const JEvtCategory& category) const override final
    {
#define RETURN_IF_DIFFERENT_RANGE(A, B)     \
      if (A != B) { return (A.getLowerLimit() == B.getLowerLimit() ?  \
                            A.getUpperLimit() <  B.getUpperLimit() :  \
                            A.getLowerLimit() <  B.getLowerLimit() ); }
      
      if (this->primaryType != category.getPrimaryType()) {
        return this->primaryType < category.getPrimaryType();
      }
      
      const JMuonBundleCategory* p = dynamic_cast<const JMuonBundleCategory*>(&category);
            
      if (p == NULL) { return false; }
 
      RETURN_IF_DIFFERENT_RANGE(this->Erange,    p->Erange);
      RETURN_IF_DIFFERENT_RANGE(this->coszRange, p->coszRange);
      RETURN_IF_DIFFERENT_RANGE(this->Mrange,    p->Mrange);
 
      return this->Rmax < p->Rmax;
      
#undef RETURN_IF_DIFFERENT_RANGE
    }
    
    
    /**
     * Get properties of this class.
     *
     * \param  equation        equation parameters
     */
    JProperties getProperties(const JEquationParameters& equation = JEvtCategory::getEquationParameters()) override final
    {
      return JMuonBundleCategoryHelper(*this, equation);
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  equation        equation parameters
     */
    JProperties getProperties(const JEquationParameters& equation = JEvtCategory::getEquationParameters()) const override final
    {
      return JMuonBundleCategoryHelper(*this, equation);
    }
 
 
  private:
 
    JRange<double> Erange;      //!< energy              range [GeV]
    JRange<double> coszRange;   //!< cosine zenith-angle range [-]
    JRange<int>    Mrange;      //!< multiplicity        range [-]
    double         Rmax;        //!< maximum radial extent     [m]
 
    
    /**
     * Auxiliary class for I/O of muon bundle category.
     */
    struct JMuonBundleCategoryHelper :
      public JProperties
    {
      /**
       * Constructor.
       *
       * \param  cat           muon bundle category
       * \param  equation      equation parameters 
       */
      template<class JMuonBundleCategory_t>
      JMuonBundleCategoryHelper(JMuonBundleCategory_t&     cat,
                                const JEquationParameters& equation) :
        JProperties(equation, 1)
      {
        this->insert(gmake_property(cat.primaryType));
        this->insert(gmake_property(cat.Erange));
        this->insert(gmake_property(cat.coszRange));
        this->insert(gmake_property(cat.Mrange));
        this->insert(gmake_property(cat.Rmax));
      }
    };
  };
}
 
#endif