#include <TrackSrim.hh>

Inheritance diagram for Garfield::TrackSrim:

Classes
struct	cluster

Public Member Functions
	TrackSrim ()
	Constructor.

virtual	~TrackSrim ()
	Destructor.

void	SetWorkFunction (const double w)
	Set/get the W value [eV].

double	GetWorkFunction () const

void	SetFanoFactor (const double f)
	Set/get the Fano factor.

double	GetFanoFactor () const

void	SetDensity (const double density)
	Set/get the density [g/cm3] of the target medium.

double	GetDensity () const

void	SetAtomicMassNumbers (const double a, const double z)
	Set/get A and Z of the target medium.

void	GetAtomicMassMumbers (double &a, double &z) const

void	SetModel (const int m)

int	GetModel () const

void	EnableTransverseStraggling ()

void	DisableTransverseStraggling ()

void	EnableLongitudinalStraggling ()

void	DisableLongitudinalStraggling ()

void	EnablePreciseVavilov ()

void	DisablePreciseVavilov ()

void	SetTargetClusterSize (const int n)

int	GetTargetClusterSize () const

void	SetClustersMaximum (const int n)

int	SetClustersMaximum () const

bool	ReadFile (const std::string &file)

void	Print ()

void	PlotEnergyLoss ()

void	PlotRange ()

void	PlotStraggling ()

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)

Public Member Functions inherited from Garfield::Track
	Track ()
	Constructor.

virtual	~Track ()
	Destructor.

virtual void	SetParticle (const std::string &part)
	Set the type of particle.

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 ( $\beta = v/c$ ) 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

double	GetBetaGamma () const

double	GetBeta () const

double	GetGamma () const

double	GetMomentum () const

double	GetKineticEnergy () const

double	GetCharge () const
	Get the charge of the projectile.

double	GetMass () const
	Get the mass [eV / c2] of the projectile.

void	SetSensor (Sensor *s)

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)

void	DisablePlotting ()

void	EnableDebugging ()

void	DisableDebugging ()

Protected Member Functions
double	DedxEM (const double e) const

double	DedxHD (const double e) const

bool	PreciseLoss (const double step, const double estart, double &deem, double &dehd) const

bool	EstimateRange (const double ekin, const double step, double &stpmax)

bool	SmallestStep (double ekin, double de, double step, double &stpmin)

double	RndmEnergyLoss (const double ekin, const double de, const double step) const

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
bool	m_precisevavilov
	Use precise Vavilov generator.

bool	m_useTransStraggle
	Include transverse straggling.

bool	m_useLongStraggle
	Include longitudinal straggling.

bool	m_chargeset
	Charge gas been defined.

double	m_density
	Density [g/cm3].

double	m_work
	Work function [eV].

double	m_fano
	Fano factor [-].

double	m_q
	Charge of ion.

double	m_mass
	Mass of ion [MeV].

double	m_a
	A and Z of the gas.

double	m_z

int	m_maxclusters
	Maximum number of clusters allowed (infinite if 0)

std::vector< double >	m_ekin
	Energy in energy loss table [MeV].

std::vector< double >	m_emloss
	EM energy loss [MeV cm2/g].

std::vector< double >	m_hdloss
	Hadronic energy loss [MeV cm2/g].

std::vector< double >	m_range
	Projected range [cm].

std::vector< double >	m_transstraggle
	Transverse straggling [cm].

std::vector< double >	m_longstraggle
	Longitudinal straggling [cm].

unsigned int	m_currcluster
	Index of the next cluster to be returned.

unsigned int	m_model

int	m_nsize
	Targeted cluster size.

std::vector< cluster >	m_clusters

Protected Attributes inherited from Garfield::Track
std::string	m_className

double	m_q

int	m_spin

double	m_mass

double	m_energy

double	m_beta2

bool	m_isElectron

std::string	m_particleName

Sensor *	m_sensor

bool	m_isChanged

bool	m_usePlotting

ViewDrift *	m_viewer

bool	m_debug

int	m_plotId

Detailed Description

Generate tracks based on SRIM energy loss, range and straggling tables.

http://www.srim.org

Definition at line 11 of file TrackSrim.hh.

Constructor & Destructor Documentation

◆ TrackSrim()

Garfield::TrackSrim::TrackSrim ( )

Constructor.

Definition at line 122 of file TrackSrim.cc.

    : Track(),
      m_precisevavilov(false),
      m_useTransStraggle(true),
      m_useLongStraggle(false),
      m_chargeset(false),
      m_density(-1.0),
      m_work(-1.),
      m_fano(-1.),
      m_a(-1.), 
      m_z(-1.),
      m_maxclusters(-1),
      m_model(4),
      m_nsize(-1) {
  m_className = "TrackSrim";
  m_mass = -1.;
}

◆ ~TrackSrim()

virtual Garfield::TrackSrim::~TrackSrim ( )

inlinevirtual

Destructor.

Definition at line 17 of file TrackSrim.hh.

17{}

Member Function Documentation

◆ DedxEM()

double Garfield::TrackSrim::DedxEM ( const double e ) const

protected

Definition at line 513 of file TrackSrim.cc.

                                             {
  return Interpolate(e, m_ekin, m_emloss);
}

Referenced by NewTrack(), and PreciseLoss().

◆ DedxHD()

double Garfield::TrackSrim::DedxHD ( const double e ) const

protected

Definition at line 517 of file TrackSrim.cc.

                                             {
  return Interpolate(e, m_ekin, m_hdloss);
}

Referenced by NewTrack(), and PreciseLoss().

◆ DisableLongitudinalStraggling()

void Garfield::TrackSrim::DisableLongitudinalStraggling ( )

inline

Definition at line 46 of file TrackSrim.hh.

46{ m_useLongStraggle = false; }

◆ DisablePreciseVavilov()

void Garfield::TrackSrim::DisablePreciseVavilov ( )

inline

Definition at line 49 of file TrackSrim.hh.

49{ m_precisevavilov = false; }

◆ DisableTransverseStraggling()

void Garfield::TrackSrim::DisableTransverseStraggling ( )

inline

Definition at line 44 of file TrackSrim.hh.

44{ m_useTransStraggle = false; }

◆ EnableLongitudinalStraggling()

void Garfield::TrackSrim::EnableLongitudinalStraggling ( )

inline

Definition at line 45 of file TrackSrim.hh.

45{ m_useLongStraggle = true; }

◆ EnablePreciseVavilov()

void Garfield::TrackSrim::EnablePreciseVavilov ( )

inline

Definition at line 48 of file TrackSrim.hh.

48{ m_precisevavilov = true; }

◆ EnableTransverseStraggling()

void Garfield::TrackSrim::EnableTransverseStraggling ( )

inline

Definition at line 43 of file TrackSrim.hh.

43{ m_useTransStraggle = true; }

◆ EstimateRange()

bool Garfield::TrackSrim::EstimateRange	(	const double	ekin,
		const double	step,
		double &	stpmax
	)

protected

Definition at line 602 of file TrackSrim.cc.

                                              {
  // Find distance over which the ion just does not lose all its energy
  // ekin       : Kinetic energy [MeV]
  // step       : Step length as guessed [cm]
  // stpmax     : Maximum step
  // SRMDEZ
 
  const std::string hdr = m_className + "::EstimateRange: ";
  // Initial estimate
  stpmax = step;
 
  // Find the energy loss expected for the present step length.
  double st1 = step;
  double deem = 0., dehd = 0.;
  PreciseLoss(st1, ekin, deem, dehd);
  double de1 = deem + dehd;
  // Do nothing if this is ok
  if (de1 < ekin) {
    if (m_debug) std::cout << hdr << "Initial step OK.\n";
    return true;
  }
  // Find a smaller step for which the energy loss is less than EKIN.
  double st2 = 0.5 * step;
  double de2 = de1;
  const unsigned int nMaxIter = 20;
  for (unsigned int iter = 0; iter < nMaxIter; ++iter) {
    // See where we stand
    PreciseLoss(st2, ekin, deem, dehd);
    de2 = deem + dehd;
    // Below the kinetic energy: done
    if (de2 < ekin) break;
    // Not yet below the kinetic energy: new iteration.
    st1 = st2;
    de1 = de2;
    st2 *= 0.5;
  }
  if (de2 >= ekin) {
    std::cerr << hdr << "\n    Did not find a smaller step in " << nMaxIter 
              << " iterations. Abandoned.\n";
    stpmax = 0.5 * (st1 + st2);
    return false;
  }
  if (m_debug)
    printf("\tstep 1 = %g cm, de 1 = %g MeV\n\tstep 2 = %g cm, de 2 = %g MeV\n",
           st1, de1 - ekin, st2, de2 - ekin);
 
  // Now perform a bisection
  for (unsigned int iter = 0; iter < nMaxIter; ++iter) {
    // Avoid division by zero.
    if (de2 == de1) {
      if (m_debug) {
        std::cerr << hdr << "Bisection failed due to equal energy loss for "
                  << "two step sizes. Abandoned.\n";
      }
      stpmax = 0.5 * (st1 + st2);
      return false;
    }
    // Estimate step to give total energy loss.
    double st3;
    if ((fabs(de1 - ekin) < 0.01 * fabs(de2 - de1)) ||
        (fabs(de1 - ekin) > 0.99 * fabs(de2 - de1))) {
      st3 = 0.5 * (st1 + st2);
    } else {
      st3 = st1 - (st2 - st1) * (de1 - ekin) / (de2 - de1);
    }
    // See how well we are doing.
    PreciseLoss(st3, ekin, deem, dehd);
    const double de3 = deem + dehd;
    if (m_debug) printf("\tStep 1 = %g cm, dE 1 = %g MeV\n", st1, de1 - ekin);
    if (m_debug) printf("\tStep 2 = %g cm, dE 2 = %g MeV\n", st2, de2 - ekin);
    if (m_debug) printf("\tStep 3 = %g cm, dE 3 = %g MeV\n", st3, de3 - ekin);
    //  Update the estimates above and below.
    if (de3 > ekin) {
      st1 = st3;
      de1 = de3;
    } else {
      st2 = st3;
      de2 = de3;
    }
    // See whether we've converged.
    if (fabs(de3 - ekin) < 1e-3 * (fabs(de3) + fabs(ekin)) ||
        fabs(st1 - st2) < 1e-3 * (fabs(st1) + fabs(st2))) {
      stpmax = st1 - (st2 - st1) * (de1 - ekin) / (de2 - de1);
      return true;
    }
  }
  if (m_debug) {
    std::cout << hdr << "Bisection did not converge in " << nMaxIter 
              << " steps. Abandoned.\n";
  }
  stpmax = st1 - (st2 - st1) * (de1 - ekin) / (de2 - de1);
  return false;
}

Referenced by NewTrack().

◆ GetAtomicMassMumbers()

void Garfield::TrackSrim::GetAtomicMassMumbers	(	double &	a,
		double &	z
	)		const

inline

Definition at line 33 of file TrackSrim.hh.

◆ GetCluster()

bool Garfield::TrackSrim::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 1330 of file TrackSrim.cc.

                                                                           {
  if (m_debug) {
    printf("Current cluster: %d, array size: %ld", 
           m_currcluster, m_clusters.size());
  }
  // Stop if we have exhausted the list of clusters.
  if (m_currcluster >= m_clusters.size()) return false;
 
  xcls = m_clusters[m_currcluster].x;
  ycls = m_clusters[m_currcluster].y;
  zcls = m_clusters[m_currcluster].z;
  tcls = m_clusters[m_currcluster].t;
 
  n = m_clusters[m_currcluster].electrons;
  e = m_clusters[m_currcluster].ec;
  extra = m_clusters[m_currcluster].kinetic;
  // Move to next cluster
  ++m_currcluster;
  return true;
}

◆ GetDensity()

double Garfield::TrackSrim::GetDensity ( ) const

inline

Definition at line 27 of file TrackSrim.hh.

27{ return m_density; }

◆ GetFanoFactor()

double Garfield::TrackSrim::GetFanoFactor ( ) const

inline

Definition at line 24 of file TrackSrim.hh.

24{ return m_fano; }

◆ GetModel()

int Garfield::TrackSrim::GetModel ( ) const

inline

Definition at line 41 of file TrackSrim.hh.

41{ return m_model; }

◆ GetTargetClusterSize()

int Garfield::TrackSrim::GetTargetClusterSize ( ) const

inline

Definition at line 52 of file TrackSrim.hh.

52{ return m_nsize; }

◆ GetWorkFunction()

double Garfield::TrackSrim::GetWorkFunction ( ) const

inline

Definition at line 21 of file TrackSrim.hh.

21{ return m_work; }

◆ NewTrack()

bool Garfield::TrackSrim::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 698 of file TrackSrim.cc.

                                           {
 
  // Generates electrons for a SRIM track
  // SRMGEN
  const std::string hdr = m_className + "::NewTrack: ";
 
  // Verify that a sensor has been set.
  if (!m_sensor) {
    std::cerr << hdr << "\n    Sensor is not defined.\n";
    return false;
  }
 
  // Get the bounding box.
  double xmin = 0., ymin = 0., zmin = 0.;
  double xmax = 0., ymax = 0., zmax = 0.;
  if (!m_sensor->GetArea(xmin, ymin, zmin, xmax, ymax, zmax)) {
    std::cerr << hdr << "\n    Drift area is not set.\n";
    return false;
  } else if (x0 < xmin || x0 > xmax || 
             y0 < ymin || y0 > ymax || z0 < zmin || z0 > zmax) {
    std::cerr << hdr << "\n    Initial position outside bounding box.\n";
    return false;
  }
 
  // Make sure the initial position is inside an ionisable medium.
  Medium* medium = NULL;
  if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
    std::cerr << hdr << "\n    No medium at initial position.\n";
    return false;
  } else if (!medium->IsIonisable()) {
    std::cerr << hdr << "\n    Medium at initial position is not ionisable.\n";
    return false;
  }
 
  // Normalise and store the direction.
  const double normdir = sqrt(dx0 * dx0 + dy0 * dy0 + dz0 * dz0);
  double xdir = dx0;
  double ydir = dy0;
  double zdir = dz0;
  if (normdir < Small) {
    if (m_debug) {
      std::cout << hdr << "\n    Direction vector has zero norm.\n"
                << "    Initial direction is randomized.\n";
    }
    // Null vector. Sample the direction isotropically.
    RndmDirection(xdir, ydir, zdir);
  } else {
    // Normalise the direction vector.
    xdir /= normdir;
    ydir /= normdir;
    zdir /= normdir;
  }
 
  // Make sure all necessary parameters have been set.
  if (m_mass < Small) {
    std::cerr << hdr << "\n    Particle mass not set.\n";
    return false;
  } else if (!m_chargeset) {
    std::cerr << hdr << "\n    Particle charge not set.\n";
    return false;
  } else if (m_energy < Small) {
    std::cerr << hdr << "\n    Initial particle energy not set.\n";
    return false;
  } else if (m_work < Small) {
    std::cerr << hdr << "\n    Work function not set.\n";
    return false;
  } else if (m_a < Small || m_z < Small) {
    std::cerr << hdr << "\n    A and/or Z not set.\n";
    return false;
  }
  // Check the initial energy (in MeV).
  const double ekin0 = 1.e-6 * GetKineticEnergy();
  if (ekin0 < 1.e-14 * m_mass || ekin0 < 1.e-3 * m_work) {
    if (m_debug) {
      std::cout << hdr << "Initial kinetic energy E = " << ekin0
                << " MeV such that beta2 = 0 or E << W; particle stopped.\n";
    }
    return true;
  }
 
  // Get an upper limit for the track length.
  const double tracklength = 10 * Interpolate(ekin0, m_ekin, m_range);
 
  // Header of debugging output.
  if (m_debug) {
    std::cout << hdr << "Track generation with the following parameters:\n";
    const unsigned int nTable = m_ekin.size();
    printf("      Table size           %u\n", nTable);
    printf("      Particle kin. energy %g MeV\n", ekin0);
    printf("      Particle mass        %g MeV\n", 1.e-6 * m_mass);
    printf("      Particle charge      %g\n", m_q);
    printf("      Work function        %g eV\n", m_work);
    if (m_fano > 0.) {
      printf("      Fano factor          %g\n", m_fano);
    } else {
      std::cout << "      Fano factor          Not set\n";
    } 
    printf("      Long. straggling:    %d\n", m_useLongStraggle);
    printf("      Trans. straggling:   %d\n", m_useTransStraggle);
    // printf("      Vavilov generator:   %d\n", m_precisevavilov);
    printf("      Cluster size         %d\n", m_nsize);
  }
 
  // Reset the cluster count
  m_currcluster = 0;
  m_clusters.clear();
 
 
  // Initial situation: starting position
  double x = x0;
  double y = y0;
  double z = z0;
 
  // Store the energy [MeV].
  double e = ekin0;
  // Total distance covered
  double dsum = 0.0;
  // Pool of unused energy
  double epool = 0.0;
 
  // Loop generating clusters
  int iter = 0;
  while (iter < m_maxclusters || m_maxclusters < 0) {
    // Work out what the energy loss per cm, straggling and projected range are
    // at the start of the step.
    const double dedxem = DedxEM(e) * m_density;
    const double dedxhd = DedxHD(e) * m_density;
    const double prange = Interpolate(e, m_ekin, m_range);
    double strlon = Interpolate(e, m_ekin, m_longstraggle);
    double strlat = Interpolate(e, m_ekin, m_transstraggle);
 
    if (!m_useLongStraggle) strlon = 0;
    if (!m_useTransStraggle) strlat = 0;
 
    if (m_debug) {
      std::cout << hdr << "\n    Energy = " << e << " MeV,\n    dEdx em, hd = "
                << dedxem << ", " << dedxhd << " MeV/cm,\n    e-/cm = " 
                << 1.e6 * dedxem / m_work << ".\n    Straggling long/lat: " 
                << strlon << ", " << strlat << " cm\n"; 
    }
    // Find the step size for which we get approximately the target # clusters.
    double step;
    if (m_nsize > 0) {
      step = m_nsize * 1.e-6 * m_work / dedxem;
    } else {
      const double ncls = m_maxclusters > 0 ? 0.5 * m_maxclusters : 100;
      step = ekin0 / (ncls * (dedxem + dedxhd));
    }
    // Truncate if this step exceeds the length.
    bool finish = false;
    // Make an accurate integration of the energy loss over the step.
    double deem = 0., dehd = 0.;
    PreciseLoss(step, e, deem, dehd);
    // If the energy loss exceeds the particle energy, truncate step.
    double stpmax;
    if (deem + dehd > e) {
      EstimateRange(e, step, stpmax);
      step = stpmax;
      PreciseLoss(step, e, deem, dehd);
      deem = e * deem / (dehd + deem);
      dehd = e - deem;
      finish = true;
      if (m_debug) std::cout << hdr << "Finish raised. Track length reached.\n";
    } else {
      stpmax = tracklength - dsum;
    }
 
    // Ensure that this is larger than the minimum modelable step size.
    double stpmin;
    if (!SmallestStep(e, deem, step, stpmin)) {
      std::cerr << hdr << "\n    Failure computing the minimum step size."
                << "\n    Clustering abandoned.\n";
      return false;
    }
 
    double eloss;
    if (stpmin > stpmax) {
      // No way to find a suitable step size: use fixed energy loss
      if (m_debug) std::cout << hdr << "stpmin > stpmax. Deposit all energy.\n";
      eloss = deem;
      if (e - eloss - dehd < 0) eloss = e - dehd;
      finish = true;
      if (m_debug) std::cout << hdr << "Finish raised. Single deposit.\n";
    } else if (step < stpmin) {
      // If needed enlarge the step size
      if (m_debug) std::cout << hdr << "Enlarging step size.\n";
      step = stpmin;
      PreciseLoss(step, e, deem, dehd);
      if (deem + dehd > e) {
        if (m_debug) std::cout << hdr << "Excess loss. Recomputing stpmax.\n";
        EstimateRange(e, step, stpmax);
        step = stpmax;
        PreciseLoss(step, e, deem, dehd);
        deem = e * deem / (dehd + deem);
        dehd = e - deem;
        eloss = deem;
      } else {
        eloss = RndmEnergyLoss(e, deem, step);
      }
    } else {
      // Draw an actual energy loss for such a step.
      if (m_debug) std::cout << hdr << "Using existing step size.\n";
      eloss = RndmEnergyLoss(e, deem, step);
    }
    // Ensure we are neither below 0 nor above the total energy.
    if (eloss < 0) {
      if (m_debug) std::cout << hdr << "Truncating negative energy loss.\n";
      eloss = 0;
    } else if (eloss > e - dehd) {
      if (m_debug) std::cout << hdr << "Excess energy loss, using mean.\n";
      eloss = deem;
      if (e - eloss - dehd < 0) {
        eloss = e - dehd;
        finish = true;
        if (m_debug) std::cout << hdr << "Finish raised. Using mean energy.\n";
      }
    }
    if (m_debug) {
      std::cout << hdr << "\n    Step length = " << step << " cm.\n    "
                << "Mean loss =   " << deem << " MeV.\n    "
                << "Actual loss = " << eloss << " MeV.\n";
    }
 
    // Check that the cluster is in an ionisable medium and within bounding box
    if (!m_sensor->GetMedium(x, y, z, medium)) {
      if (m_debug) {
        std::cout << hdr << "No medium at position (" 
                  << x << "," << y << "," << z << ").\n";
      }
      break;
    } else if (!medium->IsIonisable()) {
      if (m_debug) {
        std::cout << hdr << "Medium at (" 
                  << x << "," << y << "," << z << ") is not ionisable.\n";
      }
      break;
    } else if (!m_sensor->IsInArea(x, y, z)) {
      if (m_debug) {
        std::cout << hdr << "Cluster at (" 
                  << x << "," << y << "," << z << ") outside bounding box.\n";
      }
      break;
    }
    // Add a cluster.
    cluster newcluster;
    newcluster.x = x;
    newcluster.y = y;
    newcluster.z = z;
    newcluster.t = t0;
    if (m_fano < Small) {
      // No fluctuations.
      newcluster.electrons = int((eloss + epool) / (1.e-6 * m_work));
      newcluster.ec = m_work * newcluster.electrons;
    } else {
      double ecl = 1.e6 * (eloss + epool);
      newcluster.electrons = 0.0;
      newcluster.ec = 0.0;
      while (true) {
        // if (newcluster.ec < 100) printf("ec = %g\n", newcluster.ec);
        const double ernd1 = RndmHeedWF(m_work, m_fano);
        if (ernd1 > ecl) break;
        newcluster.electrons++;
        newcluster.ec += ernd1;
        ecl -= ernd1;
      }
      // printf("ec = %g DONE\n", newcluster.ec);
      if (m_debug)
        std::cout << hdr << "EM + pool: " << 1.e6 * (eloss + epool) 
                  << " eV, W: " << m_work << " eV, E/w: " 
                  << (eloss + epool) / (1.0e-6 * m_work) << ", n: " 
                  << newcluster.electrons << ".\n"; 
    }
    newcluster.kinetic = e;
    epool += eloss - 1.e-6 * newcluster.ec;
    if (m_debug) {
      std::cout << hdr << "Cluster " << m_clusters.size() << "\n    at (" 
                << newcluster.x << ", " << newcluster.y << ", " << newcluster.z
                << "),\n    e = " << newcluster.ec << ",\n    n = " 
                << newcluster.electrons << ",\n    pool = " 
                << epool << " MeV.\n";
    }
    m_clusters.push_back(newcluster);
 
    // Keep track of the length and energy
    dsum += step;
    e -= eloss + dehd;
    if (finish) {
      // Stop if the flag is raised
      if (m_debug) std::cout << hdr << "Finishing flag raised.\n";
      break;
    } else if (e < ekin0 * 1.e-9) {
      // No energy left
      if (m_debug) std::cout << hdr << "Energy exhausted.\n";
      break;
    }
    // Draw scattering distances
    const double scale = sqrt(step / prange);
    const double sigt1 = RndmGaussian(0., scale * strlat);
    const double sigt2 = RndmGaussian(0., scale * strlat);
    const double sigl = RndmGaussian(0., scale * strlon);
    if (m_debug) std::cout << hdr << "sigma l, t1, t2: " 
                         << sigl << ", " << sigt1 << ", " << sigt2 << "\n";
    // Rotation angles to bring z-axis in line
    double theta, phi;
    if (xdir * xdir + zdir * zdir <= 0) {
      if (ydir < 0) {
        theta = -HalfPi;
      } else if (ydir > 0) {
        theta = +HalfPi;
      } else {
        std::cerr << hdr << "\n    Zero step length; clustering abandoned.\n";
        return false;
      }
      phi = 0;
    } else {
      phi = atan2(xdir, zdir);
      theta = atan2(ydir, sqrt(xdir * xdir + zdir * zdir));
    }
 
    // Update position
    const double cp = cos(phi);
    const double ct = cos(theta);
    const double sp = sin(phi);
    const double st = sin(theta);
    x += step * xdir + cp * sigt1 - sp * st * sigt2 + sp * ct * sigl;
    y += step * ydir + ct * sigt2 + st * sigl;
    z += step * zdir - sp * sigt1 - cp * st * sigt2 + cp * ct * sigl;
 
    // (Do not) update direction
    if (false) {
      xdir = step * xdir + cp * sigt1 - sp * st * sigt2 + sp * ct * sigl;
      ydir = step * ydir + ct * sigt2 + st * sigl;
      zdir = step * zdir - sp * sigt1 - cp * st * sigt2 + cp * ct * sigl;
      double dnorm = sqrt(xdir * xdir + ydir * ydir + zdir * zdir);
      if (dnorm <= 0) {
        std::cerr << hdr << "\n    Zero step length; clustering abandoned.\n";
        return false;
      }
      xdir = xdir / dnorm;
      ydir = ydir / dnorm;
      zdir = zdir / dnorm;
    }
    // Next cluster
    iter++;
  }
  if (iter == m_maxclusters) {
    std::cerr << hdr << "Exceeded maximum number of clusters.\n";
  }
  return true;
  // finished generating
}

◆ PlotEnergyLoss()

void Garfield::TrackSrim::PlotEnergyLoss ( )

Definition at line 367 of file TrackSrim.cc.

                               {
 
  // Make a graph for the 3 curves to plot
  double xmin = 0., xmax = 0.;
  double ymin = 0., ymax = 0.;
  const unsigned int nPoints = m_ekin.size();
  TGraph* grem = new TGraph(nPoints);
  TGraph* grhd = new TGraph(nPoints);
  TGraph* grtot = new TGraph(nPoints);
  for (unsigned int i = 0; i < nPoints; ++i) {
    const double eplot = m_ekin[i];
    if (eplot < xmin) xmin = eplot;
    if (eplot > xmax) xmax = eplot;
    const double emplot = m_emloss[i] * m_density;
    const double hdplot = m_hdloss[i] * m_density;
    const double totplot = emplot + hdplot;
    if (totplot < ymin) ymin = totplot;
    if (totplot > ymax) ymax = totplot;
    grem->SetPoint(i, eplot, emplot);
    grhd->SetPoint(i, eplot, hdplot);
    grtot->SetPoint(i, eplot, totplot);
  }
 
  // Prepare a plot frame
  TCanvas* celoss = new TCanvas();
  celoss->SetLogx();
  frame(celoss, xmin, 0.0, xmax, ymax * 1.05, "Ion energy [MeV]",
        "Energy loss [MeV/cm]");
  TLegend* legend = new TLegend(0.65, 0.75, 0.94, 0.95);
  legend->SetShadowColor(0);
  legend->SetTextFont(22);
  legend->SetTextSize(0.044);
  legend->SetFillColor(kWhite);
  legend->SetBorderSize(1);
 
  grem->SetLineColor(kBlue);
  grem->SetLineStyle(kSolid);
  grem->SetLineWidth(2.0);
  grem->SetMarkerStyle(21);
  grem->SetMarkerColor(kBlack);
  grem->Draw("SAME");
  legend->AddEntry(grem, "EM energy loss", "l");
 
  grhd->SetLineColor(kGreen + 2);
  grhd->SetLineStyle(kSolid);
  grhd->SetLineWidth(2.0);
  grhd->SetMarkerStyle(21);
  grhd->SetMarkerColor(kBlack);
  grhd->Draw("SAME");
  legend->AddEntry(grhd, "HD energy loss", "l");
 
  grtot->SetLineColor(kOrange);
  grtot->SetLineStyle(kSolid);
  grtot->SetLineWidth(2.0);
  grtot->SetMarkerStyle(21);
  grtot->SetMarkerColor(kBlack);
  grtot->Draw("SAME");
  legend->AddEntry(grtot, "Total energy loss", "l");
 
  legend->Draw();
}

◆ PlotRange()

void Garfield::TrackSrim::PlotRange ( )

Definition at line 429 of file TrackSrim.cc.

                          {
 
  // Make a graph
  double xmin = 0., xmax = 0.;
  double ymin = 0., ymax = 0.;
  const unsigned int nPoints = m_ekin.size();
  TGraph* grrange = new TGraph(nPoints);
  for (unsigned int i = 0; i < nPoints; ++i) {
    const double eplot = m_ekin[i];
    if (eplot < xmin) xmin = eplot;
    if (eplot > xmax) xmax = eplot;
    const double rangeplot = m_range[i];
    if (rangeplot < ymin) ymin = rangeplot;
    if (rangeplot > ymax) ymax = rangeplot;
    grrange->SetPoint(i, eplot, rangeplot);
  }
 
  // Prepare a plot frame
  TCanvas* crange = new TCanvas();
  crange->SetLogx();
  frame(crange, xmin, 0.0, xmax, ymax * 1.05, "Ion energy [MeV]",
        "Projected range [cm]");
 
  grrange->SetLineColor(kOrange);
  grrange->SetLineStyle(kSolid);
  grrange->SetLineWidth(2.0);
  grrange->SetMarkerStyle(21);
  grrange->SetMarkerColor(kBlack);
  grrange->Draw("SAME");
}

◆ PlotStraggling()

void Garfield::TrackSrim::PlotStraggling ( )

Definition at line 460 of file TrackSrim.cc.

                               {
 
  // Make a graph for the 2 curves to plot
  double xmin = 0., xmax = 0.;
  double ymin = 0., ymax = 0.;
  const unsigned int nPoints = m_ekin.size();
  TGraph* grlong = new TGraph(nPoints);
  TGraph* grtrans = new TGraph(nPoints);
  for (unsigned int i = 0; i < nPoints; ++i) {
    const double eplot = m_ekin[i];
    if (eplot < xmin) xmin = eplot;
    if (eplot > xmax) xmax = eplot;
    const double longplot = m_longstraggle[i];
    const double transplot = m_transstraggle[i];
    if (longplot < ymin) ymin = longplot;
    if (longplot > ymax) ymax = longplot;
    if (transplot < ymin) ymin = transplot;
    if (transplot > ymax) ymax = transplot;
    grlong->SetPoint(i, eplot, longplot);
    grtrans->SetPoint(i, eplot, transplot);
  }
 
  // Prepare a plot frame
  TCanvas* cstraggle = new TCanvas();
  cstraggle->SetLogx();
  frame(cstraggle, xmin, 0.0, xmax, ymax * 1.05, "Ion energy [MeV]",
        "Straggling [cm]");
  TLegend* legend = new TLegend(0.15, 0.75, 0.44, 0.95);
  legend->SetShadowColor(0);
  legend->SetTextFont(22);
  legend->SetTextSize(0.044);
  legend->SetFillColor(kWhite);
  legend->SetBorderSize(1);
 
  grlong->SetLineColor(kOrange);
  grlong->SetLineStyle(kSolid);
  grlong->SetLineWidth(2.0);
  grlong->SetMarkerStyle(21);
  grlong->SetMarkerColor(kBlack);
  grlong->Draw("SAME");
  legend->AddEntry(grlong, "Longitudinal", "l");
 
  grtrans->SetLineColor(kGreen + 2);
  grtrans->SetLineStyle(kSolid);
  grtrans->SetLineWidth(2.0);
  grtrans->SetMarkerStyle(21);
  grtrans->SetMarkerColor(kBlack);
  grtrans->Draw("SAME");
  legend->AddEntry(grtrans, "Transversal", "l");
 
  legend->Draw();
}

◆ PreciseLoss()

bool Garfield::TrackSrim::PreciseLoss	(	const double	step,
		const double	estart,
		double &	deem,
		double &	dehd
	)		const

protected

Definition at line 521 of file TrackSrim.cc.

                                                              {
  // SRMRKS
 
  const std::string hdr = m_className + "::PreciseLoss: ";
  // Debugging
  if (m_debug) {
    std::cout << hdr << "\n"
              << "    Initial energy: " << estart << " MeV\n"
              << "    Step: " << step << " cm\n";
  }
  // Precision aimed for.
  const double eps = 1.0e-2;
  // Number of intervals.
  unsigned int ndiv = 1;
  // Loop until precision achieved
  const unsigned int nMaxIter = 10;
  bool converged = false;
  for (unsigned int iter = 0; iter < nMaxIter; ++iter) {
    double e4 = estart;
    double e2 = estart;
    deem = 0.;
    dehd = 0.;
    // Compute rk2 and rk4 over the number of sub-divisions
    const double s = m_density * step / ndiv;
    for (unsigned int i = 0; i < ndiv; i++) {
      // rk2: initial point
      const double em21 = s * DedxEM(e2);
      const double hd21 = s * DedxHD(e2);
      const double de21 = em21 + hd21;
      // Mid-way point
      const double em22 = s * DedxEM(e2 - 0.5 * de21);
      const double hd22 = s * DedxHD(e2 - 0.5 * de21);
      // Trace the rk2 energy
      e2 -= em22 + hd22;
      // rk4: initial point
      const double em41 = s * DedxEM(e4);
      const double hd41 = s * DedxHD(e4);
      const double de41 = em41 + hd41;
      // Mid-way point
      const double em42 = s * DedxEM(e4 - 0.5 * de41);
      const double hd42 = s * DedxHD(e4 - 0.5 * de41);
      const double de42 = em42 + hd42;
      // Second mid-point estimate
      const double em43 = s * DedxEM(e4 - 0.5 * de42);
      const double hd43 = s * DedxHD(e4 - 0.5 * de42);
      const double de43 = em43 + hd43;
      // End point estimate
      const double em44 = s * DedxEM(e4 - de43);
      const double hd44 = s * DedxHD(e4 - de43);
      const double de44 = em44 + hd44;
      // Store the energy loss terms (according to rk4)
      deem += (em41 + em44) / 6. + (em42 + em43) / 3.;
      dehd += (hd41 + hd44) / 6. + (hd42 + hd43) / 3.;
      // Store the new energy computed with rk4
      e4 -= (de41 + de44) / 6. + (de42 + de43) / 3.;
    }
    if (m_debug) {
      std::cout << hdr << "\n    Iteration " << iter << " has " << ndiv 
                << " division(s). Losses:\n";
      printf("\tde4 = %12g, de2 = %12g MeV\n", estart - e2, estart - e4);
      printf("\tem4 = %12g, hd4 = %12g MeV\n", deem, dehd);
    }
    // Compare the two estimates
    if (fabs(e2 - e4) > eps * (fabs(e2) + fabs(e4) + fabs(estart))) {
      // Repeat with twice the number of steps.
      ndiv *= 2;
    } else {
      converged = true;
      break;
    }
  }
 
  if (!converged) {
    std::cerr << hdr << "No convergence achieved integrating energy loss.\n";
  } else if (m_debug) {
    std::cout << hdr << "Convergence at eps = " << eps << "\n";
  }
  return converged;
}

Referenced by EstimateRange(), and NewTrack().

◆ Print()

void Garfield::TrackSrim::Print ( )

Definition at line 345 of file TrackSrim.cc.

                      {
 
  std::cout << "TrackSrim::Print:\n    SRIM energy loss table\n\n"
            << "    Energy     EM Loss     HD loss       Range  "
            << "l straggle  t straggle\n"
            << "     [MeV]    [MeV/cm]    [MeV/cm]        [cm] "
            << "      [cm]        [cm]\n\n";
  const unsigned int nPoints = m_emloss.size();
  for (unsigned int i = 0; i < nPoints; ++i) {
    printf("%10g  %10g  %10g  %10g  %10g  %10g\n", m_ekin[i],
           m_emloss[i] * m_density, m_hdloss[i] * m_density, m_range[i],
           m_longstraggle[i], m_transstraggle[i]);
  }
  std::cout << "\n";
  printf("    Work function:  %g eV\n", m_work);
  printf("    Fano factor:    %g\n", m_fano);
  printf("    Ion charge:     %g\n", m_q);
  printf("    Mass:           %g MeV\n", 1.e-6 * m_mass);
  printf("    Density:        %g g/cm3\n", m_density);
  printf("    A, Z:           %g, %g\n", m_a, m_z);
}

◆ ReadFile()

bool Garfield::TrackSrim::ReadFile ( const std::string & file )

Definition at line 140 of file TrackSrim.cc.

                                              {
  // SRMREA
 
  const std::string hdr = m_className + "::ReadFile:\n    ";
  // Open the material list.
  std::ifstream fsrim;
  fsrim.open(file.c_str(), std::ios::in);
  if (fsrim.fail()) {
    std::cerr << hdr << "Could not open SRIM  file " << file
              << " for reading.\n    The file perhaps does not exist.\n";
    return false;
  }
  unsigned int nread = 0;
 
  // Read the header
  if (m_debug) {
    std::cout << hdr << "SRIM header records from file " << file << "\n";
  }
  const int size = 100;
  char line[size];
  while (fsrim.getline(line, 100, '\n')) {
    nread++;
    if (strstr(line, "SRIM version") != NULL) {
      if (m_debug) std::cout << "\t" << line << "\n";
    } else if (strstr(line, "Calc. date") != NULL) {
      if (m_debug) std::cout << "\t" << line << "\n";
    } else if (strstr(line, "Ion =") != NULL) {
      break;
    }
  }
 
  // Identify the ion
  char* token = NULL;
  token = strtok(line, " []=");
  token = strtok(NULL, " []=");
  token = strtok(NULL, " []=");
  // Set the ion charge.
  m_q = std::atof(token);
  m_chargeset = true;
  token = strtok(NULL, " []=");
  token = strtok(NULL, " []=");
  token = strtok(NULL, " []=");
  // Set the ion mass (convert amu to eV).
  m_mass = std::atof(token) * AtomicMassUnitElectronVolt;
 
  // Find the target density
  if (!fsrim.getline(line, 100, '\n')) {
    std::cerr << hdr << "Premature EOF looking for target density (line "
              << nread << ").\n";
    return false;
  }
  nread++;
  if (!fsrim.getline(line, 100, '\n')) {
    std::cerr << hdr << "Premature EOF looking for target density (line "
              << nread << ").\n";
    return false;
  }
  nread++;
  token = strtok(line, " ");
  token = strtok(NULL, " ");
  token = strtok(NULL, " ");
  token = strtok(NULL, " ");
  SetDensity(std::atof(token));
 
  // Check the stopping units
  while (fsrim.getline(line, 100, '\n')) {
    nread++;
    if (strstr(line, "Stopping Units") == NULL) continue;
    if (strstr(line, "Stopping Units =  MeV / (mg/cm2)") != NULL) {
      if (m_debug) {
        std::cout << hdr << "Stopping units: MeV / (mg/cm2) as expected.\n";
      }
      break;
    }
    std::cerr << hdr << "Unknown stopping units. Aborting (line " << nread
              << ").\n";
    return false;
  }
 
  // Skip to the table
  while (fsrim.getline(line, 100, '\n')) {
    nread++;
    if (strstr(line, "-----------") != NULL) break;
  }
 
  // Read the table line by line
  m_ekin.clear();
  m_emloss.clear();
  m_hdloss.clear();
  m_range.clear();
  m_transstraggle.clear();
  m_longstraggle.clear();
  unsigned int ntable = 0;
  while (fsrim.getline(line, 100, '\n')) {
    nread++;
    if (strstr(line, "-----------") != NULL) break;
    // Energy
    token = strtok(line, " ");
    m_ekin.push_back(atof(token));
    token = strtok(NULL, " ");
    if (strcmp(token, "eV") == 0) {
      m_ekin[ntable] *= 1.0e-6;
    } else if (strcmp(token, "keV") == 0) {
      m_ekin[ntable] *= 1.0e-3;
    } else if (strcmp(token, "GeV") == 0) {
      m_ekin[ntable] *= 1.0e3;
    } else if (strcmp(token, "MeV") != 0) {
      std::cerr << hdr << "Unknown energy unit " << token << "; aborting\n";
      return false;
    }
    // EM loss
    token = strtok(NULL, " ");
    m_emloss.push_back(atof(token));
    // HD loss
    token = strtok(NULL, " ");
    m_hdloss.push_back(atof(token));
    // Projected range
    token = strtok(NULL, " ");
    m_range.push_back(atof(token));
    token = strtok(NULL, " ");
    if (strcmp(token, "A") == 0) {
      m_range[ntable] *= 1.0e-8;
    } else if (strcmp(token, "um") == 0) {
      m_range[ntable] *= 1.0e-4;
    } else if (strcmp(token, "mm") == 0) {
      m_range[ntable] *= 1.0e-1;
    } else if (strcmp(token, "m") == 0) {
      m_range[ntable] *= 1.0e2;
    } else if (strcmp(token, "km") == 0) {
      m_range[ntable] *= 1.0e5;
    } else if (strcmp(token, "cm") != 0) {
      std::cerr << hdr << "Unknown distance unit " << token << "; aborting\n";
      return false;
    }
    // Longitudinal straggling
    token = strtok(NULL, " ");
    m_longstraggle.push_back(atof(token));
    token = strtok(NULL, " ");
    if (strcmp(token, "A") == 0) {
      m_longstraggle[ntable] *= 1.0e-8;
    } else if (strcmp(token, "um") == 0) {
      m_longstraggle[ntable] *= 1.0e-4;
    } else if (strcmp(token, "mm") == 0) {
      m_longstraggle[ntable] *= 1.0e-1;
    } else if (strcmp(token, "m") == 0) {
      m_longstraggle[ntable] *= 1.0e2;
    } else if (strcmp(token, "km") == 0) {
      m_longstraggle[ntable] *= 1.0e5;
    } else if (strcmp(token, "cm") != 0) {
      std::cerr << hdr << "Unknown distance unit " << token << "; aborting\n";
      return false;
    }
    // Transverse straggling
    token = strtok(NULL, " ");
    m_transstraggle.push_back(atof(token));
    token = strtok(NULL, " ");
    if (strcmp(token, "A") == 0) {
      m_transstraggle[ntable] *= 1.0e-8;
    } else if (strcmp(token, "um") == 0) {
      m_transstraggle[ntable] *= 1.0e-4;
    } else if (strcmp(token, "mm") == 0) {
      m_transstraggle[ntable] *= 1.0e-1;
    } else if (strcmp(token, "m") == 0) {
      m_transstraggle[ntable] *= 1.0e2;
    } else if (strcmp(token, "km") == 0) {
      m_transstraggle[ntable] *= 1.0e5;
    } else if (strcmp(token, "cm") != 0) {
      std::cerr << hdr << "Unknown distance unit " << token << "; aborting\n";
      return false;
    }
 
    // Increment table line counter
    ++ntable;
  }
 
  // Find the scaling factor and convert to MeV/cm
  double scale = -1.;
  while (fsrim.getline(line, 100, '\n')) {
    nread++;
    if (strstr(line, "=============") != NULL) {
      break;
    } else if (strstr(line, "MeV / (mg/cm2)") != NULL ||
               strstr(line, "MeV/(mg/cm2)") != NULL) {
      token = strtok(line, " ");
      scale = std::atof(token);
    }
  }
  if (scale < 0) {
    std::cerr << hdr << "Did not find stopping unit scaling; aborting.\n";
    return false;
  }
  scale *= 1.e3;
  for (unsigned int i = 0; i < ntable; ++i) {
    m_emloss[i] *= scale;
    m_hdloss[i] *= scale;
  }
 
  // Seems to have worked
  if (m_debug) {
    std::cout << hdr << "Successfully read " << file << "(" << nread
              << " lines).\n";
  }
  return true;
}

◆ RndmEnergyLoss()

double Garfield::TrackSrim::RndmEnergyLoss	(	const double	ekin,
		const double	de,
		const double	step
	)		const

protected

Definition at line 1226 of file TrackSrim.cc.

                                                          {
  //   RNDDE  - Generates a random energy loss.
  //   VARIABLES : EKIN       : Kinetic energy [MeV]
  //            DE         : Mean energy loss over the step [MeV]
  //            STEP       : Step length [cm]
  //            BETA2      : Velocity-squared
  //            GAMMA      : Projectile gamma
  //            EMAX       : Maximum energy transfer per collision [MeV]
  //            XI         : Rutherford term [MeV]
  //            FCONST     : Proportionality constant
  //            EMASS      : Electron mass [MeV]
  // (Last changed on 26/10/07.)
 
  const std::string hdr = "TrackSrim::RndmEnergyLoss: ";
  // Check correctness.
  if (ekin <= 0 || de <= 0 || step <= 0) {
    std::cerr << hdr << "\n    Input parameters not valid.\n    Ekin = " 
              << ekin << " MeV, dE = " << de << " MeV, step length = " 
              << step << " cm.\n";
    return 0.;
  } else if (m_mass <= 0 || fabs(m_q) <= 0) {
    std::cerr << hdr << "\n    Track parameters not valid.\n    Mass = "
              << m_mass << " MeV, charge = " << m_q << ".\n";
    return 0.;
  } else if (m_a <= 0 || m_z <= 0 || m_density <= 0) {
    std::cerr << hdr << "\n    Material parameters not valid.\n A = " << m_a 
              << ", Z = " << m_z << ", density = " << m_density << " g/cm3.\n";
    return 0.;
  }
  // Basic kinematic parameters
  const double rkin = 1.e6 * ekin / m_mass;
  const double gamma = 1. + rkin;
  const double beta2 = rkin > 1.e-5 ? 1. - 1. / (gamma * gamma) : 2. * rkin;
 
  // Compute maximum energy transfer
  const double rm = ElectronMass / m_mass;
  const double emax = 2 * ElectronMass * 1.e-6 * beta2 * gamma * gamma /
                      (1. + 2 * gamma * rm + rm * rm);
  // Compute the Rutherford term
  const double fconst = 0.1534;
  const double xi = fconst * m_q * m_q * m_z * m_density * step / (m_a * beta2);
  // Compute the scaling parameter
  const double rkappa = xi / emax;
  // Debugging output.
  if (m_debug) {
    PrintSettings(hdr, de, step, ekin, beta2, gamma, m_a, m_z, m_density, m_q,
                  m_mass, emax, xi, rkappa);
  }
  double rndde = de;
  if (m_model <= 0 || m_model > 4) {
    // No fluctuations.
    if (m_debug) std::cout << hdr << "Fixed energy loss.\n";
  } else if (m_model == 1) {
    // Landau distribution
    if (m_debug) std::cout << hdr << "Landau imposed.\n";
    const double xlmean = -(log(rkappa) + beta2 + 1. - Gamma);
    rndde += xi * (RndmLandau() - xlmean);
  } else if (m_model == 2) {
    // Vavilov distribution, ensure we are in range.
    if (m_debug) std::cout << hdr << "Vavilov imposed.\n";
    if (rkappa > 0.01 && rkappa < 12) {
      const double xvav = RndmVavilov(rkappa, beta2);
      rndde += xi * (xvav + log(rkappa) + beta2 + (1 - Gamma));
    }
  } else if (m_model == 3) {
    // Gaussian model
    if (m_debug) std::cout << hdr << "Gaussian imposed.\n";
    rndde += RndmGaussian(0., sqrt(xi * emax * (1 - 0.5 * beta2)));
  } else if (rkappa < 0.05) {
    // Combined model: for low kappa, use the landau distribution.
    if (m_debug) std::cout << hdr << "Landau automatic.\n";
    const double xlmean = -(log(rkappa) + beta2 + (1 - Gamma));
    const double par[] = {0.50884,    1.26116, 0.0346688,  1.46314,
                          0.15088e-2, 1.00324, -0.13049e-3};
    const double xlmax = par[0] + par[1] * xlmean + par[2] * xlmean * xlmean +
                         par[6] * xlmean * xlmean * xlmean +
                         (par[3] + xlmean * par[4]) * exp(par[5] * xlmean);
    double xlan = RndmLandau();
    for (unsigned int iter = 0; iter < 100; ++iter) {
      if (xlan < xlmax) break;
      xlan = RndmLandau();
    }
    rndde += xi * (xlan - xlmean);
  } else if (rkappa < 5) {
    // For medium kappa, use the Vavilov distribution, precise
    // } else if (m_precisevavilov && rkappa < 5) {
    //   printf("Vavilov slow automatic\n");
    //   rndde = de+xi*(rndvvl(rkappa,beta2) + log(xi/emax)+beta2+(1-Gamma));
    //   // ... or fast.
    if (m_debug) std::cout << hdr << "Vavilov fast automatic.\n";
    const double xvav = RndmVavilov(rkappa, beta2);
    rndde += xi * (xvav + log(rkappa) + beta2 + (1 - Gamma));
  } else {
    // And for large kappa, use the Gaussian values.
    if (m_debug) std::cout << hdr << "Gaussian automatic.\n";
    rndde = RndmGaussian(de, sqrt(xi * emax * (1 - 0.5 * beta2)));
  }
  // Debugging output
  if (m_debug)
    std::cout << hdr << "Energy loss generated = " << rndde << " MeV.\n";
  return rndde;
}

Referenced by NewTrack().

◆ SetAtomicMassNumbers()

void Garfield::TrackSrim::SetAtomicMassNumbers	(	const double	a,
		const double	z
	)

inline

Set/get A and Z of the target medium.

Definition at line 29 of file TrackSrim.hh.

◆ SetClustersMaximum() [1/2]

int Garfield::TrackSrim::SetClustersMaximum ( ) const

inline

Definition at line 55 of file TrackSrim.hh.

55{ return m_maxclusters; }

◆ SetClustersMaximum() [2/2]

void Garfield::TrackSrim::SetClustersMaximum ( const int n )

inline

Definition at line 54 of file TrackSrim.hh.

54{ m_maxclusters = n; }

◆ SetDensity()

void Garfield::TrackSrim::SetDensity ( const double density )

inline

Set/get the density [g/cm3] of the target medium.

Definition at line 26 of file TrackSrim.hh.

26{ m_density = density; }

Referenced by ReadFile().

◆ SetFanoFactor()

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

inline

Set/get the Fano factor.

Definition at line 23 of file TrackSrim.hh.

23{ m_fano = f; }

◆ SetModel()

void Garfield::TrackSrim::SetModel ( const int m )

inline

Set/get the fluctuation model (0 = none, 1 = Landau, 2 = Vavilov, 3 = Gaussian, 4 = Combined)

Definition at line 40 of file TrackSrim.hh.

40{ m_model = m; }

◆ SetTargetClusterSize()

void Garfield::TrackSrim::SetTargetClusterSize ( const int n )

inline

Definition at line 51 of file TrackSrim.hh.

51{ m_nsize = n; }

◆ SetWorkFunction()

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

inline

Set/get the W value [eV].

Definition at line 20 of file TrackSrim.hh.

20{ m_work = w; }

◆ SmallestStep()

bool Garfield::TrackSrim::SmallestStep	(	double	ekin,
		double	de,
		double	step,
		double &	stpmin
	)

protected

Definition at line 1052 of file TrackSrim.cc.

                                             {
  // Determines the smallest step size for which there is little
  // or no risk of finding negative energy fluctuations.
  // SRMMST
 
  const std::string hdr = m_className + "::SmallestStep: ";
  const double expmax = 30;
 
  // By default, assume the step is right.
  stpmin = step;
  // Check correctness.
  if (ekin <= 0 || de <= 0 || step <= 0) {
    std::cerr << hdr << "\n    Input parameters not valid.\n    Ekin = "
              << ekin << " MeV, dE = " << de << " MeV, step length = " 
              << step << " cm.\n";
    return false;
  } else if (m_mass <= 0 || fabs(m_q) <= 0) {
    std::cerr << hdr << "\n    Track parameters not valid.\n    Mass = "
              << m_mass << " eV, charge = " << m_q << ".\n";
    return false;
  } else if (m_a <= 0 || m_z <= 0 || m_density <= 0) {
    std::cerr << hdr << "\n    Gas parameters not valid.\n    A = " << m_a
              << ", Z = " << m_z  << " density = " << m_density << " g/cm3.\n";
    return false;
  }
 
  // Basic kinematic parameters
  const double rkin = 1.e6 * ekin / m_mass;
  const double gamma = 1. + rkin;
  const double beta2 = rkin > 1.e-5 ? 1. - 1. / (gamma * gamma) : 2. * rkin;
 
  // Compute maximum energy transfer [MeV]
  const double rm = ElectronMass / m_mass;
  const double emax = 2 * ElectronMass * 1.e-6 * beta2 * gamma * gamma /
                      (1. + 2 * gamma * rm + rm * rm);
  // Compute the Rutherford term
  const double fconst = 0.1534;
  double xi = fconst * m_q * m_q * m_z * m_density * step / (m_a * beta2);
  // Compute the scaling parameter
  double rkappa = xi / emax;
  // Step size and energy loss
  double denow = de;
  double stpnow = step;
  const unsigned int nMaxIter = 10;
  for (unsigned int iter = 0; iter < nMaxIter; ++iter) {
    bool retry = false;
    // Debugging output.
    if (m_debug) {
      PrintSettings(hdr, denow, stpnow, ekin, beta2, gamma, m_a, m_z, m_density,
                    m_q, m_mass, emax, xi, rkappa);
    }
    double xinew = xi;
    double rknew = rkappa;
    if (m_model <= 0 || m_model > 4) {
      // No fluctuations: any step is permitted
      stpmin = stpnow;
    } else if (m_model == 1) {
      // Landau distribution
      const double xlmin = -3.;
      const double exponent = -xlmin - 1. + Gamma - beta2 - denow / xi;
      const double rklim = exponent < -expmax ? 0. : exp(exponent);
      stpmin = stpnow * (rklim / rkappa);
      if (m_debug) {
        std::cout << hdr << "Landau distribution is imposed.\n    kappa_min = "
                  << rklim << ", d_min = " << stpmin << " cm.\n";
      }
    } else if (m_model == 2) {
      // Vavilov distribution, ensure we're in range.
      const double xlmin = StepVavilov(rkappa);
      const double exponent = -xlmin - 1. + Gamma - beta2 - denow / xi;
      const double rklim = exponent < -expmax ? 0. : exp(exponent);
      stpmin = stpnow * (rklim / rkappa);
      xinew = fconst * m_q * m_q * m_z * m_density * stpmin / (m_a * beta2);
      rknew = xinew / emax;
      if (m_debug) {
        std::cout << hdr << "Vavilov distribution is imposed.\n    kappa_min = "
                  << rklim << ", d_min = " << stpmin << " cm\n    kappa_new = "
                  << rknew << ", xi_new = " << xinew << " MeV.\n";
      }
      if (stpmin > stpnow * 1.1) {
        if (m_debug) std::cout << hdr << "Step size increase. New pass.\n";
        retry = true;
      }
    } else if (m_model == 3) {
      // Gaussian model
      stpmin = stpnow * 16 * xi * emax * (1 - 0.5 * beta2) / (denow * denow);
      if (m_debug) {
        std::cout << hdr << "Gaussian distribution is imposed.\n";
        printf("\td_min = %g cm.\n\tsigma/mu_old = %g, sigma/mu_min = %g\n",
               stpmin, sqrt(xi * emax * (1 - 0.5 * beta2)) / de,
               sqrt((fconst * m_q * m_q * m_z * m_density * stpmin /
                     (m_a * beta2)) *
                    emax * (1 - 0.5 * beta2)) /
                   (stpmin * denow / stpnow));
      }
    } else if (rkappa < 0.05) {
      // Combined model: for low kappa, use the Landau distribution.
      const double xlmin = -3.;
      const double exponent = -xlmin - 1. + Gamma - beta2 - denow / xi;
      const double rklim = exponent < -expmax ? 0. : exp(exponent);
      stpmin = stpnow * (rklim / rkappa);
      xinew = fconst * m_q * m_q * m_z * m_density * stpmin / (m_a * beta2);
      rknew = xinew / emax;
      if (m_debug) {
        std::cout << hdr << "Landau distribution automatic.\n    kappa_min = "
                  << rklim << ", d_min = " << stpmin << " cm.\n";
      }
      if (rknew > 0.05 || stpmin > stpnow * 1.1) {
        retry = true;
        if (m_debug) {
          std::cout << hdr << "Model change or step increase. New pass.\n";
        }
      }
    } else if (rkappa < 5) {
      // For medium kappa, use the Vavilov distribution
      const double xlmin = StepVavilov(rkappa);
      const double exponent = -xlmin - 1. + Gamma - beta2 - denow / xi;
      const double rklim = exponent < -expmax ? 0. : exp(exponent);
      stpmin = stpnow * (rklim / rkappa);
      xinew = fconst * m_q * m_q * m_z * m_density * stpmin / (m_a * beta2);
      rknew = xinew / emax;
      if (m_debug) {
        std::cout << hdr << "Vavilov distribution automatic.\n    kappa_min = "
                  << rklim << ", d_min = " << stpmin << " cm\n    kappa_new = "
                  << rknew << ", xi_new = " << xinew << " MeV.\n";
      }
      if (rknew > 5 || stpmin > stpnow * 1.1) {
        retry = true;
        if (m_debug) {
          std::cout << hdr << "Model change or step increase. New pass.\n";
        }
      }
    } else {
      // And for large kappa, use the Gaussian values.
      stpmin = stpnow * 16 * xi * emax * (1 - 0.5 * beta2) / (denow * denow);
      if (m_debug) {
        std::cout << hdr << "Gaussian distribution automatic.\n";
        printf("\td_min = %g cm.\n\tsigma/mu_old = %g, sigma/mu_min = %g\n",
               stpmin, sqrt(xi * emax * (1 - 0.5 * beta2)) / de,
               sqrt((fconst * m_q * m_q * m_z * m_density * stpmin /
                     (m_a * beta2)) *
                    emax * (1 - 0.5 * beta2)) /
                   (stpmin * denow / stpnow));
      }
    }
    // See whether we should do another pass.
    if (stpnow > stpmin) {
      if (m_debug) {
        std::cout << hdr << "Step size ok, minimum: " << stpmin << " cm\n";
      }
      break;
    }
    if (!retry) {
      if (m_debug) {
        std::cerr << hdr << "\nStep size must be increased to " 
                  << stpmin << "cm.\n";
      }
      break;
    }
    // New iteration
    rkappa = rknew;
    xi = xinew;
    denow *= stpmin / stpnow;
    stpnow = stpmin;
    if (m_debug) std::cout << hdr << "Iteration " << iter << "\n";
    if (iter == nMaxIter - 1) {
      // Need interation, but ran out of tries
      std::cerr << hdr << "\n    No convergence reached on step size.\n";
    }
  }
  return true;
}

Referenced by NewTrack().

Member Data Documentation

◆ m_a

double Garfield::TrackSrim::m_a

protected

A and Z of the gas.

Definition at line 90 of file TrackSrim.hh.

Referenced by GetAtomicMassMumbers(), NewTrack(), Print(), RndmEnergyLoss(), SetAtomicMassNumbers(), and SmallestStep().

◆ m_chargeset

bool Garfield::TrackSrim::m_chargeset

protected

Charge gas been defined.

Definition at line 78 of file TrackSrim.hh.

Referenced by NewTrack(), and ReadFile().

◆ m_clusters

std::vector<cluster> Garfield::TrackSrim::m_clusters

protected

Definition at line 121 of file TrackSrim.hh.

Referenced by GetCluster(), and NewTrack().

◆ m_currcluster

unsigned int Garfield::TrackSrim::m_currcluster

protected

Index of the next cluster to be returned.

Definition at line 109 of file TrackSrim.hh.

Referenced by GetCluster(), and NewTrack().

◆ m_density

double Garfield::TrackSrim::m_density

protected

Density [g/cm3].

Definition at line 80 of file TrackSrim.hh.

Referenced by GetDensity(), NewTrack(), PlotEnergyLoss(), PreciseLoss(), Print(), RndmEnergyLoss(), SetDensity(), and SmallestStep().

◆ m_ekin

std::vector<double> Garfield::TrackSrim::m_ekin

protected

Energy in energy loss table [MeV].

Definition at line 96 of file TrackSrim.hh.

Referenced by DedxEM(), DedxHD(), NewTrack(), PlotEnergyLoss(), PlotRange(), PlotStraggling(), Print(), and ReadFile().

◆ m_emloss

std::vector<double> Garfield::TrackSrim::m_emloss

protected

EM energy loss [MeV cm2/g].

Definition at line 98 of file TrackSrim.hh.

Referenced by DedxEM(), PlotEnergyLoss(), Print(), and ReadFile().

◆ m_fano

double Garfield::TrackSrim::m_fano

protected

Fano factor [-].

Definition at line 84 of file TrackSrim.hh.

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

◆ m_hdloss

std::vector<double> Garfield::TrackSrim::m_hdloss

protected

Hadronic energy loss [MeV cm2/g].

Definition at line 100 of file TrackSrim.hh.

Referenced by DedxHD(), PlotEnergyLoss(), Print(), and ReadFile().

◆ m_longstraggle

std::vector<double> Garfield::TrackSrim::m_longstraggle

protected

Longitudinal straggling [cm].

Definition at line 106 of file TrackSrim.hh.

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

◆ m_mass

double Garfield::TrackSrim::m_mass

protected

Mass of ion [MeV].

Definition at line 88 of file TrackSrim.hh.

Referenced by NewTrack(), Print(), ReadFile(), RndmEnergyLoss(), SmallestStep(), and TrackSrim().

◆ m_maxclusters

int Garfield::TrackSrim::m_maxclusters

protected

Maximum number of clusters allowed (infinite if 0)

Definition at line 94 of file TrackSrim.hh.

Referenced by NewTrack(), and SetClustersMaximum().

◆ m_model

unsigned int Garfield::TrackSrim::m_model

protected

Fluctuation model (0 = none, 1 = Landau, 2 = Vavilov, 3 = Gaussian, 4 = Combined)

Definition at line 112 of file TrackSrim.hh.

Referenced by GetModel(), RndmEnergyLoss(), SetModel(), and SmallestStep().

◆ m_nsize

int Garfield::TrackSrim::m_nsize

protected

Targeted cluster size.

Definition at line 114 of file TrackSrim.hh.

Referenced by GetTargetClusterSize(), NewTrack(), and SetTargetClusterSize().

◆ m_precisevavilov

bool Garfield::TrackSrim::m_precisevavilov

protected

Use precise Vavilov generator.

Definition at line 71 of file TrackSrim.hh.

Referenced by DisablePreciseVavilov(), and EnablePreciseVavilov().

◆ m_q

double Garfield::TrackSrim::m_q

protected

Charge of ion.

Definition at line 86 of file TrackSrim.hh.

Referenced by NewTrack(), Print(), ReadFile(), RndmEnergyLoss(), and SmallestStep().

◆ m_range

std::vector<double> Garfield::TrackSrim::m_range

protected

Projected range [cm].

Definition at line 102 of file TrackSrim.hh.

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

◆ m_transstraggle

std::vector<double> Garfield::TrackSrim::m_transstraggle

protected

Transverse straggling [cm].

Definition at line 104 of file TrackSrim.hh.

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

◆ m_useLongStraggle

bool Garfield::TrackSrim::m_useLongStraggle

protected

Include longitudinal straggling.

Definition at line 75 of file TrackSrim.hh.

Referenced by DisableLongitudinalStraggling(), EnableLongitudinalStraggling(), and NewTrack().

◆ m_useTransStraggle

bool Garfield::TrackSrim::m_useTransStraggle

protected

Include transverse straggling.

Definition at line 73 of file TrackSrim.hh.

Referenced by DisableTransverseStraggling(), EnableTransverseStraggling(), and NewTrack().

◆ m_work

double Garfield::TrackSrim::m_work

protected

Work function [eV].

Definition at line 82 of file TrackSrim.hh.

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

◆ m_z

double Garfield::TrackSrim::m_z

protected

Definition at line 91 of file TrackSrim.hh.

Referenced by GetAtomicMassMumbers(), NewTrack(), Print(), RndmEnergyLoss(), SetAtomicMassNumbers(), and SmallestStep().

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

Classes

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ TrackSrim()

◆ ~TrackSrim()

Member Function Documentation

◆ DedxEM()

◆ DedxHD()

◆ DisableLongitudinalStraggling()

◆ DisablePreciseVavilov()

◆ DisableTransverseStraggling()

◆ EnableLongitudinalStraggling()

◆ EnablePreciseVavilov()

◆ EnableTransverseStraggling()

◆ EstimateRange()

◆ GetAtomicMassMumbers()

◆ GetCluster()

◆ GetDensity()

◆ GetFanoFactor()

◆ GetModel()

◆ GetTargetClusterSize()

◆ GetWorkFunction()

◆ NewTrack()

◆ PlotEnergyLoss()

◆ PlotRange()

◆ PlotStraggling()

◆ PreciseLoss()

◆ Print()

◆ ReadFile()

◆ RndmEnergyLoss()

◆ SetAtomicMassNumbers()

◆ SetClustersMaximum() [1/2]

◆ SetClustersMaximum() [2/2]

◆ SetDensity()

◆ SetFanoFactor()

◆ SetModel()

◆ SetTargetClusterSize()

◆ SetWorkFunction()

◆ SmallestStep()

Member Data Documentation

◆ m_a

◆ m_chargeset

◆ m_clusters

◆ m_currcluster

◆ m_density

◆ m_ekin

◆ m_emloss

◆ m_fano

◆ m_hdloss

◆ m_longstraggle

◆ m_mass

◆ m_maxclusters

◆ m_model

◆ m_nsize

◆ m_precisevavilov

◆ m_q

◆ m_range

◆ m_transstraggle

◆ m_useLongStraggle

◆ m_useTransStraggle

◆ m_work

◆ m_z