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 MonteCarlo 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
halflives of energy level in fs

bool success

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
electronequivalent ionization energy eV

vector<double> time
time since capture fs

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

int n
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; andsuccess
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; andsuccess
istrue
if it can be interpreted as a doubleprecision number (either in decimal or scientific notation). Acli
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 doubleprecision number (representing the neutron separation energy) OR it is a recognized isotopic symbol like74Ge
. 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 totrue
if a Weisskopf multipolarity abbreviation is recognized, likew(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 totrue
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 totrue
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 whitespacesplit 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 constantacceleration 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 sizenlev
in keV;taus
is an array of level lifetimes of sizenlev
in fs;A
is the mass number;mtrand
is the random number generator object—in this casestd::mt19937
. The function returns acri
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 sizenlev
in keV;taus
is an array of level lifetimes of sizenlev
in fs;A
is the mass number;mtrand
is the random number generator object—in this casestd::mt19937
. The function returns acri
struct.

double geDecay(double v, double M, double Egam, mt19937 *rand)
germanium func. to return the Energy after the midstop 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; andmtrand
is the random number generator object—in this casestd::mt19937
. The function returns the energy in eV after the midstep 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 casestd::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 stoppingpower 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 stoppingpower 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 stoppingpower 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 sizenlev
in keV;taus
is an array of level lifetimes of sizenlev
in fs;A
is the mass number;mtrand
is the random number generator object—in this casestd::mt19937
. The function returns acri
struct.

double siDecay(double v, double M, double Egam, mt19937 *rand)
silicon func. to return the Energy after the midstop 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; andmtrand
is the random number generator object—in this casestd::mt19937
. The function returns the energy in eV after the midstep 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 casestd::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 stoppingpower 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 stoppingpower 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 stoppingpower 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 sizenlev
in keV;taus
is an array of level lifetimes of sizenlev
in fs;A
is the mass number;mtrand
is the random number generator object—in this casestd::mt19937
. The function returns acri
struct.

double arDecay(double v, double M, double Egam, mt19937 *rand)
argon func. to return the Energy after the midstop 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; andmtrand
is the random number generator object—in this casestd::mt19937
. The function returns the energy in eV after the midstep 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 casestd::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 stoppingpower 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 stoppingpower 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 stoppingpower 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 sizenlev
in keV;taus
is an array of level lifetimes of sizenlev
in fs;A
is the mass number;mtrand
is the random number generator object—in this casestd::mt19937
. The function returns acri
struct.

double neDecay(double v, double M, double Egam, mt19937 *rand)
neon func. to return the Energy after the midstop 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; andmtrand
is the random number generator object—in this casestd::mt19937
. The function returns the energy in eV after the midstep 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 casestd::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 stoppingpower 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 stoppingpower 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 stoppingpower 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 Poissondistributed integer random number with the Knuth method.
lambda
is the parameter of the Poisson distribution;rand
is a random generator object—specifically thestd::mt19937
. The return value is a single integer realization of a Poisson random variable.

int poissonAtkinson(double lambda, mt19937 *rand)
Generate a Poissondistributed 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 thestd::mt19937
. The return value is a single integer realization of a Poisson random variable.

int poisson(double lambda, mt19937 *rand)
Generate a Poissondistributed 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 thestd::mt19937
. The return value is a single integer realization of a Poisson random variable.

int poissonFano(double lambda, double F, mt19937 *rand)
Generate a sudoPoissondistributed integer random number that sets the width based on the socalled “Fano Factor” which is an empirical factor developed in the 50’s for widening a sortof 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 thestd::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 valuex[0]
which is the input recoil energy in eV. The arraypar
needs exactly 6 parameters to construct the standard Lindhard ionization yield model. These arepar[0]
the prefactor to the first energy term in the Lindhard model (unitless);par[1]
the exponent of the first energy term;par[2]
the prefactor 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 inputx
is really an array with a single valuex[0]
which is the input recoil energy in eV. The arraypar
needs exactly 6 parameters to construct the standard Lindhard ionization yield model. These arepar[0]
the prefactor to the first energy term in the Lindhard model (unitless);par[1]
the exponent of the first energy term;par[2]
the prefactor 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—specificallystd::mt19937
.

double lindhard_si_k(vector<double> &x, vector<double> &par)
silicon Lindhard function for energy
x[0]
. The inputx
is really an array with a single valuex[0]
which is the input recoil energy in eV. The arraypar
needs exactly 6 parameters to construct the standard Lindhard ionization yield model. These arepar[0]
the prefactor to the first energy term in the Lindhard model (unitless);par[1]
the exponent of the first energy term;par[2]
the prefactor 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—specificallystd::mt19937
.

double lindhard_ar_k(vector<double> &x, vector<double> &par)
argon Lindhard function for energy
x[0]
. The inputx
is really an array with a single valuex[0]
which is the input recoil energy in eV. The arraypar
needs exactly 6 parameters to construct the standard Lindhard ionization yield model. These arepar[0]
the prefactor to the first energy term in the Lindhard model (unitless);par[1]
the exponent of the first energy term;par[2]
the prefactor 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—specificallystd::mt19937
.

double lindhard_ne_k(vector<double> &x, vector<double> &par)
neon Lindhard function for energy
x[0]
. The inputx
is really an array with a single valuex[0]
which is the input recoil energy in eV. The arraypar
needs exactly 6 parameters to construct the standard Lindhard ionization yield model. These arepar[0]
the prefactor to the first energy term in the Lindhard model (unitless);par[1]
the exponent of the first energy term;par[2]
the prefactor 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—specificallystd::mt19937
.
weisskopf.h
In this header is contained a prototype for obtaining the Weisskopf decaytime 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 exampleE1
means electric dipole. NOTE: Multipoles beyondM3
are not currently supported and that case will return the slowest multipole that is available. The return value is a doubleprecision 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 hardcoded in the library.
Functions

double getRecoilEnergy(std::string isotope = "70Ge")
get the recoil energy of the groundstate 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 aTTree
object which is aROOT
internal object for storing data;nr
is the number of recoils in a particular cascade realization;cri
is the full recoil object thecri
struct. The function returns a bool that indicates success by atrue
value.