#ifndef __JAANET__JOSCFLUX__ #define __JAANET__JOSCFLUX__ #include "km3net-dataformat/offline/Evt.hh" #include "km3net-dataformat/offline/Trk.hh" #include "flux/Flux.hh" #include "JLang/JClonable.hh" #include "JAAnet/JFlux.hh" #include "JAAnet/JDiffuseFlux.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JAAnet/JEvtWeightFactorHelper.hh" #include "JOscProb/JOscProb.hh" #include "JOscProb/JOscChannel.hh" #include "JOscProb/JOscProbHelper.hh" /** * \author bjung */ namespace JAANET {} namespace JPP { using namespace JAANET; } namespace JAANET { using JLANG::JClonable; using JOSCPROB::JOscProb; using JOSCPROB::JOscProbHelper; /** * Implementation of oscillated neutrino flux. */ struct JOscFlux : public JClonable, public JDiffuseFluxHelper, public JOscProbHelper { /** * Default constructor. */ JOscFlux() {} /** * Constructor. * * \param diffuseFlux diffuse flux object * \param oscProb oscillation probability interface */ JOscFlux(const JDiffuseFlux& diffuseFlux, const JOscProb& 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 / neutrino.dir.len(); 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::getOscProb(channel, neutrino.E, costh) ); } } return flux; } /** * 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