create a own CAL product "level"

Author: nicHoch

.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

stixcore/calibration/ecc_post_fit.py

@@ -1,27 +1,27 @@
 """
-    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
+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
 
 """
 
@@ -49,43 +49,42 @@ def open_fits_tables(fits_path):
     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))]
+    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
+    Pix = Pix + (DETECTOR_ID - 1) * 12
     obs = [0] * Nbin
     erg_c = data.ERG_center  # data.field(0)
-    obs = data.field(3+Pix)
+    obs = data.field(3 + Pix)
     nbin = Nbin
-    Nbin = [Nbin/NRebin]
+    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)
+    yerr = NSigma * np.sqrt(obs)
     return erg_c, obs, yerr
 
 
 def line(x, slope, intercept):
     """a line"""
-    return (slope*x + intercept)
+    return slope * x + intercept
 
 
 def poly(x, degree, slope, intercept):
     """a line"""
-    return (degree*x*x + slope*x + intercept)
+    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)))
+    return amp * np.exp(-((x - cen) ** 2) / (2 * wid**2))
 
 
 def Fit_Ba_Lines_Robust(erg_c, obs, PLOT_VERBOSE=1, LogScale=1):
@@ -97,66 +96,68 @@ def Fit_Ba_Lines_Robust(erg_c, obs, PLOT_VERBOSE=1, LogScale=1):
         + 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.)).nonzero()
-    y = (obs[pipo])
-    x = (erg_c[pipo])
-    x = np.array(x, dtype='float64')
-    y = np.array(y, dtype='float64')
+    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.) < 1.)):
-        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']
+    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 FOLLWING SECTION IS QUIET ROBUST BUT MIGHT OVERESTIMATE THE ENERGY RESOLUTION AT 32 keV
-    # THIS ALLOWS TO FORCE THE BASE LINE TO ADJSUT in 40-45 keV range while a simple Gaussian
+    # 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')
+    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., slope=0, intercept=0.)
+    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.) < 1.) & (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']
+    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
+    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
+    DETs = np.arange(32) + 1
     LARGEs = [0, 1, 2, 3, 4, 5, 6, 7]
     SMALLs = [8, 9, 10, 11]
-    PIXELs = LARGEs+SMALLs
+    PIXELs = LARGEs + SMALLs
 
     # Prep Pandas DataFrame
     df = pd.DataFrame()  # dict will be define later during concactenation with appropriated keys
@@ -166,63 +167,63 @@ def ecc_post_fit(erg_file, para_file, livetime):
     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]}
+            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 concatening the DataFrame row by row, this slows down dramatically progressively
+            # 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)
+    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
+    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)
+    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
+    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]
+    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]
+    df.loc[idx, "Gain_Cor"] = df["Gain_ECC"][idx]
+    df.loc[idx, "Offset_Cor"] = df["Offset_ECC"][idx]
 
     return df, idx

stixcore/calibration/elut_manager.py

@@ -10,7 +10,7 @@
 from stixcore.util.logging import get_logger
 from stixcore.util.singleton import Singleton
 
-__all__ = ['ELUTManager']
+__all__ = ["ELUTManager"]
 
 ELUT_DATA_DIR = Path(__file__).parent.parent / "config" / "data" / "common" / "elut"
 
@@ -42,10 +42,10 @@ def _load_index(self):
             logger.info(f"Loaded {len(self.elut_index)} ELUT entries from index")
 
         except FileNotFoundError:
-            logger.warning(f'No ELUT index found at: {self.elut_index_file}')
+            logger.warning(f"No ELUT index found at: {self.elut_index_file}")
             self.elut_index = []
         except Exception as e:
-            logger.error(f'Error loading ELUT index: {e}')
+            logger.error(f"Error loading ELUT index: {e}")
             self.elut_index = []
 
     def _find_elut_file(self, date):
@@ -173,7 +173,7 @@ def cache_size(self):
 
 
 # Create singleton instance
-if 'pytest' in sys.modules:
-    ELUTManager.instance = ELUTManager(test_data.elut if hasattr(test_data, 'elut') else None)
+if "pytest" in sys.modules:
+    ELUTManager.instance = ELUTManager(test_data.elut if hasattr(test_data, "elut") else None)
 else:
     ELUTManager.instance = ELUTManager()

stixcore/calibration/tests/test_elut_manager.py

@@ -34,7 +34,7 @@ def test_elut_index_loading(self):
         """Test that ELUT index is loaded correctly."""
         manager = ELUTManager()
 
-        assert hasattr(manager, 'elut_index')
+        assert hasattr(manager, "elut_index")
         assert len(manager.elut_index) > 0
 
         # Check that index contains expected structure
@@ -64,7 +64,7 @@ def test_find_elut_file_known_date(self):
 
         assert elut_file is not None
         assert isinstance(elut_file, Path)
-        assert elut_file.suffix == '.csv'
+        assert elut_file.suffix == ".csv"
 
         # Verify the file exists
         elut_path = manager.data_root / elut_file
@@ -178,12 +178,13 @@ def test_error_handling_nonexistent_file(self):
         manager = ELUTManager()
 
         from stixpy.calibration.detector import get_sci_channels
+
         sci_channels = get_sci_channels(datetime(2022, 8, 1))
 
         with pytest.raises(FileNotFoundError, match="ELUT file not found"):
             manager.read_elut("nonexistent_file.csv", sci_channels)
 
     def test_instance_attribute_exists(self):
         """Test that singleton instance is accessible."""
-        assert hasattr(ELUTManager, 'instance')
+        assert hasattr(ELUTManager, "instance")
         assert isinstance(ELUTManager.instance, ELUTManager)

stixcore/ecc/manager.py

@@ -208,6 +208,7 @@ def create_context(self, date=None):
                 Ignore_Min_Gain_Det_Pix_List=json.loads(
                     config.get("DEFAULT", "Ignore_Min_Gain_Det_Pix_List", fallback="[]")
                 ),
+                Name=configuration_name,
             )
 
             logger.info(f"Read config from in: {temp_dir / 'post_ecc.ini'}")