JNeutrinoInteractionCategory.hh

Go to the documentation of this file.

#ifndef __JAANET__JNEUTRINOINTERACTIONCATEGORY__
#define __JAANET__JNEUTRINOINTERACTIONCATEGORY__
 
#include <utility>
#include <iostream>
 
#include "km3net-dataformat/definitions/w2list_gseagen.hh"
 
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"
 
#include "JLang/JClonable.hh"
#include "JLang/JParameter.hh"
#include "JLang/JPredicate.hh"
 
#include "JTools/JRange.hh"
 
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JParticleTypes.hh"
#include "JAAnet/JGENIETypes.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JEvtCategory.hh"
 
 
/**
 * \file
 *
 * Classes and methods for defining neutrino interaction categories.
 * \author bjung
 */
namespace JAANET {}
namespace JPP { using namespace JAANET; }
 
namespace JAANET {
 
  using JLANG::JClonable;
  using JLANG::JParameter;
 
  using JTOOLS::JRange;
 
 
  static const char* const NEUTRAL_CURRENT = "NC"; //!< Neutral current header string identifier
  static const char* const CHARGED_CURRENT = "CC"; //!< Charged current header string identifier
 
  
  static const JRange<double> DEFAULT_NEUTRINO_ENERGY_RANGE              = JRange<double>(0, 1e11);   //!< Default neutrino energy range              [GeV]
  static const JRange<double> DEFAULT_NEUTRINO_COSINE_ZENITH_ANGLE_RANGE = JRange<double>(-1.0, 1.0); //!< Default neutrino cosine zenith angle range [-]
  
 
  /**
   * Auxiliary function to check if given PDG code corresponds to a neutrino.
   *
   * \param  type        PDG type
   * \return             true if given PDG type corresponds to a neutrino; else false
   */
  inline bool is_neutrino_primary(const int type)
  {
    const int absType = abs(type);
 
    return (absType == TRACK_TYPE_NUE  ||
            absType == TRACK_TYPE_NUMU ||
            absType == TRACK_TYPE_NUTAU);
  }
  
  
  /**
   * Auxiliary function to retrieve the GENIE interaction type corresponding to a given header.
   *
   * \param  header              MC header
   * \return                     GENIE interaction type
   */
  inline JInteractionTypeGENIE_t getInteractionType(const JHead& header)
  {
    using namespace std;
    using namespace JPP;
    
    if (is_gseagen(header)) {
      
      for (vector<JAANET::physics>::const_iterator i = header.physics.cbegin(); i != header.physics.cend(); ++i) {
        
        istringstream iss(i->buffer);
        
        for (string buffer; iss >> buffer ;) {
          
          if (buffer        == string(NEUTRAL_CURRENT)) {
            return JInteractionTypeGENIE_t::WEAK_NEUTRAL_CURRENT;
          } else if (buffer == string(CHARGED_CURRENT)) {
            return JInteractionTypeGENIE_t::WEAK_CHARGED_CURRENT;
          }
        }
      }
    }
 
    return JInteractionTypeGENIE_t::UNDEFINED;
  }
 
 
  /**
   * Auxiliary function to retrieve the target atomic number and mass number.
   *
   * \param  header              MC header
   * \return                     pair of atomic number (Z) and mass number (A)
   */
  inline std::pair<JParameter<int>,
                   JParameter<int> > getTargetInfo(const JHead& header)
  {
    using namespace std;
    using namespace JPP;
 
    JParameter<int> Z;
    JParameter<int> A;
    
    if (!header.target.buffer.empty()) {
      
      istringstream iss(header.target.buffer);
 
      if (! (iss.peek() == (int) 'A'  &&
             iss.ignore() && iss >> A &&
             iss.peek() == (int) 'Z'  &&
             iss.ignore() && iss >> Z &&
             Z > 0 && A >= Z) ) {
        THROW(JValueOutOfRange, "getTargetInfo(): Invalid target header-field: " << header.target.buffer << endl);
      }
    }
    
    return make_pair(Z,A);
  }
   
 
  /**
   * Class for neutrino interaction categories.
   */
  struct JNeutrinoInteractionCategory final :
    public JClonable<JEvtCategory, JNeutrinoInteractionCategory>
  {
    typedef JTOOLS::JRange<double>   JRange_t;
    
    
    /**
     * Default constructor.
     */
    JNeutrinoInteractionCategory() :
      interactionType(),
      scatteringType (),
      targetZ        (),
      targetA        (),
      Erange         (DEFAULT_NEUTRINO_ENERGY_RANGE),
      coszRange      (DEFAULT_NEUTRINO_COSINE_ZENITH_ANGLE_RANGE)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  neutrinoType       neutrino PDG type
     */
    JNeutrinoInteractionCategory(const int neutrinoType) :
      interactionType(),
      scatteringType (),
      targetZ        (),
      targetA        (),
      Erange         (DEFAULT_NEUTRINO_ENERGY_RANGE),
      coszRange      (DEFAULT_NEUTRINO_COSINE_ZENITH_ANGLE_RANGE)
    {
      this->primaryType = neutrinoType;
      
      check_validity();
    }
 
 
    /**
     * Constructor.
     *
     * \param  neutrinoType       neutrino PDG type
     * \param  Erange             neutrino energy range
     * \param  coszRange          zenith-angle    range
     */
    JNeutrinoInteractionCategory(const int       neutrinoType,
                                 const JRange_t& Erange,
                                 const JRange_t& coszRange) :
      interactionType(),
      scatteringType (),
      targetZ        (),
      targetA        (),
      Erange         (Erange),
      coszRange      (coszRange)
    {
      this->primaryType = neutrinoType;
      
      check_validity();
    }
    
 
    /**
     * Constructor.
     *
     * \param  neutrinoType       neutrino PDG type
     * \param  interactionType    interaction type
     * \param  scatteringType     scattering  type
     * \param  targetZ            target atomic number
     * \param  targetA            target mass   number
     * \param  Erange             neutrino energy range
     * \param  coszRange          zenith-angle    range
     */
    JNeutrinoInteractionCategory(const int                     neutrinoType,
                                 const JInteractionTypeGENIE_t interactionType,
                                 const JScatteringTypeGENIE_t  scatteringType,
                                 const int                     targetZ,
                                 const int                     targetA,
                                 const JRange_t&               Erange,
                                 const JRange_t&               coszRange) :
      interactionType((int) interactionType),
      scatteringType ((int) scatteringType),
      targetZ        (targetZ),
      targetA        (targetA),
      Erange         (Erange),
      coszRange      (coszRange)
    {
      this->primaryType = neutrinoType;      
      
      check_validity();
    }
 
 
    /**
     * Constructor.
     *
     * \param  header             MC header
     * \param  primaryType        primary PDG type
     */
    JNeutrinoInteractionCategory(const JHead& header,
                                 const int    primaryType) :
      JNeutrinoInteractionCategory()
    {
      using namespace std;
      using namespace JPP;
 
      if (is_gseagen(header)) {
 
        vector<JAANET::flux>::const_iterator i = find_if(header.flux.cbegin(), header.flux.cend(),
                                                         make_predicate(&JAANET::flux::type, primaryType));
          
        if (primaryType == header.primary.type || i != header.flux.cend()) {
          this->primaryType = primaryType;
        } else {
          THROW(JValueOutOfRange, "JNeutrinoInteractionCategory::JNeutrinoInteractionCategory(): Given header does not match with given primary type " << primaryType << ".");
        }
        
        const JInteractionTypeGENIE_t intType = getInteractionType(header);
 
        if (intType != JInteractionTypeGENIE_t::UNDEFINED) {
          interactionType = (int) intType;
        }
        
        const pair<JParameter<int>,
                   JParameter<int> > targetInfo = getTargetInfo(header);
        
        targetZ   = targetInfo.first;
        targetA   = targetInfo.second; 
 
        if (header.cut_nu.E.is_valid() &&
            header.cut_nu.E != cut().E) {
          Erange    = header.cut_nu.E;
        }
 
        if (header.cut_nu.cosT.is_valid() &&
            header.cut_nu.cosT != cut().cosT) {
          coszRange = header.cut_nu.cosT;
        }
      }
      
      check_validity();
    }
    
    
    /**
     * Constructor.
     *
     * N.B.: In this method, the neutrino PDG type is taken from `JAANET::JHead::primary::type`.
     *
     * \param  header             MC header
     */
    JNeutrinoInteractionCategory(const JHead& header) :
      JNeutrinoInteractionCategory(header, header.primary.type)
    {}
 
 
    /**
     * Check if neutrino interaction categories is valid.
     *
     * \return                    true if neutrino interaction category is valid; else false
     */
    inline bool is_valid() const override final
    {
      if (!is_neutrino_primary(this->primaryType)) {
        return false;
      }
 
      if (  interactionType.isDefined() &&
          !(interactionType.getValue() > (int) JInteractionTypeGENIE_t::INTERACTION_TYPES_START &&
            interactionType.getValue() < (int) JInteractionTypeGENIE_t::INTERACTION_TYPES_END) ) {
        return false;
      }
 
      if (   scatteringType.isDefined() &&
          !((scatteringType > (int) JScatteringTypeGENIE_t::NEUTRINO_SCATTERING_TYPES_START    &&
             scatteringType < (int) JScatteringTypeGENIE_t::NEUTRINO_SCATTERING_TYPES_END)     ||
            (scatteringType > (int) JScatteringTypeGENIE_t::DARK_MATTER_SCATTERING_TYPES_START &&
             scatteringType < (int) JScatteringTypeGENIE_t::DARK_MATTER_SCATTERING_TYPES_END))  ) {
        return false;
      }
 
      if (targetZ.isDefined() && targetZ.getValue() <= 0) {
        return false;
      }
 
      if (  targetA.isDefined() &&
          (!targetZ.isDefined() || targetA.getValue() < targetZ.getValue())) {
        return false;
      }
 
      return Erange.is_valid() && coszRange.is_valid();
    }
 
 
    /**
     * 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_gseagen(header)) {
        
        const JNeutrinoInteractionCategory cat(header);
 
        return (this->primaryType     == cat.primaryType     &&
                this->interactionType == cat.interactionType &&
                this->scatteringType  == cat.scatteringType  &&
                this->targetZ         == cat.targetZ         &&
                this->targetA         == cat.targetA         &&
                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_neutrino(primary)) { return false; }
        
      if (interactionType.isDefined()                      &&
          event.w2list.size() > W2LIST_GSEAGEN_CC          &&
          interactionType.getValue() != event.w2list[W2LIST_GSEAGEN_CC]) {
        return false;
      }
 
      if (scatteringType.isDefined()                       &&
          event.w2list.size() > W2LIST_GSEAGEN_ICHAN       &&
          scatteringType.getValue() != event.w2list[W2LIST_GSEAGEN_ICHAN]) {
        return false;
      }
 
      if (targetZ.isDefined()                              &&
          event.w2list.size() > W2LIST_GSEAGEN_TARGETZ     &&
          targetZ.getValue() != event.w2list[W2LIST_GSEAGEN_TARGETZ]) {
        return false;
      }
 
      if (targetA != 0                                     &&
          event.w2list.size() > W2LIST_GSEAGEN_TARGETA     &&
          targetA.getValue() != event.w2list[W2LIST_GSEAGEN_TARGETA]) {
        return false;
      }
 
      return Erange(primary.E) && coszRange(-primary.dir.z);
    }
 
 
    /**
     * Less-than method.
     *
     * Note: The energy and cosine zenith-angle 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(A, B)       \
      if (A != B) { return A < B; }
      
      RETURN_IF_DIFFERENT(this->primaryType, category.getPrimaryType());
      
      const JNeutrinoInteractionCategory* p = dynamic_cast<const JNeutrinoInteractionCategory*>(&category);
      
      if (p == NULL) { return false; }
      
      RETURN_IF_DIFFERENT(this->interactionType, p->interactionType);
      RETURN_IF_DIFFERENT(this->scatteringType,  p->scatteringType);
      RETURN_IF_DIFFERENT(this->targetZ,         p->targetZ);
      RETURN_IF_DIFFERENT(this->targetA,         p->targetA);
      
      if (p->Erange != this->Erange) {
        return (this->Erange.getLowerLimit() == p->Erange.getLowerLimit() ?
                this->Erange.getUpperLimit() <  p->Erange.getUpperLimit() :
                this->Erange.getLowerLimit() <  p->Erange.getLowerLimit());
      } else {
        return (this->coszRange.getLowerLimit() == p->coszRange.getLowerLimit() ?
                this->coszRange.getUpperLimit() <  p->coszRange.getUpperLimit() :
                this->coszRange.getLowerLimit() <  p->coszRange.getLowerLimit());
      }
#undef RETURN_IF_DIFFERENT            
    }
    
    
    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtCategory::getEquationParameters()) override final
    {
      return JNeutrinoInteractionCategoryHelper(*this, eqpars);
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtCategory::getEquationParameters()) const override final
    {
      return JNeutrinoInteractionCategoryHelper(*this, eqpars);
    }
 
    
  private:
 
    JParameter<int> interactionType;   //!< interaction  type
    JParameter<int> scatteringType;    //!< scattering   type
    JParameter<int> targetZ;           //!< target atomic number
    JParameter<int> targetA;           //!< target mass   number
    JRange_t        Erange;            //!< neutrino energy range
    JRange_t        coszRange;         //!< cosine   zenith range
    
    
    /**
     * Auxiliary class for I/O of neutrino interaction category.
     */
    struct JNeutrinoInteractionCategoryHelper :
      public JProperties
    {
      /**
       * Constructor.
       *
       * \param  cat           neutrino interaction category
       * \param  eqpars        equation parameters 
       */
      template<class JNeutrinoInteractionCategory_t>
      JNeutrinoInteractionCategoryHelper(JNeutrinoInteractionCategory_t& cat,
                                         const JEquationParameters&      eqpars) :
        JProperties(eqpars, 1)
      {
        (*this)["neutrinoType"] = cat.primaryType;
        
        this->insert(gmake_property(cat.interactionType));
        this->insert(gmake_property(cat.scatteringType));
        this->insert(gmake_property(cat.targetZ));
        this->insert(gmake_property(cat.targetA));
        this->insert(gmake_property(cat.Erange));
        this->insert(gmake_property(cat.coszRange));
      }
    };
  };
}
 
#endif