#include <HeedDeltaElectron.h>

Inheritance diagram for Heed::HeedDeltaElectron:

Public Member Functions
	HeedDeltaElectron ()=default
	Default constructor.

	HeedDeltaElectron (manip_absvol primvol, const point &pt, const vec &vel, vfloat ftime, long fparent_particle_number, HeedFieldMap fieldmap, bool fs_print_listing=false)
	Constructor.

virtual	~HeedDeltaElectron ()
	Destructor.

HeedDeltaElectron *	copy () const override
	Clone the particle.

void	print (std::ostream &file, int l) const override
	Print-out.

Public Member Functions inherited from Heed::eparticle
	eparticle ()=default
	Default constructor.

	eparticle (manip_absvol primvol, const point &pt, const vec &vel, vfloat time, particle_def fpardef, HeedFieldMap *fieldmap)
	Constructor using velocity vector.

virtual	~eparticle ()
	Destructor.

eparticle *	copy () const override
	Clone the particle.

void	print (std::ostream &file, int l) const override
	Print-out.

Public Member Functions inherited from Heed::mparticle
	mparticle ()=default
	Default constructor.

	mparticle (manip_absvol *primvol, const point &pt, const vec &vel, vfloat ftime, double fmass)
	Constructor, $\gamma - 1$ calculated from the from velocity vector.

virtual	~mparticle ()
	Destructor.

double	kinetic_energy () const
	Get the current kinetic energy.

void	print (std::ostream &file, int l) const override
	Print-out.

mparticle *	copy () const override
	Clone the particle.

Public Member Functions inherited from Heed::gparticle
	gparticle ()=default
	Default constructor.

	gparticle (manip_absvol *primvol, const point &pt, const vec &vel, vfloat time)
	Constructor.

virtual	~gparticle ()
	Destructor.

virtual void	fly (std::vector< gparticle * > &secondaries)
	Transport the particle.

virtual void	fly (std::vector< gparticle * > &secondaries, const bool one_step)

void	set_step_limits (const vfloat fmax_range, const vfloat frad_for_straight, const vfloat fmax_straight_arange, const vfloat fmax_circ_arange)
	Set limits/parameters for trajectory steps.

const vec &	position () const
	Get the current position of the particle.

vfloat	time () const
	Get the current time of the particle.

const vec &	direction () const
	Get the current direction of the particle.

bool	alive () const

virtual void	print (std::ostream &file, int l) const
	Print-out.

virtual gparticle *	copy () const
	Clone the particle.

Public Member Functions inherited from Heed::particle_type
	particle_type ()=default

	particle_type (particle_def *f)

	particle_type (const char *name, int s=0)

int	operator== (const particle_type &f)

int	operator!= (const particle_type &f)

void	print_notation (std::ostream &file) const

Public Attributes
std::vector< HeedCondElectron >	conduction_electrons

std::vector< HeedCondElectron >	conduction_ions

long	parent_particle_number

Public Attributes inherited from Heed::particle_type
particle_def *	pardef = nullptr

Protected Member Functions
void	physics_mrange (double &fmrange) override

void	physics_after_new_speed (std::vector< gparticle * > &secondaries) override
	Apply any other processes (turn the trajectory, kill the particle, ...).

Protected Member Functions inherited from Heed::eparticle
int	force (const point &pt, vec &f, vec &f_perp, vfloat &mrange) override
	Calculate force components.

Protected Member Functions inherited from Heed::mparticle
void	step (std::vector< gparticle * > &secondaries) override

void	curvature (bool &curved, vec &frelcen, vfloat &fmrange, vfloat prec) override

virtual int	force (const point &pt, vec &f, vec &f_perp, vfloat &mrange)

Protected Member Functions inherited from Heed::gparticle
virtual void	step (std::vector< gparticle * > &secondaries)

virtual void	change_vol ()
	Move from one volume to another.

virtual void	curvature (bool &curved, vec &frelcen, vfloat &fmrange, vfloat prec)

virtual void	physics_after_new_speed (std::vector< gparticle * > &)
	Apply any other processes (turn the trajectory, kill the particle, ...).

virtual void	physics (std::vector< gparticle * > &)
	Apply any other processes (turn the trajectory, kill the particle, ...).

virtual void	physics_mrange (double &fmrange)

virtual stvpoint	calc_step_to_bord ()
	Determine next position.

stvpoint	switch_new_vol ()
	Generate next position in new volume.

Additional Inherited Members
Protected Attributes inherited from Heed::eparticle
HeedFieldMap *	m_fieldMap = nullptr
	Pointer to field map.

Protected Attributes inherited from Heed::mparticle
double	m_mass = 0.
	Mass (not mass * speed_of_light^2)

double	m_curr_ekin = 0.
	Current kinetic energy.

double	m_orig_ekin = 0.
	Original kinetic energy.

double	m_prev_ekin = 0.
	Previous kinetic energy.

double	m_curr_gamma_1 = 0.
	Current $\gamma - 1$ .

double	m_orig_gamma_1 = 0.
	Original $\gamma - 1$ .

double	m_prev_gamma_1 = 0.
	Previous $\gamma - 1$ .

Protected Attributes inherited from Heed::gparticle
bool	m_alive = false
	Status flag whether the particle is active.

long	m_nstep = 0
	Step number.

long	m_nzero_step = 0
	Number of previous steps with zero range (including this step).

stvpoint	m_origin
	Original point.

double	m_total_range_from_origin = 0.
	Range from origin to current position.

stvpoint	m_prevpos
	Previous point.

stvpoint	m_currpos
	Current point.

stvpoint	m_nextpos
	Next point.

Static Protected Attributes inherited from Heed::gparticle
static constexpr long	m_max_qzero_step = 100
	Max. number of zero-steps allowed.

Detailed Description

Definition of delta-electron which can be traced through the geometry. 2003, I. Smirnov

Definition at line 15 of file HeedDeltaElectron.h.

Constructor & Destructor Documentation

◆ HeedDeltaElectron() [1/2]

Heed::HeedDeltaElectron::HeedDeltaElectron ( )

default

Default constructor.

Referenced by copy().

◆ HeedDeltaElectron() [2/2]

Heed::HeedDeltaElectron::HeedDeltaElectron	(	manip_absvol *	primvol,
		const point &	pt,
		const vec &	vel,
		vfloat	ftime,
		long	fparent_particle_number,
		HeedFieldMap *	fieldmap,
		bool	fs_print_listing = `false`
	)

Constructor.

Definition at line 48 of file HeedDeltaElectron.cpp.

    : eparticle(primvol, pt, vel, ftime, &electron_def, fieldmap),
      parent_particle_number(fparent_particle_number),
      m_particle_number(last_particle_number++),
      m_print_listing(fprint_listing) {
  mfunname("HeedDeltaElectron::HeedDeltaElectron(...)");
}

◆ ~HeedDeltaElectron()

virtual Heed::HeedDeltaElectron::~HeedDeltaElectron ( )

inlinevirtual

Destructor.

Definition at line 24 of file HeedDeltaElectron.h.

24{}

Member Function Documentation

◆ copy()

HeedDeltaElectron * Heed::HeedDeltaElectron::copy ( ) const

inlineoverridevirtual

Clone the particle.

Reimplemented from Heed::gparticle.

Definition at line 26 of file HeedDeltaElectron.h.

                                           {
    return new HeedDeltaElectron(*this);
  }

◆ physics_after_new_speed()

void Heed::HeedDeltaElectron::physics_after_new_speed ( std::vector< gparticle * > & )

overrideprotectedvirtual

Apply any other processes (turn the trajectory, kill the particle, ...).

Reimplemented from Heed::gparticle.

Definition at line 138 of file HeedDeltaElectron.cpp.

                              {
  mfunname("void HeedDeltaElectron::physics_after_new_speed()");
  if (m_print_listing) {
    mcout << "HeedDeltaElectron::physics_after_new_speed\n";
    Iprint2n(mcout, m_currpos.prange, m_curr_ekin);
  }
  check_econd11(vecerror, != 0, mcerr);
  if (m_currpos.prange <= 0.0) {
    if (m_curr_ekin <= 0.0) {
      // Get local volume.
      absvol* av = m_currpos.volume();
      if (av && av->s_sensitive && m_fieldMap->inside(m_currpos.ptloc)) {
        if (m_print_listing) mcout << "Convert to conduction electron.\n";
        conduction_electrons.emplace_back(
            HeedCondElectron(m_currpos.ptloc, m_currpos.time));
      }
      m_alive = false;
    }
    if (m_print_listing) mcout << "exit due to currpos.prange <= 0.0\n";
    return;
  }
  // Get local volume and convert it to a cross-section object.
  const absvol* av = m_currpos.volume();
  auto hdecs = dynamic_cast<const HeedDeltaElectronCS*>(av);
  if (!hdecs) return;
  double ek = m_curr_ekin / MeV;
  if (m_print_listing) {
    Iprintnf(mcout, ek);
    Iprint3n(mcout, m_stop_eloss, m_phys_mrange, m_currpos.prange);
  }
  // Calculate dE/dx and energy loss. Update the kinetic energy.
  double dedx;
  double Eloss = 0.;
  if (m_stop_eloss && m_phys_mrange == m_currpos.prange) {
    Eloss = m_curr_ekin;
    m_curr_ekin = 0.0;
    dedx = Eloss / m_currpos.prange / (MeV / cm);
  } else {
    EnergyMesh* emesh = hdecs->hmd->energy_mesh;
    dedx = interpolate(emesh, ek, hdecs->eLoss);
    Eloss = std::min(m_currpos.prange * dedx * MeV / cm, m_curr_ekin);
    m_total_eloss += Eloss;
    m_curr_ekin -= Eloss;
  }
  if (m_print_listing)
    Iprint3nf(mcout, m_prev_ekin / eV, m_curr_ekin / eV, Eloss / eV);
  if (m_curr_ekin <= 0.0) {
    if (m_print_listing) mcout << "m_curr_ekin <= 0.0\n";
    m_curr_ekin = 0.0;
    m_curr_gamma_1 = 0.0;
    m_currpos.speed = 0.0;
    m_alive = false;
  } else {
    const double resten = m_mass * c_squared;
    m_curr_gamma_1 = m_curr_ekin / resten;
    m_currpos.speed = c_light * lorbeta(m_curr_gamma_1);
  }
  absvol* vav = m_currpos.volume();
  if (vav && vav->s_sensitive) {
    if (m_print_listing) {
      mcout << "volume is sensitive\n";
      Iprint2nf(mcout, Eloss / eV, m_necessary_energy / eV);
    }
    if (Eloss > 0.0) ionisation(Eloss, dedx, hdecs->pairprod);
  }
  if (m_print_listing) {
    mcout << '\n';
    Iprintn(mcout, m_alive);
  }
  if (!m_alive) {
    // Done tracing the delta electron. Create the last conduction electron.
    vav = m_currpos.volume();
    if (vav && vav->s_sensitive && m_fieldMap->inside(m_currpos.ptloc)) {
      if (m_print_listing) mcout << "Last conduction electron\n";
      conduction_electrons.emplace_back(
          HeedCondElectron(m_currpos.ptloc, m_currpos.time));
    }
    return;
  }
 
  if (m_print_listing) mcout << "\nstart to rotate by low angle\n";
  double ek_restr = std::max(ek, 0.0005);
  if (m_print_listing) Iprint2nf(mcout, m_currpos.prange, m_phys_mrange);
  if (m_currpos.prange < m_phys_mrange) {
    // recalculate scatterings
    m_path_length = false;
    if (s_low_mult_scattering) {
      EnergyMesh* emesh = hdecs->hmd->energy_mesh;
      const double low_path_length = interpolate(emesh, ek_restr, hdecs->low_lambda) * cm;
      if (m_print_listing) Iprintnf(mcout, low_path_length / cm);
      m_mult_low_path_length = false;
      m_q_low_path_length = m_currpos.prange / low_path_length;
      if (m_print_listing) Iprintnf(mcout, m_q_low_path_length);
    }
  }
  if (m_print_listing) Iprintnf(mcout, m_q_low_path_length);
#ifdef RANDOM_POIS
  if (m_q_low_path_length > 0.0) {
    long random_q_low_path_length = pois(m_q_low_path_length);
    m_q_low_path_length = long(random_q_low_path_length);
    if (m_print_listing) {
      mcout << "After pois:\n";
      Iprintnf(mcout, m_q_low_path_length);
    }
  }
#endif
  if (m_q_low_path_length > 0) {
#ifdef DIRECT_LOW_IF_LITTLE
    const long max_q_low_path_length_for_direct = 5;
    if (m_q_low_path_length < max_q_low_path_length_for_direct) {
      // direct modeling
      if (m_print_listing) {
        mcout << "direct modeling of low scatterings\n";
        Iprint(mcout, m_currpos.dir);
      }
      EnergyMesh* emesh = hdecs->hmd->energy_mesh;
      const long n1r = findInterval(emesh, ek_restr);
      for (long nscat = 0; nscat < m_q_low_path_length; ++nscat) {
        if (m_print_listing) Iprintn(mcout, nscat);
        double theta_rot =
            hdecs->low_angular_points_ran[n1r].ran(SRANLUX()) * degree;
        if (m_print_listing) Iprintnf(mcout, theta_rot);
        vec dir = m_currpos.dir;
        basis temp(dir, "temp");
        vec vturn;
        vturn.random_round_vec();
        vturn = vturn * sin(theta_rot);
        vec new_dir(vturn.x, vturn.y, cos(theta_rot));
        new_dir.down(&temp);
        m_currpos.dir = new_dir;
        if (m_print_listing) Iprint(mcout, new_dir);
      }
      m_currpos.dirloc = m_currpos.dir;
      m_currpos.tid.up_absref(&m_currpos.dirloc);
    } else {
#endif
      double sigma_ctheta = hdecs->get_sigma(ek_restr, m_q_low_path_length);
      // actually it is mean(1-cos(theta)) or
      // sqrt( mean( square(1-cos(theta) ) ) ) depending on USE_MEAN_COEF
 
      if (m_print_listing) Iprintnf(mcout, sigma_ctheta);
      // Gauss (but actually exponential distribution fits better).
      // double ctheta = 1.0 - fabs(rnorm_improved() * sigma_ctheta);
      // Exponential:
      double ctheta = 0.0;
      {
#ifdef USE_MEAN_COEF
#else
      double sq2 = sqrt(2.0);
#endif
        do {
          double y = 0.0;
          do {  // in order to avoid SRANLUX() = 1
            y = SRANLUX();
            if (m_print_listing) Iprintnf(mcout, y);
          } while (y == 1.0);
#ifdef USE_MEAN_COEF
          double x = sigma_ctheta * (-log(1.0 - y));
#else
        double x = sigma_ctheta * 1.0 / sq2 * (-log(1.0 - y));
#endif
          ctheta = 1.0 - x;
          if (m_print_listing) Iprint2nf(mcout, x, ctheta);
        } while (ctheta <= -1.0);  // avoid absurd cos(theta)
        check_econd21(ctheta, < -1.0 ||, > 1.0, mcerr);
      }
      if (m_print_listing) Iprintnf(mcout, ctheta);
      double theta_rot = acos(ctheta);
      if (m_print_listing) Iprint2nf(mcout, theta_rot, theta_rot / degree);
      vec dir = m_currpos.dir;
      basis temp(dir, "temp");
      vec vturn;
      vturn.random_round_vec();
      vturn = vturn * sin(theta_rot);
      vec new_dir(vturn.x, vturn.y, cos(theta_rot));
      new_dir.down(&temp);
      m_currpos.dir = new_dir;
      m_currpos.dirloc = m_currpos.dir;
      m_currpos.tid.up_absref(&m_currpos.dirloc);
    }
#ifdef DIRECT_LOW_IF_LITTLE
  }
#endif
  if (m_path_length) {
    if (m_print_listing) {
      mcout << "\nstarting to rotate by large angle" << std::endl;
      Iprintnf(mcout, m_path_length);
    }
    EnergyMesh* emesh = hdecs->hmd->energy_mesh;
    const long n1r = findInterval(emesh, ek_restr);
    double theta_rot = hdecs->angular_points_ran[n1r].ran(SRANLUX()) * degree;
    if (m_print_listing) Iprintnf(mcout, theta_rot);
    vec dir = m_currpos.dir;
    basis temp(dir, "temp");
    vec vturn;
    vturn.random_round_vec();
    vturn = vturn * sin(theta_rot);
    vec new_dir(vturn.x, vturn.y, cos(theta_rot));
    new_dir.down(&temp);
    m_currpos.dir = new_dir;
    m_currpos.dirloc = m_currpos.dir;
    m_currpos.tid.up_absref(&m_currpos.dirloc);
  }
  if (m_print_listing) Iprint2nf(mcout, m_currpos.dir, m_currpos.dirloc);
}

◆ physics_mrange()

void Heed::HeedDeltaElectron::physics_mrange ( double & fmrange )

overrideprotectedvirtual

Reduce the maximal possible range due to continuous processes. Called from calc_step_to_bord after the call of curvature. but before considering the crossing with volumes. Therefore mrange may be reduced after this.

Reimplemented from Heed::gparticle.

Definition at line 60 of file HeedDeltaElectron.cpp.

                                                      {
  mfunname("void HeedDeltaElectron::physics_mrange(double& fmrange)");
  if (m_print_listing) mcout << "HeedDeltaElectron::physics_mrange\n";
 
  m_mult_low_path_length = false;
  m_q_low_path_length = 0.0;
  m_path_length = false;
  if (fmrange <= 0.0) return;
  if (m_curr_ekin <= 0.0) {
    fmrange = 0.0;
    return;
  }
  // Get local volume and convert it to a cross-section object.
  const absvol* av = m_currpos.volume();
  auto hdecs = dynamic_cast<const HeedDeltaElectronCS*>(av);
  if (!hdecs) return;
  if (m_print_listing) Iprintnf(mcout, fmrange);
  const double ek = m_curr_ekin / MeV;
  // Get the dE/dx at this kinetic energy.
  EnergyMesh* emesh = hdecs->hmd->energy_mesh;
  const double dedx = interpolate(emesh, ek, hdecs->eLoss);
  // Min. loss 50 eV.
  double eloss = std::max(0.1 * ek, 0.00005);
  if (eloss > ek) {
    eloss = ek;
    m_stop_eloss = true;
  } else {
    m_stop_eloss = false;
  }
  fmrange = std::min(fmrange, (eloss / dedx) * cm);
  if (m_print_listing) Iprint2nf(mcout, fmrange, ek);
  const double ek_restr = std::max(ek, 0.0005);
  if (m_print_listing) Iprintnf(mcout, ek_restr / keV);
 
  double low_path_length = 0.;  // in internal units
  if (s_low_mult_scattering) {
    low_path_length = interpolate(emesh, ek_restr, hdecs->low_lambda) * cm;
    if (m_print_listing) Iprintnf(mcout, low_path_length / cm);
    long qscat = hdecs->eesls->get_qscat();
    const double sigma_ctheta = hdecs->get_sigma(ek_restr, qscat);
    // Reduce the number of scatterings, if the angle is too large.
    if (sigma_ctheta > 0.3) qscat = long(qscat * 0.3 / sigma_ctheta);
    const double mult_low_path_length = qscat * low_path_length;
    if (m_print_listing) Iprintnf(mcout, mult_low_path_length);
    if (fmrange > mult_low_path_length) {
      fmrange = mult_low_path_length;
      m_mult_low_path_length = true;
      m_q_low_path_length = hdecs->eesls->get_qscat();
      m_stop_eloss = false;
    } else {
      m_mult_low_path_length = false;
      m_q_low_path_length = fmrange / low_path_length;
    }
    if (m_print_listing) Iprint2nf(mcout, fmrange, m_q_low_path_length);
  }
 
  if (s_high_mult_scattering) {
    const double mean_path = interpolate(emesh, ek_restr, hdecs->lambda);
    if (m_print_listing) Iprintnf(mcout, mean_path);
    const double path_length = -mean_path * cm * log(1.0 - SRANLUX());
    if (m_print_listing) Iprintnf(mcout, path_length);
    if (fmrange > path_length) {
      fmrange = path_length;
      m_path_length = true;
      m_mult_low_path_length = true;
      if (s_low_mult_scattering) {
        m_q_low_path_length = fmrange / low_path_length;
        if (m_print_listing) Iprintnf(mcout, m_q_low_path_length);
      }
      m_stop_eloss = false;
    } else {
      m_path_length = false;
    }
    if (m_print_listing) Iprintnf(mcout, fmrange);
  }
  m_phys_mrange = fmrange;
}

◆ print()

void Heed::HeedDeltaElectron::print	(	std::ostream &	file,
		int	l
	)		const

overridevirtual

Print-out.

Reimplemented from Heed::gparticle.

Definition at line 397 of file HeedDeltaElectron.cpp.

                                                           {
  if (l < 0) return;
  Ifile << "HeedDeltaElectron (l=" << l
        << "): particle_number=" << m_particle_number << "\n";
  if (l == 1) return;
  indn.n += 2;
  Ifile << "s_low_mult_scattering=" << s_low_mult_scattering
        << " s_high_mult_scattering=" << s_high_mult_scattering << '\n';
  Ifile << "phys_mrange=" << m_phys_mrange << " stop_eloss=" << m_stop_eloss
        << " mult_low_path_length=" << m_mult_low_path_length << '\n';
  Ifile << "q_low_path_length=" << m_q_low_path_length
        << " path_length=" << m_path_length
        << " necessary_energy/eV=" << m_necessary_energy / eV << '\n';
  Ifile << " parent_particle_number=" << parent_particle_number << '\n';
 
  mparticle::print(file, l - 1);
  indn.n -= 2;
}

Member Data Documentation

◆ conduction_electrons

std::vector<HeedCondElectron> Heed::HeedDeltaElectron::conduction_electrons

Definition at line 31 of file HeedDeltaElectron.h.

Referenced by physics_after_new_speed(), and Garfield::TrackHeed::TransportDeltaElectron().

◆ conduction_ions

std::vector<HeedCondElectron> Heed::HeedDeltaElectron::conduction_ions

Definition at line 32 of file HeedDeltaElectron.h.

Referenced by Garfield::TrackHeed::TransportDeltaElectron().

◆ parent_particle_number

long Heed::HeedDeltaElectron::parent_particle_number

Definition at line 34 of file HeedDeltaElectron.h.

Referenced by print().

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

Public Member Functions

Public Attributes

Protected Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ HeedDeltaElectron() [1/2]

◆ HeedDeltaElectron() [2/2]

◆ ~HeedDeltaElectron()

Member Function Documentation

◆ copy()

◆ physics_after_new_speed()

◆ physics_mrange()

◆ print()

Member Data Documentation

◆ conduction_electrons

◆ conduction_ions

◆ parent_particle_number