Adds support for amptek detector files

docs/conf.py

@@ -50,7 +50,7 @@
 plot_include_source = True
 
 # Set automodapi to generate files inside the generated directory
-# automodapi_toctreedirnm = "_build/html/api"
+automodapi_toctreedirnm = "_build/api"
 numpydoc_show_class_members = False
 # generate autosummary even if no references
 autosummary_generate = True
@@ -91,18 +91,24 @@
     ),
     "numpy": (
         "https://docs.scipy.org/doc/numpy/",
-        (None, "http://data.astropy.org/intersphinx/numpy.inv"),
+        (None, "https://numpy.org/doc/stable/objects.inv"),
     ),
     "scipy": (
         "https://docs.scipy.org/doc/scipy/reference/",
-        (None, "http://data.astropy.org/intersphinx/scipy.inv"),
+        (None, "https://docs.scipy.org/doc/scipy/objects.inv"),
     ),
     "matplotlib": (
         "https://matplotlib.org/",
-        (None, "http://data.astropy.org/intersphinx/matplotlib.inv"),
+        (None, "https://matplotlib.org/stable/objects.inv"),
+    ),
+    "astropy": (
+        "http://docs.astropy.org/en/stable/",
+        "https://docs.astropy.org/en/stable/objects.inv",
+    ),
+    "sunpy": (
+        "https://docs.sunpy.org/en/stable/",
+        "https://docs.sunpy.org/en/stable/objects.inv",
     ),
-    "astropy": ("http://docs.astropy.org/en/stable/", None),
-    "sunpy": ("https://docs.sunpy.org/en/stable/", None),
 }
 
 # -- Options for HTML output -------------------------------------------------

docs/user-guide/amptek.rst

@@ -0,0 +1,122 @@
+.. _amptek:
+
+*************************
+Amptek Reference Detector
+*************************
+
+Overview
+========
+
+For many ground-based calibration tasks, an amptek detector was used.
+To facilitate the analysis of these data, this package includes some support for opening and analyzing these calibration files.
+
+A test file is included of a measurement of a Mini-X xray tube with a gold target.
+
+.. plot::
+
+    An amptek Mini-X xray tube spectrum.
+
+    >>> import numpy as np
+    >>> import padre_meddea
+    >>> from padre_meddea.io.amptek import read_mca
+    >>> mca_file = padre_meddea._test_files_directory / "minix_20kV_15uA_sdd.mca"
+    >>> spec = read_mca(mca_file)
+    >>> spec.plot()
+
+
+The output of `~read_mca` is a Spectrum object.
+It contains a bunch of information inside its metadata which is a dictionary.
+For example, if the mca file contains calibration data, a calibration function is provided.
+And if regions of interest were defined those are also provided.
+
+.. plot::
+
+    An amptek Mini-X xray tube spectrum.
+
+    >>> import numpy as np
+    >>> import matplotlib.pyplot as plt
+    >>> import astropy.units as u
+    >>> import padre_meddea
+    >>> from padre_meddea.io.amptek import read_mca
+    >>> mca_file = padre_meddea._test_files_directory / "minix_20kV_15uA_sdd.mca"
+    >>> spec = read_mca(mca_file)
+    >>> f = spec.meta['calib']
+    >>> energy_ax = f(spec.spectral_axis.value)
+    >>> fig, ax = plt.subplots(layout="constrained")
+    >>> ax.plot(energy_ax, spec.flux)
+    >>> ax.set_xlabel("energy [keV]")
+    >>> ax.set_ylabel("Counts")
+    >>> ax.set_yscale("log")
+    >>> for this_roi in spec.meta['roi']:
+    >>>     ax.axvspan(f(this_roi.lower.value), f(this_roi.upper.value), alpha=0.2)
+
+The regions of interest are provided as SpectralRegion objects.
+They can be used to easily extract a sub spectrum to, for example, fit the gaussian to the emission line.
+
+.. plot::
+
+    An amptek Mini-X xray tube spectrum.
+
+    >>> import matplotlib.pyplot as plt
+    >>> import padre_meddea
+    >>> from padre_meddea.io.amptek import read_mca
+    >>> from specutils.manipulation import extract_region
+    >>> from specutils.fitting import estimate_line_parameters, fit_lines
+    >>> from astropy.modeling import models
+    >>> mca_file = padre_meddea._test_files_directory / "minix_20kV_15uA_sdd.mca"
+    >>> spec = read_mca(mca_file)
+    >>> this_roi = spec.meta['roi'][0]
+    >>> sub_spec = extract_region(spec, this_roi)
+    >>> params = estimate_line_parameters(sub_spec, models.Gaussian1D())
+    >>> g_init = models.Gaussian1D(amplitude=params.amplitude, mean=params.mean, stddev=params.stddev)
+    >>> g_fit = fit_lines(sub_spec, g_init)
+    >>> fig, ax = plt.subplots(layout="constrained")
+    >>> ax.plot(sub_spec.spectral_axis, sub_spec.flux)
+    >>> ax.plot(sub_spec.spectral_axis, g_fit(sub_spec.spectral_axis), label=f"{g_fit}")
+    >>> plt.legend(bbox_to_anchor =(0.65, 1.25))
+
+To identify these lines we can make use of the roentgen Python package.
+This x-ray tube makes use of a gold target so most of these lines are from gold.
+The gold lines can be found with the following code:
+
+.. code-block:: python
+
+    >>> from roentgen.lines import get_lines
+    >>> au_lines = get_lines(1 * u.keV, 20 * u.keV, "au")
+    >>> au_lines
+    <QTable length=7>
+    energy   z   symbol transition intensity
+    eV                                    
+    float64 int64  str2     str6      int64  
+    ------- ----- ------ ---------- ---------
+    2122.9    79     Au        Mα1       100
+    8493.9    79     Au         Ll         5
+    9628.0    79     Au        Lα2        11
+    9713.3    79     Au        Lα1       100
+    11442.3    79     Au        Lβ1        67
+    11584.7    79     Au        Lβ2        23
+    13381.7    79     Au        Lγ1        13
+
+
+In this case we only consider the brightest lines.
+Let's plot them over our spectrum which has been roughly calibrated.
+
+.. plot::
+
+    >>> import matplotlib.pyplot as plt
+    >>> from roentgen.lines import get_lines
+    >>> import astropy.units as u
+    >>> import padre_meddea
+    >>> from padre_meddea.io.amptek import read_mca
+    >>> mca_file = padre_meddea._test_files_directory / "minix_20kV_15uA_sdd.mca"
+    >>> spec = read_mca(mca_file)
+    >>> au_lines = get_lines(1 * u.keV, 20 * u.keV, "au")
+    >>> f = spec.meta['calib']
+    >>> energy_ax = f(spec.spectral_axis.value)
+    >>> fig, ax = plt.subplots(layout="constrained")
+    >>> ax.plot(energy_ax, spec.flux)
+    >>> for this_line in au_lines:
+    >>>     this_label = f"{this_line['energy'].to('keV'):0.2f} Au {this_line['transition']} {this_line['intensity']}"
+    >>>     ax.axvline(this_line['energy'].to('keV').value, label=this_label, color='red')
+    >>> plt.legend()
+

docs/user-guide/index.rst

@@ -14,6 +14,8 @@ For more details checkout the :ref:`reference`.
    data
    raw
    level0
-   spectrum
+   simul_spectrum
+   amptek
+   raw_spectrum_data
    customization
    logger

padre_meddea/calibration/__init__.py

@@ -1 +0,0 @@
-from .calibration import *

padre_meddea/calibration/spectrum.py

@@ -1,23 +1,85 @@
 """
-This module provides tools to analyze and manipulate spectrum data.
+This module provides tools to analyze and manipulate spectra data.
 """
 
 import numpy as np
+from numpy.polynomial import Polynomial
 
+import astropy.units as u
+from astropy.modeling import models
+
+from specutils import Spectrum1D, SpectralRegion
+from specutils.manipulation import extract_region
+from specutils.fitting import estimate_line_parameters, fit_lines
+
+from astropy.nddata import StdDevUncertainty
+from astropy.table import Table
+from astropy.timeseries import TimeSeries
 from astropy.timeseries import aggregate_downsample
 from specutils import Spectrum1D
 
 
-def phlist_to_lc(event_list, int_time):
-    """Convert an event list to a light curve.
+__all__ = ["get_calib_energy_func", "elist_to_spectrum", "elist_to_lc"]
+
+
+def get_calib_energy_func(spectrum: Spectrum1D, line_energies, rois, deg=1):
+    """Given a spectrum with known emission lines, return a function to transform between spectral axis units to energy units.
+
+    The method used here is to fit a Gaussian model to each region of interest.
+    A linear fit is then performed between the means of each Gaussian fit and the known energies.
+
+    Parameters
+    ----------
+    spectrum: Spectrum1D
+        A spectrum with known emission lines.
+    line_energies: u.Quantity
+        The energy of each line.
+    rois: ndarray
+        A list of regions of interest for each line.
+
+    Returns
+    -------
+    func:
+        A function to convert between channel space to energy space.
+
+    Examples
+    --------
+    """
+    if len(rois) != len(line_energies):
+        raise ValueError(
+            f"Number of line energies {len(line_energies)} does not match number of rois ({len(rois)})."
+        )
+    # if len(line_energies) < (deg + 1):
+    #    raise ValueError(f"Not enough values to perform fit with degree {deg}.")
+    fit_means = []
+    spectral_axis_units = spectrum.spectral_axis.unit
+    for this_roi in rois:
+        this_region = SpectralRegion(
+            this_roi[0] * spectral_axis_units, this_roi[1] * spectral_axis_units
+        )
+        sub_spectrum = extract_region(spectrum, this_region)
+        params = estimate_line_parameters(sub_spectrum, models.Gaussian1D())
+        g_init = models.Gaussian1D(
+            amplitude=params.amplitude, mean=params.mean, stddev=params.stddev
+        )
+        g_fit = fit_lines(sub_spectrum, g_init)
+        fit_means.append(g_fit.mean.value)
+    result = Polynomial.fit(fit_means, line_energies, deg=deg)  # fit a linear model
+    result_all = Polynomial.fit(fit_means, line_energies, deg=deg, full=True)
+    print(result_all)
+    return result.convert()  # return the function
+
+
+def elist_to_lc(event_list, int_time):
+    """Convert an event list to a light curve timeseries
 
     Parameters
     ----------
     event_list: An event list
 
     Returns
     -------
-    event_list
+    ts: TimeSeries
 
     Examples
     --------
@@ -26,13 +88,25 @@ def phlist_to_lc(event_list, int_time):
     return ts
 
 
-def phlist_to_spec(event_list, bins=None):
+def elist_to_spectrum(event_list: Table, bins=None):
+    """Convert an event list to a spectrum object.
+
+    Parameters
+    ----------
+    event_list: Table
+        An event list
+
+
+    Returns
+    -------
+    spectrum: Spectrum1D
+    """
     if bins is None:
         bins = np.arange(0, 2**12 - 1)
-    data, bins = np.histogram(event_list["energy"], bins=bins)
+    data, bins = np.histogram(event_list["atod"], bins=bins)
     result = Spectrum1D(
         flux=u.Quantity(data, "count"),
         spectral_axis=u.Quantity(bins, "pix"),
-        uncertainty=np.sqrt(data),
+        uncertainty=StdDevUncertainty(np.sqrt(data) * u.count),
     )
     return result