Merge pull request #2606 from desihub/coadd_normalization

Rescale individual exposures before coaddition

Author: sbailey

py/desispec/coaddition.py

@@ -37,7 +37,7 @@ exposure_qa:
  # bad POSITIONING flag
  bad_positions_mask : "BROKENFIBER|MISSINGPOSITION|BADPOSITION|POORPOSITION"
  # bad QAFIBERSTATUS flag
- bad_qafstatus_mask : "UNASSIGNED|BROKENFIBER|MISSINGPOSITION|BADPOSITION|POORPOSITION|BADFIBER|BADTRACE|BADFLAT|BADARC|MANYBADCOL|MANYREJECTED|BADAMPB|BADAMPR|BADAMPZ|BADPETALPOS|BADPETALSKY|BADPETALSTDSTAR|BADPETALFLUXCAL|BADPETALSNR|BADREADNOISE|LOWEFFTIME"
+ bad_qafstatus_mask : "UNASSIGNED|BROKENFIBER|MISSINGPOSITION|BADPOSITION|POORPOSITION|BADFIBER|BADTRACE|BADFLAT|BADARC|MANYREJECTED|BADAMPB|BADAMPR|BADAMPZ|BADPETALPOS|BADPETALSKY|BADPETALSTDSTAR|BADPETALFLUXCAL|BADPETALSNR|BADREADNOISE|LOWEFFTIME"
  # bad QAFIBERSTATUS per-petal flag
  bad_petal_mask : "BADPETALPOS|BADPETALSKY|BADPETALSTDSTAR|BADPETALFLUXCAL|BADPETALSNR|BADREADNOISE"
 

py/desispec/data/qa/qa-params.yaml

@@ -55,7 +55,7 @@ def get_fiberbitmasked_frame_arrays(frame,bitmask=None,ivar_framemask=True,retur
         Example bitmask list::
 
             bitmask = [fmsk.BROKENFIBER,fmsk.UNASSIGNED,fmsk.BADFIBER,\
-                        fmsk.BADTRACE,fmsk.MANYBADCOL, fmsk.MANYREJECTED]
+                        fmsk.BADTRACE, fmsk.MANYREJECTED]
             bitmask = get_fiberbitmask_comparison_value(kind='fluxcalib', band='brz')
             bitmask = 'fluxcalib'
             bitmask = 4128780
@@ -154,8 +154,8 @@ def get_flat_fiberbitmask_val(band):
     Return mask of bad FIBERSTATUS bits for fiberflats
     i.e. fibers with these bits set have a bad fiberflat and cannot be used
     """
-    return (fmsk.BROKENFIBER | fmsk.BRIGHTNEIGHBOR | fmsk.BADFIBER | fmsk.BADTRACE | fmsk.BADARC | \
-            fmsk.MANYBADCOL | fmsk.MANYREJECTED )
+    return (fmsk.BROKENFIBER | fmsk.BRIGHTNEIGHBOR | fmsk.BADFIBER | \
+            fmsk.BADTRACE | fmsk.BADARC | fmsk.MANYREJECTED )
 
 def get_fluxcalib_fiberbitmask_val(band):
     """
@@ -182,10 +182,9 @@ def get_all_nonamp_fiberbitmask_val():
         NEARCHARGETRAP and VARIABLETHRU are also not included since
         they are ok for some types of processing but not others.
     """
-    return (fmsk.BROKENFIBER | fmsk.MISSINGPOSITION | \
-            fmsk.BADPOSITION | fmsk.BRIGHTNEIGHBOR  | \
-            fmsk.BADFIBER | fmsk.BADTRACE | fmsk.BADARC | fmsk.BADFLAT | \
-            fmsk.MANYBADCOL | fmsk.MANYREJECTED )
+    return (fmsk.BROKENFIBER | fmsk.MISSINGPOSITION | fmsk.BADPOSITION | \
+            fmsk.BRIGHTNEIGHBOR | fmsk.BADFIBER | fmsk.BADTRACE | \
+            fmsk.BADARC | fmsk.BADFLAT | fmsk.MANYREJECTED )
 
 def get_justamps_fiberbitmask():
     """

py/desispec/fiberbitmasking.py

@@ -169,6 +169,8 @@
                     ('IN_COADD_B',                 'bool',           '', 'Was this exposure included in B coadd',        'coadd'),
                     ('IN_COADD_R',                 'bool',           '', 'Was this exposure included in R coadd',        'coadd'),
                     ('IN_COADD_Z',                 'bool',           '', 'Was this exposure included in Z coadd',        'coadd'),
+                    ('COADD_NORM',                 'f4',             '', 'Flux normalization coefficient for coadd',     'coadd'),
+                    ('SIGMA_COADD_NORM',            'f4',             '', 'Error on COADD_NORM',                         'coadd'),
                     )
 
 # dict to allow unit tests etc that augment fibermaps to look up the

py/desispec/io/fibermap.py

@@ -70,7 +70,7 @@
     - [BADTRACE,       17, "Bad trace solution"]
     - [BADFLAT,        18, "Bad fiber flat"]
     - [BADARC,         19, "Bad arc solution"]
-    - [MANYBADCOL,     20, ">10% of spectral bins are bad columns"]
+    - [VARIABLE,       20, "Object or calibration suspected to be variable"]
     - [MANYREJECTED,   21, ">10% of spectral bins rejected in extraction"]
     - [BADAMPB,        22, "Issues in the amplifier readouts of camera B make this unusable"]
     - [BADAMPR,        23, "Issues in the amplifier readouts of camera R make this unusable"]

py/desispec/maskbits.py

@@ -44,7 +44,10 @@ def parse(options=None):
             help="coadd spectra of different cameras. works only if wavelength grids are aligned")
     parser.add_argument("--onetile", action="store_true",
             help="input spectra are from a single tile")
-
+    parser.add_argument("--no-normalize", action="store_true",
+            help="multiplicative normalization factors are not computed for each target exposure prior to coaddition")
+    parser.add_argument("--norm-chi2-threshold", type=float, default=0.1,
+            help="minimum reduced chi‑squared value to trigger per‑exposure normalization; use 0 to renormalize all science targets")
 
     if options is None:
         args = parser.parse_args()
@@ -156,7 +159,8 @@ def main(args=None):
                 onetile=args.onetile)
     else :
         log.info("coadding ...")
-        coadd(spectra, cosmics_nsig=args.nsig, onetile=args.onetile)
+        coadd(spectra, cosmics_nsig=args.nsig, onetile=args.onetile,
+              no_normalize=args.no_normalize, norm_chi2_threshold=args.norm_chi2_threshold)
 
     if args.lin_step is not None :
         log.info("resampling ...")

py/desispec/scripts/coadd_spectra.py

@@ -230,6 +230,28 @@ def test_coadd(self):
         self.assertEqual(s1.flux['b'].shape[0], 1)
         self.assertIsInstance(s1.scores, Table)
 
+    def test_coadd_normalization(self):
+        """Test coaddition"""
+        nspec, nwave = 10, 300
+        bands = ('b','z')
+        s1 = self._random_spectra( ns=nspec, nw=nwave, bands=bands)
+        self.assertEqual(s1.flux['b'].shape[0], nspec)
+        self.assertIsInstance(s1.scores, Table)
+
+        #- three different targets, one with only one exposure
+        s1.fibermap['TARGETID'][:5] = 10
+        for b in bands:
+            # offset first spectrum
+            s1.flux[b][0] +=5
+            s1.ivar[b][0] /=3
+        s1.fibermap['TARGETID'][5:] = 11
+        s1.fibermap['TARGETID'][-1] = 12
+        s1.fibermap['OBJTYPE'] = 'TGT'
+        coadd(s1, onetile=True, no_normalize=False)
+        self.assertEqual(s1.flux['b'].shape[0], 3)
+        self.assertEqual(s1.exp_fibermap['COADD_NORM'][-1], 1)
+        self.assertIsInstance(s1.scores, Table)
+
     def test_coadd_exposures(self):
         """Test coadd_exposures doesn't mutate input and matches in-place coadd"""
         nspec, nwave = 3, 10

py/desispec/test/test_coadd.py

@@ -47,6 +47,7 @@ def setUpClass(cls):
         fibermap['TILEID'] = 1000 * np.ones(nspec, dtype='i4')
         fibermap['FIBERSTATUS'] = np.zeros(nspec, dtype='i4')
         fibermap['DESI_TARGET'] = np.arange(nspec, dtype='i8')
+        fibermap['OBJTYPE'] = 'TGT'
 
         cls.spectra = Spectra(bands=bands, wave=wave, flux=flux, ivar=ivar,
                               resolution_data=rdat, fibermap=fibermap)