Garfield::TrackBichsel Class Reference

#include <TrackBichsel.hh>

Inheritance diagram for Garfield::TrackBichsel:

Public Member Functions
	TrackBichsel ()
	Constructor.
virtual	~TrackBichsel ()
	Destructor.
virtual bool	NewTrack (const double x0, const double y0, const double z0, const double t0, const double dx0, const double dy0, const double dz0)
virtual bool	GetCluster (double &xcls, double &ycls, double &zcls, double &tcls, int &n, double &e, double &extra)
virtual double	GetClusterDensity ()
virtual double	GetStoppingPower ()
	Get the stopping power (mean energy loss [eV] per cm).
void	SetDataFile (const std::string &filename)
Public Member Functions inherited from Garfield::Track
	Track ()
	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.
virtual bool	NewTrack (const double x0, const double y0, const double z0, const double t0, const double dx0, const double dy0, const double dz0)=0
virtual bool	GetCluster (double &xcls, double &ycls, double &zcls, double &tcls, int &n, double &e, double &extra)=0
virtual double	GetClusterDensity ()
virtual double	GetStoppingPower ()
	Get the stopping power (mean energy loss [eV] per cm).
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)
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
bool	m_isElectron = false
std::string	m_particleName = "mu-"
Sensor *	m_sensor = nullptr
bool	m_isChanged = true
bool	m_usePlotting = false
ViewDrift *	m_viewer = nullptr
bool	m_debug = false
int	m_plotId = -1

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 14 of file TrackBichsel.hh.

Constructor & Destructor Documentation

TrackBichsel()

Garfield::TrackBichsel::TrackBichsel

(

)

Constructor.

Definition at line 17 of file TrackBichsel.cc.

    : Track(), m_speed(SpeedOfLight * m_bg / sqrt(1. + m_bg * m_bg)) {
  m_className = "TrackBichsel";
}

~TrackBichsel()

virtual Garfield::TrackBichsel::~TrackBichsel ( )

inlinevirtual

Destructor.

Definition at line 19 of file TrackBichsel.hh.

{}

Member Function Documentation

GetCluster()

bool Garfield::TrackBichsel::GetCluster ( double &  xcls, double &  ycls, double &  zcls, double &  tcls, int &  n, double &  e, double &  extra  )

virtual

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

Parameters

xcls,ycls,zcls	coordinates of the collision
tcls	time of the collision
n	number of electrons produced
e	deposited energy
extra	additional information (not always implemented)

Implements Garfield::Track.

Definition at line 96 of file TrackBichsel.cc.

                                                                              {
  if (!m_isInitialised || !m_isInMedium) return false;
 
  const double d = -log(RndmUniformPos()) / m_imfp;
  m_x += m_dx * d;
  m_y += m_dy * d;
  m_z += m_dz * d;
  m_t += d / m_speed;
 
  xcls = m_x;
  ycls = m_y;
  zcls = m_z;
  tcls = m_t;
  n = 0;
  e = 0.;
  extra = 0.;
 
  Medium* medium;
  if (!m_sensor->GetMedium(m_x, m_y, m_z, medium)) {
    m_isInMedium = false;
    if (m_debug) {
      std::cout << m_className << "::GetCluster: Particle left the medium.\n";
    }
    return false;
  }
 
  if (medium->GetName() != "Si" || !medium->IsIonisable()) {
    m_isInMedium = false;
    if (m_debug) {
      std::cout << m_className << "::GetCluster: Particle left the medium.\n";
    }
    return false;
  }
 
  const double u = m_nCdfEntries * RndmUniform();
  const int j = int(u);
  if (j == 0) {
    e = 0. + u * m_cdf[0][m_iCdf];
  } else if (j >= m_nCdfEntries) {
    e = m_cdf[m_nCdfEntries - 1][m_iCdf];
  } else {
    e = m_cdf[j - 1][m_iCdf] +
        (u - j) * (m_cdf[j][m_iCdf] - m_cdf[j - 1][m_iCdf]);
  }
 
  return true;
}

GetClusterDensity()

double Garfield::TrackBichsel::GetClusterDensity ( )

virtual

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

Reimplemented from Garfield::Track.

Definition at line 145 of file TrackBichsel.cc.

                                       {
  constexpr unsigned int nEntries = 38;
  constexpr std::array<double, nEntries> tabBg = {
      {0.316,   0.398,   0.501,   0.631,    0.794,    1.000,   1.259,   1.585,
       1.995,   2.512,   3.162,   3.981,    5.012,    6.310,   7.943,   10.000,
       12.589,  15.849,  19.953,  25.119,   31.623,   39.811,  50.119,  63.096,
       79.433,  100.000, 125.893, 158.489,  199.526,  251.189, 316.228, 398.107,
       501.187, 630.958, 794.329, 1000.000, 1258.926, 1584.894}};
  constexpr std::array<double, nEntries> tabImfp = {
      {30.32496, 21.14965, 15.06555, 11.05635, 8.43259, 6.72876, 5.63184,
       4.93252,  4.49174,  4.21786,  4.05090,  3.95186, 3.89531, 3.86471,
       3.84930,  3.84226,  3.83952,  3.83887,  3.83912, 3.83970, 3.84035,
       3.84095,  3.84147,  3.84189,  3.84223,  3.84249, 3.84269, 3.84283,
       3.84293,  3.84300,  3.84304,  3.84308,  3.84310, 3.84311, 3.84312,
       3.84313,  3.84313,  3.84314}};
 
  if (m_isChanged) m_bg = GetBetaGamma();
 
  if (m_bg < tabBg.front()) {
    if (m_debug) {
      std::cerr << m_className << "::GetClusterDensity:\n"
                << "    Bg is below the tabulated range.\n";
    }
    return tabImfp.front() * 1.e4;
  } else if (m_bg > tabBg.back()) {
    return tabImfp.back() * 1.e4;
  }
 
  // Locate the requested energy in the table.
  const auto it1 = std::upper_bound(tabBg.cbegin(), tabBg.cend(), m_bg);
  if (it1 == tabBg.cbegin()) return 1.e4 * tabImfp.front();
  const auto it0 = std::prev(it1);
  const double x0 = *it0;
  const double x1 = *it1;
  const double y0 = tabImfp[it0 - tabBg.cbegin()];
  const double y1 = tabImfp[it1 - tabBg.cbegin()];
  const double tol = 1.e-6 * (x1 - x0);
  if (fabs(m_bg - x0) < tol) return y0 * 1.e4;
  if (fabs(m_bg - x1) < tol) return y1 * 1.e4;
 
  // Log-log interpolation
  const double lnx0 = log(x0);
  const double lnx1 = log(x1);
  const double lny0 = log(y0);
  const double lny1 = log(y1);
  const double d = lny0 + (log(m_bg) - lnx0) * (lny1 - lny0) / (lnx1 - lnx0);
  return 1.e4 * exp(d);
}

Referenced by NewTrack().

GetStoppingPower()

double Garfield::TrackBichsel::GetStoppingPower ( )

virtual

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

Reimplemented from Garfield::Track.

Definition at line 194 of file TrackBichsel.cc.

                                      {
  constexpr unsigned int nEntries = 51;
  constexpr std::array<double, nEntries> tabBg = {
      {0.316,     0.398,     0.501,    0.631,     0.794,     1.000,
       1.259,     1.585,     1.995,    2.512,     3.162,     3.981,
       5.012,     6.310,     7.943,    10.000,    12.589,    15.849,
       19.953,    25.119,    31.623,   39.811,    50.119,    63.096,
       79.433,    100.000,   125.893,  158.489,   199.526,   251.189,
       316.228,   398.107,   501.187,  630.958,   794.329,   1000.000,
       1258.926,  1584.894,  1995.263, 2511.888,  3162.280,  3981.074,
       5011.875,  6309.578,  7943.287, 10000.010, 12589.260, 15848.940,
       19952.640, 25118.880, 31622.800}};
  constexpr std::array<double, nEntries> tabdEdx = {
      {2443.71800, 1731.65600, 1250.93400, 928.69920, 716.37140, 578.28850,
       490.83670,  437.33820,  406.58490,  390.95170, 385.29000, 386.12000,
       391.07730,  398.53930,  407.39420,  416.90860, 426.63010, 436.30240,
       445.78980,  455.02530,  463.97370,  472.61410, 480.92980, 488.90240,
       496.51900,  503.77130,  510.65970,  517.19570, 523.39830, 529.29120,
       534.90670,  540.27590,  545.42880,  550.39890, 555.20800, 559.88820,
       564.45780,  568.93850,  573.34700,  577.69140, 581.99010, 586.25090,
       590.47720,  594.68660,  598.86880,  603.03510, 607.18890, 611.33250,
       615.46810,  619.59740,  623.72150}};
 
  if (m_isChanged) m_bg = GetBetaGamma();
 
  if (m_bg < tabBg.front()) {
    if (m_debug) {
      std::cerr << m_className << "::GetStoppingPower:\n"
                << "    Bg is below the tabulated range.\n";
    }
    return tabdEdx.front() * 1.e4;
  } else if (m_bg > tabBg.back()) {
    return tabdEdx.back() * 1.e4;
  }
 
  // Locate the requested energy in the table.
  const auto it1 = std::upper_bound(tabBg.cbegin(), tabBg.cend(), m_bg);
  if (it1 == tabBg.cbegin()) return 1.e4 * tabdEdx.front();
  const auto it0 = std::prev(it1);
  const double x0 = *it0;
  const double x1 = *it1;
  if (m_debug) {
    std::cout << m_className << "::GetStoppingPower:\n"
              << "    Bg = " << m_bg << "\n"
              << "    Interpolating between " << x0 << " and " << x1 << "\n";
  }
  const double y0 = tabdEdx[it0 - tabBg.cbegin()];
  const double y1 = tabdEdx[it1 - tabBg.cbegin()];
  const double tol = 1.e-6 * (x1 - x0);
  if (fabs(m_bg - x0) < tol) return y0 * 1.e4;
  if (fabs(m_bg - x1) < tol) return y1 * 1.e4;
 
  // Log-log interpolation
  const double lnx0 = log(x0);
  const double lnx1 = log(x1);
  const double lny0 = log(y0);
  const double lny1 = log(y1);
  const double dedx = lny0 + (log(m_bg) - lnx0) * (lny1 - lny0) / (lnx1 - lnx0);
  return 1.e4 * exp(dedx);
}

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  )

virtual

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

Implements Garfield::Track.

Definition at line 22 of file TrackBichsel.cc.

                                              {
  // Make sure a sensor has been defined.
  if (!m_sensor) {
    std::cerr << m_className << "::NewTrack: Sensor is not defined.\n";
    m_isInMedium = false;
    return false;
  }
 
  // If not yet done, load the cross-section table from file.
  if (!m_isInitialised) {
    if (!LoadCrossSectionTable(m_datafile)) {
      std::cerr << m_className << "::NewTrack:\n"
                << "    Cross-section table could not be loaded.\n";
      return false;
    }
    m_isInitialised = true;
  }
 
  // Make sure we are inside a medium.
  Medium* medium;
  if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
    std::cerr << m_className << "::NewTrack:\n"
              << "    No medium at initial position.\n";
    m_isInMedium = false;
    return false;
  }
 
  // Check if the medium is silicon.
  if (medium->GetName() != "Si") {
    std::cerr << m_className << "::NewTrack:\n"
              << "    Medium at initial position is not silicon.\n";
    m_isInMedium = false;
    return false;
  }
 
  // Check if primary ionisation has been enabled.
  if (!medium->IsIonisable()) {
    std::cerr << m_className << "::NewTrack:\n"
              << "    Medium at initial position is not ionisable.\n";
    m_isInMedium = false;
    return false;
  }
 
  m_isInMedium = true;
  m_x = x0;
  m_y = y0;
  m_z = z0;
  m_t = t0;
 
  // Normalise the direction vector.
  const double d = sqrt(dx0 * dx0 + dy0 * dy0 + dz0 * dz0);
  if (d < Small) {
    // In case of a null vector, choose a random direction.
    RndmDirection(m_dx, m_dy, m_dz);
  } else {
    m_dx = dx0 / d;
    m_dy = dy0 / d;
    m_dz = dz0 / d;
  }
 
  // If the particle properties have changed, update the cross-section table.
  if (m_isChanged) {
    m_bg = GetBetaGamma();
    m_imfp = GetClusterDensity();
    m_speed = SpeedOfLight * GetBeta();
    SelectCrossSectionTable();
    m_isChanged = false;
  }
 
  return true;
}

SetDataFile()

void Garfield::TrackBichsel::SetDataFile ( const std::string &  filename )

inline

Definition at line 30 of file TrackBichsel.hh.

{ m_datafile = filename; }

---

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

• TrackBichsel.hh
• TrackBichsel.cc