marley::TabulatedXSec Class Reference

Computes inclusive lepton-nucleus cross sections using tabulated nuclear response functions. More...

#include <TabulatedXSec.hh>

Classes
struct	IntegralTerm
	Storage for terms in the integrals evaluated by compute_integral(). More...
struct	MultipoleLabel
	Simple struct representing a given multipole order, e.g., 2+. More...
struct	OptimizationMapKey
struct	OptimizationMapValue

Public Member Functions
	TabulatedXSec (int target_pdg, marley::Reaction::ProcessType p_type, marley::CoulombCorrector::CoulombMode mode, double delta_ias)
void	add_table (const std::string &file_name)
void	add_table (unsigned J, marley::Parity Pi, const std::string &file_name)
double	compute_integral (int pdg_a, double KEa, const MultipoleLabel &ml, double &diff_max, std::vector< IntegralTerm > *integral_terms=nullptr) const
	Helper function for integral that does the actual integration.
double	delta_ias () const
	Get the shift used to compute the effective energy transfer.
double	diff_xsec (int pdg_a, double KE_a, double omega, double cos_theta, const MultipoleLabel &ml) const
const ResponseTable &	get_table (const MultipoleLabel &ml) const
const std::map< MultipoleLabel, ResponseTable > &	get_table_map () const
double	integral (int pdg_a, double KEa) const
double	integral (int pdg_a, double KEa, const MultipoleLabel &ml, double &diff_max) const
bool	is_cc () const
	Returns true if this cross section represents a CC process or false otherwise.
void	optimize (int pdg_a, double max_KEa)
void	unoptimize ()

Protected Attributes
marley::CoulombCorrector::CoulombMode	coulomb_mode_
	Method to use for computing Coulomb corrections to the cross section.
double	delta_ias_
std::map< OptimizationMapKey, OptimizationMapValue >	optimization_map_
int	pdg_d_
marley::Reaction::ProcessType	proc_type_
	Kind of process for which the cross section will be computed.
std::map< MultipoleLabel, ResponseTable >	responses_
	Tables of nuclear responses organized by multipole.
TargetAtom	ta_
	Nuclide represented by this table of nuclear responses.

Detailed Description

Computes inclusive lepton-nucleus cross sections using tabulated nuclear response functions.

Definition at line 36 of file TabulatedXSec.hh.

Constructor & Destructor Documentation

TabulatedXSec()

marley::TabulatedXSec::TabulatedXSec ( int target_pdg, marley::Reaction::ProcessType p_type, marley::CoulombCorrector::CoulombMode mode, double delta_ias )

explicit

Definition at line 29 of file TabulatedXSec.cc.

                     : ta_( target_pdg ), proc_type_( p_type ),
  coulomb_mode_( mode ), delta_ias_( delta_ias )
{
  // Set the residue PDG code (pdg_d_) based on the input information
  int dummy_charge;
  marley::Reaction::get_residue_pdg_and_charge( p_type, target_pdg,
    pdg_d_, dummy_charge );
}

Member Function Documentation

add_table()

void marley::TabulatedXSec::add_table

(

const std::string &

file_name

)

Definition at line 40 of file TabulatedXSec.cc.

{
  const auto& fm = marley::FileManager::Instance();
  std::string full_file_name = fm.find_file( file_name );
  if ( full_file_name.empty() ) {
    throw marley::Error( "Could not open the nuclear response data file \""
      + file_name + '\"' );
  }
 
  std::ifstream in_file( full_file_name );
 
  // TODO: add error handling for file parsing
 
  // Multipole order
  unsigned J;
  in_file >> J;
 
  // Grid sizes
  unsigned num_w, num_q;
 
  // Temporary storage for grid points
  double w, q;
 
  // 1D grids
  auto wvec = std::make_shared< std::vector<double> >();
  auto qvec = std::make_shared< std::vector<double> >();
 
  in_file >> num_w;
  for ( unsigned iw = 0u; iw < num_w; ++iw ) {
    in_file >> w;
    wvec->push_back( w );
  }
 
  in_file >> num_q;
  for ( unsigned iq = 0u; iq < num_q; ++iq ) {
    in_file >> q;
    qvec->push_back( q );
  }
 
  // Natural parity is given by (-1)^J, while unnatural parity is the opposite
  marley::Parity natural;
  bool J_is_odd = (J % 2 == 1);
  if ( J_is_odd ) natural = -1;
  else natural = 1;
 
  marley::Parity unnatural = -natural;
 
  MultipoleLabel nat_ml( J, natural );
  MultipoleLabel unnat_ml( J, unnatural );
 
  auto nat_resp = std::make_shared<
    std::vector<marley::NuclearResponses> >();
 
  auto unnat_resp = std::make_shared<
    std::vector<marley::NuclearResponses> >();
 
  // Temporary storage for responses
  double rcc, rll, rcl, rtVV, rtAA, rtprime;
 
  // Flag used to swap back and forth between natural and unnatural parity
  bool nat = true;
  while ( in_file >> rcc >> rll >> rcl >> rtVV >> rtAA >> rtprime ) {
 
    if ( nat ) nat_resp->emplace_back( rcc, rll, rcl, rtVV, rtAA, rtprime );
    else unnat_resp->emplace_back( rcc, rll, rcl, rtVV, rtAA, rtprime );
 
    // Flip the parity flag
    nat = !nat;
  }
 
  // Store the finished tables as new objects in the map
  responses_[ nat_ml ] = ResponseTable( wvec, qvec, nat_resp );
  responses_[ unnat_ml ] = ResponseTable( wvec, qvec, unnat_resp );
 
}

compute_integral()

double marley::TabulatedXSec::compute_integral	(	int	pdg_a,
		double	KEa,
		const MultipoleLabel &	ml,
		double &	diff_max,
		std::vector< IntegralTerm > *	integral_terms = nullptr ) const

Helper function for integral that does the actual integration.

Parameters

[out]

integral_terms

Optional vector that will be cleared and filled with each term in the sum that leads to the value of the integral returned by the function. The terms are associated with energy transfer (omega, MeV) and scattering cosine (cos_theta, dimensionless) values in the IntegralTerm struct.

Definition at line 206 of file TabulatedXSec.cc.

{
  // Set the maximum differential cross section to zero to start
  diff_max = 0.;
 
  // If the input vector is not null, then clear it in preparation for
  // accumulating the terms of the integral
  if ( integral_terms ) integral_terms->clear();
 
  // Get the table of nuclear responses for the requested multipole
  const auto& rt = this->responses_.at( ml );
 
  // Get the vectors of grid points
  const auto& wvec = rt.w_grid();
  const auto& qvec = rt.q_grid();
 
  // TODO: reduce code duplication with this and diff_xsec
  // Look up the masses of the projectile and ejectile
  int pdg_c = marley::Reaction::get_ejectile_pdg( pdg_a, proc_type_ );
  const auto& mt = marley::MassTable::Instance();
  double ma = mt.get_particle_mass( pdg_a );
  double mc = mt.get_particle_mass( pdg_c );
 
  // Projectile total energy
  double Ea = KEa + ma;
 
  // Get the grid step sizes
  // NOTE: this assumes that the grid is regularly spaced (I take advantage
  // of this to simplify the trapezoidal rule for integration)
  // TODO: revisit this assumption
  size_t num_w = wvec.size();
  size_t num_q = qvec.size();
 
  size_t num_w_minus_one = num_w - 1;
  size_t num_q_minus_one = num_q - 1;
 
  double dw = ( wvec.back() - wvec.front() ) / num_w_minus_one;
  double dq = ( qvec.back() - qvec.front() ) / num_q_minus_one;
 
  // Loop over every grid point. Do an angular integral at every omega
  // point, and use the results to integrate over omega.
  double integ = 0.;
  for ( size_t iw = 0u; iw < num_w; ++iw ) {
 
    // Get the energy transfer at the current grid point. Note that the
    // tables of nuclear responses are reported on a grid that uses the
    // effective value (shifted by delta_ias_). We correct for this here.
    double w_eff = wvec.at( iw );
    double w = w_eff - this->delta_ias();
 
    // Use it to compute the ejectile total energy, etc.
    double Ec = Ea - w;
    // We can skip unphysical terms for which the total energy is smaller than
    // the final lepton mass
    if ( Ec < mc ) continue;
    double pc = marley_utils::real_sqrt( Ec*Ec - mc*mc );
    double pa = marley_utils::real_sqrt( Ea*Ea - ma*ma );
 
    double w_integ = 0.;
    for ( size_t iq = 0u; iq < num_q; ++iq ) {
 
      // Get the magnitude of the 3-momentum transfer at the current grid point
      double q = qvec.at( iq );
 
      // Convert it into a value for the scattering cosine
      double cos_theta = ( pa*pa + pc*pc - q*q ) / ( 2.*pa*pc );
 
      // Compute the differential cross section at this 2D grid point
      double diff = this->diff_xsec( pdg_a, KEa, w, cos_theta, ml );
 
      // If it is larger than any value encountered so far, record it as
      // the new maximum
      if ( diff > diff_max ) diff_max = diff;
 
      // Apply a Jacobian to transform from dwdcostheta to dwdq (since we're
      // integrating on a grid of q values)
      diff *= q / pa / pc;
 
      // Add this term to the integral over omega according to the trapezoid
      // rule
      if ( iq == 0u || iq == num_q_minus_one ) diff /= 2.;
      w_integ += diff;
 
      // If the user supplied a vector to record the terms in the integral,
      // then store the contribution from the current iteration before
      // continuing to the next grid point
      if ( integral_terms ) {
        double value = diff * dq * dw;
        if ( iw == 0u || iw == num_w_minus_one ) value /= 2.;
        integral_terms->emplace_back( w, cos_theta, value );
      }
 
    } // q grid points
 
    // We're done summing over q values. Scale the result to get the
    // integral over q for this w grid point.
    w_integ *= dq;
 
    // Now add this term to the integral over w according to the trapezoid rule
    if ( iw == 0u || iw == num_w_minus_one ) w_integ /= 2.;
 
    integ += w_integ;
 
  } // w grid points
 
  // We're done. Scale the integral over w by the needed prefactor
  integ *= dw;
 
  return integ;
}

References delta_ias(), marley::Reaction::get_ejectile_pdg(), marley::MassTable::Instance(), proc_type_, and responses_.

delta_ias()

double marley::TabulatedXSec::delta_ias

(

)

const

Get the shift used to compute the effective energy transfer.

Definition at line 448 of file TabulatedXSec.cc.

                                            {
  bool is_charged_current = this->is_cc();
 
  double result = 0.;
  if ( is_charged_current ) result = delta_ias_;
  return result;
}

References delta_ias_, and is_cc().

Referenced by compute_integral().

diff_xsec()

double marley::TabulatedXSec::diff_xsec	(	int	pdg_a,
		double	KE_a,
		double	omega,
		double	cos_theta,
		const MultipoleLabel &	ml ) const

Definition at line 116 of file TabulatedXSec.cc.

{
  int helicity = marley_utils::get_particle_helicity( pdg_a );
 
  // Look up the masses of the projectile and ejectile
  int pdg_c = marley::Reaction::get_ejectile_pdg( pdg_a, proc_type_ );
  const auto& mt = marley::MassTable::Instance();
  double ma = mt.get_particle_mass( pdg_a );
  double mc = mt.get_particle_mass( pdg_c );
 
  // Determine the total energies and momenta of the initial and final leptons
  double Ea = KEa + ma;
  double pa = marley_utils::real_sqrt( Ea*Ea - ma*ma );
  double Ec = Ea - omega;
  if ( Ec < mc ) return 0.;
 
  double pc = marley_utils::real_sqrt( Ec*Ec - mc*mc );
 
  // Get the magnitude of the 3-momentum transfer (q) from this information
  // and the scattering cosine. If the scattering cosine is unphysical, then
  // just return zero.
  if ( std::abs(cos_theta) > 1. ) return 0.;
  double q = marley_utils::real_sqrt( pa*pa + pc*pc - 2.*pa*pc*cos_theta );
 
  // Get the table of pre-computed nuclear responses for the requested
  // multipole
  const auto& rt = this->responses_.at( ml );
 
  // Shift the energy transfer at which the responses are evaluated. This
  // effective value of the energy transfer tries to correct for using the
  // same nuclear potential for the initial and final states despite
  // a change of nuclear charge for charged-current interactions.
  double omega_eff = omega + this->delta_ias();
 
  // Interpolate a set of nuclear responses for the given omega_eff and q values
  if ( omega_eff < rt.w_min() || omega_eff > rt.w_max()
    || q < rt.q_min() || q > rt.q_max() ) return 0.;
  auto nr = rt.interpolate( omega_eff, q );
 
  // Compute the lepton factors
  double beta = pc / Ec;
  double sin_theta2 = 1. - cos_theta*cos_theta;
  double q2 = q*q;
  double vcc = 1. + beta * cos_theta;
  double vll = vcc - 2.*Ea*Ec*sin_theta2*beta*beta/q2;
  double vcl = -2. * ( omega*vcc/q + mc*mc/Ec/q );
  double vT = 1. - beta*cos_theta + Ea*Ec*beta*beta*sin_theta2/q2;
  double vTprime = 2. * helicity * ( (Ea + Ec)*(1 - beta*cos_theta)/q
    - mc*mc/q/Ec );
  LeptonFactors lf( vcc, vll, vcl, vT, vTprime );
 
  // Compute the double-differential cross section with respect to
  // the energy transfer and the scattering cosine
  double xsec = lf * nr;
  xsec *= 2. * marley_utils::GF2 * marley_utils::Vud2 * Ec * pc;
 
  // For CC cross sections, apply the appropriate correction factor for
  // Coulomb corrections based on the active "Coulomb mode"
  bool is_charged_current = this->is_cc();
 
  if ( is_charged_current ) {
 
    // Get the mass of the atomic target
    double mb = mt.get_atomic_mass( ta_.pdg() );
 
    // Lab-frame total energy of the residue (via energy conservation)
    double Ed = omega + mb;
 
    // Mass of the residue (via 3-momentum conservation)
    double md = marley_utils::real_sqrt( Ed*Ed - q*q );
 
    // Dot product of the 4-momenta of the ejectile and residue
    double pc_dot_pd = Ec*Ed + pc*pc - pa*pc*cos_theta;
 
    // Manifestly Lorentz-invariant relative speed of particles c and d
    double beta_rel_cd = marley_utils::real_sqrt(
    std::pow(pc_dot_pd, 2) - mc*mc*md*md ) / pc_dot_pd;
 
    // Compute the Coulomb correction factor using the relative speed
    CoulombCorrector coul_corr( pdg_c, pdg_d_, coulomb_mode_ );
    double FC = coul_corr.coulomb_correction_factor( beta_rel_cd );
 
    // Apply it to the differential cross section
    xsec *= FC;
  }
 
  return xsec;
}

get_table()

const ResponseTable & marley::TabulatedXSec::get_table ( const MultipoleLabel & ml ) const

inline

Definition at line 83 of file TabulatedXSec.hh.

        { return responses_.at( ml ); }

get_table_map()

const std::map< MultipoleLabel, ResponseTable > & marley::TabulatedXSec::get_table_map ( ) const

inline

Definition at line 86 of file TabulatedXSec.hh.

              { return responses_; }

integral() [1/2]

double marley::TabulatedXSec::integral	(	int	pdg_a,
		double	KEa ) const

Definition at line 354 of file TabulatedXSec.cc.

{
  double integ = 0.;
  for ( const auto& pair : responses_ ) {
    double dummy;
    const auto& ml = pair.first;
    double integ_ml = this->integral( pdg_a, KEa, ml, dummy );
    integ += integ_ml;
  }
  return integ;
}

integral() [2/2]

double marley::TabulatedXSec::integral	(	int	pdg_a,
		double	KEa,
		const MultipoleLabel &	ml,
		double &	diff_max ) const

Definition at line 319 of file TabulatedXSec.cc.

{
  // First attempt to look up a pair of ChebyshevInterpolatingFunction
  // objects for the given projectile PDG code and multipole
  OptimizationMapKey key( pdg_a, ml );
  auto end = optimization_map_.end();
  auto iter = optimization_map_.find( key );
  // If we found one, attempt to compute the cross sections using them
  if ( iter != end ) {
    const auto& omv = iter->second;
    const auto& tot = omv.tot_xsec_;
    // If the requested kinetic energy is below threshold, then just
    // return zero
    if ( KEa < tot.x_min() ) return 0.;
    // If the kinetic energy is within the range covered by the
    // interpolating functions, then go ahead and use them
    else if ( KEa <= tot.x_max() ) {
      double tot_xsec = omv.tot_xsec_.evaluate( KEa );
      diff_max = omv.max_diff_xsec_.evaluate( KEa );
 
      // Ensure that interpolation problems will not give us a negative value
      // for the total or maximum differential cross sections
      tot_xsec = std::max( 0., tot_xsec );
      diff_max = std::max( 0., diff_max );
 
      return tot_xsec;
    }
  }
 
  // If we couldn't find a suitable interpolating function, then fall back to
  // brute-force integration
  return this->compute_integral( pdg_a, KEa, ml, diff_max );
}

is_cc()

bool marley::TabulatedXSec::is_cc ( ) const

inline

Returns true if this cross section represents a CC process or false otherwise.

Definition at line 51 of file TabulatedXSec.hh.

                                 {
        bool is_cc = ( proc_type_ == ProcType::NeutrinoCC_Discrete
         || proc_type_ == ProcType::AntiNeutrinoCC_Discrete
         || proc_type_ == ProcType::NeutrinoCC_Continuum
         || proc_type_ == ProcType::AntiNeutrinoCC_Continuum );
        return is_cc;
      }

References is_cc(), and proc_type_.

Referenced by delta_ias(), and is_cc().

optimize()

void marley::TabulatedXSec::optimize	(	int	pdg_a,
		double	max_KEa )

Definition at line 366 of file TabulatedXSec.cc.

                                                              {
  MARLEY_LOG( INFO, "physics.reaction.xsec" ) << "Optimizing CRPA cross"
    " section up to " << max_KEa << " MeV";
  // Loop over each of the multipoles
  for ( const auto& pair : responses_ ) {
    const auto& ml = pair.first;
 
    double min_KEa = 0.;
    std::function<double(double)> tot_xsec_func = [](double)
      -> double { return 0.; };
    std::function<double(double)> max_diff_xsec_func = tot_xsec_func;
 
    // Verify that the total cross section for this multipole is non-vanishing
    // for the requested maximum projectile kinetic energy KEa. If it vanishes,
    // then just skip the current multipole.
    double dummy;
    double xsec_at_max = this->integral( pdg_a, max_KEa, ml, dummy );
 
    if ( xsec_at_max > 0. ) {
      // Do a binary search to find the threshold for this multipole. Continue
      // until we've found it within the given tolerance.
      constexpr double thresh_tol = 1e-6; // MeV
      // Set up the bounds of a bracketing interval that will contain the
      // kinetic energy threshold
      double low_KEa = 0.;
      double high_KEa = max_KEa;
      do {
        // Check the total cross section at the midpoint of the current
        // bracketing interval
        double cur_KEa = (low_KEa + high_KEa) / 2.;
        double xsec = this->integral( pdg_a, cur_KEa, ml, dummy );
        // If it vanishes, move the lower bound up
        if ( xsec <= 0. ) low_KEa = cur_KEa;
        // If it doesn't, move the upper bound down
        else high_KEa = cur_KEa;
        // Continue until the bracketing interval is no larger than the
        // tolerance defined above
      } while ( std::abs(high_KEa - low_KEa) > thresh_tol );
 
      // Adopt the lower bound of the bracketing interval as the threshold
      min_KEa = low_KEa;
 
      tot_xsec_func = [&, this](double KEa)
        -> double { return this->integral( pdg_a, KEa, ml, dummy ); };
 
      max_diff_xsec_func = [&, this](double KEa) -> double {
        double max_diff;
        this->integral( pdg_a, KEa, ml, max_diff );
        return max_diff;
      };
    }
 
    MARLEY_LOG( DEBUG, "physics.reaction.xsec" ) << "Optimizing total cross section for "
      << ml.J_ << ml.Pi_ << " over KE in ["
      << min_KEa << ", " << max_KEa << "] MeV";
 
    // Now we're ready to build the Chebyshev interpolating functions for
    // this multipole. We need one for the total cross section and the
    // other one for the maximum of the differential cross section.
    ChebyshevInterpolatingFunction tot_xs_cif( tot_xsec_func, min_KEa,
      max_KEa, 64 );
    // TODO: do you want adaptive grid sizing here?
 
    MARLEY_LOG( DEBUG, "physics.reaction.xsec" ) << "Optimizing max diff for "
      << ml.J_ << ml.Pi_;
 
    ChebyshevInterpolatingFunction max_diff_xs_cif( max_diff_xsec_func,
      min_KEa, max_KEa, 64 );
 
    // Build the key and value we need for the map entry
    OptimizationMapKey key( pdg_a, ml );
    OptimizationMapValue value( tot_xs_cif, max_diff_xs_cif );
 
    optimization_map_.emplace( std::make_pair(key, value) );
  }
}

unoptimize()

void marley::TabulatedXSec::unoptimize ( )

inline

Definition at line 96 of file TabulatedXSec.hh.

{ optimization_map_.clear(); }

Member Data Documentation

coulomb_mode_

marley::CoulombCorrector::CoulombMode marley::TabulatedXSec::coulomb_mode_

protected

Method to use for computing Coulomb corrections to the cross section.

Definition at line 158 of file TabulatedXSec.hh.

delta_ias_

double marley::TabulatedXSec::delta_ias_

protected

Effective shift in the energy transfer that accounts for the mass difference between the ground state of the initial nucleus and the isobaric analog state in the daughter nucleus

Definition at line 174 of file TabulatedXSec.hh.

Referenced by delta_ias().

optimization_map_

std::map< OptimizationMapKey, OptimizationMapValue > marley::TabulatedXSec::optimization_map_

protected

Storage for interpolating functions used to optimize calculations requiring numerical integration

Definition at line 169 of file TabulatedXSec.hh.

pdg_d_

int marley::TabulatedXSec::pdg_d_

protected

PDG code of the nuclear residue produced by the reaction described by this cross section

Definition at line 162 of file TabulatedXSec.hh.

proc_type_

marley::Reaction::ProcessType marley::TabulatedXSec::proc_type_

protected

Kind of process for which the cross section will be computed.

Definition at line 155 of file TabulatedXSec.hh.

Referenced by compute_integral(), and is_cc().

responses_

std::map< MultipoleLabel, ResponseTable > marley::TabulatedXSec::responses_

protected

Tables of nuclear responses organized by multipole.

Definition at line 165 of file TabulatedXSec.hh.

Referenced by compute_integral().

ta_

TargetAtom marley::TabulatedXSec::ta_

protected

Nuclide represented by this table of nuclear responses.

Definition at line 152 of file TabulatedXSec.hh.

---

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

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