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 "flux/Flux.hh"

#include "JLang/JEquals.hh"
#include "JLang/JClonable.hh"

#include "JAAnet/JFlux.hh"
#include "JAAnet/JDiffuseFlux.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JParticleTypes.hh"
#include "JAAnet/JEvtWeightFactorHelper.hh"

#include "JOscProb/JOscChannel.hh"
#include "JOscProb/JOscProbHelper.hh"
#include "JOscProb/JOscProbInterface.hh"


/**
 * \author bjung
 */

namespace JAANET {}
namespace JPP { using namespace JAANET; }

namespace JAANET {

  using JLANG::JEquals;
  using JLANG::JClonable;

  using JOSCPROB::JOscProbHelper;
  using JOSCPROB::JOscProbInterface;
  

  /**
   * Implementation of oscillated neutrino flux.
   */ 
  struct JOscFlux :
    public JEquals<JOscFlux>,
    public JClonable<JFlux, JOscFlux>,
    public JDiffuseFluxHelper,
    public JOscProbHelper
  {
    /**
     * Default constructor.
     */
    JOscFlux()
    {}


    /**
     * Constructor.
     *
     * \param  diffuseFlux       diffuse flux object
     * \param  oscProb           oscillation probability interface
     */
    JOscFlux(const JDiffuseFlux&      diffuseFlux,
             const JOscProbInterface& oscProb)
    {
      JDiffuseFluxHelper::configure(diffuseFlux);
      JOscProbHelper::configure(oscProb);
    }
    
    
    /**
     * 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 operator()(const Evt& evt) const
    {
      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 == 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 += ( JDiffuseFluxHelper::getFactor(inType, log10(neutrino.E), costh) *
                      JOscProbHelper::getP(channel, neutrino.E, costh) );
          }
        }
        
      } else if (interactionType == WEAK_NEUTRAL_CURRENT) { // For NC events, the neutrino flavour is irrelevant

        const int Cparity = (int) getChargeParity(neutrino);

        flux += (JDiffuseFluxHelper::getFactor(Cparity * TRACK_TYPE_NUE,   log10(neutrino.E), costh) +
                 JDiffuseFluxHelper::getFactor(Cparity * TRACK_TYPE_NUMU,  log10(neutrino.E), costh) +
                 JDiffuseFluxHelper::getFactor(Cparity * TRACK_TYPE_NUTAU, log10(neutrino.E), costh));
      }
      
      return flux;
    }


    /**
     * Check if this flux is equal to given flux.
     *
     * \param  object            flux object
     * \return                   true if this flux is identical to given flux; else false
     */
    bool equals(const JOscFlux& object) const
    {
      return (static_cast<const JDiffuseFluxHelper&>(object) == static_cast<const JDiffuseFluxHelper&>(*this) &&
              static_cast<const JOscProbHelper&>(object)     == static_cast<const JOscProbHelper&>    (*this));
    }


    /**
     * Get event-weight factor for given event.
     *
     * \param  evt               event
     * \return                   event-weight factor [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFactor(const Evt& evt) const override
    {
      return (*this)(evt);
    }
  };
}

#endif