#include <HeedDeltaElectron.h>

Inheritance diagram for Heed::HeedDeltaElectron:

Public Member Functions
	HeedDeltaElectron ()
	Default constructor.

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

virtual	~HeedDeltaElectron ()
	Destructor.

virtual void	physics_mrange (double &fmrange)

virtual void	physics_after_new_speed (std::vector< gparticle * > &secondaries)

virtual HeedDeltaElectron *	copy () const

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

Public Member Functions inherited from Heed::eparticle
	eparticle ()
	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.

virtual eparticle *	copy () const

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

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

Public Member Functions inherited from Heed::mparticle
void	check_consistency () const
	Check consistency of kin_energy, gamma_1, speed, speed_of_light and mass.

virtual void	step (std::vector< gparticle * > &secondaries)

virtual void	curvature (int &fs_cf, vec &frelcen, vfloat &fmrange, vfloat prec)

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

void	new_speed ()
	Set new speed, direction and time for currpos.

	mparticle ()
	Default constructor.

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

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

virtual mparticle *	copy () const

virtual	~mparticle ()
	Destructor.

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

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

virtual	~gparticle ()
	Destructor.

virtual void	step (std::vector< gparticle * > &secondaries)

virtual void	change_vol (void)

virtual void	curvature (int &fs_cf, vec &frelcen, vfloat &fmrange, vfloat prec)

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

virtual void	physics (std::vector< gparticle * > &)

virtual void	physics_mrange (double &fmrange)

virtual stvpoint	calc_step_to_bord ()
	Produces nextpos.

stvpoint	switch_new_vol ()

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

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

virtual gparticle *	copy () const

Public Member Functions inherited from Heed::RegPassivePtr
	RegPassivePtr (void)

	RegPassivePtr (char fs_ban_del, char fs_ban_sub, char fs_ban_cop=0)

	RegPassivePtr (const RegPassivePtr &f)

RegPassivePtr &	operator= (const RegPassivePtr &f)

CountPP_ns::CountPassivePtr *	book (void) const

void	clear_pointers (void) const

virtual RegPassivePtr *	copy () const

virtual	~RegPassivePtr ()

virtual void	print (std::ostream &file, int l=1) const

void	set_s_ban_del (char fs_ban_del)

char	get_s_ban_del (void) const

void	set_s_ban_sub (char fs_ban_sub)

char	get_s_ban_sub (void) const

void	set_s_ban_cop (char fs_ban_cop)

char	get_s_ban_cop (void) const

void	set_s_allow_del_at_zero_count (char fs_allow_del_at_zero_count)

char	get_s_allow_del_at_zero_count (void) const

long	get_total_number_of_references (void) const

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

	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::mparticle
double	mass
	Mass (not mass * speed_of_light^2)

double	orig_kin_energy

double	orig_gamma_1

double	prev_kin_energy

double	prev_gamma_1

double	curr_kin_energy

double	curr_gamma_1

Public Attributes inherited from Heed::gparticle
bool	s_life

long	nstep
	Step number.

double	total_range_from_origin
	Range from origin to currpos.

long	n_zero_step
	Number of previous steps with zero range (including this step).

stvpoint	origin

stvpoint	prevpos

stvpoint	currpos

stvpoint	nextpos

vec	curr_relcen

Public Attributes inherited from Heed::particle_type
PassivePtr< particle_def >	pardef

Additional Inherited Members
Static Public Member Functions inherited from Heed::RegPassivePtr
static void	set_s_ban_del_ignore (char fs_ban_del_ignore)

static char	get_s_ban_del_ignore (void)

static void	set_s_print_adr_cpp (char fs_print_adr_cpp)

static char	get_s_print_adr_cpp (void)

Static Public Attributes inherited from Heed::gparticle
static long	max_q_zero_step = 100

Protected Attributes inherited from Heed::eparticle
HeedFieldMap *	m_fieldMap

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 ( )

inline

Default constructor.

Definition at line 18 of file HeedDeltaElectron.h.

18: eparticle(), s_print_listing(false) {}

Heed::eparticle::eparticle

eparticle()

Default constructor.

Definition: eparticle.h:18

Referenced by copy().

◆ HeedDeltaElectron() [2/2]

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

Constructor.

Definition at line 36 of file HeedDeltaElectron.cpp.

    : eparticle(primvol, pt, vel, time, &electron_def, fieldmap),
      parent_particle_number(fparent_particle_number),
      particle_number(last_particle_number++),
      s_print_listing(fs_print_listing),
      phys_mrange(0.0),
      s_stop_eloss(false),
      s_mult_low_path_length(false),
      q_low_path_length(0.0),
      s_path_length(false),
      necessary_energy(0.0),
      total_Eloss(0.) {
  mfunname("HeedDeltaElectron::HeedDeltaElectron(...)");
}

◆ ~HeedDeltaElectron()

virtual Heed::HeedDeltaElectron::~HeedDeltaElectron ( )

inlinevirtual

Destructor.

Definition at line 24 of file HeedDeltaElectron.h.

24{}

Member Function Documentation

◆ copy()

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

inlinevirtual

Reimplemented from Heed::eparticle.

Definition at line 33 of file HeedDeltaElectron.h.

                                          {
    return new HeedDeltaElectron(*this);
  }

◆ physics_after_new_speed()

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

virtual

Reimplemented from Heed::gparticle.

Definition at line 153 of file HeedDeltaElectron.cpp.

                              {
  mfunname("void HeedDeltaElectron::physics_after_new_speed()");
  if (s_print_listing) {
    mcout << "HeedDeltaElectron::physics_after_new_speed is started\n";
    Iprint2n(mcout, currpos.prange, curr_kin_energy);
  }
  check_econd11(vecerror, != 0, mcerr);
  if (currpos.prange <= 0.0) {
    if (curr_kin_energy == 0.0) {
      // Get local volume.
      absvol* av = currpos.tid.G_lavol();
      if (av && av->s_sensitive && m_fieldMap->inside(currpos.ptloc)) {
        if (s_print_listing) {
          mcout << "HeedDeltaElectron::physics_after_new_speed: \n";
          mcout << "This is converted to conduction\n";
        }
        // TODO: replace push_back by emplace_back.
        conduction_electrons.push_back(
            HeedCondElectron(currpos.ptloc, currpos.time));
      }
      s_life = false;
    }
    if (s_print_listing) mcout << "exit due to currpos.prange <= 0.0\n";
    return;
  }
  if (s_print_listing) mcout << "physics_after_new_speed: started" << std::endl;
  // Get local volume and convert it to a cross-section object.
  const absvol* av = currpos.tid.G_lavol();
  const HeedDeltaElectronCS* hdecs =
      dynamic_cast<const HeedDeltaElectronCS*>(av);
  if (!hdecs) return;
  double ek = curr_kin_energy / MeV;
  if (s_print_listing) {
    Iprintnf(mcout, ek);
    Iprint3n(mcout, s_stop_eloss, phys_mrange, currpos.prange);
  }
  // Calculate dE/dx and energy loss. Update the kinetic energy.
  double dedx;
  double Eloss = 0.;
  if (s_stop_eloss && phys_mrange == currpos.prange) {
    Eloss = curr_kin_energy;
    curr_kin_energy = 0.0;
    dedx = Eloss / currpos.prange / (MeV / cm);
  } else {
    EnergyMesh* emesh = hdecs->hmd->energy_mesh.get();
    const long n1 = findInterval(emesh, ek);
    if (s_print_listing) Iprintn(mcout, n1);
    const double e1 = emesh->get_ec(n1);
    const double e2 = emesh->get_ec(n1 + 1);
    const double y1 = hdecs->eLoss[n1];
    const double y2 = hdecs->eLoss[n1 + 1];
    dedx = y1 + (y2 - y1) / (e2 - e1) * (ek - e1);
    Eloss = (dedx * MeV / cm) * currpos.prange;
    total_Eloss += Eloss;
    curr_kin_energy -= Eloss;
  }
  if (s_print_listing)
    Iprint3nf(mcout, prev_kin_energy / eV, curr_kin_energy / eV, Eloss / eV);
  if (curr_kin_energy <= 0.0) {
    if (s_print_listing) mcout << "curr_kin_energy <= 0.0\n";
    curr_kin_energy = 0.0;
    curr_gamma_1 = 0.0;
    currpos.speed = 0.0;
    s_life = false;
  } else {
    const double resten = mass * c_squared;
    curr_gamma_1 = curr_kin_energy / resten;
    currpos.speed = c_light * lorbeta(curr_gamma_1);
  }
  absvol* vav = currpos.tid.G_lavol();
  if (vav && vav->s_sensitive) {
    if (s_print_listing) {
      mcout << "volume is sensitive\n";
      Iprint2nf(mcout, Eloss / eV, necessary_energy / eV);
    }
    if (Eloss > 0.0) ionisation(Eloss, dedx, hdecs->pairprod.get());
  }
  if (s_print_listing) {
    mcout << '\n';
    Iprintn(mcout, s_life);
  }
  if (!s_life) {
    // Done tracing the delta electron. Create the last conduction electron.
    vav = currpos.tid.G_lavol();
    if (vav && vav->s_sensitive && m_fieldMap->inside(currpos.ptloc)) {
      if (s_print_listing) mcout << "Last conduction electron\n";
      conduction_electrons.push_back(
          HeedCondElectron(currpos.ptloc, currpos.time));
    }
    return;
  }
 
  if (s_print_listing) mcout << "\nstart to rotate by low angle\n";
  double ek_restr = std::max(ek, 0.0005);
  if (s_print_listing) Iprint2nf(mcout, currpos.prange, phys_mrange);
  if (currpos.prange < phys_mrange) {
    // recalculate scatterings
    s_path_length = false;
    if (s_low_mult_scattering) {
      EnergyMesh* emesh = hdecs->hmd->energy_mesh.get();
      const long n1r = findInterval(emesh, ek_restr);
      const double e1 = emesh->get_ec(n1r);
      const double e2 = emesh->get_ec(n1r + 1);
      const double y1 = hdecs->low_lambda[n1r];
      const double y2 = hdecs->low_lambda[n1r + 1];
      // Path length in internal units.
      double low_path_length =
          (y1 + (y2 - y1) / (e2 - e1) * (ek_restr - e1)) * cm;
      if (s_print_listing) Iprintnf(mcout, low_path_length / cm);
      s_mult_low_path_length = false;
      q_low_path_length = currpos.prange / low_path_length;
      if (s_print_listing) Iprintnf(mcout, q_low_path_length);
    }
  }
  if (s_print_listing) Iprintnf(mcout, q_low_path_length);
#ifdef RANDOM_POIS
  if (q_low_path_length > 0.0) {
    int ierror = 0;
    long random_q_low_path_length = pois(q_low_path_length, ierror);
    check_econd11a(ierror, == 1,
                   " q_low_path_length=" << q_low_path_length << '\n', mcerr);
    q_low_path_length = long(random_q_low_path_length);
    if (s_print_listing) {
      mcout << "After pois:\n";
      Iprintnf(mcout, q_low_path_length);
    }
  }
#endif
  if (q_low_path_length > 0) {
#ifdef DIRECT_LOW_IF_LITTLE
    const long max_q_low_path_length_for_direct = 5;
    if (q_low_path_length < max_q_low_path_length_for_direct) {
      // direct modeling
      if (s_print_listing) {
        mcout << "direct modeling of low scatterings\n";
        Iprint(mcout, currpos.dir);
      }
      EnergyMesh* emesh = hdecs->hmd->energy_mesh.get();
      const long n1r = findInterval(emesh, ek_restr);
      for (long nscat = 0; nscat < q_low_path_length; ++nscat) {
        if (s_print_listing) Iprintn(mcout, nscat);
        double theta_rot =
            hdecs->low_angular_points_ran[n1r].ran(SRANLUX()) * degree;
        if (s_print_listing) Iprintnf(mcout, theta_rot);
        vec dir = 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);
        currpos.dir = new_dir;
        if (s_print_listing) Iprint(mcout, new_dir);
      }
      currpos.dirloc = currpos.dir;
      currpos.tid.up_absref(&currpos.dirloc);
    } else {
#endif
      double sigma_ctheta = hdecs->get_sigma(ek_restr, q_low_path_length);
      // actually it is mean(1-cos(theta)) or
      // sqrt( mean( square(1-cos(theta) ) ) ) depending on USE_MEAN_COEF
 
      if (s_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 (s_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 (s_print_listing) Iprint2nf(mcout, x, ctheta);
        } while (ctheta <= -1.0);  // avoid absurd cos(theta)
        check_econd21(ctheta, < -1.0 ||, > 1.0, mcerr);
      }
      if (s_print_listing) Iprintnf(mcout, ctheta);
      double theta_rot = acos(ctheta);
      if (s_print_listing) Iprint2nf(mcout, theta_rot, theta_rot / degree);
      vec dir = 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);
      currpos.dir = new_dir;
      currpos.dirloc = currpos.dir;
      currpos.tid.up_absref(&currpos.dirloc);
    }
#ifdef DIRECT_LOW_IF_LITTLE
  }
#endif
  if (s_path_length) {
    if (s_print_listing) {
      mcout << "\nstarting to rotate by large angle" << std::endl;
      Iprintnf(mcout, s_path_length);
    }
    EnergyMesh* emesh = hdecs->hmd->energy_mesh.get();
    const long n1r = findInterval(emesh, ek_restr);
    double theta_rot = hdecs->angular_points_ran[n1r].ran(SRANLUX()) * degree;
    if (s_print_listing) Iprintnf(mcout, theta_rot);
    vec dir = 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);
    currpos.dir = new_dir;
    currpos.dirloc = currpos.dir;
    currpos.tid.up_absref(&currpos.dirloc);
  }
  if (s_print_listing) Iprint2nf(mcout, currpos.dir, currpos.dirloc);
}

◆ physics_mrange()

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

virtual

Reimplemented from Heed::gparticle.

Definition at line 55 of file HeedDeltaElectron.cpp.

                                                      {
  mfunname("void HeedDeltaElectron::physics_mrange(double& fmrange)");
  if (s_print_listing) mcout << "void HeedDeltaElectron::physics_mrange\n";
 
  s_mult_low_path_length = false;
  q_low_path_length = 0.0;
  s_path_length = false;
  if (fmrange <= 0.0) return;
  if (curr_kin_energy == 0.0) {
    fmrange = 0.0;
    return;
  }
  // Get local volume and convert it to a cross-section object.
  const absvol* av = currpos.tid.G_lavol();
  const HeedDeltaElectronCS* hdecs =
      dynamic_cast<const HeedDeltaElectronCS*>(av);
  if (!hdecs) return;
  if (s_print_listing) Iprintnf(mcout, fmrange);
  const double ek = curr_kin_energy / MeV;
  EnergyMesh* emesh = hdecs->hmd->energy_mesh.get();
  const long n1 = findInterval(emesh, ek);
  const double e1 = emesh->get_ec(n1);
  const double e2 = emesh->get_ec(n1 + 1);
  const double loss1 = hdecs->eLoss[n1];
  const double loss2 = hdecs->eLoss[n1 + 1];
  double dedx = loss1 + (loss2 - loss1) / (e2 - e1) * (ek - e1);
  // Min. loss 50 eV.
  double eloss = std::max(0.1 * ek, 0.00005);
  if (eloss > ek) {
    eloss = ek;
    s_stop_eloss = true;
  } else {
    s_stop_eloss = false;
  }
  fmrange = std::min(fmrange, (eloss / dedx) * cm);
  if (s_print_listing) Iprint2nf(mcout, fmrange, ek);
  const double ek_restr = std::max(ek, 0.0005);
  if (s_print_listing) Iprintnf(mcout, ek_restr / keV);
 
  const long n1r = findInterval(emesh, ek_restr);
  double low_path_length = 0.;  // in internal units
  if (s_low_mult_scattering) {
    const double e1r = emesh->get_ec(n1r);
    const double e2r = emesh->get_ec(n1r + 1);
    const double y1 = hdecs->low_lambda[n1r];
    const double y2 = hdecs->low_lambda[n1r + 1];
    low_path_length = (y1 + (y2 - y1) / (e2r - e1r) * (ek_restr - e1r)) * cm;
    if (s_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);
    double mult_low_path_length = qscat * low_path_length;
    if (s_print_listing) Iprintnf(mcout, mult_low_path_length);
    if (fmrange > mult_low_path_length) {
      fmrange = mult_low_path_length;
      s_mult_low_path_length = true;
      q_low_path_length = hdecs->eesls->get_qscat();
      s_stop_eloss = false;
    } else {
      s_mult_low_path_length = false;
      q_low_path_length = fmrange / low_path_length;
    }
    if (s_print_listing) Iprint2nf(mcout, fmrange, q_low_path_length);
  }
 
  if (s_high_mult_scattering) {
    const double e1r = emesh->get_ec(n1r);
    const double e2r = emesh->get_ec(n1r + 1);
    if (s_print_listing) {
      Iprintf(mcout, currpos.pt);
      Iprintnf(mcout, n1r);
      Iprintnf(mcout, ek_restr);
      Iprint2nf(mcout, e1r, e2r);
    }
    const double y1 = hdecs->lambda[n1r];
    const double y2 = hdecs->lambda[n1r + 1];
    const double mean_path = y1 + (y2 - y1) / (e2r - e1r) * (ek_restr - e1r);
    if (s_print_listing) Iprintnf(mcout, mean_path);
    const double path_length = -mean_path * cm * log(1.0 - SRANLUX());
    if (s_print_listing) Iprintnf(mcout, path_length);
    if (fmrange > path_length) {
      fmrange = path_length;
      s_path_length = true;
      s_mult_low_path_length = true;
      if (s_low_mult_scattering) {
        q_low_path_length = fmrange / low_path_length;
        if (s_print_listing) Iprintnf(mcout, q_low_path_length);
      }
      s_stop_eloss = false;
    } else {
      s_path_length = false;
    }
    if (s_print_listing) Iprintnf(mcout, fmrange);
  }
  phys_mrange = fmrange;
}

◆ print()

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

virtual

Reimplemented from Heed::eparticle.

Definition at line 434 of file HeedDeltaElectron.cpp.

                                                           {
  if (l < 0) return;
  Ifile << "HeedDeltaElectron (l=" << l
        << "): particle_number=" << 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=" << phys_mrange << " s_stop_eloss=" << s_stop_eloss
        << " s_mult_low_path_length=" << s_mult_low_path_length << '\n';
  Ifile << "q_low_path_length=" << q_low_path_length
        << " s_path_length=" << s_path_length
        << " necessary_energy/eV=" << 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 26 of file HeedDeltaElectron.h.

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

◆ conduction_ions

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

Definition at line 27 of file HeedDeltaElectron.h.

Referenced by Garfield::TrackHeed::GetCluster(), Garfield::TrackHeed::TransportDeltaElectron(), and Garfield::TrackHeed::TransportPhoton().

◆ parent_particle_number

long Heed::HeedDeltaElectron::parent_particle_number

Definition at line 29 of file HeedDeltaElectron.h.

Referenced by print().

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

Public Member Functions

Public Attributes

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