Decay Chain

class curie.DecayChain(parent_isotope, R=None, A0=None, units='s', timestamp=True)

Radioactive Decay Chain

Uses the Bateman equations to calculate the activities and number of decays from a radioactive decay chain as a function of time, both in production and decay. Also, initial isotope activities and production rates can be fit to observed count data, or directly fit to HPGe spectra using the get_counts() function.

Parameters:

parent_isotopestr

Parent isotope in the chain.

Rarray_like, dict, str or pd.DataFrame

Production rate for each isotope in the decay chain as a function of time. If a Nx2 np.ndarray, element n gives the production rate R_n up until time t_n for the parent isotope. E.g. If the production rate of the parent is 5 for 1 hour, and 8 for 3 hours, the array will be [[5, 1], [8, 4]]. If instead time intervals are preferred to a monotonically increasing grid of timestamps, set ‘timestamp=False’. In this case the production rate array will be [[5, 1], [8, 3]]. (R=5 for 1 hour, R=8 for 3 hours).

If R is a dict, it specifies the production rate for multiple isotopes, where the keys are the isotopes and the values are type np.ndarray.

If R is a pd.DataFrame, it must have columns ‘R’ and ‘time’, and optionally ‘isotope’ if R>0 for any isotopes other than the parent. If R is a str, it must be a path to a file where the same data is provided. Supported file types are .csv, .json and .db files, where .json files must be in the ‘records’ format, and .db files must have a table named ‘R’. Also, each isotope must have the same time grid, for both timestamp=True and timestamp=False.

A0float or dict

Initial activity. If a float, the initial activity of the parent isotope. If a dict, the keys are the isotopes for which the values represent the initial activity.

unitsstr, optional

Units of time for the chain. Options are ‘ns’, ‘us’, ‘ms’, ‘s’, ‘m’, ‘h’, ‘d’, ‘y’, ‘ky’, ‘My’, ‘Gy’. Default is ‘s’.

timestampbool, optional

Determines if the ‘time’ variable in R is to be read as a timestamp-like grid, i.e. in monotonically increasing order, or as a series of time intervals. Default is True.

Examples

>>> dc = ci.DecayChain('Ra-225', R=[[1.0, 1.0], [0.5, 1.5], [2.0, 6]], units='d')
>>> print(dc.isotopes)
['225RAg', '225ACg', '221FRg', '217ATg', '213BIg', '217RNg', '209TLg', '213POg', '209PBg', '209BIg']
>>> print(dc.R_avg)
  R_avg isotope
0  1.708333  225RAg

>>> dc = ci.DecayChain('152EU', A0=3.7E3, units='h')
>>> print(ci.isotopes)
['152EUg', '152GDg', '152SMg']

Attributes:

Rpd.DataFrame

Production rate as a function of time, for each isotope in the chain. This will be modified if fit_R() is called.

A0dict

Initial activity of each isotope in the chain.

isotopeslist

List of isotopes in the decay chain.

countspd.DataFrame

Observed counts from isotopes in the decay chain, which can be used to determine the initial activities or average production rates using the fit_R() or fit_A0() functions.

R_avgpd.DataFrame

Time-averaged production rate for each isotope where R>0. This will be modified if fit_R() is called.

Methods

activity(isotope, time[, units, _R_dict, ...])	Activity of an isotope in the chain
decays(isotope, t_start, t_stop[, units, ...])	Number of decays in a given time interval
fit_A0([max_error, min_counts])	Fit the initial activity to count data
fit_R([max_error, min_counts])	Fit the production rate to count data
get_counts(spectra, EoB[, peak_data])	Retrieves the number of measured decays
plot([time, max_plot, max_label, ...])	Plot the activities in the decay chain

activity(isotope, time, units=None, _R_dict=None, _A_dict=None)

Activity of an isotope in the chain

Computes the activity of a given isotope in the decay chain at a given time. Units of activity are in Bq. Units of time must be either the units for the DecayChain (default ‘s’), or specified by the units keyword.

Parameters:

isotopestr

Isotope for which the activity is calculated.

timearray_like

Time to calculate the activity. Units of time must be the same as the decay chain, or be given by units. Note that if R!=0, time=0 is defined as the end of production time. Else, if A0!=0, time=0 is defined as the time at which the specified activities equaled A0. t<0 is not allowed.

unitsstr, optional

Units of time, if different from the units of the decay chain.

Returns:

activitynp.ndarray

Activity of the given isotope in Bq.

Examples

>>> dc = ci.DecayChain('152EU', A0=3.7E3, units='h')
>>> print(dc.activity('152EU', time=0))
3700.0
>>> print(dc.activity('152EU', time=13.537, units='y'))
1849.999906346199

decays(isotope, t_start, t_stop, units=None, _A_dict=None)

Number of decays in a given time interval

Computes the number of decays from a given isotope in the decay chain in the time interal t_start to t_stop. The units of t_start and t_stop must be either the same units as the decay chain, or be specified by the units keyword.

Parameters:

isotopestr

Isotope for which the number of decays is calculated.

t_startarray_like

Time of the start of the interval.

t_stoparray_like

Time of the end of the interval.

unitsstr, optional

Units of time, if different from the units of the decay chain.

Returns:

decaysnp.ndarray

Number of decays

Examples

>>> dc = ci.DecayChain('152EU', A0=3.7E3, units='h')
>>> print(dc.decays('152EU', t_start=1, t_stop=2))
13319883.293399204
>>> print(dc.decays('152EU', t_start=50, t_stop=50.1, units='y'))
900151618.5228329

fit_A0(max_error=0.4, min_counts=1)

Fit the initial activity to count data

Fits a scalar multiplier to the initial activity for each isotope specified in self.A0. The fit minimizes to the number of measured decays (self.counts) as a function of time, rather than the activity, because the activity at each time point may be sensitive to the shape of the decay curve.

Parameters:

max_errorfloat, optional

The maximum relative error of a count datum to include in the fit. E.g. 0.25=25% (ony points with less than 25% error will be shown). Default, 0.4 (40%).

min_countsfloat or int, optional

The minimum number of counts (decays) for a datum in self.counts to be included in the fit. Default, 1.

Returns:

isotopeslist

List of isotopes where A0>0. Same indices as fit. (i.e. isotope[0] corresponds to fit[0] and cov[0][0].)

fitnp.ndarray

The initial activity for each isotope where A0>0.

covnp.ndarray

Covariance matrix on the fit.

Examples

>>> sp = ci.Spectrum('eu_calib_7cm.Spe')
>>> sp.isotopes = ['152EU']

>>> dc = ci.DecayChain('152EU', A0=3.7E4, units='d')
>>> dc.get_counts([sp], EoB='01/01/2016 08:39:08')
>>> print(dc.fit_A0())
(['152EUg'], array([6501.93665952]), array([[42425.53832341]]))

fit_R(max_error=0.4, min_counts=1)

Fit the production rate to count data

Fits a scalar multiplier to the production rate (as a function of time) for each isotope specified in self.R. The fit minimizes to the number of measured decays (self.counts) as a function of time, rather than the activity, because the activity at each time point may be sensitive to the shape of the decay curve.

Parameters:

max_errorfloat, optional

The maximum relative error of a count datum to include in the fit. E.g. 0.25=25% (ony points with less than 25% error will be shown). Default, 0.4 (40%).

min_countsfloat or int, optional

The minimum number of counts (decays) for a datum in self.counts to be included in the fit. Default, 1.

Returns:

isotopeslist

List of isotopes where R>0. Same indices as fit. (i.e. isotope[0] corresponds to fit[0] and cov[0][0].)

fitnp.ndarray

The fitted time-averaged production rate for each isotope where R>0.

covnp.ndarray

Covariance matrix on the fit.

Examples

>>> sp = ci.Spectrum('eu_calib_7cm.Spe')
>>> sp.isotopes = ['152EU']

>>> dc = ci.DecayChain('152EU', R=[[3E5, 36.0]], units='d')
>>> dc.get_counts([sp], EoB='01/01/2016 08:39:08')
>>> print(dc.fit_R())
(array(['152EUg'], dtype=object), array([1291584.51735774]), array([[1.67412376e+09]]))

get_counts(spectra, EoB, peak_data=None)

Retrieves the number of measured decays

Takes the number of measured decays from one of the following: a list of spectra, a file with peak data, or a pandas DataFrame with peak data.

Parameters:

spectralist or str

List of ci.Spectrum objects, or str of spectra filenames. If list of str, peak_data must be specified. In this case the filenames must be an exact match of the filenames in peak_data. If spectra is a str, it is assumed to be a regex match for the filenames in peak_data.

EoBstr or datetime.datetime

Date/time of end-of-bombardment (t=0). Must be a datetime object or a string in the format ‘%m/%d/%Y %H:%M:%S’. This is used to calculate the decay time for the count.

peak_datastr or pd.DataFrame, optional

Either a file path to a file that was created using ci.Spectrum.saveas() or a DataFrame with the same structure as ci.Spectrum.peaks.

Examples

>>> sp = ci.Spectrum('eu_calib_7cm.Spe')
>>> sp.isotopes = ['152EU']
>>> sp.saveas('test_spec.json')

>>> dc = ci.DecayChain('152EU', A0=3.7E3, units='h')
>>> dc.get_counts([sp], EoB='01/01/2016 08:39:08')
>>> dc.get_counts(['eu_calib_7cm.Spe'], EoB='01/01/2016 08:39:08', peak_data='test_spec.json')
>>> print(dc.counts)

plot(time=None, max_plot=10, max_label=10, max_plot_error=0.4, max_plot_chi2=10, **kwargs)

Plot the activities in the decay chain

Plots the activities as a function of time for all radioactive isotopes in the decay chain. Can plot along a specified time grid, else the time will be inferred from the half-life of the parent isotope, or any count information given to self.counts.

Parameters:

timearray_like, optional

Time grid along which to plot. Units must be the same as the decay chain.

max_plotint, optional

Maximum number of isotope activities to plot in the decay chain. Default, 10.

max_labelint, optional

Maximum number of isotope activities to label in the legend. Default, 10.

max_plot_errorfloat, optional

The maximum relative error of a count point to include on the plot. E.g. 0.25=25% (ony points with less than 25% error will be shown). Default, 0.4.

max_plot_chi2float or int, optional

Maximum chi^2 of a count point to include on the plot. Only points with a chi^2 less than this value will be shown. Default, 10.

Other Parameters:

**kwargs

Optional keyword arguments for plotting. See the plotting section of the curie API for a complete list of kwargs.

Examples

>>> sp = ci.Spectrum('eu_calib_7cm.Spe')
>>> sp.isotopes = ['152EU']

>>> dc = ci.DecayChain('152EU', R=[[3E5, 36.0]], units='d')
>>> dc.get_counts([sp], EoB='01/01/2016 08:39:08')
>>> dc.fit_R()
>>> dc.plot()

>>> dc = ci.DecayChain('99MO', A0=350E6, units='d')
>>> dc.plot()