Garfield::TrackBichsel Class Reference

#include <TrackBichsel.hh>

Inheritance diagram for Garfield::TrackBichsel:

Classes
struct	Cluster

Public Member Functions
	TrackBichsel ()
	Default constructor.
	TrackBichsel (Sensor *sensor)
	Constructor.
virtual	~TrackBichsel ()
	Destructor.
bool	NewTrack (const double x0, const double y0, const double z0, const double t0, const double dx0, const double dy0, const double dz0) override
bool	GetCluster (double &xc, double &yc, double &zc, double &tc, int &nc, double &ec, double &extra) override
const std::vector< Cluster > &	GetClusters () const
double	GetClusterDensity () override
double	GetStoppingPower () override
	Get the stopping power (mean energy loss [eV] per cm).
bool	Initialise ()
bool	ComputeCrossSection ()
Public Member Functions inherited from Garfield::Track
	Track ()=delete
	Default constructor.
	Track (const std::string &name)
	Constructor.
virtual	~Track ()
	Destructor.
virtual void	SetParticle (const std::string &part)
void	SetEnergy (const double e)
	Set the particle energy.
void	SetBetaGamma (const double bg)
	Set the relative momentum of the particle.
void	SetBeta (const double beta)
	Set the speed ( ) of the particle.
void	SetGamma (const double gamma)
	Set the Lorentz factor of the particle.
void	SetMomentum (const double p)
	Set the particle momentum.
void	SetKineticEnergy (const double ekin)
	Set the kinetic energy of the particle.
double	GetEnergy () const
	Return the particle energy.
double	GetBetaGamma () const
	Return the of the projectile.
double	GetBeta () const
	Return the speed ( ) of the projectile.
double	GetGamma () const
	Return the Lorentz factor of the projectile.
double	GetMomentum () const
	Return the particle momentum.
double	GetKineticEnergy () const
	Return the kinetic energy of the projectile.
double	GetCharge () const
	Get the charge of the projectile.
double	GetMass () const
	Get the mass [eV / c2] of the projectile.
void	SetSensor (Sensor *s)
	Set the sensor through which to transport the particle.
void	EnablePlotting (ViewDrift *viewer)
	Switch on plotting.
void	DisablePlotting ()
	Switch off plotting.
void	EnableDebugging ()
	Switch on debugging messages.
void	DisableDebugging ()
	Switch off debugging messages.

Additional Inherited Members
Protected Member Functions inherited from Garfield::Track
void	PlotNewTrack (const double x0, const double y0, const double z0)
void	PlotCluster (const double x0, const double y0, const double z0)
Static Protected Member Functions inherited from Garfield::Track
static std::array< double, 3 >	StepBfield (const double dt, const double qoverm, const double vmag, double bx, double by, double bz, std::array< double, 3 > &dir)
Protected Attributes inherited from Garfield::Track
std::string	m_className = "Track"
double	m_q = -1.
int	m_spin = 1
double	m_mass
double	m_energy = 0.
double	m_beta2 = 1.
bool	m_isElectron = false
std::string	m_particleName = "mu-"
Sensor *	m_sensor = nullptr
bool	m_isChanged = true
ViewDrift *	m_viewer = nullptr
bool	m_debug = false
size_t	m_plotId = 0

Detailed Description

Generate tracks using differential cross-sections for silicon computed by Hans Bichsel. References:

• H. Bichsel, Rev. Mod. Phys. 60 (1988), 663-699
• https://faculty.washington.edu/hbichsel/

Definition at line 17 of file TrackBichsel.hh.

Constructor & Destructor Documentation

TrackBichsel() [1/2]

Garfield::TrackBichsel::TrackBichsel ( )

inline

Default constructor.

Definition at line 25 of file TrackBichsel.hh.

: TrackBichsel(nullptr) {}

Referenced by TrackBichsel().

TrackBichsel() [2/2]

Garfield::TrackBichsel::TrackBichsel

(

Sensor *

sensor

)

Constructor.

Definition at line 26 of file TrackBichsel.cc.

                                         : Track("Bichsel") {
  m_sensor = sensor;
  Initialise();
}

~TrackBichsel()

virtual Garfield::TrackBichsel::~TrackBichsel ( )

inlinevirtual

Destructor.

Definition at line 29 of file TrackBichsel.hh.

{}

Member Function Documentation

ComputeCrossSection()

bool Garfield::TrackBichsel::ComputeCrossSection

(

)

Definition at line 189 of file TrackBichsel.cc.

                                       {
 
  if (!m_initialised) {
    std::cerr << m_className << "::ComputeCrossSection: Not initialised.\n";
    return false;
  }
 
  const double bg = GetBetaGamma();
  if (m_debug) {
    std::cerr << m_className << "::ComputeCrossSection:\n"
              << "    Calculating differential cross-section for bg = "
              << bg << ".\n";
  }
  const double gamma = sqrt(bg * bg + 1.);
  // Prefactors in Bhabha formula.
  const double g1 = (gamma - 1.) * (gamma - 1.) / (gamma * gamma);
  const double g2 = (2. * gamma  - 1.) / (gamma * gamma);
 
  const double ek = GetKineticEnergy();
  const double rm = ElectronMass / m_mass;
  // Maximum energy transfer.
  double emax = 2 * ElectronMass * bg * bg;
  if (m_isElectron) {
    emax = 0.5 * ek;
  } else {
    emax = 2 * ElectronMass * bg * bg / (1. + 2 * gamma * rm + rm * rm);
  }
  if (m_debug) std::printf("    Max. energy transfer: %12.4f eV\n", emax);
  const double betaSq = bg * bg / (1. + bg * bg);
  m_speed = SpeedOfLight * sqrt(betaSq);
  constexpr size_t nTerms = 3;
  std::array<double, nTerms + 1> m0;
  std::array<double, nTerms + 1> m1;
  std::array<double, nTerms + 1> m2; 
  m0.fill(0.);
  m1.fill(0.);
  m2.fill(0.);
 
  std::array<std::array<double, NEnergyBins>, nTerms + 1> cs;
  std::vector<double> cdf(NEnergyBins, 0.);
 
  const auto betaSqOverEmax = betaSq / emax;
  const auto twoMeBetaSq = 2 * ElectronMass * betaSq;
 
  size_t jmax = 0; 
  for (size_t j = 0; j < NEnergyBins; ++j) {
    ++jmax;
    if (m_E[j] > emax) break;
    const auto e2 = m_E[j] * m_E[j];
    double q1 = RydbergEnergy;
    // CCS-33, 39 & 47
    if (m_E[j] < 11.9) {
      q1 = m_k1[j] * m_k1[j] * RydbergEnergy;
    } else if (m_E[j] < 100.) {
      constexpr double k1 = 0.025;
      q1 = k1 * k1 * RydbergEnergy;
    }
    const auto qmin = e2 / twoMeBetaSq;
    if (m_E[j] < 11.9 && q1 <= qmin) {
      cs[0][j] = 0.;
    } else {
      cs[0][j] = m_E[j] * m_dfdE[j] * log(q1 / qmin);
    }
    // Fano Eq. (47).
    auto b1 = 1. - betaSq * m_eps1[j];
    if (b1 == 0.) b1 = 1.e-20;
    const auto b2 = betaSq * m_eps2[j];
    const auto g = m_E[j] * m_dfdE[j] * (-0.5) * log(b1 * b1 + b2 * b2);
    auto theta = atan(b2 / b1);
    if (theta < 0.) theta += Pi;
    const auto epsSq = m_eps1[j] * m_eps1[j] + m_eps2[j] * m_eps2[j];
    const auto h = m_conv * e2 * (betaSq - m_eps1[j] / epsSq) * theta;
    cs[1][j] = g + h;
    // The integral was over d(lnK) rather than d(lnQ)
    cs[2][j] = 2 * m_int[j];
    if (m_isElectron) {
      // Uehling Eq. (9).
      const double u = m_E[j] / (ek - m_E[j]);
      const double v = m_E[j] / ek;
      cs[2][j] *= (1. + u * u + g1 * v * v - g2 * u);
    } else {
      // Uehling Eq. (2).
      cs[2][j] *= (1. - m_E[j] * betaSqOverEmax);
    } 
    cs[3][j] = 0.;
    cs[nTerms][j] = 0.;
    const double dE = m_E[j + 1] - m_E[j];
    for (size_t k = 0; k < nTerms; ++k) {
      m0[k] += cs[k][j] * dE / e2;
      m1[k] += cs[k][j] * dE / m_E[j];
      m2[k] += cs[k][j] * dE;
      cs[nTerms][j] += cs[k][j];
    } 
    m0[nTerms] += cs[nTerms][j] * dE / e2;
    m1[nTerms] += cs[nTerms][j] * dE / m_E[j]; 
    m2[nTerms] += cs[nTerms][j] * dE;
    cs[nTerms][j] /= e2;
    if (!m_debug) continue;
    if ((j + 1) % 10 != 0) continue;
    std::printf(" %4zu %9.1f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f\n", 
                j + 1, m_E[j], m_dfdE[j], g, h, cs[0][j], cs[1][j], cs[2][j]);
  }
  if (jmax < NEnergyBins) cdf.resize(jmax);
  if (m_debug) {
    printf("  M0: %15.4f %15.4f %15.4f %15.4f\n", m0[0], m0[1], m0[2], m0[3]);
    printf("  M1: %15.4f %15.4f %15.4f %15.4f\n", m1[0], m1[1], m1[2], m1[3]);
    printf("  M2: %15.4f %15.4f %15.4f %15.4f\n", m2[0], m2[1], m2[2], m2[3]);
  }
  // Calculate the cumulative differential cross-section.
  cdf[0] = 0.5 * m_E[0] * cs[nTerms][0];
  for (size_t j = 1; j < jmax; ++j) {
    const double dE = m_E[j] - m_E[j - 1];
    cdf[j] = cdf[j - 1] + 0.5 * (cs[nTerms][j - 1] + cs[nTerms][j]) * dE;
  }
 
  constexpr double ary = BohrRadius * RydbergEnergy;
  constexpr double prefactor = 8 * Pi * ary * ary / ElectronMass;
  const double dec = m_q * m_q * m_density * prefactor / betaSq;
  m_imfp = m0.back() * dec;
  m_dEdx = m1.back() * dec;
  if (m_debug) {
    std::printf("    M0 = %12.4f cm-1      ... inverse mean free path\n", 
                m_imfp);
    std::printf("    M1 = %12.4f keV/cm    ... dE/dx\n", m_dEdx * 1.e-3);
    std::printf("    M2 = %12.4f keV2/cm\n", m2.back() * dec * 1.e-6);
  }
  // Calculate the residual cross-section.
  double integral = cdf.back();
  if (emax > m_E.back()) {
    const double e1 = m_E.back();
    const double rm0 = (1. - 720. * betaSqOverEmax) * (
      (1. / e1 - 1. / emax) + 
      2 * (1. / (e1 * e1) - 1. / (emax * emax))) -
      betaSqOverEmax * log(emax / e1);
    if (m_debug) std::printf("    Residual M0 = %15.5f\n", rm0 * 14. * dec);
    m_imfp += rm0 * 14. * dec;
    m0.back() += rm0;
    integral += 14. * rm0;
  }
  const double scale = 1. / integral;
  for (size_t j = 0; j < jmax; ++j) {
    cdf[j] *= scale;
    if (!m_debug) continue;
    if ((j + 1) % 20 != 0) continue;
    std::printf(" %4zu %9.1f %15.6f\n", j + 1, m_E[j], cdf[j]); 
  } 
  m_tab.fill(0.);
  for (size_t i = 0; i < NCdfBins; ++i) {
    constexpr double step = 1. / NCdfBins;
    const double x = (i + 1) * step;
    // Interpolate.
    m_tab[i] = 0.;
    const auto it1 = std::upper_bound(cdf.cbegin(), cdf.cend(), x);
    if (it1 == cdf.cbegin()) {
      m_tab[i] = m_E.front();
      continue;
    }
    const auto it0 = std::prev(it1);
    const double x0 = *it0;
    const double x1 = *it1;
    const double y0 = m_E[it0 - cdf.cbegin()];
    const double y1 = m_E[it1 - cdf.cbegin()];
    const double f0 = (x - x0) / (x1 - x0);
    const double f1 = 1. - f0;
    m_tab[i] = f0 * y0 + f1 * y1;
  } 
  return true;
}

Referenced by GetClusterDensity(), GetStoppingPower(), and NewTrack().

GetCluster()

bool Garfield::TrackBichsel::GetCluster ( double & xc, double & yc, double & zc, double & tc, int & nc, double & ec, double & extra )

overridevirtual

Get the next "cluster" (ionising collision of the charged particle).

Parameters

xc,yc,zc	coordinates of the collision
tc	time of the collision
nc	number of electrons produced
ec	deposited energy
extra	additional information (not always implemented)

Implements Garfield::Track.

Definition at line 451 of file TrackBichsel.cc.

                                                                              {
  xc = yc = zc = tc = ec = extra = 0.;
  ne = 0;
  if (m_clusters.empty()) return false;
  // Increment the cluster index.
  if (m_cluster < m_clusters.size()) {
    ++m_cluster;
  } else if (m_cluster > m_clusters.size()) {
    m_cluster = 0;
  } 
  if (m_cluster >= m_clusters.size()) return false;
 
  xc = m_clusters[m_cluster].x;
  yc = m_clusters[m_cluster].y;
  zc = m_clusters[m_cluster].z;
  tc = m_clusters[m_cluster].t;
  ec = m_clusters[m_cluster].energy;
  return true;
}

GetClusterDensity()

double Garfield::TrackBichsel::GetClusterDensity ( )

overridevirtual

Get the cluster density (number of ionizing collisions per cm or inverse mean free path for ionization).

Reimplemented from Garfield::Track.

Definition at line 472 of file TrackBichsel.cc.

                                       {
 
  if (m_isChanged) {
    if (!ComputeCrossSection()) return 0.;
    m_isChanged = false;
  }
  return m_imfp;
}

GetClusters()

const std::vector< Cluster > & Garfield::TrackBichsel::GetClusters ( ) const

inline

Definition at line 36 of file TrackBichsel.hh.

{ return m_clusters; }

GetStoppingPower()

double Garfield::TrackBichsel::GetStoppingPower ( )

overridevirtual

Get the stopping power (mean energy loss [eV] per cm).

Reimplemented from Garfield::Track.

Definition at line 481 of file TrackBichsel.cc.

                                      {
 
  if (m_isChanged) {
    if (!ComputeCrossSection()) return 0.;
    m_isChanged = false;
  }
  return m_dEdx;
}

Initialise()

bool Garfield::TrackBichsel::Initialise

(

)

Definition at line 31 of file TrackBichsel.cc.

                              {
 
  std::cout << m_className << "::Initialize:\n";
  // Reset the tables.
  m_E.fill(0.);
  m_dfdE.fill(0.);
  m_eps1.fill(0.);
  m_eps2.fill(0.);
  m_int.fill(0.);
  m_k1.fill(0.);
 
  MediumSilicon si;
  m_density = si.GetNumberDensity();
  // Conversion from loss function to oscillator strength density.
  m_conv = ElectronMass / (2 * Pi2 * FineStructureConstant * pow(HbarC, 3) * m_density);
  // Number of bins for each factor of 2 in energy
  constexpr unsigned int n2 = 64;
  const double u = log(2.) / n2;
  const double um  = exp(u);
  constexpr double kEdge = 1839.;
  const double ken = log(kEdge / 1.5) / u;
  m_E[0] = kEdge / pow(2, ken / n2);
  if (m_debug) std::cout << "  Bin       Energy [eV]\n";
  for (size_t j = 0; j < NEnergyBins; ++j) {
    m_E[j + 1] = m_E[j] * um;
    if (!m_debug) continue;
    if ((j + 1) % 50 == 0) std::printf(" %4zu %16.8f\n", j + 1, m_E[j]);
  }
 
  // Read in the complex dielectric function (epsilon). This is used 
  // for the cross section of small momentum transfer excitations.
  std::string path = "";
  auto installdir = std::getenv("GARFIELD_INSTALL");
  if (!installdir) {
    std::cerr << "    Environment variable GARFIELD_INSTALL not set.\n";
  } else {
    path = std::string(installdir) + "/share/Garfield/Data/";
  }
  std::ifstream infile;
  std::cout << "    Reading dielectric function.\n";
  infile.open(path + "heps.tab", std::ios::in);
  if (!infile) {
    std::cerr << "    Could not open heps.tab.\n";
    return false;
  }
  bool ok = true;
  for (std::string line; std::getline(infile, line);) {
    ltrim(line);
    // Skip comments and empty lines.
    if (IsComment(line)) continue;
    std::istringstream data(line);
    size_t j;
    double energy, eps1, eps2, lossFunction;
    data >> j >> energy >> eps1 >> eps2 >> lossFunction;
    if (j < 1 || j > NEnergyBins) {
      std::cerr << "    Index out of range.\n";
      ok = false;
      break;
    }
    --j;
    m_eps1[j] = eps1;
    m_eps2[j] = eps2;
    m_dfdE[j] = lossFunction * m_conv * m_E[j];
  }
  infile.close();
  if (!ok) {
    std::cerr << "    Error reading dielectric function.\n";
    return false;
  }
 
  // Read in a table of the generalized oscillator strength density 
  // integrated over the momentum transfer K, corresponding to the 
  // function A(E) in Equation (2.11) in (Bichsel, 1988).
  // These values are used for calculating the cross section for 
  // longitudinal excitations of K and L shell electrons 
  // with large momentum transfer.
  std::cout << "    Reading K and L shell calculations.\n";
  infile.open(path + "macom.tab", std::ios::in);
  if (!infile) {
    std::cerr << "    Could not open macom.tab.\n";
    return false;
  }
  for (std::string line; std::getline(infile, line);) {
    ltrim(line);
    if (IsComment(line)) continue;
    std::istringstream data(line);
    size_t j;
    double energy, val;
    data >> j >> energy >> val;
    if (j < 1 || j > NEnergyBins) {
      std::cerr << "    Index out of range.\n";
      ok = false;
      break;
    }
    --j;
    m_int[j] = val;
  }
  infile.close();
  if (!ok) {
    std::cerr << "    Error reading K/L shell data.\n";
    return false;
  }
 
  // Read in a table of the generalized oscillator strength density 
  // for M shell electrons integrated over the momentum transfer K.
  // Based on equations in the appendix in Emerson et al., 
  // Phys. Rev. B 7 (1973), 1798 (DOI 10.1103/PhysRevB.7.1798)
  std::cout << "    Reading M shell calculations.\n";
  infile.open(path + "emerc.tab", std::ios::in);
  if (!infile) {
    std::cerr << "    Could not open emerc.tab.\n";
    return false;
  }
  for (std::string line; std::getline(infile, line);) {
    ltrim(line);
    if (IsComment(line)) continue;
    std::istringstream data(line);
    size_t j;
    double energy, a, k1;
    data >> j >> energy >> a >> k1;
    if (j < 1 || j > NEnergyBins) {
      std::cerr << "    Index out of range.\n";
      ok = false;
      break;
    }
    --j;
    m_int[j] = a;
    m_k1[j] = k1;
    if (!m_debug) continue;
    if (j < 19 || j > 30) continue;
    std::printf(" %4zu %11.2f %11.2f %12.6f %12.6f\n", 
                j + 1, m_E[j], energy, m_int[j], m_k1[j]); 
  }
  infile.close();
  if (!ok) {
    std::cerr << "    Error reading M shell data.\n";
    return false;
  }
 
  double s0 = 0.;
  double s1 = 0.;
  double avI = 0.;
  for (size_t j = 0; j < NEnergyBins; ++j) {
    const auto logE = log(m_E[j]);
    const double dE = m_E[j + 1] - m_E[j];
    s0 += m_dfdE[j] * dE;
    s1 += m_dfdE[j] * dE * m_E[j];
    avI += m_dfdE[j] * dE * logE;
  }
  if (m_debug) {
    std::printf("    S0  = %15.5f\n", s0);
    std::printf("    S1  = %15.5f\n", s1);
    std::printf("    lnI = %15.5f\n", avI);
  }
  m_initialised = true;
  return true;
}

Referenced by TrackBichsel().

NewTrack()

bool Garfield::TrackBichsel::NewTrack ( const double x0, const double y0, const double z0, const double t0, const double dx0, const double dy0, const double dz0 )

overridevirtual

Calculate a new track starting from (x0, y0, z0) at time t0 in direction (dx0, dy0, dz0).

Implements Garfield::Track.

Definition at line 358 of file TrackBichsel.cc.

                                              {
 
  // Reset the list of clusters.
  m_clusters.clear();
  m_cluster = 0;
 
  // Make sure a sensor has been defined.
  if (!m_sensor) {
    std::cerr << m_className << "::NewTrack: Sensor is not defined.\n";
    return false;
  }
 
  // If not yet done, compute the cross-section table.
  if (m_isChanged) {
    if (!ComputeCrossSection()) {
      std::cerr << m_className << "::NewTrack:\n"
                << "    Could not calculate cross-section table.\n";
      return false;
    }
    m_isChanged = false;
  }
 
  // Make sure we are inside a medium.
  Medium* medium = m_sensor->GetMedium(x0, y0, z0);
  if (!medium || !medium->IsIonisable()) {
    std::cerr << m_className << "::NewTrack:\n"
              << "    No ionisable medium at initial position.\n";
    return false;
  }
 
  // Check if the medium is silicon.
  if (medium->GetName() != "Si") {
    std::cerr << m_className << "::NewTrack: Medium is not silicon.\n";
    return false;
  }
 
  // Normalise the direction vector.
  double dx = dx0, dy = dy0, dz = dz0;
  const double d = sqrt(dx0 * dx0 + dy0 * dy0 + dz0 * dz0);
  if (d < Small) {
    // In case of a null vector, choose a random direction.
    RndmDirection(dx, dy, dz);
  } else {
    dx = dx0 / d;
    dy = dy0 / d;
    dz = dz0 / d;
  }
  const double dt = m_speed > 0. ? 1. / m_speed : 0.;
 
  double x = x0;
  double y = y0;
  double z = z0;
  double t = t0;
  const double mfp = 1. / m_imfp;
  double ekin = GetKineticEnergy();
  while (ekin > 0.) {
    const double step = -log(RndmUniformPos()) * mfp;
    x += dx * step;
    y += dy * step;
    z += dz * step;
    t += dt * step;
 
    medium = m_sensor->GetMedium(x, y, z);
    if (!medium || !medium->IsIonisable() || medium->GetName() != "Si") {
      if (m_debug) {
        std::cout << m_className << "::NewTrack: Particle left the medium.\n";
      }
      break;
    }
 
    Cluster cluster;
    cluster.x = x;
    cluster.y = y;
    cluster.z = z;
    cluster.t = t;
    const double u = NCdfBins * RndmUniform();
    const size_t j = static_cast<size_t>(std::floor(u));
    if (j == 0) {
      cluster.energy = u * m_tab.front();
    } else if (j >= NCdfBins) {
      cluster.energy = m_tab.back();
    } else {
      cluster.energy = m_tab[j - 1] + (u - j) * (m_tab[j] - m_tab[j - 1]);
    }
    ekin -= cluster.energy;
    m_clusters.push_back(std::move(cluster));
  }
  m_cluster = m_clusters.size() + 2; 
  return true;
}

---

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

• TrackBichsel.hh
• TrackBichsel.cc