3. The libraries of nrCascadeSim

The functionality of the executables of nrCascadeSim such as realizeCascades is constructed from a set of library functions whose code is stored in the libncap.so file and whose prototypes are stored in several header files. Below we list the different header files with the internal functions, data structures, and their uses.

cascadeProd.h

This header has the prototypes for the functions that generate the statistical realizations of each cascade by Monte-Carlo simulation. The structures and functions related to our data structures are prototyped below:

struct cli

A structure to hold input information on each cascade. the “Cascade Level Info” struct

Public Members

bool success: Has the structure been set successfully and written correctly?

int n: number of levels in this cascade

int cid: cascade id

double frac: probability of observing this particular cascade

double Sn: neutron separation energy in MeV

int A: mass number of capturing isotope

vector<double> Elev: array of energy levels in keV

vector<double> taus: half-lives of energy level in fs

struct cri

A structure to hold output information on each cascade. the “Cascade Recoil Info” struct

Public Members

int n: number of levels in the cascade that gave this recoil

int cid: cascade id for the cascade that gave this recoil

vector<double> E: recoil energy at beginning of recoil step eV

vector<double> delE: energy deposited in step eV

vector<int> I: ionization created in each step in n e/h pairs

vector<double> Ei: electron-equivalent ionization energy eV

vector<double> time: time since capture fs

vector<double> Eg: gamma energy emitted in each step MeV

There are also some utility functions that are used for reading the cascade input files, and they are prototyped below.

vector<cli> readCascadeDistributionFile(int &n, string file, bool &success): function to read in the cascade file with n lines. n is the number of cascade rows; file is the full path to the file; and success is a flag that indicates successful reading, true for a success.

double interpretDbl(string in, bool &success): function to interpret doubles from input using regex. in is the input string to interpret; and success is true if it can be interpreted as a double-precision number (either in decimal or scientific notation). A cli object is returned.

double interpretSn(string in, bool &success): function to read Sn from input correctly using regex. in is the input string to interpret; success will be true if it can be read as a double-precision number (representing the neutron separation energy) OR it is a recognized isotopic symbol like 74Ge. A double representing the string is returned.

double interpretWeisskopf(string in, double Egam, double A, bool &success): function to convert Weisskopf abbreviations using regex. in is the input string to interpret; success will be set to true if a Weisskopf multipolarity abbreviation is recognized, like w(E1) for an electric dipole. Egam is the energy of the emitted gamma; A is the mass number. A double giving the transition lifetime in fs is returned.

vector<double> interpretElevVector(int &n, string in, bool &success): function to read E levels from input correctly. n is the number of elements in the vector to be interpreted from the input; in is the input string to interpret; success will be set to true if a vector of doubles in the correct format (surrounded by brackets [] and separated by whitespace) is detected.

vector<double> interpretTauVector(int n, string in, double A, vector<double> &Elev, bool &success): function to read lifetimes from input correctly. n is the number of elements in the vector to be interpreted from the input; in is the input string to interpret; success will be set to true if a vector of doubles or symbolic Weisskopof indicators in the correct format (surrounded by brackets and separated by whitespace) is detected.

vector<string> vsplit(string in): function for splitting strings. in is the input string to interpret. A vector of strings is returned represented the whitespace-split version of the original.

The functions below provide the functionality to calculate various details of the atom/ion trajectories for the supported elements: germanium, silicon, argon, neon. At this time there are separate functions for each of the supported elements; this is meant to be unified in the future in order to support a wider range of elements. For now we always use constant-acceleration S2 stopping. S2 refers to the parameter from the Lindhard paper [Lindhard1963].

vector<cri> Cascade(int n, int cid, double Sn, int nlev, vector<double> &Elev, vector<double> &taus, double A, mt19937 *mtrand): Generic function for realizing n cascades, WARNING currently only works with Si, and Ge cid is the cascade id; Sn is the neutron separation energy in MeV; nlev is the number of levels of the cascade; Elev is an array of level energys of size nlev in keV; taus is an array of level lifetimes of size nlev in fs; A is the mass number; mtrand is the random number generator object—in this case std::mt19937. The function returns a cri struct.

vector<cri> geCascade(int n, int cid, double Sn, int nlev, vector<double> &Elev, vector<double> &taus, double A, mt19937 *mtrand): A function to simulate n cascade realizations for a germanium isotope. cid is the cascade id; Sn is the neutron separation energy in MeV; nlev is the number of levels of the cascade; Elev is an array of level energys of size nlev in keV; taus is an array of level lifetimes of size nlev in fs; A is the mass number; mtrand is the random number generator object—in this case std::mt19937. The function returns a cri struct.

double geDecay(double v, double M, double Egam, mt19937 *rand): germanium func. to return the Energy after the mid-stop decay. v is the velocity in units of the speed of light (c); M is the mass of the slowing object in GeV, Egam is the energy of the emitted gamma ray in MeV; and mtrand is the random number generator object—in this case std::mt19937. The function returns the energy in eV after the mid-step decay kick.

vector<double> geStop(double E, double M, double tau, mt19937 *rand): germanium func. to return the velocity at a random stopping time. E is the initial atom energy in eV; M is the mass of the slowing object in GeV; tau is the (random) decay time of the current intermediate nuclear level; mtrand is the random generator object—in this case std::mt19937. The function returns the velocity at the random decay time in units of the speed of light (c).

double rgeS2(double E, double M, double t): germanium func. constant stopping-power stopping. E is the initial energy of the stopping atom in eV; M is the mass of the slowing object in GeV; t is the time (in fs) at which to return the distance traveled. Function returns distance in nm

double vgeS2(double E, double M, double t): germanium func. constant stopping-power stopping. E is the initial energy of the stopping atom in eV; M is the mass of the slowing object in GeV; t is the time (in fs) at which to return the atoms velocity. Function returns velocity in units of the speed of light (c).

double vgeS2func(vector<double> &x, vector<double> &par): germanium func. constant stopping-power stopping. x[0] is the time (in fs) at which to return the atoms velocity; par[0] is the atoms initial energy in eV; `par[1] is the atoms mass in GeV. Function returns velocity in units of the speed of light (c).

vector<cri> siCascade(int n, int cid, double Sn, int nlev, vector<double> &Elev, vector<double> &taus, double A, mt19937 *mtrand): A function to simulate n cascade realizations for a silicon isotope. cid is the cascade id; Sn is the neutron separation energy in MeV; nlev is the number of levels of the cascade; Elev is an array of level energys of size nlev in keV; taus is an array of level lifetimes of size nlev in fs; A is the mass number; mtrand is the random number generator object—in this case std::mt19937. The function returns a cri struct.

double siDecay(double v, double M, double Egam, mt19937 *rand): silicon func. to return the Energy after the mid-stop decay. v is the velocity in units of the speed of light (c); M is the mass of the slowing object in GeV, Egam is the energy of the emitted gamma ray in MeV; and mtrand is the random number generator object—in this case std::mt19937. The function returns the energy in eV after the mid-step decay kick.

vector<double> siStop(double E, double M, double tau, mt19937 *rand): silicon func. to return the velocity at a random stopping time. E is the initial atom energy in eV; M is the mass of the slowing object in GeV; tau is the (random) decay time of the current intermediate nuclear level; mtrand is the random generator object—in this case std::mt19937. The function returns the velocity at the random decay time in units of the speed of light (c).

double rsiS2(double E, double M, double t): silicon func. constant stopping-power stopping. E is the initial energy of the stopping atom in eV; M is the mass of the slowing object in GeV; t is the time (in fs) at which to return the distance traveled. Function returns distance in nm

double vsiS2(double E, double M, double t): silicon func. constant stopping-power stopping. E is the initial energy of the stopping atom in eV; M is the mass of the slowing object in GeV; t is the time (in fs) at which to return the atoms velocity. Function returns velocity in units of the speed of light (c).

double vsiS2func(vector<double> &x, vector<double> &par): silicon func. constant stopping-power stopping. x[0] is the time (in fs) at which to return the atoms velocity; par[0] is the atoms initial energy in eV; `par[1] is the atoms mass in GeV. Function returns velocity in units of the speed of light (c).

vector<cri> arCascade(int n, int cid, double Sn, int nlev, vector<double> &Elev, vector<double> &taus, double A, mt19937 *mtrand): A function to simulate n cascade realizations for a argon isotope. cid is the cascade id; Sn is the neutron separation energy in MeV; nlev is the number of levels of the cascade; Elev is an array of level energys of size nlev in keV; taus is an array of level lifetimes of size nlev in fs; A is the mass number; mtrand is the random number generator object—in this case std::mt19937. The function returns a cri struct.

double arDecay(double v, double M, double Egam, mt19937 *rand): argon func. to return the Energy after the mid-stop decay. v is the velocity in units of the speed of light (c); M is the mass of the slowing object in GeV, Egam is the energy of the emitted gamma ray in MeV; and mtrand is the random number generator object—in this case std::mt19937. The function returns the energy in eV after the mid-step decay kick.

vector<double> arStop(double E, double M, double tau, mt19937 *rand): argon func. to return the velocity at a random stopping time. E is the initial atom energy in eV; M is the mass of the slowing object in GeV; tau is the (random) decay time of the current intermediate nuclear level; mtrand is the random generator object—in this case std::mt19937. The function returns the velocity at the random decay time in units of the speed of light (c).

double rarS2(double E, double M, double t): argon func. constant stopping-power stopping. E is the initial energy of the stopping atom in eV; M is the mass of the slowing object in GeV; t is the time (in fs) at which to return the distance traveled. Function returns distance in nm

double varS2(double E, double M, double t): argon func. constant stopping-power stopping. E is the initial energy of the stopping atom in eV; M is the mass of the slowing object in GeV; t is the time (in fs) at which to return the atoms velocity. Function returns velocity in units of the speed of light (c).

double varS2func(vector<double> &x, vector<double> &par): argon func. constant stopping-power stopping. x[0] is the time (in fs) at which to return the atoms velocity; par[0] is the atoms initial energy in eV; `par[1] is the atoms mass in GeV. Function returns velocity in units of the speed of light (c).

vector<cri> neCascade(int n, int cid, double Sn, int nlev, vector<double> &Elev, vector<double> &taus, double A, mt19937 *mtrand): A function to simulate n cascade realizations for a neon isotope. cid is the cascade id; Sn is the neutron separation energy in MeV; nlev is the number of levels of the cascade; Elev is an array of level energys of size nlev in keV; taus is an array of level lifetimes of size nlev in fs; A is the mass number; mtrand is the random number generator object—in this case std::mt19937. The function returns a cri struct.

double neDecay(double v, double M, double Egam, mt19937 *rand): neon func. to return the Energy after the mid-stop decay. v is the velocity in units of the speed of light (c); M is the mass of the slowing object in GeV, Egam is the energy of the emitted gamma ray in MeV; and mtrand is the random number generator object—in this case std::mt19937. The function returns the energy in eV after the mid-step decay kick.

vector<double> neStop(double E, double M, double tau, mt19937 *rand): neon func. to return the velocity at a random stopping time. E is the initial atom energy in eV; M is the mass of the slowing object in GeV; tau is the (random) decay time of the current intermediate nuclear level; mtrand is the random generator object—in this case std::mt19937. The function returns the velocity at the random decay time in units of the speed of light (c).

double rneS2(double E, double M, double t): neon func. constant stopping-power stopping. E is the initial energy of the stopping atom in eV; M is the mass of the slowing object in GeV; t is the time (in fs) at which to return the distance traveled. Function returns distance in nm

double vneS2(double E, double M, double t): neon func. constant stopping-power stopping. E is the initial energy of the stopping atom in eV; M is the mass of the slowing object in GeV; t is the time (in fs) at which to return the atoms velocity. Function returns velocity in units of the speed of light (c).

double vneS2func(vector<double> &x, vector<double> &par): neon func. constant stopping-power stopping. x[0] is the time (in fs) at which to return the atoms velocity; par[0] is the atoms initial energy in eV; `par[1] is the atoms mass in GeV. Function returns velocity in units of the speed of light (c).

edepmath.h

In this header is contained prototypes for supporting mathematical functions. Most of the content are functions to assist with drawing numbers from specific probability distributions.

Functions

int poissonKnuth(double lambda, mt19937 *rand): Generate a Poisson-distributed integer random number with the Knuth method. lambda is the parameter of the Poisson distribution; rand is a random generator object—specifically the std::mt19937. The return value is a single integer realization of a Poisson random variable.

int poissonAtkinson(double lambda, mt19937 *rand): Generate a Poisson-distributed integer random number with the Atkinson method. DANGEROUS: NEVER USE for lambda<30. lambda is the parameter of the Poisson distribution; rand is a random generator object—specifically the std::mt19937. The return value is a single integer realization of a Poisson random variable.

int poisson(double lambda, mt19937 *rand): Generate a Poisson-distributed integer random number that chooses between the Knuth and Atkinson method based on performance. Uses the Atkinson function when lambda<30 and the Knuth method otherwise. lambda is the parameter of the Poisson distribution; rand is a random generator object—specifically the std::mt19937. The return value is a single integer realization of a Poisson random variable.

int poissonFano(double lambda, double F, mt19937 *rand): Generate a sudo-Poisson-distributed integer random number that sets the width based on the so-called “Fano Factor” which is an empirical factor developed in the 50’s for widening a sort-of Poisson distrubtion. A Fano Factor of 1 corresponds to the standard Poisson distrubtion. lambda is the parameter of the Poisson distribution; F is the dimensionless Fano Factor typically a number less than 1 (0.13 for germanium, 0.115 for silicon are nominal values) for Electron Recoils; rand is a random generator object—specifically the std::mt19937. The return value is a single integer realization of a Poisson random variable.

lindhard.h

In this header is contained prototypes for functions to furnish simple representations of the Lindhard ionization model [Lindhard1963]. They generally help return the ionization yield fraction given at a particular starting energy (in eV). There are also specified functions to return the ionization for an atom slowing down from one starting energy to another (as would happen in one step of the cascade). Again, as in cascadeProd.h there are separate functions for each isotope currently and this is intended to be unified in the future.

Functions

double lindhard(vector<double> &x, vector<double> &par): generic lindhard function (do not use). The input x is really an array with a single value x[0] which is the input recoil energy in eV. The array par needs exactly 6 parameters to construct the standard Lindhard ionization yield model. These are par[0] the pre-factor to the first energy term in the Lindhard model (unitless); par[1] the exponent of the first energy term; par[2] the pre-factor of thee second energy term; par[3] the exponent of the second energy term; par[4] the atomic number Z; par[5] the mass number A. The return value is the ionization yield (the ratio between the deposited ionization energy and the deposited total energy, < 1.0).

double lindhard_ge_k(vector<double> &x, vector<double> &par): germanium Lindhard function for energy x[0]. The input x is really an array with a single value x[0] which is the input recoil energy in eV. The array par needs exactly 6 parameters to construct the standard Lindhard ionization yield model. These are par[0] the pre-factor to the first energy term in the Lindhard model (unitless); par[1] the exponent of the first energy term; par[2] the pre-factor of thee second energy term; par[3] the exponent of the second energy term; par[4] the atomic number Z; par[5] the mass number A. The return value is the ionization yield (the ratio between the deposited ionization energy and the deposited total energy, < 1.0).

vector<double> geIonizationInRange_k(double E0, double E1, double k, mt19937 *rand): germanium Lindhard ionization in an energy range. E0 is the initial energy of the atom in eV**; E1 is the final energy of the recoiling atom in eV; rand is a random number generator object—specifically std::mt19937.

double lindhard_si_k(vector<double> &x, vector<double> &par): silicon Lindhard function for energy x[0]. The input x is really an array with a single value x[0] which is the input recoil energy in eV. The array par needs exactly 6 parameters to construct the standard Lindhard ionization yield model. These are par[0] the pre-factor to the first energy term in the Lindhard model (unitless); par[1] the exponent of the first energy term; par[2] the pre-factor of thee second energy term; par[3] the exponent of the second energy term; par[4] the atomic number Z; par[5] the mass number A. The return value is the ionization yield (the ratio between the deposited ionization energy and the deposited total energy, < 1.0).

vector<double> siIonizationInRange_k(double E0, double E1, double k, mt19937 *rand): silicon Lindhard ionization in an energy range. E0 is the initial energy of the atom in eV**; E1 is the final energy of the recoiling atom in eV; rand is a random number generator object—specifically std::mt19937.

double lindhard_ar_k(vector<double> &x, vector<double> &par): argon Lindhard function for energy x[0]. The input x is really an array with a single value x[0] which is the input recoil energy in eV. The array par needs exactly 6 parameters to construct the standard Lindhard ionization yield model. These are par[0] the pre-factor to the first energy term in the Lindhard model (unitless); par[1] the exponent of the first energy term; par[2] the pre-factor of thee second energy term; par[3] the exponent of the second energy term; par[4] the atomic number Z; par[5] the mass number A. The return value is the ionization yield (the ratio between the deposited ionization energy and the deposited total energy, < 1.0).

vector<double> arIonizationInRange_k(double E0, double E1, double k, mt19937 *rand): argon Lindhard ionization in an energy range. E0 is the initial energy of the atom in eV**; E1 is the final energy of the recoiling atom in eV; rand is a random number generator object—specifically std::mt19937.

double lindhard_ne_k(vector<double> &x, vector<double> &par): neon Lindhard function for energy x[0]. The input x is really an array with a single value x[0] which is the input recoil energy in eV. The array par needs exactly 6 parameters to construct the standard Lindhard ionization yield model. These are par[0] the pre-factor to the first energy term in the Lindhard model (unitless); par[1] the exponent of the first energy term; par[2] the pre-factor of thee second energy term; par[3] the exponent of the second energy term; par[4] the atomic number Z; par[5] the mass number A. The return value is the ionization yield (the ratio between the deposited ionization energy and the deposited total energy, < 1.0).

vector<double> neIonizationInRange_k(double E0, double E1, double k, mt19937 *rand): neon Lindhard ionization in an energy range. E0 is the initial energy of the atom in eV**; E1 is the final energy of the recoiling atom in eV; rand is a random number generator object—specifically std::mt19937.

weisskopf.h

In this header is contained a prototype for obtaining the Weisskopf decay-time estimate [Weisskopf1951] for a gamma decay of a certain energy (in MeV) and certain multipolarity (like M1, E1, etc.).

Functions

double we(double Egam, double A, std::string transition = "E1"): return the Weisskopf estimated lifetime. Egam is the emitted gamma energy in MeV; A is the mass number of the decaying nucleus; transition is a string that represents the multipolarity of the transition, for example E1 means electric dipole. NOTE: Multipoles beyond M3 are not currently supported and that case will return the slowest multipole that is available. The return value is a double-precision number in fs represented the estimated lifetime of the transition.

isotope_info.h

In this header is contained prototypes for getting various isotope information. In the future this should be replaced with a more robust API to a database to get all of this information. For now, the information needed is hard-coded in the library.

Functions

double getRecoilEnergy(std::string isotope = "70Ge"): get the recoil energy of the ground-state transition for isotope (default 70Ge). Return value is energy in eV.

double getMass(std::string isotope = "70Ge"): get the mass for isotope (default 70Ge). Return value is in GeV.

double getDelta(std::string isotope = "70Ge"): get the mass deficit for isotope (default 70Ge). Return value is in MeV.

double getN(std::string isotope = "70Ge"): get the neutron number for isotope (default 70Ge). Return value is an integer cast to double.

double getZ(std::string isotope = "70Ge"): get the proton (atomic) number for isotope (default 70Ge). Return value is an integer cast to double.

double getSn(std::string isotope = "70Ge"): get the neutron separation energy MeV for isotope (default 70Ge). Return value is in MeV**.

void listStuff(): Print all available information.

rootUtil.h

In this header is contained prototypes for interfacing with the ROOT [ROOT1997] system. This is only for the writing of the output file.

Defines

ROOT_UTILITIES

Functions

bool addToNRTTree(TTree *t, int nr, vector<cri> &recoil_lists, cli cascade_spec): Add data to the output TTree (ROOT data class). t is a pointer to a TTree object which is a ROOT internal object for storing data; nr is the number of recoils in a particular cascade realization; cri is the full recoil object the cri struct. The function returns a bool that indicates success by a true value.