ExitChannel.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
 
// standard library includes
#include <memory>
#include <vector>
 
// MARLEY includes
#include "marley/ChebyshevInterpolatingFunction.hh"
#include "marley/Fragment.hh"
#include "marley/Generator.hh"
#include "marley/IteratorToPointerMember.hh"
#include "marley/Level.hh"
#include "marley/MassTable.hh"
#include "marley/Parity.hh"
#include "marley/marley_utils.hh"
 
namespace marley {
 
  class ExitChannel {
 
    public:
 
      ExitChannel() = default;
 
      ExitChannel(int pdgi, int qi, double Exi, int twoJi, marley::Parity Pi,
        double rho_i, marley::StructureDatabase& sdb ) : pdgi_( pdgi ),
        qi_( qi ), Exi_( Exi ), twoJi_( twoJi ), Pi_( Pi ), sdb_( &sdb )
      {
        one_over_two_pi_rho_i_ = std::pow( 2. * marley_utils::pi * rho_i, -1 );
      }
 
      virtual ~ExitChannel() = default;
 
      virtual bool is_continuum() const = 0;
 
      virtual bool emits_fragment() const = 0;
 
      virtual void do_decay(double& Exf, int& two_Jf,
        marley::Parity& Pf, const marley::Particle& compound_nucleus,
        marley::Particle& emitted_particle, marley::Particle& residual_nucleus,
        marley::Generator& gen) const = 0;
 
      template<typename It> static inline
        marley::IteratorToPointerMember<It, double> make_width_iterator(It it)
      {
        return marley::IteratorToPointerMember<It,
          double>(it, &marley::ExitChannel::width_);
      }
 
      inline double width() const { return width_; }
 
      virtual int emitted_particle_pdg() const = 0;
 
      virtual int final_nucleus_pdg() const = 0;
 
    protected:
 
      virtual void compute_total_width() = 0;
 
      virtual void prepare_products(const marley::Particle& compound_nucleus,
        marley::Particle& emitted_particle, marley::Particle& residual_nucleus,
        double Exf, marley::Generator& gen) const;
 
      int pdgi_;
 
      int qi_;
 
      double Exi_;
 
      int twoJi_;
 
      marley::Parity Pi_;
 
      double one_over_two_pi_rho_i_;
 
      double width_;
 
      marley::StructureDatabase* sdb_;
  };
 
  class DiscreteExitChannel : virtual public ExitChannel {
    public:
 
      DiscreteExitChannel( const marley::Level& flev ) : final_level_( flev ) {}
 
      virtual void do_decay(double& Ex, int& two_J,
        marley::Parity& Pi, const marley::Particle& compound_nucleus,
        marley::Particle& emitted_particle, marley::Particle& residual_nucleus,
        marley::Generator& gen) const final override;
 
      inline virtual bool is_continuum() const final override { return false; }
 
      inline const marley::Level& get_final_level() const
        { return final_level_; }
 
    protected:
 
      const marley::Level& final_level_;
  };
 
  class FragmentExitChannel : virtual public ExitChannel {
    public:
 
      FragmentExitChannel(const marley::Fragment& fragment)
        : fragment_pdg_( fragment.get_pid() ) {}
 
      virtual int emitted_particle_pdg() const final override
        { return fragment_pdg_; }
 
      virtual bool emits_fragment() const final override
        { return true; }
 
      virtual int final_nucleus_pdg() const final override;
 
    protected:
 
      double max_Exf() const;
 
      int fragment_pdg_;
  };
 
  class GammaExitChannel : virtual public ExitChannel {
    public:
 
      GammaExitChannel() {}
 
      virtual int emitted_particle_pdg() const final override
        { return marley_utils::PHOTON; }
 
      virtual bool emits_fragment() const final override
        { return false; }
 
      inline virtual int final_nucleus_pdg() const final override
        { return pdgi_; }
 
    protected:
 
      // Returns the gamma-ray energy corresponding to a particular
      // final nuclear excitation energy
      // @param Exf Final nuclear excitation energy (MeV)
      // @return Energy of the gamma-ray emitted in this exit channel (MeV)
      double gamma_energy( double Exf ) const;
 
      marley::GammaStrengthFunctionModel::TransitionType get_transition_type(
        int mpol, marley::Parity Pf ) const;
  };
 
  class ContinuumExitChannel : virtual public ExitChannel
  {
    public:
 
      ContinuumExitChannel(double Ec_min, int lmax) : E_c_min_( Ec_min ),
        l_max_( lmax ) {}
 
      virtual void compute_total_width() final override;
 
      virtual void do_decay(double& Ex, int& two_J,
        marley::Parity& Pi, const marley::Particle& compound_nucleus,
        marley::Particle& emitted_particle, marley::Particle& residual_nucleus,
        marley::Generator& gen) const final override;
 
      virtual double differential_width( double Exf,
        bool store_jpi_widths = false ) const = 0;
 
      inline virtual bool is_continuum() const final override { return true; }
 
      inline void set_skip_jpi_sampling(bool skip_it) const
        { skip_jpi_sampling_ = skip_it; }
 
      struct SpinParityWidth {
 
        SpinParityWidth(int twoJ, marley::Parity p, double w)
          : twoJf(twoJ), Pf(p), diff_width(w) {}
 
        int twoJf; 
        marley::Parity Pf; 
        double diff_width; 
      };
 
      double sample_Exf(marley::Generator& gen) const;
 
      void sample_spin_parity(double Exf, int& two_Jf, marley::Parity& Pf,
        marley::Generator& gen) const;
 
      inline double E_c_min() const { return E_c_min_; }
 
      virtual double E_c_max() const = 0;
 
    protected:
 
      double E_c_min_;
 
      int l_max_;
 
      mutable std::vector<SpinParityWidth> jpi_widths_table_;
 
      mutable bool skip_jpi_sampling_ = false;
 
      mutable std::unique_ptr<marley::ChebyshevInterpolatingFunction> Exf_cdf_;
  };
 
  class FragmentDiscreteExitChannel : public DiscreteExitChannel,
    public FragmentExitChannel
  {
    public:
 
      FragmentDiscreteExitChannel(int pdgi, int qi, double Exi, int twoJi,
        marley::Parity Pi, double rho_i, marley::StructureDatabase& sdb,
        const marley::Level& flev, const marley::Fragment& frag)
        : ExitChannel( pdgi, qi, Exi, twoJi, Pi, rho_i, sdb ),
        DiscreteExitChannel( flev ), FragmentExitChannel( frag )
      {
        this->compute_total_width();
      }
 
      virtual void compute_total_width() final override;
  };
 
  class GammaDiscreteExitChannel : public DiscreteExitChannel,
    public GammaExitChannel
  {
    public:
 
      GammaDiscreteExitChannel(int pdgi, int qi, double Exi, int twoJi,
        marley::Parity Pi, double rho_i, marley::StructureDatabase& sdb,
        const marley::Level& flev) : ExitChannel( pdgi, qi, Exi, twoJi, Pi,
        rho_i, sdb ), DiscreteExitChannel( flev ), GammaExitChannel()
      {
        this->compute_total_width();
      }
 
      virtual void compute_total_width() final override;
  };
 
 
  class FragmentContinuumExitChannel : public ContinuumExitChannel,
    public FragmentExitChannel
  {
    public:
 
      FragmentContinuumExitChannel(int pdgi, int qi, double Exi, int twoJi,
        marley::Parity Pi, double rho_i, marley::StructureDatabase& sdb,
        double Ec_min, const marley::Fragment& frag)
        : ExitChannel( pdgi, qi, Exi, twoJi, Pi, rho_i, sdb ),
        ContinuumExitChannel( Ec_min, sdb.get_fragment_l_max() ),
        FragmentExitChannel( frag )
      {
        this->compute_total_width();
      }
 
      virtual double differential_width( double Exf,
        bool store_jpi_widths = false ) const final override;
 
      inline virtual double E_c_max() const final override
        { return this->max_Exf(); }
  };
 
  class GammaContinuumExitChannel : public ContinuumExitChannel,
    public GammaExitChannel
  {
    public:
 
      GammaContinuumExitChannel(int pdgi, int qi, double Exi, int twoJi,
        marley::Parity Pi, double rho_i, marley::StructureDatabase& sdb,
        double Ec_min) : ExitChannel( pdgi, qi, Exi, twoJi, Pi, rho_i, sdb ),
        ContinuumExitChannel( Ec_min, sdb.get_gamma_l_max() ),
        GammaExitChannel()
      {
        this->compute_total_width();
      }
 
      virtual double differential_width( double Exf,
        bool store_jpi_widths = false ) const final override;
 
      inline virtual double E_c_max() const final override
        { return Exi_; }
  };
}