Heed::ExAtomPhotoAbsCS Class Reference

Atomic photo-absorption with excitation. More...

#include <PhotoAbsCS.h>

Inheritance diagram for Heed::ExAtomPhotoAbsCS:

Public Member Functions
virtual double	get_threshold (int nshell) const
	Get the ionization threshold for a given shell.
virtual double	get_ACS (double energy) const
virtual double	get_integral_ACS (double energy1, double energy2) const
	Integrated photo-absorption cross-section overa given interval.
virtual double	get_ACS (int nshell, double energy) const
	Sub-shell photo-absorption cross-section [Mbarn] at a given energy [MeV].
virtual double	get_integral_ACS (int nshell, double energy1, double energy2) const
	Integrated sub-shell photo-absorption cross-section.
virtual double	get_ICS (double energy) const
virtual double	get_integral_ICS (double energy1, double energy2) const
	Integrated photo-ionization cross-section over a given interval.
virtual double	get_ICS (int nshell, double energy) const
	Sub-shell photo-ionization cross-section at a given energy.
virtual double	get_integral_ICS (int nshell, double energy1, double energy2) const
	Integrated sub-shell photo-ionization cross-section.
virtual int	get_main_shell_number (int nshell) const
void	replace_shells_by_average (double fwidth, double fstep, long fmax_q_step)
virtual void	print (std::ostream &file, int l) const
virtual ExAtomPhotoAbsCS *	copy () const
	ExAtomPhotoAbsCS ()
	Default constructor.
	ExAtomPhotoAbsCS (int fZ, const std::string &fthreshold_file_name, const std::string &fsimple_table_file_name, const std::string &fname="none", double fminimal_threshold=0.0)
	ExAtomPhotoAbsCS (int fZ, const std::string &fname, const std::string &fBT_file_name, int id, double fminimal_threshold=0.0)
	ExAtomPhotoAbsCS (int fZ, const std::string &fname, const std::string &fFitBT_file_name, int id, int s_no_scale, double fminimal_threshold=0.0)
	ExAtomPhotoAbsCS (int fZ, const std::string &fname, const std::string &fFitBT_file_name, const std::string &fsimple_table_file_name, double emax_repl, int id, double fminimal_threshold=0.0)
virtual	~ExAtomPhotoAbsCS ()
	Destructor.
Public Member Functions inherited from Heed::AtomPhotoAbsCS
	AtomPhotoAbsCS ()
	Default constructor.
int	get_Z () const
	Get the atomic number.
unsigned int	get_qshell () const
	Get the number of shells.
virtual double	get_threshold (int nshell) const =0
	Get the ionization threshold for a given shell.
virtual double	get_I_min () const
	Get the lowest ionization threshold among all shells.
virtual double	get_ACS (double energy) const =0
virtual double	get_integral_ACS (double energy1, double energy2) const =0
	Integrated photo-absorption cross-section overa given interval.
virtual double	get_ACS (int nshell, double energy) const =0
	Sub-shell photo-absorption cross-section [Mbarn] at a given energy [MeV].
virtual double	get_integral_ACS (int nshell, double energy1, double energy2) const =0
	Integrated sub-shell photo-absorption cross-section.
virtual double	get_ICS (double energy) const =0
virtual double	get_TICS (double energy, double factual_minimal_threshold) const
virtual double	get_integral_ICS (double energy1, double energy2) const =0
	Integrated photo-ionization cross-section over a given interval.
virtual double	get_integral_TICS (double energy1, double energy2, double factual_minimal_threshold) const
	Integral photo-ionization cross-section with redefined threshold.
virtual double	get_ICS (int nshell, double energy) const =0
	Sub-shell photo-ionization cross-section at a given energy.
virtual double	get_TICS (int nshell, double energy, double factual_minimal_threshold) const
	Sub-shell photo-ionization cross-section with redefined threshold.
virtual double	get_integral_ICS (int nshell, double energy1, double energy2) const =0
	Integrated sub-shell photo-ionization cross-section.
virtual double	get_integral_TICS (int nshell, double energy1, double energy2, double factual_minimal_threshold) const
	Integrated sub-shell photo-ionization cross-section (redefined threshold).
virtual void	get_escape_particles (const int nshell, double energy, std::vector< double > &el_energy, std::vector< double > &ph_energy) const
virtual int	get_main_shell_number (int nshell) const =0
virtual void	remove_shell (int nshell)
	Deactivate a sub-shell. Set s_ignore_shell flag to true.
virtual void	restore_shell (int nshell)
	Activate a sub-shell. Set s_ignore_shell flag to false.
virtual void	print (std::ostream &file, int l) const
virtual AtomPhotoAbsCS *	copy () const =0
AtomicSecondaryProducts *	get_asp (int nshell)

Protected Attributes
std::string	threshold_file_name
std::string	simple_table_file_name
std::string	BT_file_name
std::vector< std::shared_ptr< PhotoAbsCS > >	m_acs
double	integ_abs_before_corr = 0.
double	integ_abs_after_corr = 0.
double	integ_ioniz_after_corr = 0.
double	height_of_excitation = 0.
	Excitation cross-section (assumed in the lowest shell).
double	exener [2] = {0., 0.}
	Boundaries of excitation.
double	minimal_threshold = 0.
Protected Attributes inherited from Heed::AtomPhotoAbsCS
std::string	name
	Name of the atom.
int	Z
	Atomic number.
int	qshell
	Number of shells.
std::vector< bool >	s_ignore_shell
std::vector< AtomicSecondaryProducts >	asp
	Sampling of relaxation products for each shell.

Static Protected Attributes
static const int	s_add_excitations_to_normalize = 1
static const int	s_scale_to_normalize_if_more = 1

Detailed Description

Atomic photo-absorption with excitation.

Definition at line 430 of file PhotoAbsCS.h.

Constructor & Destructor Documentation

ExAtomPhotoAbsCS() [1/5]

Heed::ExAtomPhotoAbsCS::ExAtomPhotoAbsCS ( )

inline

Default constructor.

Definition at line 455 of file PhotoAbsCS.h.

: AtomPhotoAbsCS() {}

Referenced by copy(), and ExAtomPhotoAbsCS().

ExAtomPhotoAbsCS() [2/5]

Heed::ExAtomPhotoAbsCS::ExAtomPhotoAbsCS	(	int	fZ,
		const std::string &	fthreshold_file_name,
		const std::string &	fsimple_table_file_name,
		const std::string &	fname = "none",
		double	fminimal_threshold = 0.0
	)

Constructor,

Parameters

fZ	atomic number
fthreshold_file_name	file from which to read name and shell energies
fsimple_table_file_name	file from which to read the cross-sections
fname	name of the atom, if "none" it is taken from fthreshold_file_name
fminimal_threshold	threshold

Definition at line 1027 of file PhotoAbsCS.cpp.

    : threshold_file_name(fthreshold_file_name),
      simple_table_file_name(fsimple_table_file_name),
      BT_file_name("none"),
      minimal_threshold(fminimal_threshold) {
  mfunnamep("ExAtomPhotoAbsCS::ExAtomPhotoAbsCS(...) const");
  check_econd11(fZ, < 1, mcerr);
  std::ifstream threshold_file(threshold_file_name.c_str());
  if (!threshold_file) {
    funnw.ehdr(mcerr);
    mcerr << "cannot open file " << threshold_file_name << std::endl;
    spexit(mcerr);
  }
  std::vector<double> thr;
  std::vector<int> Zshell;
  std::vector<double> fl;
  std::vector<std::string> shell_name;
  bool foundZ = false;
  while (findmark(threshold_file, "#") == 1) {
    threshold_file >> Z;
    if (Z != fZ) continue;
    threshold_file >> qshell;
    check_econd21(qshell, < 1 ||, > 10000, mcerr);
    s_ignore_shell.resize(qshell, false);
    // Iprintn(mcout, qshell);
    thr.resize(qshell, 0.0);
    Zshell.resize(qshell, 0);
    fl.resize(qshell, 0.0);
    shell_name.resize(qshell);
    m_acs.resize(qshell);
    asp.resize(qshell);
    std::string temp_name;
    threshold_file >> temp_name;
    name = fname == "none" ? temp_name : fname;
    int sZshell = 0;
    for (int nshell = 0; nshell < qshell; nshell++) {
      threshold_file >> thr[nshell];
      check_econd11(thr[nshell], <= 0.0, mcerr);
      thr[nshell] *= 1.0e-6;
      threshold_file >> Zshell[nshell];
      check_econd11(Zshell[nshell], <= 0, mcerr);
      sZshell += Zshell[nshell];
      threshold_file >> fl[nshell];
      findmark(threshold_file, "!");
      threshold_file >> shell_name[nshell];
    }
    check_econd12(sZshell, !=, Z, mcerr);
    // currently the minimal shell is the last,
    // but to avoid this assumption, we check all.
    int n_min = 0;
    double st = DBL_MAX;
    for (int nshell = 0; nshell < qshell; nshell++) {
      if (thr[nshell] < st) {
        n_min = nshell;
        st = thr[nshell];
      }
    }
    for (int nshell = 0; nshell < qshell; nshell++) {
      if (fl[nshell] <= 0) continue;
      check_econd12(nshell, ==, n_min, mcerr);
      std::vector<double> felectron_energy;
      std::vector<double> fphoton_energy;
      fphoton_energy.push_back(thr[nshell] - thr[n_min]);
      asp[nshell].add_channel(fl[nshell], felectron_energy, fphoton_energy);
    }
    foundZ = true;
    break;
  }
  if (!foundZ) {
    funnw.ehdr(mcerr);
    mcerr << "there is no element Z=" << fZ << " in file "
          << threshold_file_name << std::endl;
    spexit(mcerr);
  }
  // Here it reads the PACS as an one shell curve:
  SimpleTablePhotoAbsCS stpacs(name, Z, 0.0, fsimple_table_file_name);
  const std::vector<double>& ener = stpacs.get_arr_ener();
  const std::vector<double>& CS = stpacs.get_arr_CS();
  std::vector<double> left_CS = CS;  // used in sequencial algorithm
  const long qe = ener.size();
  // here cs is saved
  std::vector<std::vector<double> > SCS(qshell, std::vector<double>(qe, 0.));
  int nct = qshell - 1;   // "current" threshold index
  unsigned long nce = 0;  // "current" energy index
  // We ignore values below the lowest threshold.
  // It is not clear whether it is right, perhaps this is one of possible ways
  if (ener[0] < thr[qshell - 1]) {
    for (long ne = 0; ne < qe && (ener[ne] < thr[qshell - 1] ||
                                  (ne > 1 && CS[ne - 1] <= CS[ne - 2]));
         ne++) {
      if (ne > 0) left_CS[ne - 1] = 0.0;
      nce = ne;
    }
  }
  int s_more;
  int nt2 = 0;  // < nt1
  int s_spes = 0;
  // Actually this is a loop by the group of thresholds
  do {
    // Find all thresholds which are less then the current energy
    int nt;
    // sign that there are thresholds more than the current energy
    s_more = 0;
    for (nt = nct; nt >= 0; nt--) {
      if (s_spes == 0) {
        if (thr[nt] > ener[nce]) {
          s_more = 1;
          break;
        }
      } else {
        if (thr[nt] > ener[nce + 1]) {
          s_more = 1;
          break;
        }
      }
    }
    // nt is now index of the next threshold or -1, if the thresholds are
    // made up.
    int nt1 = nct;
    int nce_next = ener.size();
    nt2 = nt + 1;
    // mcout<<"nt="<<nt<<" nt1="<<nt1<<" nt2="<<nt2<<" s_more="<<s_more<<'\n';
    if (s_more == 1) {
      // if(nt >= 0)  // so if there are other larger thresholds,
      //{        // we should check how far we can pass at this step
      unsigned long ne;
      // finding energy larger than the next threshold
      for (ne = nce; ne < ener.size(); ne++) {
        if (thr[nt] <= ener[ne]) {
          nce_next = ne;
          s_spes = 0;
          break;
        }
        // At the following condition energy could be less than threshold,
        // but this point will anyway be associated with the next shell
        // corresponding to this threshold.
        // This is related to not precise measurement of cross section
        // and not precise knowledge of shell energies.
        // Occurence of this condition is marked by s_spes = 1.
        // At the next passing of this loop the thresholds are compared with
        // the next energy.
        if (ne > 1 && ne < ener.size() - 1 && ne > nce + 2 &&
            thr[nt] <= ener[ne + 1] &&
            (thr[nt] - ener[ne]) / (ener[ne + 1] - ener[ne]) < 0.1 &&
            CS[ne] > CS[ne - 1]) {
          // mcout<<"special condition is satisf.\n";
          nce_next = ne;
          s_spes = 1;
          break;
        }
      }
      if (ne == ener.size())  // threshold is larger then energy mesh
        s_more = 0;           // to finish the loop
    }
    // Iprintn(mcout, nce_next);
    // Iprintn(mcout, ener[nce_next-1]);
    int qt = nt1 - nt2 + 1;  // quantity of the new thresholds
    // Iprintn(mcout, qt);
 
    // Calculate sum of Z.
    int s = 0;
    for (nt = 0; nt < qt; nt++) {
      const int nshell = nct - nt;
      s += Zshell[nshell];
    }
    // Weights according to charges
    std::vector<double> w(qt);
    for (nt = 0; nt < qt; nt++) {
      const int nshell = nct - nt;
      w[nt] = double(Zshell[nshell]) / s;
    }
    double save_left_CS = left_CS[nce_next - 1];
    for (long ne = 0; ne < nce_next; ne++) {
      for (nt = 0; nt < qt; nt++) {
        int nshell = nct - nt;
        SCS[nshell][ne] = left_CS[ne] * w[nt];
      }
      left_CS[ne] = 0.0;
    }
    for (unsigned long ne = nce_next; ne < ener.size(); ne++) {
      double extrap_CS =
          save_left_CS * pow(ener[nce_next - 1], 2.75) / pow(ener[ne], 2.75);
      if (extrap_CS > left_CS[ne]) extrap_CS = left_CS[ne];
      for (nt = 0; nt < qt; nt++) {
        int nshell = nct - nt;
        SCS[nshell][ne] = extrap_CS * w[nt];
      }
      left_CS[ne] -= extrap_CS;
    }
    nce = nce_next;
    nct = nt2 - 1;
  } while (s_more != 0);
  // now nt2 will be index of last filled shell
  // Now to fill the shells which are absent in the input table.
  // They will be initialized phenomenologically, based on the sum rule.
  for (int ns = 0; ns < nt2; ns++) {
    m_acs[ns].reset(new PhenoPhotoAbsCS(shell_name[ns], Zshell[ns], thr[ns]));
  }
  // Initialization of input shells.
  for (int ns = nt2; ns < qshell; ns++) {
    auto adr = new SimpleTablePhotoAbsCS(shell_name[ns], Zshell[ns], thr[ns], ener, SCS[ns]);
    adr->remove_leading_zeros();
    m_acs[ns].reset(adr);
  }
  height_of_excitation = 0.0;
  exener[0] = exener[1] = 0.0;
  double integ = get_integral_ACS(0.0, DBL_MAX);
  // Iprintn(mcout, integ);
  integ_abs_before_corr = integ;
  double pred_integ = Thomas_sum_rule_const_Mb * Z;
  // Iprintn(mcout, pred_integ);
  if (pred_integ > integ) {
    if (s_add_excitations_to_normalize == 1) {
      const double threshold = m_acs[qshell - 1]->get_threshold();
      // add excitation
      exener[0] = low_boundary_of_excitations * threshold;
      exener[1] = threshold;
      height_of_excitation = (pred_integ - integ) / (exener[1] - exener[0]);
      if (minimal_threshold > 0.0) {
        if (minimal_threshold > threshold) {
          // currently the minimal shell is the last one
          exener[0] += minimal_threshold - threshold;
          exener[1] += minimal_threshold - threshold;
        }
      }
    }
  } else if (pred_integ < integ) {
    if (s_scale_to_normalize_if_more == 1) {
      const double fact = pred_integ / integ;
      for (int nshell = 0; nshell < qshell; ++nshell) {
        m_acs[nshell]->scale(fact);
      }
    }
  }
  integ_abs_after_corr = get_integral_ACS(0.0, DBL_MAX);
  integ_ioniz_after_corr = get_integral_ICS(0.0, DBL_MAX);
}

ExAtomPhotoAbsCS() [3/5]

Heed::ExAtomPhotoAbsCS::ExAtomPhotoAbsCS	(	int	fZ,
		const std::string &	fname,
		const std::string &	fBT_file_name,
		int	id,
		double	fminimal_threshold = 0.0
	)

Constructor, shells from Band and Thragzkovskaya.

Parameters

fZ	atomic number
fname	name of the atom
fBT_file_name	file with shell names and energies
id	1 - old files without fluorescence rate 2 - new files with fluorescence rate other values - error
fminimal_threshold	threshold

Definition at line 1269 of file PhotoAbsCS.cpp.

    : threshold_file_name("none"),
      simple_table_file_name("none"),
      BT_file_name(fBT_file_name),
      minimal_threshold(fminimal_threshold) {
  mfunnamep(
      "ExAtomPhotoAbsCS::ExAtomPhotoAbsCS(int fZ, const std::string& fname, "
      "const std::string& fBT_file_name, int id, double fminimal_threshold)");
  check_econd11(fZ, < 1, mcerr);
  check_econd21(id, < 1 ||, > 2, mcerr);
 
  name = fname;
  std::ifstream BT_file(BT_file_name.c_str());
  if (!BT_file) {
    funnw.ehdr(mcerr);
    mcerr << "cannot open file " << BT_file_name << std::endl;
    spexit(mcerr);
  }
  std::vector<double> thresh;
  std::vector<double> fl;
  Z = fZ;
  int i = findmark(BT_file, "NUCLEAR CHARGE =");
  check_econd11a(i, != 1, "wrong file format", mcerr);
  int Z_from_file;
  BT_file >> Z_from_file;
  check_econd12(Z_from_file, !=, Z, mcerr);
  qshell = 0;
  while ((i = findmark(BT_file, "Z =")) == 1) {
    BT_file >> i;
    check_econd11(i, != Z, mcerr);
    std::string shellname;
    BT_file >> shellname;
    // Iprintn(mcout, shellname);
    i = findmark(BT_file, "$");
    check_econd11(i, != 1, mcerr);
    long qen;
    BT_file >> qen;
    check_econd11(qen, <= 0, mcerr);
    std::vector<double> fener(qen, 0.0);
    std::vector<double> fcs(qen, 0.0);
    double thr = 0.0;
    BT_file >> thr;
    check_econd11(thr, <= 0, mcerr);
    thr *= 1.0e-3;  // pass from keV to MeV
    if (id == 2) {
      thresh.push_back(thr);
      fl.resize(fl.size() + 1);
      BT_file >> fl[qshell];
      check_econd21(fl[qshell], < 0.0 ||, > 1.0, mcerr);
      // Iprintn(mcout, fl[qshell]);
    }
    long nen;
    for (nen = 0; nen < qen; nen++) {
      BT_file >> fener[nen] >> fcs[nen];
      check_econd11(fener[nen], <= 0.0, mcerr);
      check_econd11(fcs[nen], < 0.0, mcerr);
      fener[nen] *= 1.0e-3;  // pass from keV to MeV
    }
    qshell++;
    m_acs.resize(qshell);
    m_acs.back().reset(new SimpleTablePhotoAbsCS(shellname, 0, thr, fener, fcs));
  }
  if (id == 2) {
    // a copy of similar thing from subroutine above
    int n_min = 0;
    double st = DBL_MAX;
    for (int nshell = 0; nshell < qshell; ++nshell) {
      // currently the minimal shell is the last,
      // but to avoid this assumption we check all
      if (thresh[nshell] < st) {
        n_min = nshell;
        st = thresh[nshell];
      }
    }
    asp.resize(qshell);
    for (int nshell = 0; nshell < qshell; ++nshell) {
      if (fl[nshell] > 0) {
        check_econd12(nshell, ==, n_min, mcerr);
        std::vector<double> felectron_energy;
        std::vector<double> fphoton_energy;
        fphoton_energy.push_back(thresh[nshell] - thresh[n_min]);
        asp[nshell].add_channel(fl[nshell], felectron_energy, fphoton_energy);
      }
    }
  }
 
  check_econd11(qshell, <= 0, mcerr);
  s_ignore_shell.resize(qshell, false);
  height_of_excitation = 0.0;
  exener[0] = exener[1] = 0.0;
  double integ = get_integral_ACS(0.0, DBL_MAX);
  // Iprintn(mcout, integ);
  integ_abs_before_corr = integ;
  double pred_integ = Thomas_sum_rule_const_Mb * Z;
  // Iprintn(mcout, pred_integ);
  if (pred_integ > integ) {
    if (s_add_excitations_to_normalize == 1) {
      const double thr = m_acs[qshell - 1]->get_threshold();
      // add excitation
      exener[0] = low_boundary_of_excitations * thr;
      exener[1] = thr;
      height_of_excitation = (pred_integ - integ) / (exener[1] - exener[0]);
      if (minimal_threshold > 0.0) {
        if (minimal_threshold > thr) {
          // currently the minimal shell is the last one
          exener[0] += minimal_threshold - thr;
          exener[1] += minimal_threshold - thr;
        }
      }
    }
  } else {
    if (s_scale_to_normalize_if_more == 1) {
      const double fact = pred_integ / integ;
      for (int nshell = 0; nshell < qshell; ++nshell) {
        m_acs[nshell]->scale(fact);
      }
    }
  }
  integ_abs_after_corr = get_integral_ACS(0.0, DBL_MAX);
  integ_ioniz_after_corr = get_integral_ICS(0.0, DBL_MAX);
}

ExAtomPhotoAbsCS() [4/5]

Heed::ExAtomPhotoAbsCS::ExAtomPhotoAbsCS	(	int	fZ,
		const std::string &	fname,
		const std::string &	fFitBT_file_name,
		int	id,
		int	s_no_scale,
		double	fminimal_threshold = 0.0
	)

Constructor, shells and fit parameters from Band and Thragzkovskaya.

Parameters

fZ	atomic number
fname	name of the atom
fFitBT_file_name	file with shell names, energies, and fit parameters
id	1 - old files without fluorescence rate 2 - new files with fluorescence rate other values - error
s_no_scale	scaling is not done, needs for next (?)
fminimal_threshold	threshold

Definition at line 1395 of file PhotoAbsCS.cpp.

    : threshold_file_name("none"),
      simple_table_file_name("none"),
      BT_file_name(fFitBT_file_name),
      minimal_threshold(fminimal_threshold) {
  mfunnamep(
      "ExAtomPhotoAbsCS::ExAtomPhotoAbsCS(int fZ, const std::string& fname, "
      "const "
      "std::string& fFitBT_file_name, int id, int id1, double "
      "fminimal_threshold)");
  check_econd11(fZ, < 1, mcerr);
  check_econd21(id, < 1 ||, > 2, mcerr);
  Z = fZ;
  name = fname;
  std::ifstream BT_file(fFitBT_file_name.c_str());
  if (!BT_file) {
    funnw.ehdr(mcerr);
    mcerr << "cannot open file " << BT_file_name << std::endl;
    spexit(mcerr);
  }
  std::vector<double> thresh;
  std::vector<double> fl;
  int i = 0;
  while ((i = findmark(BT_file, "$")) == 1) {
    long iZ;
    BT_file >> iZ;
    if (iZ != Z) continue;
    BT_file >> qshell;
    // Iprintn(mcout, qshell);
    check_econd11(qshell, <= 0, mcerr);
    check_econd11(qshell, > 1000, mcerr);
    m_acs.resize(qshell);
    if (id == 2) {
      thresh.resize(qshell);
      fl.resize(qshell);
    }
    for (int nshell = 0; nshell < qshell; ++nshell) {
#ifdef READ_FILE_WITH_PRINCIPAL_NUMBERS
      int n_princ = 0;
#endif
      int l;
      double threshold;
      double E0;
      double yw;
      double ya;
      double P;
      double sigma;
      if (BT_file.eof()) {
        mcerr << "unexpected end of file " << BT_file_name << '\n';
        spexit(mcerr);
      }
      if (!BT_file.good()) {
        mcerr << "bad format of file " << BT_file_name << '\n';
        spexit(mcerr);
      }
#ifdef READ_FILE_WITH_PRINCIPAL_NUMBERS
      BT_file >> n_princ;
      if (!BT_file.good()) {
        mcerr << "bad format of file " << BT_file_name << '\n';
        spexit(mcerr);
      }
      check_econd21(n_princ, < 0 ||, > 10, mcerr);
#endif
      BT_file >> l >> threshold >> E0 >> sigma >> ya >> P >> yw;
      check_econd11(l, < 0, mcerr);
      check_econd11(l, > 20, mcerr);
      threshold *= 1.0e-6;
      E0 *= 1.0e-6;
 
      check_econd11a(threshold, <= 2.0e-6,
                     "n_princ=" << n_princ << " l=" << l << '\n', mcerr);
      check_econd11(E0, <= 0, mcerr);
      double flu = 0.0;
      if (id == 2) {
        if (BT_file.eof()) {
          mcerr << "unexpected end of file " << BT_file_name << '\n';
          spexit(mcerr);
        }
        if (!BT_file.good()) {
          mcerr << "bad format of file " << BT_file_name << '\n';
          spexit(mcerr);
        }
        BT_file >> flu;
        check_econd11(flu, < 0.0, mcerr);
        check_econd11(flu, > 1.0, mcerr);
        thresh[nshell] = threshold;
        fl[nshell] = flu;
      }
#ifdef READ_FILE_WITH_PRINCIPAL_NUMBERS
      // necessary for generation escape products
      std::string shellname(std::to_string(n_princ) + " shell number " +
                            std::to_string(nshell));
#else
      std::string shellname("shell number " + std::to_string(nshell));
#endif
      m_acs[nshell].reset(new SimpleTablePhotoAbsCS(shellname, 0, threshold, 
            l, E0, yw, ya, P, sigma));
      // Iprintn(mcout, nshell);
      // Iprint3n(mcout, l, threshold, E0);
      // Iprint4n(mcout, yw, ya, P, sigma);
      // acs[nshell]->print(mcout, 5);
    }
    goto mark1;
  }
  funnw.ehdr(mcerr);
  mcerr << "there is no element Z=" << fZ << " in file " << fFitBT_file_name
        << std::endl;
  spexit(mcerr);
mark1:
  if (id == 2) {
    // a copy of similar thing from subroutine above
    int n_min = 0;
    double st = DBL_MAX;
    for (int nshell = 0; nshell < qshell; ++nshell) {
      // currently the minimal shell is the last,
      // but to avoid this assumption we check all
      if (thresh[nshell] < st) {
        n_min = nshell;
        st = thresh[nshell];
      }
    }
    asp.resize(qshell);
    for (int nshell = 0; nshell < qshell; ++nshell) {
      if (fl[nshell] > 0) {
        check_econd12(nshell, ==, n_min, mcerr);
        std::vector<double> felectron_energy;
        std::vector<double> fphoton_energy;
        fphoton_energy.push_back(thresh[nshell] - thresh[n_min]);
        asp[nshell].add_channel(fl[nshell], felectron_energy, fphoton_energy);
      }
    }
  }
 
  check_econd11(qshell, <= 0, mcerr);
  s_ignore_shell.resize(qshell, false);
  height_of_excitation = 0.0;
  exener[0] = exener[1] = 0.0;
  double integ = get_integral_ACS(0.0, DBL_MAX);
  // Iprintn(mcout, integ);
  integ_abs_before_corr = integ;
  double pred_integ = Thomas_sum_rule_const_Mb * Z;
  // Iprintn(mcout, pred_integ);
  if (pred_integ > integ) {
    if (s_add_excitations_to_normalize == 1) {
      const double thr = m_acs[qshell - 1]->get_threshold();
      // add excitation
      exener[0] = low_boundary_of_excitations * thr;
      exener[1] = thr;
      height_of_excitation = (pred_integ - integ) / (exener[1] - exener[0]);
      if (minimal_threshold > 0.0) {
        if (minimal_threshold > thr) {
          // currently the minimal shell is the last one
          exener[0] += minimal_threshold - thr;
          exener[1] += minimal_threshold - thr;
        }
      }
    }
  } else {
    if (s_scale_to_normalize_if_more == 1 && s_no_scale == 0) {
      const double fact = pred_integ / integ;
      for (int nshell = 0; nshell < qshell; ++nshell) {
        m_acs[nshell]->scale(fact);
      }
    }
  }
  integ_abs_after_corr = get_integral_ACS(0.0, DBL_MAX);
  integ_ioniz_after_corr = get_integral_ICS(0.0, DBL_MAX);
}

ExAtomPhotoAbsCS() [5/5]

Heed::ExAtomPhotoAbsCS::ExAtomPhotoAbsCS	(	int	fZ,
		const std::string &	fname,
		const std::string &	fFitBT_file_name,
		const std::string &	fsimple_table_file_name,
		double	emax_repl,
		int	id,
		double	fminimal_threshold = 0.0
	)

Constructor, combination of Band and Thragzkovskaya fit and Henke tables. Initialize BT fit and replaces the part of the first shell from threshold taken from BT- fit to emax_repl by values from the table.

Parameters

fZ	atomic number
fname	name of the atom
fFitBT_file_name	file with shell names, energies, and fit parameters
fsimple_table_file_name	file with cross-section table
emax_repl	energy up to which to use the cross-section table
id	1 - old files without fluorescence rate 2 - new files with fluorescence rate other values - error
fminimal_threshold	threshold

Definition at line 1567 of file PhotoAbsCS.cpp.

                               {
  mfunname("ExAtomPhotoAbsCS::ExAtomPhotoAbsCS(...)");
  Z = fZ;
  name = fname;
  int s_no_scale = 1;
  *this = ExAtomPhotoAbsCS(fZ, fname, fFitBT_file_name, id, s_no_scale,
                           fminimal_threshold);
 
  height_of_excitation = 0.0;
  exener[0] = exener[1] = 0.0;
 
  double thrmin = DBL_MAX;
  long nsmin = -1;
  // Look for minimal shell (usually the last).
  for (long ns = 0; ns < qshell; ++ns) {
    if (thrmin > m_acs[ns]->get_threshold()) {
      nsmin = ns;
      thrmin = m_acs[ns]->get_threshold();
    }
  }
  check_econd11(nsmin, < 0, mcerr);
  check_econd11(nsmin, != qshell - 1, mcerr);
  
  PhotoAbsCS* apacs = m_acs[nsmin].get();
  auto first_shell = dynamic_cast<SimpleTablePhotoAbsCS*>(apacs);
  check_econd11(first_shell, == nullptr, mcerr);
 
  SimpleTablePhotoAbsCS stpacs(name, Z, 0.0, fsimple_table_file_name);
  stpacs.remove_leading_tiny(1.0e-10);
 
  // Merging shells:
  SimpleTablePhotoAbsCS* merged = new SimpleTablePhotoAbsCS(*first_shell, stpacs, emax_repl); 
  m_acs[nsmin].reset(merged);
 
  s_ignore_shell.resize(qshell, false);
  height_of_excitation = 0.0;
  exener[0] = exener[1] = 0.0;
  double integ = get_integral_ACS(0.0, DBL_MAX);
  integ_abs_before_corr = integ;
  double pred_integ = Thomas_sum_rule_const_Mb * Z;
  if (pred_integ > integ) {
    if (s_add_excitations_to_normalize == 1) {
      const double thr = m_acs[qshell - 1]->get_threshold();
      // add excitation
      exener[0] = low_boundary_of_excitations * thr;
      exener[1] = thr;
      height_of_excitation = (pred_integ - integ) / (exener[1] - exener[0]);
      if (minimal_threshold > 0.0) {
        if (minimal_threshold > thr) {
          // currently the minimal shell is the last one
          exener[0] += minimal_threshold - thr;
          exener[1] += minimal_threshold - thr;
        }
      }
    }
  } else {
    if (s_scale_to_normalize_if_more == 1) {
      const double fact = pred_integ / integ;
      for (int nshell = 0; nshell < qshell; ++nshell) {
        m_acs[nshell]->scale(fact);
      }
    }
  }
  integ_abs_after_corr = get_integral_ACS(0.0, DBL_MAX);
  integ_ioniz_after_corr = get_integral_ICS(0.0, DBL_MAX);
}

~ExAtomPhotoAbsCS()

virtual Heed::ExAtomPhotoAbsCS::~ExAtomPhotoAbsCS ( )

inlinevirtual

Destructor.

Definition at line 538 of file PhotoAbsCS.h.

{}

Member Function Documentation

copy()

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

inlinevirtual

Implements Heed::AtomPhotoAbsCS.

Definition at line 452 of file PhotoAbsCS.h.

{ return new ExAtomPhotoAbsCS(*this); }

get_ACS() [1/2]

double Heed::ExAtomPhotoAbsCS::get_ACS ( double  energy ) const

virtual

Photo-absorption cross-section [Mbarn] at a given energy [MeV]. The photo-absorption cross-section can include excitation.

Implements Heed::AtomPhotoAbsCS.

Definition at line 1704 of file PhotoAbsCS.cpp.

                                                    {
  mfunname("double ExAtomPhotoAbsCS::get_ACS(double energy) const");
  double s = 0.0;
  for (int n = 0; n < qshell; ++n) {
    if (s_ignore_shell[n]) continue;
    double shift = 0.0;
    const double t = m_acs[n]->get_threshold();
    if (minimal_threshold > 0.0) {
      if (t < minimal_threshold) shift = minimal_threshold - t;
    }
    s += m_acs[n]->get_CS(energy - shift);
  }
  if (energy >= exener[0] && energy <= exener[1]) s += height_of_excitation;
  return s;
}

get_ACS() [2/2]

double Heed::ExAtomPhotoAbsCS::get_ACS ( int  nshell, double  energy  ) const

virtual

Sub-shell photo-absorption cross-section [Mbarn] at a given energy [MeV].

Implements Heed::AtomPhotoAbsCS.

Definition at line 1738 of file PhotoAbsCS.cpp.

                                                                {
  mfunname("double ExAtomPhotoAbsCS::get_ACS(int nshell, double energy)");
  check_econd21(nshell, < 0 ||, > qshell, mcerr);
  if (s_ignore_shell[nshell]) return 0.;
  double shift = 0.0;
  const double t = m_acs[nshell]->get_threshold();
  if (minimal_threshold > 0.0) {
    if (t < minimal_threshold) shift = minimal_threshold - t;
  }
  double s = m_acs[nshell]->get_CS(energy - shift);
  if (nshell == qshell - 1 && energy >= exener[0] && energy <= exener[1]) {
    s += height_of_excitation;
  }
  return s;
}

get_ICS() [1/2]

double Heed::ExAtomPhotoAbsCS::get_ICS ( double  energy ) const

virtual

Photo-ionization cross-section [Mbarn] at a given energy [MeV]. The photo-ionization cross-section does not include excitation.

Implements Heed::AtomPhotoAbsCS.

Definition at line 1648 of file PhotoAbsCS.cpp.

                                                    {
  mfunname("double ExAtomPhotoAbsCS::get_ACS(double energy) const");
  double s = 0.0;
  for (int n = 0; n < qshell; ++n) {
    if (s_ignore_shell[n]) continue;
    double shift = 0.0;
    const double t = m_acs[n]->get_threshold();
    if (minimal_threshold > 0.0) {
      if (t < minimal_threshold) shift = minimal_threshold - t;
    }
    s += m_acs[n]->get_CS(energy - shift);
  }
  return s;
}

get_ICS() [2/2]

double Heed::ExAtomPhotoAbsCS::get_ICS ( int  nshell, double  energy  ) const

virtual

Sub-shell photo-ionization cross-section at a given energy.

Implements Heed::AtomPhotoAbsCS.

Definition at line 1679 of file PhotoAbsCS.cpp.

                                                                {
  mfunname("double ExAtomPhotoAbsCS::get_ICS(int nshell, double energy) const");
  check_econd21(nshell, < 0 ||, > qshell, mcerr);
  if (s_ignore_shell[nshell]) return 0.;
  double shift = 0.0;
  const double t = m_acs[nshell]->get_threshold();
  if (minimal_threshold > 0.0) {
    if (t < minimal_threshold) shift = minimal_threshold - t;
  }
  return m_acs[nshell]->get_CS(energy - shift);
}

get_integral_ACS() [1/2]

double Heed::ExAtomPhotoAbsCS::get_integral_ACS ( double  energy1, double  energy2  ) const

virtual

Integrated photo-absorption cross-section overa given interval.

Implements Heed::AtomPhotoAbsCS.

Definition at line 1720 of file PhotoAbsCS.cpp.

                                                                {
  mfunname("double ExAtomPhotoAbsCS::get_integral_ACS(...) const");
  double s = 0.0;
  for (int n = 0; n < qshell; ++n) {
    if (s_ignore_shell[n]) continue;
    double shift = 0.0;
    const double t = m_acs[n]->get_threshold();
    if (minimal_threshold > 0.0) {
      if (t < minimal_threshold) shift = minimal_threshold - t;
    }
    s += m_acs[n]->get_integral_CS(energy1 - shift, energy2 - shift);
  }
  double b[2] = {std::max(exener[0], energy1), std::min(exener[1], energy2)};
  if (b[1] >= b[0]) s += height_of_excitation * (b[1] - b[0]);
  return s;
}

Referenced by ExAtomPhotoAbsCS(), and print().

get_integral_ACS() [2/2]

double Heed::ExAtomPhotoAbsCS::get_integral_ACS ( int  nshell, double  energy1, double  energy2  ) const

virtual

Integrated sub-shell photo-absorption cross-section.

Implements Heed::AtomPhotoAbsCS.

Definition at line 1754 of file PhotoAbsCS.cpp.

                                                                {
  mfunname("double ExAtomPhotoAbsCS::get_integral_ACS(int nshell, ...) const");
  check_econd21(nshell, < 0 ||, > qshell, mcerr);
  if (s_ignore_shell[nshell]) return 0.;
  double shift = 0.0;
  const double t = m_acs[nshell]->get_threshold();
  if (minimal_threshold > 0.0) {
    if (t < minimal_threshold) shift = minimal_threshold - t;
  }
  double s = m_acs[nshell]->get_integral_CS(energy1 - shift, energy2 - shift);
  if (nshell == qshell - 1) {
    double b[2] = {std::max(exener[0], energy1), std::min(exener[1], energy2)};
    if (b[1] >= b[0]) s += height_of_excitation * (b[1] - b[0]);
  }
  return s;
}

get_integral_ICS() [1/2]

double Heed::ExAtomPhotoAbsCS::get_integral_ICS ( double  energy1, double  energy2  ) const

virtual

Integrated photo-ionization cross-section over a given interval.

Implements Heed::AtomPhotoAbsCS.

Definition at line 1663 of file PhotoAbsCS.cpp.

                                                                {
  mfunname("double ExAtomPhotoAbsCS::get_integral_ICS(double energy) const");
  double s = 0.0;
  for (int n = 0; n < qshell; ++n) {
    if (s_ignore_shell[n]) continue;
    double shift = 0.0;
    const double t = m_acs[n]->get_threshold();
    if (minimal_threshold > 0.0) {
      if (t < minimal_threshold) shift = minimal_threshold - t;
    }
    s += m_acs[n]->get_integral_CS(energy1 - shift, energy2 - shift);
  }
  return s;
}

Referenced by ExAtomPhotoAbsCS(), and print().

get_integral_ICS() [2/2]

double Heed::ExAtomPhotoAbsCS::get_integral_ICS ( int  nshell, double  energy1, double  energy2  ) const

virtual

Integrated sub-shell photo-ionization cross-section.

Implements Heed::AtomPhotoAbsCS.

Definition at line 1691 of file PhotoAbsCS.cpp.

                                                                {
  mfunname("double ExAtomPhotoAbsCS::get_integral_ICS(int nshell, ...) const");
  check_econd21(nshell, < 0 ||, > qshell, mcerr);
  if (s_ignore_shell[nshell]) return 0.;
  double shift = 0.0;
  const double t = m_acs[nshell]->get_threshold();
  if (minimal_threshold > 0.0) {
    if (t < minimal_threshold) shift = minimal_threshold - t;
  }
  return m_acs[nshell]->get_integral_CS(energy1 - shift, energy2 - shift);
}

get_main_shell_number()

virtual int Heed::ExAtomPhotoAbsCS::get_main_shell_number ( int  nshell ) const

inlinevirtual

Return the shell number (1, 2, ...) for a given index. The number is taken from the shell name. If the shell number cannot be determined, the function returns -1.

Implements Heed::AtomPhotoAbsCS.

Definition at line 445 of file PhotoAbsCS.h.

                                                      {
    return m_acs[nshell]->get_number();
  }

get_threshold()

double Heed::ExAtomPhotoAbsCS::get_threshold ( int  nshell ) const

virtual

Get the ionization threshold for a given shell.

Implements Heed::AtomPhotoAbsCS.

Definition at line 1638 of file PhotoAbsCS.cpp.

                                                       {
  mfunname("double ExAtomPhotoAbsCS::get_threshold(int nshell) const");
  check_econd21(nshell, < 0 ||, > qshell, mcerr);
  double r = m_acs[nshell]->get_threshold();
  if (minimal_threshold > 0.0) {
    if (r < minimal_threshold) r = minimal_threshold;
  }
  return r;
}

Referenced by ExAtomPhotoAbsCS().

print()

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

virtual

Reimplemented from Heed::AtomPhotoAbsCS.

Definition at line 1772 of file PhotoAbsCS.cpp.

                                                          {
  if (l <= 0) return;
  Ifile << "ExAtomPhotoAbsCS(l=" << l << "): name=" << name << " Z = " << Z
        << " qshell = " << qshell << std::endl;
  indn.n += 2;
  Ifile << "threshold_file_name=" << threshold_file_name << '\n';
  Ifile << "simple_table_file_name=" << simple_table_file_name << '\n';
  Ifile << "BT_file_name=" << BT_file_name << std::endl;
  Ifile << "Thomas_sum_rule_const_Mb * Z = " << Thomas_sum_rule_const_Mb* Z
        << '\n';
  Ifile << "integ_abs_before_corr        = " << integ_abs_before_corr << '\n';
  Ifile << "integ_abs_after_corr         = " << integ_abs_after_corr << '\n';
  Ifile << "integ_ioniz_after_corr       = " << integ_ioniz_after_corr << '\n';
  Ifile << "height_of_excitation=" << height_of_excitation
        << " exener=" << exener[0] << ' ' << exener[1] << '\n';
  Iprintn(file, minimal_threshold);
  Ifile << "integrals by shells:\n";
  Ifile << "nshell, int(acs), int(ics)\n";
  for (long n = 0; n < qshell; n++) {
    double ainteg = get_integral_ACS(n, 0.0, DBL_MAX);
    double iinteg = get_integral_ICS(n, 0.0, DBL_MAX);
    Ifile << n << "    " << ainteg << "    " << iinteg << '\n';
  }
 
  if (l > 1) {
    l--;
    indn.n += 2;
    for (long n = 0; n < qshell; ++n) {
      Ifile << "nshell=" << n << std::endl;
      m_acs[n]->print(file, l);
    }
    AtomPhotoAbsCS::print(file, l);
    indn.n -= 2;
  }
  indn.n -= 2;
}

replace_shells_by_average()

void Heed::ExAtomPhotoAbsCS::replace_shells_by_average	(	double	fwidth,
		double	fstep,
		long	fmax_q_step
	)

Definition at line 1809 of file PhotoAbsCS.cpp.

                                                                   {
  mfunname("void ExAtomPhotoAbsCS::replace_shells_by_average(...)");
  for (long n = 0; n < qshell; n++) {
    if (!m_acs[n]) continue;
    PhotoAbsCS* a =
        new AveragePhotoAbsCS(m_acs[n].get(), fwidth, fstep, fmax_q_step);
    m_acs[n].reset(a);
  }
}

Member Data Documentation

BT_file_name

std::string Heed::ExAtomPhotoAbsCS::BT_file_name

protected

Definition at line 543 of file PhotoAbsCS.h.

Referenced by ExAtomPhotoAbsCS(), and print().

exener

double Heed::ExAtomPhotoAbsCS::exener[2] = {0., 0.}

protected

Boundaries of excitation.

Definition at line 555 of file PhotoAbsCS.h.

Referenced by ExAtomPhotoAbsCS(), get_ACS(), get_integral_ACS(), and print().

height_of_excitation

double Heed::ExAtomPhotoAbsCS::height_of_excitation = 0.

protected

Excitation cross-section (assumed in the lowest shell).

Definition at line 553 of file PhotoAbsCS.h.

Referenced by ExAtomPhotoAbsCS(), get_ACS(), get_integral_ACS(), and print().

integ_abs_after_corr

double Heed::ExAtomPhotoAbsCS::integ_abs_after_corr = 0.

protected

Definition at line 550 of file PhotoAbsCS.h.

Referenced by ExAtomPhotoAbsCS(), and print().

integ_abs_before_corr

double Heed::ExAtomPhotoAbsCS::integ_abs_before_corr = 0.

protected

Definition at line 549 of file PhotoAbsCS.h.

Referenced by ExAtomPhotoAbsCS(), and print().

integ_ioniz_after_corr

double Heed::ExAtomPhotoAbsCS::integ_ioniz_after_corr = 0.

protected

Definition at line 551 of file PhotoAbsCS.h.

Referenced by ExAtomPhotoAbsCS(), and print().

m_acs

std::vector<std::shared_ptr<PhotoAbsCS> > Heed::ExAtomPhotoAbsCS::m_acs

protected

Ionization cross-section (the name acs is misleading). Excitations are added separately as height_of_excitation.

Definition at line 546 of file PhotoAbsCS.h.

Referenced by ExAtomPhotoAbsCS(), get_ACS(), get_ICS(), get_integral_ACS(), get_integral_ICS(), get_main_shell_number(), get_threshold(), print(), and replace_shells_by_average().

minimal_threshold

double Heed::ExAtomPhotoAbsCS::minimal_threshold = 0.

protected

Definition at line 557 of file PhotoAbsCS.h.

Referenced by ExAtomPhotoAbsCS(), get_ACS(), get_ICS(), get_integral_ACS(), get_integral_ICS(), get_threshold(), and print().

s_add_excitations_to_normalize

const int Heed::ExAtomPhotoAbsCS::s_add_excitations_to_normalize = 1

staticprotected

Flag whether to add excitations. If 0 excitations will not be added (useful for debugging and for checking the effect produced by adding excitations). For real work, this variable should always be set to 1.

Definition at line 569 of file PhotoAbsCS.h.

Referenced by ExAtomPhotoAbsCS().

s_scale_to_normalize_if_more

const int Heed::ExAtomPhotoAbsCS::s_scale_to_normalize_if_more = 1

staticprotected

Definition at line 570 of file PhotoAbsCS.h.

Referenced by ExAtomPhotoAbsCS().

simple_table_file_name

std::string Heed::ExAtomPhotoAbsCS::simple_table_file_name

protected

Definition at line 542 of file PhotoAbsCS.h.

Referenced by print().

threshold_file_name

std::string Heed::ExAtomPhotoAbsCS::threshold_file_name

protected

Definition at line 541 of file PhotoAbsCS.h.

Referenced by ExAtomPhotoAbsCS(), and print().

---

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

• PhotoAbsCS.h
• PhotoAbsCS.cpp