Heed::HeedDeltaElectronCS Class Reference

#include <HeedDeltaElectronCS.h>

Inheritance diagram for Heed::HeedDeltaElectronCS:

Public Member Functions
	HeedDeltaElectronCS ()
	Default constructor.
	HeedDeltaElectronCS (HeedMatterDef *fhmd, ElElasticScat *fees, ElElasticScatLowSigma *feesls, PairProd *fpairprod, int fsruth=2, double fmlambda=0.001 *4.0e-3, double fmthetac=0.1)
	Constructor.
double	get_sigma (double energy, double nscat) const
void	print (std::ostream &file, int l) const
HeedDeltaElectronCS *	copy () const

Public Attributes
HeedMatterDef *	hmd = nullptr
ElElasticScat *	ees = nullptr
ElElasticScatLowSigma *	eesls = nullptr
PairProd *	pairprod = nullptr
std::vector< double >	beta
	Table of velocities.
std::vector< double >	momentum
	Table of momenta [MeV/c].
std::vector< double >	eLoss
double	mlambda
std::vector< double >	lambda
std::vector< double >	low_lambda
int	sruth
double	mthetac
std::vector< double >	angular_mesh_c
	Angular mesh, centers, angles in degrees.
std::vector< PointsRan >	angular_points_ran
std::vector< PointsRan >	low_angular_points_ran
std::vector< double >	mean_coef_low_sigma

Static Public Attributes
static constexpr long	q_angular_mesh = 50
static constexpr double	low_cut_angle_deg = 20.

Detailed Description

Cross sections and various parameters for delta-electron transport. 2003, I. Smirnov

Definition at line 21 of file HeedDeltaElectronCS.h.

Constructor & Destructor Documentation

HeedDeltaElectronCS() [1/2]

Heed::HeedDeltaElectronCS::HeedDeltaElectronCS

(

)

Default constructor.

Referenced by copy().

HeedDeltaElectronCS() [2/2]

Heed::HeedDeltaElectronCS::HeedDeltaElectronCS	(	HeedMatterDef *	fhmd,
		ElElasticScat *	fees,
		ElElasticScatLowSigma *	feesls,
		PairProd *	fpairprod,
		int	fsruth = 2,
		double	fmlambda = 0.001 * 4.0e-3,
		double	fmthetac = 0.1
	)

Constructor.

Definition at line 24 of file HeedDeltaElectronCS.cpp.

    : hmd(fhmd),
      ees(fees),
      eesls(feesls),
      pairprod(fpairprod),
      mlambda(fmlambda),
      sruth(fsruth),
      mthetac(fmthetac) {
  mfunname("HeedDeltaElectronCS::HeedDeltaElectronCS(...)");
 
  const long qe = hmd->energy_mesh->get_q();
  eLoss.resize(qe, 0.);
  beta.resize(qe, 0.);
  momentum.resize(qe, 0.);
  std::vector<double> beta2(qe, 0.);
  std::vector<double> momentum2(qe, 0.);
 
  const double rho = hmd->matter->density();
  const double zmean = hmd->matter->Z_mean();
  const double amean = hmd->matter->A_mean();
  const double ZA = zmean / amean;
  const double I_eff = 15.8 * eV * zmean;
  double smax = 0.0;
  for (long ne = 0; ne < qe; ne++) {
    const double ec = hmd->energy_mesh->get_ec(ne) * MeV;
    const double gamma_1 = ec / electron_mass_c2;
    beta[ne] = lorbeta(gamma_1);
    beta2[ne] = beta[ne] * beta[ne];
    const double en = electron_mass_c2 + ec;
    momentum2[ne] =
        (en * en - electron_mass_c2 * electron_mass_c2) / (MeV * MeV);
    momentum[ne] = sqrt(momentum2[ne]);
    const double dedx = e_cont_enloss(ZA, I_eff, rho, ec, DBL_MAX, -1);
    if (smax < dedx) smax = dedx;
    eLoss[ne] = dedx / (MeV / cm);
  }
  smax = smax / (MeV / cm);
  // looking for minimal (maximal by value but negative) derivative
  double deriv_min = 0.0;
  long n_deriv_min = 0;
  for (long ne = 1; ne < qe; ne++) {
    double d =
        (eLoss[ne] * beta[ne] - eLoss[ne - 1] * beta[ne - 1]) /
        (hmd->energy_mesh->get_ec(ne) - hmd->energy_mesh->get_ec(ne - 1));
    if (deriv_min > d) {
      deriv_min = d;
      n_deriv_min = ne - 1;
    }
  }
  for (long ne = 0; ne < n_deriv_min - 1; ne++) {
    eLoss[ne] = (eLoss[n_deriv_min - 1] * beta[n_deriv_min - 1] +
                 deriv_min * (hmd->energy_mesh->get_ec(ne) -
                              hmd->energy_mesh->get_ec(n_deriv_min - 1))) /
                beta[ne];
  }
 
  lambda.resize(qe);
  if (sruth == 2) {
    angular_mesh_c.resize(q_angular_mesh);
    angular_mesh_c[0] = 0.0;
    const double amin = 0.3;
    const double amax = 180.0;
    const double rk = pow(amax / amin, (1.0 / (q_angular_mesh - 2.)));
    double ar = amin;
    angular_mesh_c[1] = ar;
    for (long n = 2; n < q_angular_mesh; n++) {
      angular_mesh_c[n] = angular_mesh_c[n - 1] * rk;
    }
    angular_mesh_c[q_angular_mesh - 1] = 180.0;
    angular_points_ran.resize(qe);
    low_angular_points_ran.resize(qe);
    low_lambda.resize(qe);
    const long qes = eesls->get_ees()->get_qe();
#ifdef USE_MEAN_COEF
    mean_coef_low_sigma.resize(qes);
#else
    coef_low_sigma.resize(qes);
#endif
    for (long ne = 0; ne < qes; ne++) {
      long qat = hmd->matter->qatom();
      double s = 0.0;
      for (long nat = 0; nat < qat; nat++) {
#ifdef USE_MEAN_COEF
        s += eesls->get_mean_coef(hmd->matter->atom(nat)->Z(), ne) *
             hmd->matter->weight_quan(nat);
#else
        s += eesls->get_coef(hmd->matter->atom(nat)->Z(), ne) *
             hmd->matter->weight_quan(nat);
#endif
      }
#ifdef USE_MEAN_COEF
      mean_coef_low_sigma[ne] = s;
#else
      coef_low_sigma[ne] = s;
#endif
    }
  }
 
  for (long ne = 0; ne < qe; ne++) {
    double rr;
    double ek = hmd->energy_mesh->get_ec(ne) * 1000.0;
    if (ek <= 10.0) {
      rr = 1.0e-3 * amean / (gram / mole) / zmean * 3.872e-3 * pow(ek, 1.492);
    } else {
      rr = 1.0e-3 * 6.97e-3 * pow(ek, 1.6);
    }
    rr = rr / (rho / (gram / cm3));
    rr = rr * 0.1;
    // Iprintn(mcout, rr);
    double cor = 1.0;
    {
      // b-k*(x-a)**2 = 0  =>  x= a +- sqrt(b/k)
      // k = b / (x - a)**2
      double a = 2.5;
      double b = 4;
      // k=1.0/4.0
      double x = 0.0;
      double k = b / ((x - a) * (x - a));
      x = ek * 1000.0;
      double r = b - k * (x - a) * (x - a);
      if (r < 0.0)
        r = 1;
      else
        r = r + 1;
      cor = r;
    }
    if (sruth == 1) {
      lambda[ne] = mlambda / (rho / (gram / cm3));
      if (lambda[ne] < rr) lambda[ne] = rr;
      lambda[ne] = lambda[ne] * cor;
      // Calculate the minimum angle for restriction of field by atomic shell
      double mT = 2.0 * asin(1.0 / (2.0 * momentum[ne] * zmean * 5.07e2));
      // Throw out too slow interactions. They do not have any influence.
      if (mT < mthetac) mT = mthetac;
      // Calculate the cut angle due to mean free part
      double A = hmd->Rutherford_const / cor / (momentum2[ne] * beta2[ne]) /
                 pow(5.07e10, 2);
      double B = lambda[ne] * A;
      B = sqrt(B / (B + 1.0));
      // Threshold turn angle
      double thetac = 2.0 * asin(B);
 
      // If it too little, reset it. It will lead to increasing
      // of lambda and decreasing of calculation time.
      if (thetac < mT) {
        thetac = mT;
        B = mT;  // B is double precision
        double r = sin(0.5 * B);
        lambda[ne] = 1 / A * 2.0 * r * r / (1 + cos(B));
        // r=cos(TetacBdel(nen,nm))
        // lamBdel=A*(1.0+r)/(1.0-r)
        // lamBdel=1.0/lamBdel
        // lamBdel=(p2*bet2*sin(TetacBdel/2.0)**2) / A
      }
      // const double CosThetac12 = cos(0.5 * thetac);
      // const dobule SinThetac12 = sin(0.5 * thetac);
      // Flag that scattering is spherical.
      // const bool sisfera = thetac > 1.5 ? true : false;
    } else if (sruth == 0) {
      // Gauss formula
      const double msig = sqrt(2.0) * 13.6 / (beta[ne] * momentum[ne]);
      double x = mthetac / msig;
      x = x * x;
      x = x * hmd->radiation_length * cor;
      lambda[ne] = mlambda / (rho / (gram / cm3));
      if (lambda[ne] < rr) lambda[ne] = rr;
      lambda[ne] = lambda[ne] * cor;
      if (lambda[ne] < x) lambda[ne] = x;
    } else if (sruth == 2) {
      // Cross section for material per one av. atom, in cm^2/rad.
      std::vector<double> smat(q_angular_mesh, 0.);
      const double energy = hmd->energy_mesh->get_ec(ne);
      const long qat = hmd->matter->qatom();
      for (long nan = 0; nan < q_angular_mesh; nan++) {
        const double angle = angular_mesh_c[nan] * degree;
        double s = 0.0;
        for (long nat = 0; nat < qat; nat++) {
          const int z = hmd->matter->atom(nat)->Z();
          const double w = hmd->matter->weight_quan(nat);
          s += ees->get_CS(z, energy, angle) * w;
        }
        s = s * 1.0E-16;
        // s = s * 1.0E-16 * C1_MEV_CM * C1_MEV_CM;
        // Angstrem**2 -> cm**2
        // cm**2 -> MeV**-2
        smat[nan] = s * twopi * sin(angle);  // sr -> dtheta
      }
      angular_points_ran[ne] =
          PointsRan(angular_mesh_c, smat, low_cut_angle_deg, 180);
      low_angular_points_ran[ne] =
          PointsRan(angular_mesh_c, smat, 0., low_cut_angle_deg);
      const double coef =
          Avogadro * rho / (gram / cm3) / (amean / (gram / mole));
      lambda[ne] =
          1. / (angular_points_ran[ne].get_integ_active() * degree * coef);
      low_lambda[ne] =
          1. / (low_angular_points_ran[ne].get_integ_active() * degree * coef);
    }
  }
}

Member Function Documentation

copy()

HeedDeltaElectronCS * Heed::HeedDeltaElectronCS::copy ( ) const

inline

Definition at line 35 of file HeedDeltaElectronCS.h.

{ return new HeedDeltaElectronCS(*this); }

get_sigma()

double Heed::HeedDeltaElectronCS::get_sigma	(	double	energy,
		double	nscat
	)		const

Definition at line 229 of file HeedDeltaElectronCS.cpp.

                                                                       {
  mfunname("double HeedDeltaElectronCS::get_sigma(...)");
  check_econd11(nscat, < 0, mcerr);
  // check_econd21(nscat , < 0 || , > eesls->get_qscat() , mcerr);
  // ^ not compatible with Poisson
  const long qe = ees->get_qe();
  double energyKeV = energy * 1000.0;
  energyKeV = std::max(energyKeV, ees->get_energy_mesh(0));
  energyKeV = std::min(energyKeV, ees->get_energy_mesh(qe - 1));
  long n1 = 0;
  long n2 = qe - 1;
  while (n2 - n1 > 1) {
    const long n3 = n1 + (n2 - n1) / 2;
    if (energyKeV < ees->get_energy_mesh(n3))
      n2 = n3;
    else
      n1 = n3;
  }
  const double e1 = ees->get_energy_mesh(n1);
  const double e2 = ees->get_energy_mesh(n2);
#ifdef USE_MEAN_COEF
  const double v1 = nscat * mean_coef_low_sigma[n1];
  const double v2 = nscat * mean_coef_low_sigma[n2];
#else
  const double v1 = nscat * coef_low_sigma[n1];
  const double v2 = nscat * coef_low_sigma[n2];
#endif
  return v1 + (v2 - v1) / (e2 - e1) * (energyKeV - e1);
}

print()

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

Definition at line 259 of file HeedDeltaElectronCS.cpp.

                                                             {
  if (l <= 0) return;
  Ifile << "HeedDeltaElectronCS(l=" << l << "):";
  long qe = hmd->energy_mesh->get_q();
  // Iprintn(mcout, qe);
  // mcout<<std::endl;
  Iprintn(file, mlambda);
  Iprintn(file, mthetac);
  Iprintn(file, sruth);
  Ifile << "         get_ec,        beta,      momentum,    eLoss,    lambda,  "
           " low_lambda:" << std::endl;
  indn.n += 2;
  for (long ne = 0; ne < qe; ne++) {
    Ifile << std::setw(3) << ne << ' ' << std::setw(12)
          << hmd->energy_mesh->get_ec(ne) << ' ' << std::setw(12) << beta[ne]
          << ' ' << std::setw(12) << momentum[ne] << ' ' << std::setw(12)
          << eLoss[ne] << ' ' << std::setw(12) << lambda[ne] << ' '
          << std::setw(12) << low_lambda[ne] << '\n';
  }
  indn.n -= 2;
  Ifile << "na, angular_mesh_c:" << std::endl;
  indn.n += 2;
  for (long na = 0; na < q_angular_mesh; na++) {
    Ifile << na << ' ' << std::setw(12) << angular_mesh_c[na] << '\n';
  }
  indn.n -= 2;
  Iprintn(file, eesls->get_ees()->get_qe());
  indn.n += 2;
#ifdef USE_MEAN_COEF
  Ifile << "ne, energy_mesh(ne), mean_coef_low_sigma:" << std::endl;
  for (long ne = 0; ne < eesls->get_ees()->get_qe(); ne++) {
    Ifile << std::setw(3) << ne << ' ' << std::setw(12)
          << eesls->get_ees()->get_energy_mesh(ne) << " KeV " << std::setw(12)
          << mean_coef_low_sigma[ne] << '\n';
  }
#else
  Ifile << "ne, energy_mesh(ne), coef_low_sigma:" << std::endl;
  for (long ne = 0; ne < eesls->get_ees()->get_qe(); ne++) {
    Ifile << std::setw(3) << ne << ' ' << std::setw(12)
          << eesls->get_ees()->get_energy_mesh(ne) << " KeV " << std::setw(12)
          << coef_low_sigma[ne] << '\n';
  }
#endif
  indn.n -= 2;
}

Member Data Documentation

angular_mesh_c

std::vector<double> Heed::HeedDeltaElectronCS::angular_mesh_c

Angular mesh, centers, angles in degrees.

Definition at line 80 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

angular_points_ran

std::vector<PointsRan> Heed::HeedDeltaElectronCS::angular_points_ran

Definition at line 83 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS().

beta

std::vector<double> Heed::HeedDeltaElectronCS::beta

Table of velocities.

Definition at line 45 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

ees

ElElasticScat* Heed::HeedDeltaElectronCS::ees = nullptr

Definition at line 41 of file HeedDeltaElectronCS.h.

Referenced by get_sigma(), and HeedDeltaElectronCS().

eesls

ElElasticScatLowSigma* Heed::HeedDeltaElectronCS::eesls = nullptr

Definition at line 42 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

eLoss

std::vector<double> Heed::HeedDeltaElectronCS::eLoss

Energy losses [MeV/cm] according to very advanced formula with Bethe-Bloch and density effect as in GEANT3, corresponding to centers of intervals of the common mesh.

Definition at line 52 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

hmd

HeedMatterDef* Heed::HeedDeltaElectronCS::hmd = nullptr

Definition at line 40 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

lambda

std::vector<double> Heed::HeedDeltaElectronCS::lambda

Path length [cm] at the energy values of the mesh. For sruth == 2 mean free path.

Definition at line 59 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

low_angular_points_ran

std::vector<PointsRan> Heed::HeedDeltaElectronCS::low_angular_points_ran

Definition at line 85 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS().

low_cut_angle_deg

constexpr double Heed::HeedDeltaElectronCS::low_cut_angle_deg = 20.

staticconstexpr

Definition at line 38 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS().

low_lambda

std::vector<double> Heed::HeedDeltaElectronCS::low_lambda

Path length for low angle scatterings. This is without multiplication used for coef_low_sigma (cm).

Definition at line 62 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

mean_coef_low_sigma

std::vector<double> Heed::HeedDeltaElectronCS::mean_coef_low_sigma

Definition at line 89 of file HeedDeltaElectronCS.h.

Referenced by get_sigma(), HeedDeltaElectronCS(), and print().

mlambda

double Heed::HeedDeltaElectronCS::mlambda

Mminimum mean length of range, multiplied by density. sm*gr/sm**3 = gr/sm**2

Definition at line 56 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

momentum

std::vector<double> Heed::HeedDeltaElectronCS::momentum

Table of momenta [MeV/c].

Definition at line 47 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

mthetac

double Heed::HeedDeltaElectronCS::mthetac

Minimum threshold turn angle. For Rutherford: the interactions with less angle will not be taken into account. The actual threshold angle can be larger. The second restriction is going from restriction of atomic shell. The third one is from mlamBdel. For usual multiple scattering: Assuming that sigma = mTetacBdel the path lengt is calculated. If mlamBdel/density is less then the last is used.

Definition at line 77 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

pairprod

PairProd* Heed::HeedDeltaElectronCS::pairprod = nullptr

Definition at line 43 of file HeedDeltaElectronCS.h.

q_angular_mesh

constexpr long Heed::HeedDeltaElectronCS::q_angular_mesh = 50

staticconstexpr

Definition at line 37 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

sruth

int Heed::HeedDeltaElectronCS::sruth

Formula to use. 0 - usual multiple scattering formula 1 - Rutherford cross-section 2 - precise theory?

Definition at line 68 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

---

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

• HeedDeltaElectronCS.h
• HeedDeltaElectronCS.cpp