<jupyter_start><jupyter_text>Gamma spectroscopy (of irradiated samples at HISKP Bonn cyclotron)This notebook contains a step-by-step procedure on how to get the activity of each isotope in an irradiated sample.    - For a given detector the following spectra have to be recorded with the **same setup / fixed geometry** :      1. Spectrum of the irradiated sample     2. Spectrum of peaks with well-known energies for an **energy calibration**     3. Spectrum of a calibrated source with well-known activity for an **efficiency calibration**        Spectra 2 & 3 can be the same spectrum if the activity-calibrated source has enough lines e.g. 152Eu. For a change in the detector setup, *each of the above spectra must be measured again*.      - The spectroscopy is performed by the following steps :      1. Fitting the peaks of well-known energy and do an **energy calibration**    2. Fitting peaks of well-known activity of a calibrated source and do an **efficiency calibration**     3. Fitting the energy-calibrated spectrum of the irradiated sample in order to identify the prominent peaks:      - Fit spectrum and try finding isotopes from an isotope library file      - Identify remaining isotopes      - Use efficiency calibration to determine actual activity from integrated counts for each identified isotope<jupyter_code># imports
import os
import yaml
import logging
import numpy as np  # c-like, vectorized arrays
import irrad_spectroscopy as isp
import irrad_spectroscopy.spectroscopy as sp  # import main spectroscopy methods
import matplotlib # plotting library
matplotlib.use('qt5agg')  # set plotting backend
import matplotlib.pyplot as plt  # shortcut to plotting functions
from collections import OrderedDict  # dict with chronological key order

# import general utils
from irrad_spectroscopy.spec_utils import get_measurement_time, get_isotope_info, source_to_dict, select_peaks

# plotting
from irrad_spectroscopy.plot_utils import plot_spectrum

# interactive plotting
%matplotlib notebook

# figsize
plt.rcParams['figure.figsize'] = [9.5, 5.5]<jupyter_output><empty_output><jupyter_text>Set paths and load filesDefine locations and files. Give the names of each required file as they are named in their respective location. **All recorded spectra are expected to be in a `.txt` file and the meta data in a `.mcd` file with the same name.**<jupyter_code># name
SAMPLE_NAME = 'example_sample'

# paths
package_path = isp.package_path

example_path = os.path.join(package_path, 'examples')

data_path = os.path.join(example_path, 'example_data')

# Create output subfolder where all output data and plots are stored
output_path = os.path.join(example_path, 'output')
if not os.path.exists(output_path):
    os.makedirs(output_path)

# load files, all have leading _
try:
    _spectrum_sample = np.loadtxt(os.path.join(data_path, 'example_sample.txt'), unpack=True)
    
    _spectrum_bkg = np.loadtxt(os.path.join(data_path, 'background.txt'), unpack=True)

    _t_sample = get_measurement_time(os.path.join(data_path, 'example_sample.txt'))
    
    _t_bkg = get_measurement_time(os.path.join(data_path, 'background.txt'))

    _spectrum_energy_calibration = np.loadtxt(os.path.join(data_path, 'Eu152_point.txt'), unpack=True)
    
    # use same spectrum for energy and efficiency calibration; 152Eu spectrum
    _spectrum_efficiency_calibration = _spectrum_energy_calibration
    
    with open(os.path.join(data_path, '152Eu_peaks.yaml'), 'r') as energy_peaks:
        _peaks_energy_calibration = yaml.safe_load(energy_peaks)
        
    with open(os.path.join(data_path, 'background_peaks.yaml'), 'r') as background_peaks:
        _peaks_bkg = yaml.safe_load(background_peaks)
        
    with open(os.path.join(data_path, '152Eu_point_source.yaml'), 'r') as source_specs:
        _calibrated_source_specs = yaml.safe_load(source_specs)

except IOError:
    logging.exception('Failed to load some of the data!')<jupyter_output><empty_output><jupyter_text>Energy calibration with well-known spectrumFit well-known spectrum and make a calibration from fitted channel numbers and known peak energies.<jupyter_code># fit peaks for energy calibration
peaks_energy_cal, bkg_energy_cal = sp.fit_spectrum(counts=_spectrum_energy_calibration[1],
                                                   expected_peaks=_peaks_energy_calibration['channel'])
# do energy_calibration
energy_calibration = sp.do_energy_calibration(peaks_energy_cal, _peaks_energy_calibration['energy'])
# plot everything
plot_spectrum(counts=_spectrum_energy_calibration[1],
              peaks=peaks_energy_cal,
              bkg=bkg_energy_cal,
              plot_calib=energy_calibration,
              title='Energy calibration',
              output_plot=os.path.join(output_path, '%s_energy_calibration.pdf' % SAMPLE_NAME))<jupyter_output><empty_output><jupyter_text>Efficiency calibration of detector with calibrated sourceDetermination by comparing calibrated activity with integrated counts in spectrum. The current acitivity of the source is calculated by the time passed between the date of the calibration and date of measurement which can be accessed by the keys `datetime_calibration` and `datetime_measurement` of the `_calibrated_source_specs` dictionary.<jupyter_code># fit peaks for efficiency calibration
peaks_efficiency_cal, bkg_efficiency_cal = sp.fit_spectrum(counts=_spectrum_efficiency_calibration[1],
                                                           energy_cal=energy_calibration['func'],
                                                           expected_peaks=source_to_dict(_calibrated_source_specs))
# do energy_calibration
efficiency_calibration = sp.do_efficiency_calibration(observed_peaks=peaks_efficiency_cal,
                                                      source_specs=_calibrated_source_specs)
# plot everything
plot_spectrum(counts=_spectrum_efficiency_calibration[1],
              peaks=peaks_efficiency_cal,
              bkg=bkg_efficiency_cal,
              energy_cal=energy_calibration['func'],
              plot_calib=efficiency_calibration,
              title='Efficiency calibration',
              output_plot=os.path.join(output_path, '%s_efficiency_calibration.pdf' % SAMPLE_NAME))<jupyter_output><empty_output><jupyter_text>Identify peaks in sample spectrum  - Optionally compare background peaks to sample peaks in order to identify sample peaks / growth of background  - Optionally subtract background from sample spectrum   - Optionally fit global background in spectrum with energy calibration  - Try fitting all peaks from the isotope library as expected peaks in spectrum with energy calibration  - Fit additional n_peaks peaks which have not been identified from the isotope library    - Identify peaks with help of http://nucleardata.nuclear.lu.se/toi/radSearch.asp data bank, the isotope table PDF in `static` and http://www.oecd-nea.org/janisweb/ cross sections<jupyter_code># checking all peaks from library and background first
peaks_sample, bkg_sample = sp.fit_spectrum(counts=_spectrum_sample[1],
                                           energy_cal=energy_calibration['func'],
                                           efficiency_cal=efficiency_calibration['func'],
                                           t_spec=_t_sample,
                                           expected_accuracy=energy_calibration['accuracy'])
# plot everything
plot_spectrum(counts=_spectrum_sample[1],
              peaks=peaks_sample,
              bkg=bkg_sample,
              energy_cal=energy_calibration['func'],
              title='Gamma spectrum of sample',
              output_plot=os.path.join(output_path, '%s_sample_spectrum.pdf' % SAMPLE_NAME))

# plot backgroud and check which lines are unique to the sample
plt.plot(energy_calibration['func'](_spectrum_bkg[0]), _spectrum_bkg[1], label='Background spectrum')
_ = plt.legend()<jupyter_output><empty_output><jupyter_text>Isotope selectionSelect the isotopes which are unique to the spectrum of the sample<jupyter_code># select the isotopes unique to the sample
selected_peaks = select_peaks(selection=['65_Zn', '48_V', '7_Be'], peaks=peaks_sample)<jupyter_output><empty_output><jupyter_text>Activity determinationSumming the probability of each peak as well as the activity to calculate the summed and scaled to 100% activity for each isotope<jupyter_code>sample_activities = sp.get_activity(selected_peaks)
for s in sample_activities:
    print('{}:\t({}+-{}) Bq'.format(s, sample_activities[s]['nominal'], sample_activities[s]['sigma']))<jupyter_output><empty_output><jupyter_text>Dose calculationCalculates the dose in uSv or dose rate in uSv/h isotope-wise and the summed total dose of the all isotopes e.g. the sample<jupyter_code>sample_doses = sp.get_dose(selected_peaks, distance=50, time=2000, material='air')
print('({}+-{}) {}'.format(sample_doses['nominal'], sample_doses['sigma'], sample_doses['unit']))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07590280873829693, 'fraction_numerical': 0.005572893446277307, 'mean_word_length': 4.1186537364517966, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '331033', 'license': 'MIT', 'n_tokens_mistral': 2610, 'n_tokens_neox': 2435, 'n_words': 738, 'path': 'irrad_spectroscopy/examples/example_spectroscopy.ipynb', 'repo_name': 'leloup314/irrad_spectroscopy', 'size': '8972'}