tickets/PIPE2D-680: Add a config to change the order of merging and 1-D sky subtraction

python/pfs/drp/stella/mergeArms.py

@@ -32,6 +32,8 @@ class MergeArmsConfig(Config):
     """Configuration for MergeArmsTask"""
     wavelength = ConfigField(dtype=WavelengthSamplingConfig, doc="Wavelength configuration")
     doSubtractSky1d = Field(dtype=bool, default=True, doc="Do 1D sky subtraction?")
+    doSubtractSky1dBeforeMerge = Field(dtype=bool, default=False,
+                                       doc="Do 1D sky subtraction before merging arms?")
     subtractSky1d = ConfigurableField(target=SubtractSky1dTask, doc="1d sky subtraction")
     doBarycentricCorr = Field(dtype=bool, default=True, doc="Do barycentric correction?")
     mask = ListField(dtype=str, default=["NO_DATA", "CR", "INTRP", "SAT"],
@@ -93,12 +95,25 @@ def runDataRef(self, expSpecRefList):
         lsf = [[None for dataRef in specRefList] for specRefList in expSpecRefList]
         pfsConfig = expSpecRefList[0][0].get("pfsConfig")
         if self.config.doSubtractSky1d:
-            sky1d = self.subtractSky1d.run(sum(spectra, []), pfsConfig, sum(lsf, []))
-            expSpecRefList[0][0].put(sky1d, "sky1d")
+            if self.config.doSubtractSky1dBeforeMerge:
+                sky1d = self.subtractSky1d.run(sum(spectra, []), pfsConfig, sum(lsf, []))
+                expSpecRefList[0][0].put(sky1d, "sky1d")
+            else:
+                sky1d = None
 
         spectrographs = [self.mergeSpectra(ss) for ss in spectra]  # Merge in wavelength
         merged = PfsMerged.fromMerge(spectrographs, metadata=getPfsVersions())  # Merge across spectrographs
+
+        if self.config.doSubtractSky1d:
+            if sky1d is None:
+                assert not self.config.doSubtractSky1dBeforeMerge
+
+                lsf = [None for m in range(len(merged))] # total hack!
+                sky1d = self.subtractSky1d.estimateSkyFromMerged(merged, pfsConfig, lsf)
+                expSpecRefList[0][0].put(sky1d, "sky1d")
+
         expSpecRefList[0][0].put(merged, "pfsMerged")
+
         return Struct(spectra=merged, pfsConfig=pfsConfig)
 
     def mergeSpectra(self, spectraList):

python/pfs/drp/stella/subtractSky1d.py

@@ -165,6 +165,39 @@ def subtractSkySpectrum(self, spectrum, lsf, fiberId, pfsConfig, sky1d):
         spectrum.mask[np.array(sky.masks)] |= bitmask
         spectrum.covariance[0] += sky.variance
 
+    def estimateSkyFromMerged(self, merged, pfsConfig, lsfList):
+        """Measure and subtract the sky from the merged 1D spectra
+
+        Parameters
+        ----------
+        merged : `pfs.datamodel.PfsMerged`
+            List of merged spectra from which to subtract the sky,
+            including sky spectra
+        pfsConfig : `pfs.datamodel.PfsConfig`
+            Configuration of the top-end, for identifying sky fibers.
+        lsfList : iterable of LSF (type TBD)
+            List of line-spread functions.
+
+        Returns
+        -------
+        sky1d : `pfs.drp.stella.FocalPlaneFunction`
+            1D sky model.
+        """
+        if not np.any(pfsConfig.targetType == TargetType.SKY):
+            raise RuntimeError("No sky fibers found")
+
+        if self.debugInfo.plotSkyFluxes:
+            self.plotSkyFibers([merged], pfsConfig, "Sky flux")
+
+        sky1d = self.measureSky([merged], pfsConfig, [lsfList])
+
+        self.subtractSkySpectra(merged, lsfList, pfsConfig, sky1d)
+
+        if self.debugInfo.plotSkyResiduals:
+            self.plotSkyFibers([merged], pfsConfig, "Sky residuals")
+
+        return sky1d
+
     def plotSkyFibers(self, spectraList, pfsConfig, title):
         """Plot spectra from sky fibers