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.

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(**kwargs)	Fit the initial activity to count data
fit_R(**kwargs)	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

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', '205TLg']
>>> print(dc.R_avg)
  R_avg isotope
0  1.708333  225RAg

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

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

property diagnostics

Fit diagnostics from the most recent fit_R or fit_A0 call

Read-only pd.DataFrame with one row per fitted isotope (fit is ‘fit_R’ or ‘fit_A0’; the fit is joint, so chi2, dof, n_points, n_dropped and scale_factor repeat on every row): columns chi2 (reduced, of the unscaled fit), dof, n_points (counts the fit used), n_dropped (counts removed by the max_error/min_counts filters), converged, model (’’ - the chain fit has no model selection), scale_factor (uncertainty inflation applied; 1.0 = none), flags (comma-joined, e.g. ‘unmoved’, ‘at_bound:R’, ‘singular_cov’), message (summary and warning text) and isotope. Empty (with the full schema) before any fit; rebuilt on each fit_R()/fit_A0() call. Accessing it never triggers a fit.

fit_A0(**kwargs)

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.

Keyword arguments other than cov are fit_config keys: they merge into dc.fit_config (like Spectrum.fit_config) and persist for subsequent fits.

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.

max_chi2float, optional

Exclude counts whose originating peak fit had chi2/dof above this bound (needs the chi2 column get_counts carries over from the peak data; counts without one are unaffected). Default None (off).

exclude_lineslist, optional

Named lines to exclude from the fit: entries are an isotope name (‘152EU’) or an (isotope, energy_keV) tuple. The energy matches the closest of that isotope’s line energies in the counts, accepted within +-0.5 keV; an entry matching nothing excludes nothing and warns, naming the nearest candidate. Default None (off).

time_rangedict, optional

Per-isotope count-time window, {isotope: (t_lo, t_hi)} in chain units; None for an open bound. Counts starting outside the window are excluded. Default None (off).

corrfloat, optional

Opt-in uniform-correlation mode for counts that carry only total uncertainties: the fraction (0-1) of each point’s variance treated as fully correlated across points (or within corr_group groups). Default None (off).

corr_groupstr, optional

Name of a counts column defining the correlation groups for corr. Default None (global).

norm_fracfloat, optional

Fraction of each line’s intensity variance treated as common to all lines of the isotope (the decay-scheme normalization). Default 1.0 (fully common - conservative until the intensity data carry the normalization uncertainty separately).

scale_factorbool, optional

If True (default), the fitted covariance is inflated by chi-square per degree of freedom when that exceeds 1 (mutually inconsistent count data); it is never deflated.

covarray_like, optional

Full covariance matrix of the count data, overriding the assembled one. For advanced use (kwarg only - not a fit_config key).

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')
>>> itp, A0, cov = dc.fit_A0()
>>> print(itp[0], A0[0])
152EUg 138582.75831764835

fit_R(**kwargs)

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.

Keyword arguments other than p0 and cov are fit_config keys: they merge into dc.fit_config (like Spectrum.fit_config) and persist for subsequent fits.

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.

max_chi2float, optional

Exclude counts whose originating peak fit had chi2/dof above this bound (needs the chi2 column get_counts carries over from the peak data; counts without one are unaffected). Default None (off).

exclude_lineslist, optional

Named lines to exclude from the fit: entries are an isotope name (‘152EU’) or an (isotope, energy_keV) tuple. The energy matches the closest of that isotope’s line energies in the counts, accepted within +-0.5 keV; an entry matching nothing excludes nothing and warns, naming the nearest candidate. Default None (off).

time_rangedict, optional

Per-isotope count-time window, {isotope: (t_lo, t_hi)} in chain units; None for an open bound (e.g. {'196AUm2': (None, 60.0)}). Counts starting outside the window are excluded. Default None (off).

unc_R_floorfloat, optional

Relative floor on the returned production-rate uncertainty: unc_R >= floor x R, per isotope, applied after the fit (announced). E.g. 0.05 keeps every unc_R at or above 5%. Default None (off).

corrfloat, optional

Opt-in uniform-correlation mode for counts that carry only total uncertainties: the fraction (0-1) of each point’s variance treated as fully correlated across points (or within corr_group groups). Default None (off).

corr_groupstr, optional

Name of a counts column defining the correlation groups for corr. Default None (global).

norm_fracfloat, optional

Fraction of each line’s intensity variance treated as common to all lines of the isotope (the decay-scheme normalization). Default 1.0 (fully common - conservative until the intensity data carry the normalization uncertainty separately).

scale_factorbool, optional

If True (default), the fitted covariance is inflated by chi-square per degree of freedom when that exceeds 1 (mutually inconsistent count data); it is never deflated.

p0dict, optional

Starting estimate override, {isotope: R_estimate} in production-rate units: seeds the fit at the given rate instead of the data-derived estimate (kwarg only - not a fit_config key). Default None.

covarray_like, optional

Full covariance matrix of the count data, overriding the assembled one. For advanced use (kwarg only - not a fit_config key).

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())
(['152EUg'], array([27527059.31414273]), array([[2.03699956e+11]]))

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 draw as a filled marker. E.g. 0.25=25%. Points above the threshold are drawn as open grey markers rather than hidden. Default, 0.4.

max_plot_chi2float or int, optional

Maximum chi^2 of a count point to draw as a filled marker. Points above the threshold are drawn as open grey markers rather than hidden. 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()