marley::CoulombCorrector Class Reference

Computes Coulomb correction factors for neutrino-nucleus differential cross sections. More...

#include <CoulombCorrector.hh>

Public Types
enum class	CoulombMode {   NO_CORRECTION , FERMI_FUNCTION , EMA , MEMA ,   FERMI_AND_EMA , FERMI_AND_MEMA }
	Enumerated type used to set the method for handling Coulomb corrections for CC nuclear reactions. More...

Public Member Functions
	CoulombCorrector (int pdg_c, int pdg_d, CoulombMode mode=CoulombMode::FERMI_AND_MEMA)
double	coulomb_correction_factor (double beta_rel_cd) const
CoulombMode	coulomb_mode () const
double	ema_factor (double beta_rel_cd, bool &ok, bool modified_ema) const
double	fermi_function (double beta_c) const
	Compute the Fermi function
void	set_coulomb_mode (CoulombMode mode)
	Set the method for handling Coulomb corrections for this reaction.

Static Public Member Functions
static CoulombMode	coulomb_mode_from_string (const std::string &str)
	Convert a string to a CoulombMode value.
static std::string	string_from_coulomb_mode (CoulombMode mode)
	Convert a CoulombMode value to a string.

Protected Attributes
int	Af_
	Residue mass number.
CoulombMode	coulomb_mode_
	The method to use when computing Coulomb corrections.
double	mc_
	Ejectile mass.
int	pdg_c_
	Ejectile PDG code.
int	Zf_
	Residue atomic number.

Static Protected Attributes
static std::map< CoulombMode, std::string >	coulomb_mode_string_map_

Detailed Description

Computes Coulomb correction factors for neutrino-nucleus differential cross sections.

Definition at line 26 of file CoulombCorrector.hh.

Member Enumeration Documentation

CoulombMode

enum class marley::CoulombCorrector::CoulombMode

strong

Enumerated type used to set the method for handling Coulomb corrections for CC nuclear reactions.

Definition at line 32 of file CoulombCorrector.hh.

                             { NO_CORRECTION, FERMI_FUNCTION, EMA, MEMA,
        FERMI_AND_EMA, FERMI_AND_MEMA };

Constructor & Destructor Documentation

CoulombCorrector()

marley::CoulombCorrector::CoulombCorrector	(	int	pdg_c,
		int	pdg_d,
		CoulombMode	mode = CoulombMode::FERMI_AND_MEMA )

Parameters

pdg_c	Ejectile PDG code
pdg_d	Residue PDG code

Definition at line 62 of file CoulombCorrector.cc.

                                : pdg_c_( pdg_c ), coulomb_mode_( mode )
{
  // Extract the mass number and proton number from the residue PDG code
  Af_ = marley_utils::get_particle_A( pdg_d );
  Zf_ = marley_utils::get_particle_Z( pdg_d );
 
  // Get the ejectile mass from the mass table
  const marley::MassTable& mt = marley::MassTable::Instance();
  mc_ = mt.get_particle_mass( pdg_c_ );
 
  // Do a sanity check that the ejectile is actually a charged (anti)lepton. If
  // it's not, then set the Coulomb mode to NO_CORRECTION to prevent misuse
  // of this object.
  if ( !is_charged_lepton_or_antilepton(pdg_c_) ) {
    coulomb_mode_ = CoulombMode::NO_CORRECTION;
  }
}

References Af_, coulomb_mode_, marley::MassTable::get_particle_mass(), marley::MassTable::Instance(), mc_, pdg_c_, and Zf_.

Member Function Documentation

coulomb_correction_factor()

double marley::CoulombCorrector::coulomb_correction_factor

(

double

beta_rel_cd

)

const

Computes an approximate correction factor to account for effects of the Coulomb potential when calculating cross sections

Parameters

beta_rel_cd

The relative speed of the final particles c and d (dimensionless)

Definition at line 117 of file CoulombCorrector.cc.

{
  // Don't do anything if Coulomb corrections are switched off
  if ( coulomb_mode_ == CoulombMode::NO_CORRECTION ) {
    MARLEY_LOG( TRACE, "physics.coulomb" ) << "Coulomb correction factor = 1"
      " (NO_CORRECTION mode)";
    return 1.;
  }
 
  // Fermi function approach to the Coulomb correction
  double fermi_func = fermi_function( beta_rel_cd );
 
  // Unconditionally return the value of the Fermi function if the user
  // has configured things this way
  if ( coulomb_mode_ == CoulombMode::FERMI_FUNCTION ) {
    MARLEY_LOG( TRACE, "physics.coulomb" ) << "Coulomb correction factor = "
      << fermi_func << " (Fermi function, beta_rel = " << beta_rel_cd << ")";
    return fermi_func;
  }
 
  bool use_mema = false;
  if ( coulomb_mode_ == CoulombMode::MEMA
    || coulomb_mode_ == CoulombMode::FERMI_AND_MEMA )
  {
    use_mema = true;
  }
 
  // Effective momentum approximation for the Coulomb correction
  bool EMA_ok = false;
  double factor_EMA = ema_factor( beta_rel_cd, EMA_ok, use_mema );
 
  if ( coulomb_mode_ == CoulombMode::EMA
    || coulomb_mode_ == CoulombMode::MEMA )
  {
    if ( EMA_ok ) return factor_EMA;
    else {
      std::string model_name( "EMA" );
      if ( coulomb_mode_ == CoulombMode::MEMA ) model_name = "MEMA";
      throw marley::Error( "Invalid " + model_name + " factor encountered"
        " in marley::CoulombCorrector::coulomb_correction_factor()" );
    }
  }
 
  if ( coulomb_mode_ != CoulombMode::FERMI_AND_EMA
    && coulomb_mode_ != CoulombMode::FERMI_AND_MEMA )
  {
    throw marley::Error( "Unrecognized Coulomb correction mode encountered"
      " in marley::CoulombCorrector::coulomb_correction_factor()" );
  }
 
  // If we've gotten this far, then we're interpolating between the Fermi
  // function and the (M)EMA correction factor. If the (modified) effective
  // momentum approximation is invalid because subtracting off the Coulomb
  // potential brings the reaction below threshold, then just use the Fermi
  // function
  if ( !EMA_ok ) return fermi_func;
 
  // Otherwise, choose the approach that yields the smaller correction (i.e.,
  // the correction factor that is closest to unity).
  double diff_Fermi = std::abs( fermi_func - 1. );
  double diff_EMA = std::abs( factor_EMA - 1. );
 
  double result = ( diff_Fermi < diff_EMA ) ? fermi_func : factor_EMA;
  MARLEY_LOG( TRACE, "physics.coulomb" ) << "Coulomb correction factor = "
    << result << " (beta_rel = " << beta_rel_cd
    << ", Fermi = " << fermi_func << ", EMA = " << factor_EMA << ")";
  return result;
}

References coulomb_mode_, ema_factor(), and fermi_function().

coulomb_mode()

CoulombMode marley::CoulombCorrector::coulomb_mode ( ) const

inline

Return the method used by this reaction for handling Coulomb corrections

Definition at line 67 of file CoulombCorrector.hh.

        { return coulomb_mode_; }

References coulomb_mode_.

coulomb_mode_from_string()

CMode marley::CoulombCorrector::coulomb_mode_from_string ( const std::string & str )

static

Convert a string to a CoulombMode value.

Definition at line 248 of file CoulombCorrector.cc.

{
  for ( const auto& pair : coulomb_mode_string_map_ ) {
    if ( str == pair.second ) return pair.first;
  }
  throw marley::Error( "The string \"" + str + "\" was not recognized"
    " as a valid Coloumb mode setting" );
}

References coulomb_mode_string_map_.

ema_factor()

double marley::CoulombCorrector::ema_factor	(	double	beta_rel_cd,
		bool &	ok,
		bool	modified_ema ) const

Computes a Coulomb correction factor according to the effective momentum approximation. See J. Engel, Phys. Rev. C 57, 2004 (1998)

Parameters

	beta_rel_cd	The relative speed of the final particles c and d (dimensionless)
[out]	ok	Flag that is set to false if subtracting the Coulomb potential pulls the event below threshold
	modified_ema	If true, the modified EMA correction factor will be returned instead of that specified by the original EMA

Definition at line 188 of file CoulombCorrector.cc.

{
  // If particle c has a positive PDG code, then it is a negatively-charged
  // lepton
  bool minus_c = ( pdg_c_ > 0 );
 
  // Approximate nuclear radius (MeV^(-1))
  double R_nuc = nuclear_radius_natural_units( Af_ );
 
  // Approximate Coulomb potential
  double Vc = ( -3. * Zf_ * marley_utils::alpha ) / ( 2. * R_nuc );
  // Adjust if needed for a final-state charged antilepton
  if ( !minus_c ) Vc *= -1;
 
  // Like the Fermi function, this approximation uses a static nuclear Coulomb
  // potential (a sphere at the origin). Typically nuclear recoil is neglected,
  // allowing one to compute the effective lepton momentum in the lab frame. In
  // MARLEY's case, we do this by calculating it in the rest frame of the final
  // nucleus ("FNR" frame). We already have the relative speed of the final
  // nucleus and outgoing lepton, so this is easy.
  double gamma_rel_cd = std::pow( 1. - std::pow(beta_rel_cd, 2), -0.5 );
 
  // Check for numerical errors from the square root
  if ( !std::isfinite(gamma_rel_cd) ) {
    MARLEY_LOG( WARN, "physics.coulomb" ) << "Invalid beta_rel = "
      << beta_rel_cd
      << " encountered in marley::CoulombCorrector::ema_factor()";
  }
 
  // Total energy of the outgoing lepton in the FNR frame
  double E_c_FNR = gamma_rel_cd * mc_;
 
  // Lepton momentum in FNR frame
  double p_c_FNR = beta_rel_cd * E_c_FNR;
 
  // Effective FNR frame total energy
  double E_c_FNR_eff = E_c_FNR - Vc;
 
  // If subtracting off the Coulomb potential drops the effective energy
  // below the lepton mass, then the expression for the effective momentum
  // will give an imaginary value. Signal this by setting the "ok" flag to
  // false.
  ok = ( E_c_FNR_eff >= mc_ );
 
  // Effective momentum in FNR frame
  double p_c_FNR_eff = marley_utils::real_sqrt(
    std::pow(E_c_FNR_eff, 2) - mc_*mc_ );
 
  // Coulomb correction factor for the original EMA
  double F_EMA = std::pow( p_c_FNR_eff / p_c_FNR, 2 );
 
  // Coulomb correction factor for the modified EMA
  double F_MEMA = ( p_c_FNR_eff * E_c_FNR_eff ) / ( p_c_FNR * E_c_FNR );
 
  if ( modified_ema ) return F_MEMA;
  else return F_EMA;
}

References Af_, mc_, pdg_c_, and Zf_.

Referenced by coulomb_correction_factor().

fermi_function()

double marley::CoulombCorrector::fermi_function

(

double

beta_c

)

const

Compute the Fermi function

Parameters

beta_c

Dimensionless speed of the ejectile

Definition at line 87 of file CoulombCorrector.cc.

                                                                   {
 
  // If the PDG code for particle c is positive, then it is a
  // negatively-charged lepton.
  bool c_minus = ( pdg_c_ > 0 );
 
  // Lorentz factor gamma for particle c
  double gamma_c = std::pow( 1. - beta_c*beta_c, -marley_utils::ONE_HALF );
 
  double s = std::sqrt( 1. - std::pow(marley_utils::alpha * Zf_, 2) );
 
  // Estimate the nuclear radius (in MeV^(-1))
  double rho = nuclear_radius_natural_units( Af_ );
 
  // Sommerfeld parameter
  double eta = marley_utils::alpha * Zf_ / beta_c;
 
  // Adjust the value of eta if the light product from this reaction is
  // an antilepton
  if ( !c_minus ) eta *= -1;
 
  // Complex variable for the gamma function
  std::complex<double> a( s, eta );
  double b = std::tgamma( 1 + 2*s );
 
  return 2. * ( 1. + s ) * std::pow( 2.*beta_c*gamma_c*rho*mc_, 2.*s - 2. )
    * std::exp( marley_utils::pi*eta ) * std::norm( marley_utils::gamma(a) )
    / std::pow( b, 2 );
}

References Af_, mc_, pdg_c_, and Zf_.

Referenced by coulomb_correction_factor().

set_coulomb_mode()

void marley::CoulombCorrector::set_coulomb_mode

(

CoulombMode

mode

)

Set the method for handling Coulomb corrections for this reaction.

Definition at line 269 of file CoulombCorrector.cc.

                                                                {
  if ( is_charged_lepton_or_antilepton(pdg_c_) ) coulomb_mode_ = mode;
  else coulomb_mode_ = CoulombMode::NO_CORRECTION;
}

References coulomb_mode_, and pdg_c_.

string_from_coulomb_mode()

std::string marley::CoulombCorrector::string_from_coulomb_mode ( CoulombMode mode )

static

Convert a CoulombMode value to a string.

Definition at line 259 of file CoulombCorrector.cc.

                                                                       {
  auto it = coulomb_mode_string_map_.find( mode );
  if ( it != coulomb_mode_string_map_.end() ) return it->second;
  else throw marley::Error( "Unrecognized CoulombMode value encountered in"
    " marley::CoulombCorrector::string_from_coulomb_mode()" );
}

References coulomb_mode_string_map_.

Member Data Documentation

Af_

int marley::CoulombCorrector::Af_

protected

Residue mass number.

Definition at line 89 of file CoulombCorrector.hh.

Referenced by CoulombCorrector(), ema_factor(), and fermi_function().

coulomb_mode_

CoulombMode marley::CoulombCorrector::coulomb_mode_

protected

The method to use when computing Coulomb corrections.

Definition at line 92 of file CoulombCorrector.hh.

Referenced by CoulombCorrector(), coulomb_correction_factor(), coulomb_mode(), and set_coulomb_mode().

coulomb_mode_string_map_

std::map< CMode, std::string > marley::CoulombCorrector::coulomb_mode_string_map_

staticprotected

Initial value:

=
{
  { CMode::NO_CORRECTION, "none" },
  { CMode::FERMI_FUNCTION, "Fermi" },
  { CMode::EMA, "EMA" },
  { CMode::MEMA, "MEMA" },
  { CMode::FERMI_AND_EMA, "Fermi-EMA" },
  { CMode::FERMI_AND_MEMA, "Fermi-MEMA" },
}

Helper map used by the methods to convert a CoulombMode value to and from a std::string

Definition at line 28 of file CoulombCorrector.hh.

                             { NO_CORRECTION, FERMI_FUNCTION, EMA, MEMA,
        FERMI_AND_EMA, FERMI_AND_MEMA };

      CoulombCorrector( int pdg_c, int pdg_d,

Referenced by coulomb_mode_from_string(), and string_from_coulomb_mode().

mc_

double marley::CoulombCorrector::mc_

protected

Ejectile mass.

Definition at line 86 of file CoulombCorrector.hh.

Referenced by CoulombCorrector(), ema_factor(), and fermi_function().

pdg_c_

int marley::CoulombCorrector::pdg_c_

protected

Ejectile PDG code.

Definition at line 85 of file CoulombCorrector.hh.

Referenced by CoulombCorrector(), ema_factor(), fermi_function(), and set_coulomb_mode().

Zf_

int marley::CoulombCorrector::Zf_

protected

Residue atomic number.

Definition at line 88 of file CoulombCorrector.hh.

Referenced by CoulombCorrector(), ema_factor(), and fermi_function().

---

The documentation for this class was generated from the following files:

• /home/sjg/Desktop/marley/include/marley/CoulombCorrector.hh
• /home/sjg/Desktop/marley/src/CoulombCorrector.cc