CAL spectra : daily energy calibration

.github/workflows/end2end_run.yml

@@ -43,7 +43,7 @@ jobs:
         python -m pip install -e .[dev]
         if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
     - name: show environment
-      run: python -m pip freeze 
+      run: python -m pip freeze
     - name: end2end test with pytest
       run: |
         pytest -p no:warnings --doctest-rst -m "end2end" --pyargs stixcore --remote-data=any

docs/pipelineconfiguration.rst

@@ -175,37 +175,6 @@ without provided arguments the default values from ``stixcore.ini`` are used
     0 5 * * * cd /home/stixcore/STIXCore/ && /home/stixcore/STIXCore/venv/bin/python /home/stixcore/STIXCore/stixcore/processing/publish.py --update_rid_lut
 
 
-Run the pipeline monitor
-************************
-
-The event based pipeline (observing incoming telemetry files) gets stuck from time to time. There is a process observing the number of open to process files. If the number of open files is constantly increasing over a longer period a notification mail is send out:
-
-.. code-block::
-
-    # run pipeline monitor task to check for pipeline not stuck
-    0 */3 * * * cd /home/stixcore/STIXCore/ && /home/stixcore/STIXCore/venv/bin/python /home/stixcore/STIXCore/stixcore/processing/pipeline_monitor.py -s /home/stixcore/monitor_status.json
-
-
-In case of a pipeline stuck restart the event based processing pipeline.
-
-.. code-block::
-
-    # stop the system.d process
-    sudo systemctl stop stix-pipeline.service
-
-    # wait 20sec so that all open sockets also gets closed
-    # start the process again
-
-    sudo systemctl start stix-pipeline.service
-
-In order to process all tm files that have not been processed so fare the config parameter start_with_unprocessed should be set to true:
-
-.. code-block::
-
-    [Pipeline]
-    start_with_unprocessed = True
-
-
 Run the 'daily' pipeline
 ************************
 

pyproject.toml

@@ -30,6 +30,7 @@ dependencies = [
     "numpy==2.3.4",
     "sunpy~=6.0.0",
     "watchdog~=6.0",
+    "lmfit==1.3.4",
     "stixdcpy==3.0",
     "stixpy>=0.2.0rc1,<0.3"
 ]

stixcore/__init__.py

@@ -3,4 +3,9 @@
 from .version import version as __version__
 from .version_conf import __version_conf__
 
+try:
+    from .version_conf import __version_conf__
+except ImportError:
+    __version_conf__ = "unknown"
+
 logger = get_logger(__name__)

stixcore/calibration/ecc_post_fit.py

@@ -0,0 +1,229 @@
+"""
+Author: O. Limousin, CEA
+Date: Oct 23, 2024
+This script to:
+
+#%%%%%%%%%%%%%% FOR fit of ECC_ULTRA_FINE %%%%%%%%%%%%%%
+../00_CALIBRATION_MONITORING_ULTRA_FINE/02_MONITOR_ECC_SPECTRA ==> XX_erg.fits
+
+../00_CALIBRATION_MONITORING_ULTRA_FINE/01_MONITOR_ECC_PARA/ ==> ECC_para_XXXX.fits
+
+    + read the log (.xlsx file)
+    + open the ECC calibrated spectra XX_erg.fits
+    + fit 31 and 81 keV lines (Fit_Ba_robust) (fit right hand side, and baseline)
+    + register the results in a Pandas Dataframe
+        + date, Run Number, Fit Goodness flags,DET, PIX,
+            P31, err_P31, dE31, err_dE31, H31
+            P81, err_P81, dE81, err_dE81, H81
+    + compute gain/offset corrections
+    + fill the dataframe with new values
+        + include errors
+        + include a Flag (Flag31 and Flag81) which is True if fit was OK
+            if was not OK, the ECC value in NOT corrected
+    + store the results in pkl
+    + generate update ECC_para files to recalibrate uncalibrated files
+
+"""
+
+import numpy as np
+import pandas as pd
+from lmfit import Model
+
+from astropy.io import fits
+
+from stixcore.util.logging import get_logger
+
+logger = get_logger(__name__)
+
+
+def open_fits_tables(fits_path):
+    # Get the data from .fits
+    data = []
+    with fits.open(fits_path, memmap=False) as hdul:
+        header = [hdul[0].header]
+        for i in range(1, len(hdul)):
+            header.append(hdul[i].header)
+            data.append(hdul[i].data)
+
+    # Get the fields of each data table and sort the data in lists
+    data_names = [data_i.columns.names for data_i in data]
+    data_format = [data_i.columns.formats for data_i in data]
+    data_unit = [data_i.columns.units for data_i in data]
+    data_list = [[list(data[i][data_names[i][j]]) for j in range(len(data_names[i]))] for i in range(len(data))]
+    return header, data, data_list, data_names, data_format, data_unit
+
+
+def Read_fits_STIX_One_Pixel(data, PIX=0, DETECTOR_ID=1, Nbin=2024, NRebin=1, NSigma=1):
+    data = data[0]
+    Pix = PIX
+    Pix = Pix + (DETECTOR_ID - 1) * 12
+    obs = [0] * Nbin
+    erg_c = data.ERG_center  # data.field(0)
+    obs = data.field(3 + Pix)
+    nbin = Nbin
+    Nbin = [Nbin / NRebin]
+    # TODO: check if we really need to rebin with congrid here
+    # erg_c = congrid(erg_c[:nbin], Nbin)
+    # obs = congrid(obs[:nbin], Nbin)
+    erg_c = erg_c[:nbin]
+    obs = obs[:nbin]
+
+    yerr = NSigma * np.sqrt(obs)
+    return erg_c, obs, yerr
+
+
+def line(x, slope, intercept):
+    """a line"""
+    return slope * x + intercept
+
+
+def poly(x, degree, slope, intercept):
+    """a line"""
+    return degree * x * x + slope * x + intercept
+
+
+def gaussian(x, amp, cen, wid):
+    #    """1-d gaussian: gaussian(x, amp, cen, wid)"""
+    return amp * np.exp(-((x - cen) ** 2) / (2 * wid**2))
+
+
+def Fit_Ba_Lines_Robust(erg_c, obs, PLOT_VERBOSE=1, LogScale=1):
+    """
+    OL, oct 22, 2024
+    Robust fit procedure to adjust peak position post ECC
+    The idea is the exact same as previous function but to return
+        + 0's in case the fit fails
+        + flag to say if the fit was successful or not
+    """
+    # Select Energy range for 81 keV Ba-133 line
+    pipo = ((erg_c > 80.5) & (erg_c < 90.0)).nonzero()
+    y = obs[pipo]
+    x = erg_c[pipo]
+    x = np.array(x, dtype="float64")
+    y = np.array(y, dtype="float64")
+    mod = Model(gaussian) + Model(line)
+    pars = mod.make_params(amp=10, cen=81, wid=0.5, slope=0, intercept=0)
+    result = mod.fit(y, pars, x=x)
+
+    if (
+        (result.params["wid"].stderr is not None)
+        & (result.params["cen"].stderr is not None)
+        & (np.abs(result.best_values["cen"] - 81.0) < 1.0)
+    ):
+        dE81 = result.best_values["wid"] * 2.35
+        err_dE81 = result.params["wid"].stderr * 2.35
+        P81 = result.best_values["cen"]
+        err_P81 = result.params["cen"].stderr
+        H81 = result.best_values["amp"]
+        Goodness_Flag_81 = True
+    else:
+        dE81, err_dE81, P81, err_P81, H81 = 0, 0, 0, 0, 0
+        Goodness_Flag_81 = False
+
+    # Select Energy range for 81 keV Ba-133 line
+    # THE FOLLOWING SECTION IS QUIET ROBUST BUT MIGHT OVERESTIMATE THE ENERGY RESOLUTION AT 32 keV
+    # THIS ALLOWS TO FORCE THE BASE LINE TO ADJUST in 40-45 keV range while a simple Gaussian
+    # is used to fit the right hand side of the 32 keV Line
+    pipo = (((erg_c > 30.2) & (erg_c < 33.0)) | ((erg_c > 40) & (erg_c < 45))).nonzero()
+    y = obs[pipo]
+    x = erg_c[pipo]
+    x = np.array(x, dtype="float64")
+    y = np.array(y, dtype="float64")
+
+    mod = Model(gaussian, prefix="g1_") + Model(poly)
+    pars = mod.make_params(g1_amp=10, g1_cen=30.6, g1_wid=0.4, degree=0.0, slope=0, intercept=0.0)
+
+    result = mod.fit(y, pars, x=x)
+
+    if (
+        (result.params["g1_wid"].stderr is not None)
+        & (result.params["g1_cen"].stderr is not None)
+        & ((np.abs(result.best_values["g1_cen"] - 31.0) < 1.0) & (Goodness_Flag_81))
+    ):
+        dE31 = result.best_values["g1_wid"] * 2.35
+        err_dE31 = result.params["g1_wid"].stderr * 2.35
+        P31 = result.best_values["g1_cen"]
+        err_P31 = result.params["g1_cen"].stderr
+        H31 = result.best_values["g1_amp"]
+        Goodness_Flag_31 = True
+    else:
+        dE31, err_dE31, P31, err_P31, H31 = 0, 0, 0, 0, 0
+        Goodness_Flag_31 = False
+
+    return P31, P81, dE31, dE81, err_P31, err_P81, err_dE31, err_dE81, H31, H81, Goodness_Flag_31, Goodness_Flag_81
+
+
+def ecc_post_fit(erg_file, para_file, livetime):
+    DETs = np.arange(32) + 1
+    LARGEs = [0, 1, 2, 3, 4, 5, 6, 7]
+    SMALLs = [8, 9, 10, 11]
+    PIXELs = LARGEs + SMALLs
+
+    # Prep Pandas DataFrame
+    df = pd.DataFrame()  # dict will be define later during concactenation with appropriated keys
+    accumulator = []
+
+    # Proceed to fit individually each pixel spectrum in the list of files
+    data = open_fits_tables(erg_file)[1]  # only need data
+    for DET in DETs:
+        for PIX in PIXELs:  # or in  LARGEs, SMALLs
+            erg_c, obs, yerr = Read_fits_STIX_One_Pixel(data, PIX=PIX, DETECTOR_ID=DET, Nbin=2024, NRebin=1, NSigma=1)
+
+            P31, P81, dE31, dE81, err_P31, err_P81, err_dE31, err_dE81, H31, H81, Goodness_Flag31, Goodness_Flag81 = (
+                Fit_Ba_Lines_Robust(erg_c, obs / livetime, PLOT_VERBOSE=0, LogScale=0)
+            )
+
+            dict = {
+                "DET": [DET],
+                "PIX": [PIX],
+                "P31": [P31],
+                "err_P31": [err_P31],
+                "dE31": [dE31],
+                "err_dE31": [err_dE31],
+                "Flag31": [Goodness_Flag31],
+                "P81": [P81],
+                "err_P81": [err_P81],
+                "dE81": [dE81],
+                "err_dE81": [err_dE81],
+                "Flag81": [Goodness_Flag81],
+            }
+            logger.debug(dict)
+
+            # NB: this is faster to append a list of dict and create DataFrame at the end
+            # than concatenating the DataFrame row by row, this slows down dramatically progressively
+            accumulator.append(pd.DataFrame(dict))
+
+    df = pd.concat(accumulator)
+    df = df.reset_index(drop=True)
+
+    # 3- gain and offset correction factors of ECC pre-calibrated data
+    G_prime = (df["P81"] - df["P31"]) / (80.9979 - (30.6254 * 33.8 + 30.9731 * 62.4) / (62.4 + 33.8))
+    O_prime = df["P31"] - G_prime * (30.6254 * 33.8 + 30.9731 * 62.4) / (62.4 + 33.8)
+
+    # 4- add correction factors to the DataFrame
+    df["Gain_Prime"] = G_prime
+    df["Offset_Prime"] = O_prime
+
+    # check, Run number and pixel number prior to assign Gain and Offset for further correction
+    data_para = open_fits_tables(para_file)[1]  # only need data
+
+    df["Gain_ECC"] = np.float32(data_para[0].gain)
+    df["Offset_ECC"] = np.float32(data_para[0].off)
+    df["goc"] = np.float32(data_para[0].goc)
+
+    # 7 - Compute corrected Gain and Offset and fill df
+    # 7.2 - Now assign the corrected Gain and Offset values
+    # except when ECC works better
+    df["Gain_Cor"] = 0.0
+    df["Offset_Cor"] = 0.0
+
+    # apply correction to gain and offset when fit is ok
+    idx = df["Flag31"] & df["Flag81"]
+    df.loc[idx, "Gain_Cor"] = df["Gain_ECC"][idx] * df["Gain_Prime"][idx]
+    df.loc[idx, "Offset_Cor"] = df["Gain_ECC"][idx] * df["Offset_Prime"][idx] + df["Offset_ECC"][idx]
+    # otherwise keep uncorrected ECC Values
+    idx = ~idx
+    df.loc[idx, "Gain_Cor"] = df["Gain_ECC"][idx]
+    df.loc[idx, "Offset_Cor"] = df["Offset_ECC"][idx]
+
+    return df, idx