Garfield::TrackTrim Class Reference

#include <TrackTrim.hh>

Inheritance diagram for Garfield::TrackTrim:

Classes
struct	Cluster

Public Member Functions
	TrackTrim ()
	Default constructor.
	TrackTrim (Sensor *sensor)
	Constructor.
virtual	~TrackTrim ()
	Destructor.
bool	ReadFile (const std::string &file, const unsigned int nIons=0, const unsigned int nSkip=0)
	Load data from an EXYZ.txt file.
void	Print ()
	Print a summary of the available TRIM data.
void	SetCharge (const double q)
	Set the projectile charge [-].
void	SetWorkFunction (const double w)
	Set the W value [eV].
double	GetWorkFunction () const
	Get the W value [eV].
void	SetFanoFactor (const double f)
	Set the Fano factor.
void	UnsetFanoFactor ()
	Use the default Fano factor.
double	GetFanoFactor () const
	Get the Fano factor.
void	EnableStepSizeLimit (const double dmax)
	Set a max. distance between subsequent steps.
void	DisableStepSizeLimit ()
	Do not reduce the steps with respect to the TRIM input file (default).
void	EnableEnergyLossLimit (const double emax)
	Set a max. energy loss per cluster.
void	DisableEnergyLossLimit ()
void	SetParticle (const std::string &part) override
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
Public Member Functions inherited from Garfield::Track
	Track ()=delete
	Default constructor.
	Track (const std::string &name)
	Constructor.
virtual	~Track ()
	Destructor.
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 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.

Protected Member Functions
void	AddIon (const std::vector< float > &x, const std::vector< float > &y, const std::vector< float > &z, const std::vector< float > &dedx, const std::vector< float > &ekin)
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
double	m_work = -1.
	Work function [eV] of the target.
bool	m_fset = false
	Has the Fano factor been set?
double	m_fano = -1.
	Fano factor [-] of the target.
double	m_ekin = 0.
	Projectile energy [eV].
std::vector< std::vector< std::array< float, 6 > > >	m_ions
	List of tracks.
size_t	m_ion = 0
	Index of the current track.
std::vector< Cluster >	m_clusters
	Clusters on the current track.
size_t	m_cluster = 0
	Index of the next cluster to be returned.
double	m_maxStepSize = -1.
	Step size limit.
double	m_maxLossPerStep = -1.
	Energy loss limit per step.
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

Additional Inherited Members
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)

Detailed Description

Generate tracks based on TRIM output files.

• http://www.srim.org

Definition at line 14 of file TrackTrim.hh.

Constructor & Destructor Documentation

TrackTrim() [1/2]

Garfield::TrackTrim::TrackTrim ( )

inline

Default constructor.

Definition at line 17 of file TrackTrim.hh.

: TrackTrim(nullptr) {}

Referenced by TrackTrim().

TrackTrim() [2/2]

Garfield::TrackTrim::TrackTrim

(

Sensor *

sensor

)

Constructor.

Definition at line 90 of file TrackTrim.cc.

                                   : Track("Trim") {
  m_sensor = sensor; 
  m_q = 1.;
}

~TrackTrim()

virtual Garfield::TrackTrim::~TrackTrim ( )

inlinevirtual

Destructor.

Definition at line 21 of file TrackTrim.hh.

{}

Member Function Documentation

AddIon()

void Garfield::TrackTrim::AddIon ( const std::vector< float > & x, const std::vector< float > & y, const std::vector< float > & z, const std::vector< float > & dedx, const std::vector< float > & ekin )

protected

Definition at line 194 of file TrackTrim.cc.

                                                     {
 
  const size_t nPoints = x.size();
  if (nPoints < 2) return;
  std::vector<std::array<float, 6> > path;
  for (size_t i = 0; i < nPoints - 1; ++i) {
    const float dx = x[i + 1] - x[i];
    const float dy = y[i + 1] - y[i];
    const float dz = z[i + 1] - z[i];
    const float dmag = sqrt(dx * dx + dy * dy + dz * dz);
    const float scale = dmag > 0. ? 1. / dmag : 1.;
    float eloss = 0.;
    if (i == 0 && dedx[i] > 10. * dedx[i + 1]) {
      eloss = dmag * dedx[i + 1];
    } else { 
      eloss = dmag * dedx[i];
    }
    const float dekin = ekin[i] - ekin[i + 1];
    if (dekin > 0.) eloss = std::min(eloss, dekin); 
    path.push_back({dx * scale, dy * scale, dz * scale, dmag, eloss, ekin[i]});
  }
  m_ions.push_back(std::move(path));
}

Referenced by ReadFile().

DisableEnergyLossLimit()

void Garfield::TrackTrim::DisableEnergyLossLimit ( )

inline

Do not limit the energy loss per cluster (default). Each segment in the TRIM input file corresponds to one cluster.

Definition at line 53 of file TrackTrim.hh.

{ m_maxLossPerStep = -1.; }

DisableStepSizeLimit()

void Garfield::TrackTrim::DisableStepSizeLimit ( )

inline

Do not reduce the steps with respect to the TRIM input file (default).

Definition at line 48 of file TrackTrim.hh.

{ m_maxStepSize = -1.; }

EnableEnergyLossLimit()

void Garfield::TrackTrim::EnableEnergyLossLimit ( const double emax )

inline

Set a max. energy loss per cluster.

Definition at line 50 of file TrackTrim.hh.

{ m_maxLossPerStep = emax; }

EnableStepSizeLimit()

void Garfield::TrackTrim::EnableStepSizeLimit ( const double dmax )

inline

Set a max. distance between subsequent steps.

Definition at line 46 of file TrackTrim.hh.

{ m_maxStepSize = dmax; }

GetCluster()

bool Garfield::TrackTrim::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 393 of file TrackTrim.cc.

                                                                           {
  if (m_debug) {
    std::cout << m_className << "::GetCluster: Cluster " << m_cluster
              << " of " << m_clusters.size() << "\n";
  }
  // Stop if we have exhausted the list of clusters.
  if (m_cluster >= m_clusters.size()) return false;
 
  const auto& cluster = m_clusters[m_cluster];
  xcls = cluster.x;
  ycls = cluster.y;
  zcls = cluster.z;
  tcls = cluster.t;
 
  n = cluster.n;
  e = cluster.energy;
  extra = cluster.ekin;
  // Move to the next cluster.
  ++m_cluster;
  return true;
}

GetClusters()

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

inline

Definition at line 70 of file TrackTrim.hh.

                                                {
    return m_clusters;
  }

GetFanoFactor()

double Garfield::TrackTrim::GetFanoFactor ( ) const

inline

Get the Fano factor.

Definition at line 43 of file TrackTrim.hh.

{ return m_fano; }

GetWorkFunction()

double Garfield::TrackTrim::GetWorkFunction ( ) const

inline

Get the W value [eV].

Definition at line 34 of file TrackTrim.hh.

{ return m_work; }

NewTrack()

bool Garfield::TrackTrim::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 243 of file TrackTrim.cc.

                                                                           {
 
  // TRMGEN - Generates TRIM clusters
 
  if (m_ions.empty()) {
    std::cerr << m_className << "::NewTrack: No TRIM data present.\n";
    return false;
  }
 
  if (m_ion >= m_ions.size()) {
    // Rewind.
    std::cout << m_className << "::NewTrack: Rewinding.\n";
    m_ion = 0;
  }
 
  // Verify that a sensor has been set.
  if (!m_sensor) {
    std::cerr << m_className << "::NewTrack: Sensor is not defined.\n";
    return false;
  }
 
  // Set the initial direction.
  std::array<double, 3> v = {dx0, dy0, dz0};
  const double d0 = sqrt(dx0 * dx0 + dy0 * dy0 + dz0 * dz0);
  if (d0 < Small) {
    if (m_debug) {
      std::cout << m_className << "::NewTrack: Sampling initial direction.\n";
    }
    // Null vector. Sample the direction isotropically.
    RndmDirection(v[0], v[1], v[2]);
  } else {
    const double scale = 1. / d0;
    for (size_t i = 0; i < 3; ++i) v[i] *= scale;
  }
 
  Medium* medium = m_sensor->GetMedium(x0, y0, z0);
  if (!medium) {
    std::cerr << m_className << "::NewTrack: No medium at initial point.\n";
    return false;
  }
  double w = m_work > 0. ? m_work : medium->GetW();
  // Warn if the W value is not defined.
  if (w < Small) {
    std::cerr << m_className << "::NewTrack: "
              << "Work function at initial point is not defined.\n";
  }
  double fano = m_fset ? m_fano : medium->GetFanoFactor();
  fano = std::max(fano, 0.);
 
  // Charge-over-mass ratio.
  const double qoverm = m_mass > 0. ? m_q / m_mass : 0.;
 
  // Plot.
  if (m_viewer) PlotNewTrack(x0, y0, z0);
 
  // Reset the cluster count.
  m_cluster = 0;
  m_clusters.clear();
 
  // Pool of unused energy
  double epool = 0.0;
 
  const double ekin0 = GetKineticEnergy();
 
  const auto& path = m_ions[m_ion];
  const size_t nPoints = path.size();
  std::array<double, 3> x = {x0, y0, z0};
  double t = t0;
  for (size_t i = 0; i < nPoints; ++i) {
    // Skip points with kinetic energy below the initial one set by the user.
    double ekin = path[i][5];
    if (ekin > ekin0) continue;
    // Rotate the direction vector.
    if (i == 0) {
      Rotate({1, 0, 0}, {path[i][0], path[i][1], path[i][2]}, v);
    } else {
      Rotate({path[i - 1][0], path[i - 1][1], path[i - 1][2]},
             {path[i][0], path[i][1], path[i][2]}, v);
    }
    // Get the step length and energy loss.
    double dmag = path[i][3];
    double eloss = path[i][4];
    std::array<double, 3> d = {dmag * v[0], dmag * v[1], dmag * v[2]};
    // Subdivide the step if necessary.
    unsigned int nSteps = 1;
    if (m_maxStepSize > 0. && dmag > m_maxStepSize) {
      nSteps = std::ceil(dmag / m_maxStepSize);
    }
    if (m_maxLossPerStep > 0. && eloss > nSteps * m_maxLossPerStep) {
      nSteps = std::ceil(eloss / m_maxLossPerStep);
    }
    if (nSteps > 1) {
      const double scale = 1. / nSteps;
      for (size_t j = 0; j < 3; ++j) d[j] *= scale;
      dmag *= scale;
      eloss *= scale;
    }
    // Compute the particle velocity.
    const double vmag = Speed(ekin, m_mass);
    // Compute the timestep.
    const double dt = vmag > 0. ? dmag / vmag : 0.;
    for (unsigned int j = 0; j < nSteps; ++j) {
      Cluster cluster;
      if (m_sensor->HasMagneticField()) {
        double bx = 0., by = 0., bz = 0.;
        int status = 0;
        m_sensor->MagneticField(x[0], x[1], x[2], bx, by, bz, status);
        d = StepBfield(dt, qoverm, vmag, bx, by, bz, v);
      }
      cluster.x = x[0] + d[0];
      cluster.y = x[1] + d[1];
      cluster.z = x[2] + d[2];
      cluster.t = t + dt;
      x = {cluster.x, cluster.y, cluster.z};
      t = cluster.t;
      // Is this point inside an ionisable medium?
      medium = m_sensor->GetMedium(cluster.x, cluster.y, cluster.z);
      if (!medium || !medium->IsIonisable()) continue;
      if (m_work < Small) w = medium->GetW();
      if (w < Small) continue;
      if (fano < Small) {
        // No fluctuations.
        cluster.n = int((eloss + epool) / w);
        cluster.energy = w * cluster.n;
      } else {
        double ec = eloss + epool;
        cluster.n = 0;
        cluster.energy = 0.0;
        while (true) {
          const double er = RndmHeedWF(w, fano);
          if (er > ec) break;
          cluster.n++;
          cluster.energy += er;
          ec -= er;
        }
      }
      // TODO
      cluster.ekin = i < nPoints - 1 ? path[i + 1][5] : ekin;
      epool += eloss - cluster.energy;
      if (cluster.n == 0) continue;
      m_clusters.push_back(std::move(cluster));
      if (m_viewer) PlotCluster(cluster.x, cluster.y, cluster.z);
    } 
  }
  // Move to the next ion in the list.
  ++m_ion;
  return true;
}

Print()

void Garfield::TrackTrim::Print

(

)

Print a summary of the available TRIM data.

Definition at line 222 of file TrackTrim.cc.

                      {
  std::cout << m_className << "::Print:\n";
  if (m_ions.empty()) {
    std::cerr << "    No TRIM data present.\n";
    return;
  }
  std::cout << "    Projectile: " << m_particleName << ", "
            << m_ekin * 1.e-3 << " keV\n"
            << "    Number of tracks: " << m_ions.size() << "\n";
  if (m_work > 0.) {
    std::cout << "    Work function: " << m_work << " eV\n";
  } else {
    std::cout << "    Work function: Automatic\n";
  }
  if (m_fset) {
    std::cout << "    Fano factor: " << m_fano << "\n";
  } else {
    std::cout << "    Fano factor: Automatic\n";
  }
}

ReadFile()

bool Garfield::TrackTrim::ReadFile	(	const std::string &	file,
		const unsigned int	nIons = 0,
		const unsigned int	nSkip = 0 )

Load data from an EXYZ.txt file.

Definition at line 99 of file TrackTrim.cc.

                                                                             {
 
  // TRMREE - Reads the TRIM EXYZ file.
 
  // Reset.
  m_ekin = 0.;
  m_ions.clear();
  m_ion = 0;
  m_clusters.clear();
  m_cluster = 0;
 
  std::ifstream infile(filename);
  if (!infile) {
    std::cerr << m_className << "::ReadFile:\n"
              << "    Unable to open the EXYZ file (" << filename << ").\n";
    return false;
  }
 
  constexpr double Angstrom = 1.e-8;
  unsigned int nRead = 0;
  unsigned int ionNumber = 0;
  bool header = true;
  double mass = 0.;
  std::vector<float> x;
  std::vector<float> y;
  std::vector<float> z;
  std::vector<float> dedx;
  std::vector<float> ekin;
  for (std::string line; std::getline(infile, line);) {
    if (line.find("------- ") != std::string::npos) {
      // Reached the end of the header.
      header = false;
      continue;
    } else if (header) {
      if (line.find("Ion Data: ") != std::string::npos) {
        // Read the next line.
        std::getline(infile, line);
        auto words = tokenize(line);
        if (words.size() >= 3) {
          m_particleName = words[0];
          mass = std::stod(words[1]);
          auto pos = words[2].find("keV");
          if (pos != std::string::npos) {
            m_ekin = 1.e3 * std::stod(words[2].substr(0, pos));
          }
        } 
      } 
      // Otherwise, skip the header.
      continue;
    } 
    auto words = tokenize(line);
    if (words.size() < 6) {
      std::cerr << m_className << "::ReadFile: Unexpected line:\n"
                << line << "\n";
      continue;
    }
    if (ionNumber != std::stoul(words[0])) {
      // New ion.
      if (ionNumber > 0) {
        if (nRead >= nSkip) AddIon(x, y, z, dedx, ekin);
        x.clear();
        y.clear();
        z.clear();
        dedx.clear(); 
        ekin.clear();
        ++nRead;
        // Stop if we are done reading the requested number of ions.
        if (nIons > 0 && m_ions.size() >= nIons) break;
      }
      ionNumber = std::stoi(words[0]);
    }
    if (nRead < nSkip) continue;
    // Convert coordinates to cm.
    x.push_back(std::stof(words[2]) * Angstrom);
    y.push_back(std::stof(words[3]) * Angstrom);
    z.push_back(std::stof(words[4]) * Angstrom);
    // Convert stopping power from eV/A to eV/cm.
    dedx.push_back(std::stof(words[5]) / Angstrom);
    // Convert ion energy from keV to eV.
    ekin.push_back(std::stof(words[1]) * 1.e3);
  }
  infile.close();
  AddIon(x, y, z, dedx, ekin);
  std::cout << m_className << "::ReadFile: Read energy vs position for " 
            << m_ions.size() << " ions.\n";
  if (m_ekin > 0. && mass > 0.) {
    std::cout << "    Initial kinetic energy set to " 
              << m_ekin * 1.e-3 << " keV. Mass number: " << mass << ".\n";
    m_mass = AtomicMassUnitElectronVolt * mass;
    SetKineticEnergy(m_ekin);
  }
  return true;
}

SetCharge()

void Garfield::TrackTrim::SetCharge ( const double q )

inline

Set the projectile charge [-].

Definition at line 30 of file TrackTrim.hh.

{ if (fabs(q) > 0.) m_q = q; }

SetFanoFactor()

void Garfield::TrackTrim::SetFanoFactor ( const double f )

inline

Set the Fano factor.

Definition at line 36 of file TrackTrim.hh.

                                     {
    m_fano = f;
    m_fset = true;
  }

SetParticle()

void Garfield::TrackTrim::SetParticle ( const std::string & part )

overridevirtual

Set the type of charged particle.

• electron,e-
• positron,e+
• muon,mu-
• mu+
• pion,pi-
• pi+
• kaon,K-
• K+
• proton,p
• anti-proton,p-bar
• deuteron,d
• alpha

Reimplemented from Garfield::Track.

Definition at line 95 of file TrackTrim.cc.

                                              {
  std::cerr << m_className << "::SetParticle: Not applicable.\n";
}

SetWorkFunction()

void Garfield::TrackTrim::SetWorkFunction ( const double w )

inline

Set the W value [eV].

Definition at line 32 of file TrackTrim.hh.

{ m_work = w; }

UnsetFanoFactor()

void Garfield::TrackTrim::UnsetFanoFactor ( )

inline

Use the default Fano factor.

Definition at line 41 of file TrackTrim.hh.

{ m_fset = false; }

Member Data Documentation

m_cluster

size_t Garfield::TrackTrim::m_cluster = 0

protected

Index of the next cluster to be returned.

Definition at line 92 of file TrackTrim.hh.

Referenced by GetCluster(), NewTrack(), and ReadFile().

m_clusters

std::vector<Cluster> Garfield::TrackTrim::m_clusters

protected

Clusters on the current track.

Definition at line 90 of file TrackTrim.hh.

Referenced by GetCluster(), GetClusters(), NewTrack(), and ReadFile().

m_ekin

double Garfield::TrackTrim::m_ekin = 0.

protected

Projectile energy [eV].

Definition at line 83 of file TrackTrim.hh.

Referenced by Print(), and ReadFile().

m_fano

double Garfield::TrackTrim::m_fano = -1.

protected

Fano factor [-] of the target.

Definition at line 80 of file TrackTrim.hh.

Referenced by GetFanoFactor(), NewTrack(), Print(), and SetFanoFactor().

m_fset

bool Garfield::TrackTrim::m_fset = false

protected

Has the Fano factor been set?

Definition at line 78 of file TrackTrim.hh.

Referenced by NewTrack(), Print(), SetFanoFactor(), and UnsetFanoFactor().

m_ion

size_t Garfield::TrackTrim::m_ion = 0

protected

Index of the current track.

Definition at line 87 of file TrackTrim.hh.

Referenced by NewTrack(), and ReadFile().

m_ions

std::vector<std::vector<std::array<float, 6> > > Garfield::TrackTrim::m_ions

protected

List of tracks.

Definition at line 85 of file TrackTrim.hh.

Referenced by AddIon(), NewTrack(), Print(), and ReadFile().

m_maxLossPerStep

double Garfield::TrackTrim::m_maxLossPerStep = -1.

protected

Energy loss limit per step.

Definition at line 97 of file TrackTrim.hh.

Referenced by DisableEnergyLossLimit(), EnableEnergyLossLimit(), and NewTrack().

m_maxStepSize

double Garfield::TrackTrim::m_maxStepSize = -1.

protected

Step size limit.

Definition at line 95 of file TrackTrim.hh.

Referenced by DisableStepSizeLimit(), EnableStepSizeLimit(), and NewTrack().

m_work

double Garfield::TrackTrim::m_work = -1.

protected

Work function [eV] of the target.

Definition at line 76 of file TrackTrim.hh.

Referenced by GetWorkFunction(), NewTrack(), Print(), and SetWorkFunction().

---

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

• TrackTrim.hh
• TrackTrim.cc