From 0354ea2a67595651830246570cdd30f0afd296fd Mon Sep 17 00:00:00 2001
From: Masato Onodera <masato.onodera@gmail.com>
Date: Tue, 3 Dec 2024 18:42:36 -1000
Subject: [PATCH] Update data loader to the latest PFS datamodel (as of Nov
 2024)

In the last 2yrs since the data loader for Subaru-pfsObject is
created, significant changes has been made for the datamodel.
Although it's still actively developed, I'd like to initiate
updating the data loader.

Unfortunately, at this moment, the actual pfsObject file is still
proprietary, but I would like to share an example file once it's
publicly available (openuse operation will start in early 2025).

Reference:
- PFS datamodel: https://github.com/Subaru-PFS/datamodel/blob/master/datamodel.txt
---
 .../io/default_loaders/subaru_pfs_spec.py     | 93 ++++++++++++-------
 1 file changed, 58 insertions(+), 35 deletions(-)

diff --git a/specutils/io/default_loaders/subaru_pfs_spec.py b/specutils/io/default_loaders/subaru_pfs_spec.py
index 8e142426d..b138f94f5 100644
--- a/specutils/io/default_loaders/subaru_pfs_spec.py
+++ b/specutils/io/default_loaders/subaru_pfs_spec.py
@@ -3,27 +3,26 @@
 
 https://github.com/Subaru-PFS/datamodel/blob/master/datamodel.txt
 """
+
 import os
 import re
 
-from astropy.units import Unit
-from astropy.nddata import StdDevUncertainty
-
 import numpy as np
+from astropy.nddata import StdDevUncertainty
+from astropy.units import Unit
 
 from ...spectra import Spectrum1D
-from ..registers import data_loader
 from ..parsing_utils import _fits_identify_by_name, read_fileobj_or_hdulist
+from ..registers import data_loader
 
-__all__ = ['identify_pfs_spec', 'pfs_spec_loader']
+__all__ = ["identify_pfs_spec", "pfs_spec_loader"]
 
 # This RE matches the file name pattern defined in Subaru-PFS' datamodel.txt :
-# "pfsObject-%05d-%s-%3d-%08x-%02d-0x%08x.fits" % (tract, patch, catId, objId,
-#                                                  nVisit % 100, pfsVisitHash)
-_spec_pattern = re.compile(r'pfsObject-(?P<tract>\d{5})-(?P<patch>.{3})-'
-                           r'(?P<catId>\d{3})-(?P<objId>\d{8})-'
-                           r'(?P<nVisit>\d{2})-(?P<pfsVisitHash>0x\w{8})'
-                           r'\.fits')
+# "pfsObject-%05d-%05d-%s-%016x-%03d-0x%016x.fits" % (catId, tract, patch, objId,
+#                                                     nVisit % 1000, pfsVisitHash)
+_spec_pattern = re.compile(
+    r"pfsObject-(?P<catId>\d{5})-(?P<tract>\d{5})-(?P<patch>.{3})-(?P<objId>\d{16})-(?P<nVisit>\d{3})-(?P<pfsVisitHash>0x\w{16})\.fits"
+)
 
 
 def identify_pfs_spec(origin, *args, **kwargs):
@@ -35,7 +34,9 @@ def identify_pfs_spec(origin, *args, **kwargs):
 
 
 @data_loader(
-    label="Subaru-pfsObject", identifier=identify_pfs_spec, extensions=['fits'],
+    label="Subaru-pfsObject",
+    identifier=identify_pfs_spec,
+    extensions=["fits"],
     priority=10,
 )
 def pfs_spec_loader(file_obj, **kwargs):
@@ -64,28 +65,50 @@ def pfs_spec_loader(file_obj, **kwargs):
 
     with read_fileobj_or_hdulist(file_obj, **kwargs) as hdulist:
         header = hdulist[0].header
-        meta = {'header': header,
-                'tract': m['tract'],
-                'patch': m['patch'],
-                'catId': m['catId'],
-                'objId': m['objId'],
-                'nVisit': m['nVisit'],
-                'pfsVisitHash': m['pfsVisitHash']}
-
-        # spectrum is in HDU 2
-        data = hdulist[2].data['flux']
-        unit = Unit('nJy')
-
-        error = hdulist[2].data['fluxVariance']
+        meta = {
+            "header": header,
+            "tract": m["tract"],
+            "patch": m["patch"],
+            "catId": m["catId"],
+            "objId": m["objId"],
+            "nVisit": m["nVisit"],
+            "pfsVisitHash": m["pfsVisitHash"],
+        }
+
+        # PFS datamodel: https://github.com/Subaru-PFS/datamodel/blob/master/datamodel.txt
+        #
+        # HDU #8 FLUXTABLE   Binary table  [FITS BINARY TABLE] NOBS*NROW
+        # Columns for:
+        #   wavelength in units of nm (vacuum)          [64-bit FLOAT]
+        #   intensity in units of nJy                   [FLOAT]
+        #   intensity error in same units as intensity  [FLOAT]
+        #   mask                                        [32-bit INT]
+        #
+        # In the datamodel, the FLUXTABLE is the 8th HDU, but in fact it's in the 10th HDU with EXTNAME of FLUX_TABLE.
+        #
+        # NOTE: No backward compatibility is guaranteed for the FLUXTABLE format.
+        data = hdulist[10].data["flux"]
+        unit = Unit("nJy")
+
+        error = hdulist[10].data["error"]
         uncertainty = StdDevUncertainty(np.sqrt(error))
 
-        wave = hdulist[2].data['lambda']
-        wave_unit = Unit('nm')
-
-        mask = hdulist[2].data['mask'] != 0
-
-    return Spectrum1D(flux=data * unit,
-                      spectral_axis=wave * wave_unit,
-                      uncertainty=uncertainty,
-                      meta=meta,
-                      mask=mask)
+        wave = hdulist[10].data["wavelength"]
+        wave_unit = Unit("nm")
+
+        # NOTE: REFLINE mask is the 15th bit in the mask bits, and can be ignored.
+        # TODO: Mask condition needs to be refined once more information becomes available.
+        mask = ~np.logical_or(
+            # mask==0: good data
+            hdulist[10].data["mask"] == 0,
+            # mask==2**REFLINE (32768) can be ignored
+            hdulist[10].data["mask"] == 2 ** (hdulist[10].header["MP_REFLINE"]),
+        )
+
+    return Spectrum1D(
+        flux=data * unit,
+        spectral_axis=wave * wave_unit,
+        uncertainty=uncertainty,
+        meta=meta,
+        mask=mask,
+    )