NeutrinoSource.hh

//
// This file is part of MARLEY (Model of Argon Reaction Low Energy Yields)
//
// MARLEY is free software: you can redistribute it and/or modify it under the
// terms of version 3 of the GNU General Public License as published by the
// Free Software Foundation.
//
// For the full text of the license please see COPYING or
// visit http://opensource.org/licenses/GPL-3.0
//
// Please respect the MCnet academic usage guidelines. See GUIDELINES
// or visit https://www.montecarlonet.org/GUIDELINES for details.
 
#pragma once
#include <functional>
#include <set>
 
#include "marley/marley_utils.hh"
#include "marley/InterpolationGrid.hh"
 
namespace marley {
 
  class Generator;

  class NeutrinoSource {
    public:

      NeutrinoSource(int particle_id);
 
      virtual ~NeutrinoSource() = default;

      virtual double get_Emax() const = 0;

      virtual double get_Emin() const = 0;

      inline virtual int get_pid() const;

      virtual double pdf(double E) const = 0;

      static inline bool pdg_is_allowed(const int pdg);

      virtual double sample_incident_neutrino(int& pdg,
        marley::Generator& gen) const;
 
    protected:
 
      int pid_; 
 
    private:
      static const std::set<int> pids_;
  };

  class MonoNeutrinoSource : public NeutrinoSource {
    public:
      inline MonoNeutrinoSource(int particle_id =
        marley_utils::ELECTRON_NEUTRINO, double E = 10.);
 
      inline virtual double get_Emax() const override;
 
      inline virtual double get_Emin() const override;
 
      inline virtual double pdf(double E) const override;
 
    protected:
      double energy_; 
  };

  class FermiDiracNeutrinoSource : public NeutrinoSource {
    public:

      FermiDiracNeutrinoSource(int particle_id
        = marley_utils::ELECTRON_NEUTRINO, double Emin = 0.,
        double Emax = 50., double temp = 3.5, double eta = 0.);
 
      inline virtual double get_Emin() const override;
 
      inline virtual double get_Emax() const override;
 
      virtual double pdf(double E) const override;
 
    protected:
 
      double Emin_; 
      double Emax_; 
      double temperature_;  
      double eta_; 
      double C_; 
  };

  class AlphaFitNeutrinoSource : public NeutrinoSource {
    public:
      AlphaFitNeutrinoSource(int particle_id
        = marley_utils::ELECTRON_NEUTRINO, double Emin = 0.,
        double Emax = 50., double Emean = 13., double alpha = 2.);
 
      inline virtual double get_Emin() const override;
      inline virtual double get_Emax() const override;
 
      virtual double pdf(double E) const override;
 
    protected:
      double Emin_; 
      double Emax_; 
      double Emean_;
      double alpha_; 
      double C_; 
  };

  class BetaFitNeutrinoSource : public NeutrinoSource {
    public:
      BetaFitNeutrinoSource(int particle_id
        = marley_utils::ELECTRON_NEUTRINO, double Emin = 0.,
        double Emax = 50., double Emean = 13., double beta = 4.5);
 
      inline virtual double get_Emin() const override;
      inline virtual double get_Emax() const override;
 
      virtual double pdf(double E) const override;
 
    protected:
      double Emin_; 
      double Emax_; 
      double Emean_;
      double alpha_; 
      double C_; 
  };

  class FunctionNeutrinoSource : public NeutrinoSource {
    public:
      FunctionNeutrinoSource(int particle_id
        = marley_utils::ELECTRON_NEUTRINO, double Emin = 0., double Emax = 50.,
        std::function<double(double)> prob_dens_func
        = [](double) -> double { return 1.; });
 
      inline virtual double get_Emax() const override;
 
      inline virtual double get_Emin() const override;
 
      inline virtual double pdf(double E) const override;
 
    private:
      double Emin_; 
      double Emax_; 
      std::function<double(double)> probability_density_;
  };

  class DecayAtRestNeutrinoSource : public NeutrinoSource {
    public:
      DecayAtRestNeutrinoSource(int particle_id
        = marley_utils::ELECTRON_NEUTRINO);
 
      inline virtual double get_Emax() const override;
 
      inline virtual double get_Emin() const override;
 
      virtual double pdf(double E) const override;
 
    private:
      // Muon mass stuff (m_mu^(-4) pre-computed for speed)
      static constexpr double m_mu_ = marley_utils::m_mu
        * marley_utils::micro_amu; // MeV
      static constexpr double m_mu_to_the_minus_four_
        = 1. / (m_mu_ * m_mu_ * m_mu_ * m_mu_);
      static constexpr double Emin_ = 0.; // MeV
      static constexpr double Emax_ = m_mu_ / 2.; // MeV
  };

  class GridNeutrinoSource : public NeutrinoSource {
    public:
      using Grid = InterpolationGrid<double>;
      using Method = Grid::InterpolationMethod;

      inline GridNeutrinoSource(const Grid& g, int particle_id
        = marley_utils::ELECTRON_NEUTRINO);

      inline GridNeutrinoSource(const std::vector<double>& Es,
        const std::vector<double>& PDs, int particle_id
        = marley_utils::ELECTRON_NEUTRINO, Method method
        = Method::LinearLinear);
 
      inline virtual double get_Emax() const override;
 
      inline virtual double get_Emin() const override;
 
      inline virtual double pdf(double E) const override;
 
    protected:
      Grid grid_;
 
    private:
      void check_for_errors();
  };
 
  // Inline function definitions
  inline int NeutrinoSource::get_pid() const { return pid_; }
  inline bool NeutrinoSource::pdg_is_allowed(const int pdg)
    { return (pids_.count(pdg) > 0); }
 
  inline MonoNeutrinoSource::MonoNeutrinoSource(int particle_id, double E)
    : NeutrinoSource(particle_id), energy_(E) {}
  inline double MonoNeutrinoSource::get_Emax() const { return energy_; }
  inline double MonoNeutrinoSource::get_Emin() const { return energy_; }
  inline double MonoNeutrinoSource::pdf(double E) const
    { if (energy_ == E) return 1.; else return 0.; }
 
  inline double FermiDiracNeutrinoSource::get_Emax() const { return Emax_; }
  inline double FermiDiracNeutrinoSource::get_Emin() const { return Emin_; }
 
  inline double AlphaFitNeutrinoSource::get_Emax() const { return Emax_; }
  inline double AlphaFitNeutrinoSource::get_Emin() const { return Emin_; }
 
  inline double BetaFitNeutrinoSource::get_Emax() const { return Emax_; }
  inline double BetaFitNeutrinoSource::get_Emin() const { return Emin_; }
 
  inline double FunctionNeutrinoSource::get_Emax() const { return Emax_; }
  inline double FunctionNeutrinoSource::get_Emin() const { return Emin_; }
  inline double FunctionNeutrinoSource::pdf(double E) const {
    if (E < Emin_ || E > Emax_) return 0.;
    else return probability_density_(E);
  }
 
  inline double DecayAtRestNeutrinoSource::get_Emax() const { return Emax_; }
  inline double DecayAtRestNeutrinoSource::get_Emin() const { return Emin_; }
 
  inline double GridNeutrinoSource::get_Emax() const
    { return grid_.back().first; }
  inline double GridNeutrinoSource::get_Emin() const
    { return grid_.front().first; }
  inline double GridNeutrinoSource::pdf(double E) const
    { return grid_.interpolate(E); }
 
  inline GridNeutrinoSource::GridNeutrinoSource(const Grid& g, int particle_id)
    : NeutrinoSource(particle_id), grid_(g) { check_for_errors(); }
 
  inline GridNeutrinoSource::GridNeutrinoSource(const std::vector<double>& Es,
    const std::vector<double>& prob_densities, int particle_id, Method method)
    : NeutrinoSource(particle_id), grid_(Es, prob_densities, method)
    { check_for_errors(); }
}