bethe_bloch.cpp

Go to the documentation of this file.

#include <stdlib.h>
#include <limits.h>
#include <cmath>
#include "wcpplib/ioniz/bethe_bloch.h"
#include "wcpplib/math/lorgamma.h"
#include "wcpplib/stream/prstream.h"
#include "wcpplib/clhep_units/WPhysicalConstants.h"
/*
2002, I. Smirnov
*/
 
namespace Heed {
 
double Bethe_Bloch_energy_loss(const double ratio_Z_to_A, const double I_eff,
                               const double beta, const double z) {
 
  const double beta2 = beta * beta;
  const double gamma2 = pow(lorgamma_1(beta) + 1.0, 2);
  const double coef1 =
      4 * M_PI * pow(classic_electr_radius, 2) * electron_mass_c2 * Avogadro;
  const double coef2 = z * z * ratio_Z_to_A / beta2;
  const double sum =
      log(2. * electron_mass_c2 * beta2 * gamma2 / I_eff) - beta2;
  return coef1 * coef2 * sum;
}
 
double Bethe_Bloch_energy_loss_gamma_1(const double ratio_Z_to_A,
                                       const double I_eff, const double gamma_1,
                                       const double z) {
 
  const double beta = lorbeta(gamma_1);
  const double beta2 = beta * beta;
  const double gamma = gamma_1 + 1.0;
  const double gamma2 = gamma * gamma;
  const double coef1 =
      4 * M_PI * pow(classic_electr_radius, 2) * electron_mass_c2 *
      Avogadro;  // should be 0.3071 according to PDG
  const double coef2 = z * z * ratio_Z_to_A / beta2;
  const double sum =
      log(2. * electron_mass_c2 * beta2 * gamma2 / I_eff) - beta2;
  return coef1 * coef2 * sum;
}
 
double Bethe_Bloch_restricted_energy_loss_gamma_1(
    const double ratio_Z_to_A, const double I_eff, const double m,
    const double gamma_1, const double ecut, const double z) {
 
  const double beta = lorbeta(gamma_1);
  const double beta2 = beta * beta;
  const double gamma = gamma_1 + 1.0;
  const double gamma2 = gamma * gamma;
  const double coef1 =
      2 * M_PI * pow(classic_electr_radius, 2) * electron_mass_c2 *
      Avogadro;  // should be 0.3071 according to PDG
  const double coef2 = z * z * ratio_Z_to_A / beta2;
  const double mrat = electron_mass_c2 / (m * c_squared);
  const double emax = 2.0 * electron_mass_c2 * beta2 * gamma2 /
                      (1.0 + 2.0 * gamma * mrat + mrat * mrat);
  double sum = 0.;
  if (ecut >= emax) {
    sum = log(2.0 * electron_mass_c2 * beta2 * gamma2 * emax / pow(I_eff, 2)) -
          2.0 * beta2;
  } else {
    sum = log(2.0 * electron_mass_c2 * beta2 * gamma2 * ecut / pow(I_eff, 2)) -
          beta2 * (1.0 + ecut / emax);
  }
  return coef1 * coef2 * sum;
}
 
}