KoningDelarocheOpticalModel.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 <cmath>
#include <complex>
#include <vector>
 
#include "marley/DecayScheme.hh"
#include "marley/MassTable.hh"
#include "marley/OpticalModel.hh"
#include "marley/marley_utils.hh"
 
namespace marley {
 
  class KoningDelarocheOpticalModel : public OpticalModel {
 
    public:
 
      KoningDelarocheOpticalModel(int Z, int A, double step_size
        = DEFAULT_NUMEROV_STEP_SIZE_);
 
      virtual std::complex<double> optical_model_potential(double r,
        double fragment_KE_lab, int fragment_pdg, int two_j, int l, int two_s,
        int target_charge = 0) override;
 
      virtual double transmission_coefficient(double total_KE_CM,
        int fragment_pdg, int two_j, int l, int two_s, int target_charge = 0)
        override;
 
      virtual double total_cross_section(double fragment_KE_lab,
        int fragment_pdg, int two_s, size_t l_max, int target_charge = 0)
        override;
 
    private:
 
      double total_CM_frame_KE_;
      double fragment_mass_;
      double fragment_KE_lab_;
      double CM_frame_momentum_squared_;
 
      // Storage for the mass of the nuclear target represented by this
      // optical model potential. Its value will be adjusted based on
      // the ionization state of the target.
      double target_mass_;
 
      // Helper function for computing optical model transmission coefficients
      // and cross sections
      std::complex<double> s_matrix_element(int fragment_pdg, int two_j,
        int l, int two_s);
 
      // Helper functions for computing the optical model potential
      void calculate_om_parameters(int fragment_pdg, int two_j, int l,
        int two_s);
 
      // Compute the optical model potential at radius r
      std::complex<double> omp(double r) const;
 
      // Computes the optical model potential minus the Coulomb potential at
      // radius r
      std::complex<double> omp_minus_Vc(double r) const;
 
      // Woods-Saxon shape
      double f(double r, double R, double a) const;
 
      // Partial derivative with respect to r of the Woods-Saxon shape
      double dfdr(double r, double R, double a) const;
 
      // Coulomb potential for a point particle with charge q*e interacting
      // with a uniformly charged sphere with radius R and charge Q*e
      double Vc(double r, double R, int Q, int q) const;
 
      // Non-derivative radial Schrödinger equation terms to use for computing
      // transmission coefficients via the Numerov method
      std::complex<double> a(double r, int l);
 
      // Version of Schrodinger equation terms with the optical model potential
      // U pre-computed
      std::complex<double> a(double r, int l, std::complex<double> U) const;
 
      // Neutron parameters
      double v1n, v2n, v3n, v4n, w1n, w2n, d1n, d2n, d3n, vso1n, vso2n;
      double wso1n, wso2n, Efn, Rvn, avn, Rdn, adn, Rso_n, aso_n;
      // Proton parameters
      double v1p, v2p, v3p, v4p, w1p, w2p, d1p, d2p, d3p, vso1p, vso2p;
      double wso1p, wso2p, Efp, Vcbar_p, Rvp, avp, Rdp, adp, Rso_p, aso_p;
      // Radius for nuclear Coulomb potential
      double Rc;
 
      // Mass of a charged pion
      static constexpr double mpiplus = 139.57018; // MeV
      // Squared pion Compton wavelength
      static constexpr double lambda_piplus2 = (marley_utils::hbar_c
        / mpiplus) * (marley_utils::hbar_c / mpiplus); // fm
 
      // Temporary storage for optical model calculations
      double Rv, av, Rd, ad, Rso, aso; // Geometrical parameters
      double Vv, Wv, Wd, Vso, Wso; // Energy-dependent terms in the potential
      double spin_orbit_eigenvalue; // Eigenvalue of the spin-orbit operator
      int z; // Fragment atomic number
 
      // Threshold for abs(U - Vc) used to find a suitable matching radius for
      // computing transmission coefficients.
      static constexpr double MATCHING_RADIUS_THRESHOLD = 1e-3;
 
      double step_size_ = DEFAULT_NUMEROV_STEP_SIZE_;
 
      static constexpr double DEFAULT_NUMEROV_STEP_SIZE_ = 0.1;
 
      // More helper functions
      void calculate_kinematic_variables(double KE_tot_CM, int fragment_pdg);
      void update_target_mass(int target_charge);
  };
 
}