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
virtual void	print (std::ostream &file, int l) const
virtual HeedDeltaElectronCS *	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 Attributes
PassivePtr< HeedMatterDef >	hmd
PassivePtr< ElElasticScat >	ees
PassivePtr< ElElasticScatLowSigma >	eesls
PassivePtr< PairProd >	pairprod
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 const long	q_angular_mesh = 50
static const double	low_cut_angle_deg = 20.0

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)

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

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

inlinevirtual

Reimplemented from Heed::RegPassivePtr.

Reimplemented in Garfield::HeedChamber.

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 230 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);
}

Referenced by Heed::HeedDeltaElectron::physics_after_new_speed(), and Heed::HeedDeltaElectron::physics_mrange().

print()

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

virtual

Reimplemented from Heed::RegPassivePtr.

Definition at line 260 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 82 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

angular_points_ran

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

Definition at line 85 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and Heed::HeedDeltaElectron::physics_after_new_speed().

beta

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

Table of velocities.

Definition at line 47 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

ees

PassivePtr<ElElasticScat> Heed::HeedDeltaElectronCS::ees

Definition at line 43 of file HeedDeltaElectronCS.h.

Referenced by get_sigma(), and HeedDeltaElectronCS().

eesls

PassivePtr<ElElasticScatLowSigma> Heed::HeedDeltaElectronCS::eesls

Definition at line 44 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), Heed::HeedDeltaElectron::physics_mrange(), 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 54 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), Heed::HeedDeltaElectron::physics_after_new_speed(), Heed::HeedDeltaElectron::physics_mrange(), and print().

hmd

PassivePtr<HeedMatterDef> Heed::HeedDeltaElectronCS::hmd

Definition at line 42 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), Heed::HeedPhoton::physics(), Heed::HeedDeltaElectron::physics_after_new_speed(), Heed::HeedPhoton::physics_after_new_speed(), Heed::HeedDeltaElectron::physics_mrange(), 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 61 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), Heed::HeedDeltaElectron::physics_mrange(), and print().

low_angular_points_ran

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

Definition at line 87 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and Heed::HeedDeltaElectron::physics_after_new_speed().

low_cut_angle_deg

const double Heed::HeedDeltaElectronCS::low_cut_angle_deg = 20.0

static

Definition at line 40 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 64 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), Heed::HeedDeltaElectron::physics_after_new_speed(), Heed::HeedDeltaElectron::physics_mrange(), and print().

mean_coef_low_sigma

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

Definition at line 91 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 58 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

momentum

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

Table of momenta [MeV/c].

Definition at line 49 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 79 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

pairprod

PassivePtr<PairProd> Heed::HeedDeltaElectronCS::pairprod

Definition at line 45 of file HeedDeltaElectronCS.h.

Referenced by Heed::HeedDeltaElectron::physics_after_new_speed().

q_angular_mesh

const long Heed::HeedDeltaElectronCS::q_angular_mesh = 50

static

Definition at line 39 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 70 of file HeedDeltaElectronCS.h.

Referenced by HeedDeltaElectronCS(), and print().

---

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

• HeedDeltaElectronCS.h
• HeedDeltaElectronCS.cpp