JOscFlux.hh

Go to the documentation of this file.

#ifndef __JAANET__JOSCFLUX__
#define __JAANET__JOSCFLUX__
 
#include "km3net-dataformat/definitions/w2list_gseagen.hh"
 
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"
 
#include "JLang/JClonable.hh"
 
#include "JAAnet/JGENIETypes.hh"
#include "JAAnet/JParticleTypes.hh"
#include "JAAnet/JAAnetToolkit.hh"
 
#include "JOscProb/JOscChannel.hh"
#include "JOscProb/JOscProbHelper.hh"
#include "JOscProb/JOscProbInterface.hh"
 
#include "JAAnet/JFlux.hh"
#include "JAAnet/JDiffuseFlux.hh"
#include "JAAnet/JDiffuseFluxHelper.hh"
 
 
/**
 * \author bjung
 */
 
namespace JAANET {}
namespace JPP { using namespace JAANET; }
 
namespace JAANET {
 
  using JLANG::JClonable;
  
  using JOSCPROB::JOscProbHelper;
  using JOSCPROB::JOscProbInterface;
  
 
  /**
   * Implementation of oscillated neutrino flux.
   */ 
  struct JOscFlux :
    public JClonable<JFlux, JOscFlux>
  {
    /**
     * Constructor.
     *
     * \param  diffuseFlux       diffuse flux function object
     * \param  oscProb           oscillation probability function object
     */
    JOscFlux(const JDiffuseFluxHelper& diffuseFlux,
             const JOscProbHelper&     oscProb) :
      P(oscProb),
      F(diffuseFlux)
    {
      check_validity();
    }
 
 
    /**
     * Get oscillation probability calculator.
     *
     * \return                   oscillation probability calculator
     */
    const JOscProbHelper& getOscProb() const
    {
      return P;
    }
 
 
    /**
     * Set oscillation probability calculator.
     *
     * \param  oscProb         oscillation probability calculator.
     */
    void setOscProb(const JOscProbHelper oscProb)
    {
      using namespace JPP;
      
      if (oscProb) {
        this->P = oscProb;
      } else {
        THROW(JNullPointerException, "JOscFlux::setOscProb(): Given oscillation probability calculator is invalid.");
      }
    }
 
 
    /**
     * Get diffuse flux function.
     *
     * \return                 diffuse flux function
     */
    const JDiffuseFluxHelper& getDiffuseFlux() const
    {
      return F;
    }
 
 
    /**
     * Set diffuse flux function.
     *
     * \param  flux            diffuse flux function
     */
    void setDiffuseFlux(const JDiffuseFluxHelper flux)
    {
      using namespace std;
      using namespace JPP;
      
      if (flux && flux->is_valid()) {
        this->F = flux;
      } else if (flux) {
        THROW(JValueOutOfRange,      "JOscFlux::setFlux(): Given flux function is invalid: " << endl << *flux);
      } else {
        THROW(JNullPointerException, "JOscFlux::setFlux(): Given flux function is invalid.");   
      }
    }
    
    
    /**
     * Check whether this oscillated neutrino flux object is valid.
     *
     * \return                 true if valid; else false
     */
    bool is_valid() const override final
    {
      return (P && F && F->is_valid());
    }
    
    
    /**
     * Get flux for given event.
     *
     * Note that in this evaluation the zenith-angle is defined\n
     * with respect to the line of sight (i.e. a neutrino pointing straight at you\n
     * from the center of the Earth has \f$ cos(\theta) = -1.0 \f$).
     *
     * \param  evt               event
     * \return                   flux  \f$ \left[\mathrm{GeV}^{-1} \, \mathrm{m}^{-2} \, \mathrm{sr}^{-1} \, \mathrm{s}^{-1}\right] \f$
     */
    double getFactor(const Evt& evt) const override final    
    {
      using namespace JPP;
 
      double flux = 0.0;
 
      const Trk& neutrino = get_neutrino(evt);
      const double costh  = -neutrino.dir.z;
 
      const int interactionType = evt.w2list[W2LIST_GSEAGEN_CC];
 
      if        (interactionType == (int) JInteractionTypeGENIE_t::WEAK_CHARGED_CURRENT) {
 
        for (int i = 0; i != NUMBER_OF_OSCCHANNELS; ++i) {
 
          const JOscChannel& channel = getOscChannel[i];
 
          const int inType  = ((int) channel.Cparity) * ((int) channel.in);
          const int outType = ((int) channel.Cparity) * ((int) channel.out);
        
          if (outType == neutrino.type) {
          
            flux += ( F.getFlux(inType,  log10(neutrino.E), costh) *
                      P.getP   (channel, neutrino.E,        costh) );
          }
        }
        
      } else if (interactionType == (int) JInteractionTypeGENIE_t::WEAK_NEUTRAL_CURRENT) { // For NC events, the neutrino flavour is irrelevant
 
        const int Cparity = (int) getChargeParity(neutrino);
 
        flux += (F.getFlux(Cparity * TRACK_TYPE_NUE,   log10(neutrino.E), costh) +
                 F.getFlux(Cparity * TRACK_TYPE_NUMU,  log10(neutrino.E), costh) +
                 F.getFlux(Cparity * TRACK_TYPE_NUTAU, log10(neutrino.E), costh));
      }
      
      return flux;
    }
    
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) override final
    {
      return JOscFluxHelper(*this, eqpars);
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) const override final
    {
      return JOscFluxHelper(*this, eqpars);      
    }
 
    
  private:
 
    /**
     * Auxiliary class for I/O of oscillated flux.
     */
    struct JOscFluxHelper :
      public JProperties
    {
      /**
       * Constructor.
       *
       * \param  oscflux       oscillated flux
       * \param  eqpars        equation parameters 
       */
      template<class JOscFlux_t>
      JOscFluxHelper(JOscFlux_t&                oscflux,
                     const JEquationParameters& eqpars) :
        JProperties(eqpars, 1)
      {
        (*this)[JEvtWeightFactor::getTypeKey()] = "oscillated neutrino flux";   
        
        if (oscflux.F) {
          
          JProperties properties = oscflux.F.getProperties(eqpars);
          
          (*this)["flux"] = properties;
        }
        
        if (oscflux.P) {
 
          JProperties properties = oscflux.P.getProperties(eqpars);
          
          (*this)["oscprob"] = properties;
        }
      }
    };
 
    
    JOscProbHelper     P; //!< Oscillation probability calculator
    JDiffuseFluxHelper F; //!< Diffuse flux function
  };
}
 
#endif