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 std::shared_ptr< HepMC3::GenParticle >&
        compound_nucleus, std::shared_ptr< HepMC3::GenParticle >&
        emitted_particle, std::shared_ptr< HepMC3::GenParticle >&
        residual_nucleus, int& qIon, 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 std::shared_ptr< HepMC3::GenParticle >& compound_nucleus,
        std::shared_ptr< HepMC3::GenParticle >& emitted_particle,
        std::shared_ptr< HepMC3::GenParticle >& residual_nucleus,
        double Exf, int& qf, 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 std::shared_ptr< HepMC3::GenParticle >&
        compound_nucleus, std::shared_ptr< HepMC3::GenParticle >&
        emitted_particle, std::shared_ptr< HepMC3::GenParticle >&
        residual_nucleus, int& qIon, 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 std::shared_ptr< HepMC3::GenParticle >&
        compound_nucleus, std::shared_ptr< HepMC3::GenParticle >&
        emitted_particle, std::shared_ptr< HepMC3::GenParticle >&
        residual_nucleus, int& qIon, 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 ) {}
 
        virtual ~SpinParityWidth() {}
 
        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;

      inline const SpinParityWidth* get_last_sampled_spw() const
        { return last_sampled_spw_; }

      inline const std::vector< std::unique_ptr< SpinParityWidth > >&
        get_spw_table() const { return jpi_widths_table_; }
 
    protected:

      double E_c_min_;

      int l_max_;

      mutable std::vector< std::unique_ptr< SpinParityWidth > >
        jpi_widths_table_;

      mutable bool skip_jpi_sampling_ = false;

      mutable std::unique_ptr<marley::ChebyshevInterpolatingFunction> Exf_cdf_;

      mutable SpinParityWidth* last_sampled_spw_ = nullptr;

      void clear_jpi_widths() const;
  };

  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:
 
      struct FragmentSpinParityWidth :
        public ContinuumExitChannel::SpinParityWidth
      {
        FragmentSpinParityWidth( int twoJ, marley::Parity p, double w,
          int two_j_fr, int ell )
          : SpinParityWidth( twoJ, p, w ), two_j_frag( two_j_fr ),
          orb_l( ell ) {}
 
        virtual ~FragmentSpinParityWidth() {}
 
        int two_j_frag; 
        int orb_l; 
      };

      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:
 
      struct GammaSpinParityWidth :
        public ContinuumExitChannel::SpinParityWidth
      {
        GammaSpinParityWidth( int twoJ, marley::Parity p, double w, int mpol )
          : SpinParityWidth( twoJ, p, w ), multipolarity( mpol ) {}
 
        int multipolarity; 
      };

      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_; }
  };
}