{ "cells": [ { "cell_type": "markdown", "id": "b6c7b01d", "metadata": {}, "source": [ "# The `Spectrum` Class\n", "\n", "The {py:obj}`Spectrum ` class, defined in the {py:obj}`snf_simulations.spec ` module, is the core data class used in the package to represent the antineutrino spectrum emitted by spent nuclear fuel (SNF).\n", "\n", "Each `Spectrum` is a histogram, with energy bins as the x-axis and the corresponding antineutrino flux values, and uncertainties, as the y-axis. This class provides methods for manipulating and analyzing the spectrum data, such as equalising, integration and sampling.\n", "\n", "\n", "## Creating a `Spectrum`\n", "\n", "You can create a new `Spectrum` object by providing the necessary parameters: energy bins (in keV), flux values (in keV^-1) and uncertainties on the flux values. Each of these should be Numpy arrays.\n", "\n", ":::{important}\n", "When giving an array of energy values, the values are defined as the edges of the histogram bins, meaning that if you have $N$ bins, you should provide $N+1$ energy values.\n", ":::\n", "\n", "For example:" ] }, { "cell_type": "code", "execution_count": 1, "id": "72228bca", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "import numpy as np\n", "from snf_simulations.spec import Spectrum\n", "\n", "# Input arrays\n", "energy_bins = np.array([0.0, 1.0, 2.0, 3.0, 4.0])\n", "flux_values = np.array([10, 20, 15, 5])\n", "flux_errors = np.array([1, 2, 1.5, 0.5])\n", "\n", "# Create a Spectrum object\n", "spectrum = Spectrum(energy_bins, flux_values, flux_errors)\n", "print(spectrum)" ] }, { "cell_type": "markdown", "id": "c23582a1", "metadata": {}, "source": [ "### Creating a `Spectrum` for a known isotope\n", "\n", "In practice, often be easier to create a `Spectrum` object for a known isotope using the {py:func}`from_isotope ` class method, which retrieves the spectrum data from the internal database. For example:" ] }, { "cell_type": "code", "execution_count": 2, "id": "5b9eae40", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "spectrum = Spectrum.from_isotope(\"Am242\")\n", "print(spectrum)" ] }, { "cell_type": "markdown", "id": "c68e3a39", "metadata": {}, "source": [ "This uses the {py:func}`load_spectrum ` function to load the spectrum data from one of the included antineutrino spectrum data files. For isotopes that are not included in the database, this will raise a `ValueError`, and you will need to create the `Spectrum` object manually using the constructor as shown in the previous example.\n", "\n", ":::{note}\n", "In the future, isotope data will be downloaded directly from the IAEA. See [this GitHub Issue](https://github.com/ekneale/SNF-simulations/issues/16) for details.\n", ":::" ] }, { "cell_type": "markdown", "id": "76b541a5", "metadata": {}, "source": [ "## Manipulating `Spectrum` Objects\n", "\n", "The `Spectrum` class provides several methods for manipulating the spectrum data.\n", "\n", "### Equalising a `Spectrum`\n", "\n", "As Spectrum objects are histograms, they can be equalised to a different set of energy bins using the {py:func}`equalise ` method. This is useful when you want to compare spectra that have different binning, or when you want to re-bin a spectrum to match the binning of another spectrum.\n", "\n", "The {py:func}`equalise ` method takes three optional parameters: `width`, `min_energy` and `max_energy`. If `width` is provided, the method will create new energy bins with the specified width, starting from `min_energy` and ending at `max_energy`. If `width` is not provided, the method will default to bins of width 1 keV. `min_energy` and `max_energy` default to the minimum and maximum energy values of the original spectrum, respectively.\n", "\n", "For example:" ] }, { "cell_type": "code", "execution_count": 3, "id": "a3a83e04", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]\n", "[10. 12.5 17.5 18.75 16.25 12.5 7.5 5. 0. 0. ]\n", "[1. 0.90138782 1.52069063 1.54616461 1.23110723 1.13192314\n", " 0.53033009 0.5 0. 0. ]\n" ] } ], "source": [ "energy_bins = np.array([0, 2, 4, 6, 8])\n", "flux_values = np.array([10, 20, 15, 5])\n", "flux_errors = np.array([1, 2, 1.5, 0.5])\n", "\n", "spectrum = Spectrum(energy_bins, flux_values, flux_errors)\n", "\n", "# Equalise the spectrum to new bins\n", "spectrum.equalise(width=1, max_energy=10)\n", "print(spectrum.energy)\n", "print(spectrum.flux)\n", "print(spectrum.errors)" ] }, { "cell_type": "markdown", "id": "d3defb94", "metadata": {}, "source": [ "Note that any bins outside of the original energy range will be filled with zeros, and the corresponding uncertainties will also be set to zero.\n", "\n", "### Scaling a `Spectrum`\n", "\n", "Spectrum objects can be scaled by a constant factor using the `*` operator. This can be useful for normalising spectra, or to scale by an activity level. For example:" ] }, { "cell_type": "code", "execution_count": 4, "id": "f0c0c983", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[20 40 30 10]\n", "[2. 4. 3. 1.]\n" ] } ], "source": [ "energy_bins = np.array([0, 2, 4, 6, 8])\n", "flux_values = np.array([10, 20, 15, 5])\n", "flux_errors = np.array([1, 2, 1.5, 0.5])\n", "\n", "spectrum = Spectrum(energy_bins, flux_values, flux_errors)\n", "\n", "# Scale the spectrum by a factor of 2\n", "scaled_spectrum = spectrum * 2\n", "print(scaled_spectrum.flux)\n", "print(scaled_spectrum.errors)" ] }, { "cell_type": "markdown", "id": "175c83f8", "metadata": {}, "source": [ "### Integrating a `Spectrum`\n", "\n", "The {py:func}`integrate ` method can be used to calculate the total flux in a specified energy range. This is done by integrating the flux values over the specified energy range, taking into account the bin widths and uncertainties. The method takes two parameters: `lower_energy` and `upper_energy`, which define the energy range for the integration. For example:" ] }, { "cell_type": "code", "execution_count": 5, "id": "9547f1d2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total flux: 65.0\n" ] } ], "source": [ "energy_bins = np.array([0, 2, 4, 6, 8])\n", "flux_values = np.array([10, 20, 15, 5])\n", "flux_errors = np.array([1, 2, 1.5, 0.5])\n", "\n", "spectrum = Spectrum(energy_bins, flux_values, flux_errors)\n", "\n", "# Integrate the spectrum from 1 keV to 5 keV\n", "total_flux = spectrum.integrate(lower_energy=1, upper_energy=5)\n", "print(f\"Total flux: {total_flux}\")" ] }, { "cell_type": "markdown", "id": "16e8f8b8", "metadata": {}, "source": [ "### Adding `Spectrum` objects\n", "\n", "Two `Spectrum` objects can be added together using the `+` operator, which will add the flux values and combine the uncertainties in quadrature. For example:" ] }, { "cell_type": "code", "execution_count": 6, "id": "72a66822", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[15 30 30 25]\n", "[1.11803399 2.23606798 2.12132034 2.06155281]\n" ] } ], "source": [ "energy_bins1 = np.array([0, 2, 4, 6, 8])\n", "flux_values1 = np.array([10, 20, 15, 5])\n", "flux_errors1 = np.array([1, 2, 1.5, 0.5])\n", "\n", "spectrum1 = Spectrum(energy_bins1, flux_values1, flux_errors1)\n", "\n", "energy_bins2 = np.array([0, 2, 4, 6, 8])\n", "flux_values2 = np.array([5, 10, 15, 20])\n", "flux_errors2 = np.array([0.5, 1, 1.5, 2])\n", "\n", "spectrum2 = Spectrum(energy_bins2, flux_values2, flux_errors2)\n", "\n", "# Add the two spectra together\n", "combined_spectrum = spectrum1 + spectrum2\n", "print(combined_spectrum.flux)\n", "print(combined_spectrum.errors)" ] }, { "cell_type": "markdown", "id": "707d257a", "metadata": {}, "source": [ "The energy bins of the two spectra **must** be the same for the addition to work. If they are not, you can use the {py:func}`equalise ` method to re-bin one of the spectra to match the other before adding them together." ] }, { "cell_type": "markdown", "id": "6d52008a", "metadata": {}, "source": [ "### Sampling from a `Spectrum`\n", "\n", "The {py:func}`sample ` method can be used to generate random samples of antineutrino energies based on the flux values in the spectrum. This is done by treating the flux values as a probability distribution and sampling from it accordingly. The method takes one parameter: `samples`, which specifies the number of random samples to generate. For example:" ] }, { "cell_type": "code", "execution_count": 7, "id": "2061ab0c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[6.86677734 2.16792332 7.05107869 3.7394425 4.95141982 5.8755514\n", " 5.88090146 3.72116845 3.73950457 7.78831592]\n" ] } ], "source": [ "energy_bins = np.array([0, 2, 4, 6, 8])\n", "flux_values = np.array([10, 20, 15, 5])\n", "flux_errors = np.array([1, 2, 1.5, 0.5])\n", "\n", "spectrum = Spectrum(energy_bins, flux_values, flux_errors)\n", "\n", "# Generate 10 random samples\n", "samples = spectrum.sample(samples=10)\n", "print(samples)" ] }, { "cell_type": "markdown", "id": "67ba541b", "metadata": {}, "source": [ "## Saving and Loading `Spectrum` objects\n", "\n", "`Spectrum` objects can be saved to a file using the {py:func}`write_csv ` method, which saves the spectrum data in a simple text format. The method takes one parameter: `filename`, which specifies the name of the file to save the spectrum to. For example:" ] }, { "cell_type": "code", "execution_count": 8, "id": "608ba722", "metadata": {}, "outputs": [], "source": [ "energy_bins = np.array([0, 2, 4, 6, 8])\n", "flux_values = np.array([10, 20, 15, 5])\n", "flux_errors = np.array([1, 2, 1.5, 0.5])\n", "\n", "spectrum = Spectrum(energy_bins, flux_values, flux_errors)\n", "\n", "# Save the spectrum to a file\n", "spectrum.write_csv(\"spectrum.csv\")" ] }, { "cell_type": "markdown", "id": "ae7b323f", "metadata": {}, "source": [ "`Spectrum` objects can be loaded from a file using the {py:func}`from_file ` class method, which reads the spectrum data from a file in the same format as the `write_csv` method. The method takes one parameter: `filename`, which specifies the name of the file to load the spectrum from. For example:" ] }, { "cell_type": "code", "execution_count": 9, "id": "e902623f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "# Load a spectrum from a file\n", "loaded_spectrum = Spectrum.from_file(\"spectrum.csv\")\n", "print(loaded_spectrum)" ] }, { "cell_type": "code", "execution_count": null, "id": "09d09fcd", "metadata": {}, "outputs": [], "source": [ "import os\n", "os.remove(\"spectrum.csv\")" ] } ], "metadata": { "kernelspec": { "display_name": "snf-simulations (3.10.12)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 5 }