NuclearReaction.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 <map>
#include <memory>
#include <regex>
#include <string>
#include <vector>
 
#include "marley/DecayScheme.hh"
#include "marley/Event.hh"
#include "marley/Level.hh"
#include "marley/MassTable.hh"
#include "marley/MatrixElement.hh"
#include "marley/Reaction.hh"
#include "marley/StructureDatabase.hh"
 
namespace marley {
 
  class Generator;
 
  class NuclearReaction : public marley::Reaction {
 
    public:
 
      NuclearReaction(ProcessType pt, int pdg_a, int pdg_b, int pdg_c,
        int pdg_d, int q_d,
        const std::shared_ptr<std::vector<marley::MatrixElement> >& mat_els);
 
      enum class CoulombMode { NO_CORRECTION, FERMI_FUNCTION, EMA, MEMA,
        FERMI_AND_EMA, FERMI_AND_MEMA };
 
      inline virtual marley::TargetAtom atomic_target() const override final
        { return marley::TargetAtom( pdg_b_ ); }
 
      virtual marley::Event create_event(int particle_id_a,
        double KEa, marley::Generator& gen) const override;
 
      double fermi_function(double beta_c) const;
 
      double max_level_energy(double KEa) const;
 
      void set_decay_scheme(marley::DecayScheme* scheme);
 
      double threshold_kinetic_energy() const override;
 
      virtual double total_xs(int pdg_a, double KEa) const override;
 
      double diff_xs(const marley::MatrixElement& mat_el, double KEa,
        double cos_theta_c_cm) const;
 
      virtual double diff_xs(int pdg_a, double KEa,
        double cos_theta_c_cm) const override;
 
      virtual double total_xs(const marley::MatrixElement& me, double KEa,
        double& beta_c_cm, bool check_max_E_level = true) const;
 
      double coulomb_correction_factor(double beta_rel_cd) const;
 
      double ema_factor(double beta_rel_cd, bool& ok, bool modified_ema) const;
 
      double weak_nuclear_charge() const;
 
      inline const std::vector<marley::MatrixElement>& matrix_elements() const
        { return *matrix_elements_; }
 
      inline CoulombMode coulomb_mode() const
        { return coulomb_mode_; }
 
      inline void set_coulomb_mode( CoulombMode mode )
        { coulomb_mode_ = mode; }
 
      static CoulombMode coulomb_mode_from_string( const std::string& str );
 
      static std::string string_from_coulomb_mode( CoulombMode mode );
 
    protected:
 
      static std::map<CoulombMode, std::string> coulomb_mode_string_map_;
 
      virtual marley::Event make_event_object(double KEa,
        double pc_cm, double cos_theta_c_cm, double phi_c_cm, double Ec_cm,
        double Ed_cm, double E_level, int twoJ, const marley::Parity& P)
        const override;
 
      double sample_cos_theta_c_cm(const marley::MatrixElement& matrix_el,
        double beta_c_cm, marley::Generator& gen) const;
 
      double summed_xs_helper(int pdg_a, double KEa, double cos_theta_c_cm,
        std::vector<double>* level_xsecs, bool differential) const;
 
      void set_description();
 
      double md_gs_; 
 
      int Zi_; 
      int Ai_; 
      int Zf_; 
      int Af_; 
 
      int q_d_;
 
      double KEa_threshold_;
 
      CoulombMode coulomb_mode_ = CoulombMode::FERMI_AND_MEMA;
 
      std::shared_ptr< std::vector<marley::MatrixElement> > matrix_elements_;
  };
 
}