Scattered light model

doc/releases/1.14.1dev.rst

@@ -24,6 +24,8 @@ Script Changes
 
 - When making datacubes, the user can select a separate frame to use for the sky subtraction.
   In this case, it is the processed data that will be used for sky subtraction (akin to nodding).
+- A new script has been written (``chk_scattlight``) to check the generation of the scattered
+  light calibration frame model.
 
 Datamodel Changes
 -----------------

pypeit/calibrations.py

@@ -20,13 +20,15 @@
 from pypeit import alignframe
 from pypeit import flatfield
 from pypeit import edgetrace
+from pypeit import scattlight
 from pypeit import slittrace
 from pypeit import wavecalib
 from pypeit import wavetilts
 from pypeit.calibframe import CalibFrame
 from pypeit.images import buildimage
 from pypeit.metadata import PypeItMetaData
 from pypeit.core import framematch
+from pypeit.core import parse
 from pypeit.par import pypeitpar
 from pypeit.spectrographs.spectrograph import Spectrograph
 from pypeit import io
@@ -181,6 +183,7 @@ def __init__(self, fitstbl, par, spectrograph, caldir, qadir=None,
         self.msbpm = None
         self.wv_calib = None
         self.slits = None
+        self.msscattlight = None
 
         self.wavetilts = None
         self.flatimages = None
@@ -321,8 +324,8 @@ def get_arc(self):
         self.msarc = buildimage.buildimage_fromlist(self.spectrograph, self.det,
                                                     self.par['arcframe'], raw_files,
                                                     bias=self.msbias, bpm=self.msbpm,
-                                                    dark=self.msdark, calib_dir=self.calib_dir,
-                                                    setup=setup, calib_id=calib_id)
+                                                    dark=self.msdark, scattlight=self.msscattlight,
+                                                    calib_dir=self.calib_dir, setup=setup, calib_id=calib_id)
         # Save the result
         self.msarc.to_file()
         # Return it
@@ -365,6 +368,7 @@ def get_tiltimg(self):
                                                      self.par['tiltframe'], raw_files,
                                                      bias=self.msbias, bpm=self.msbpm,
                                                      dark=self.msdark, slits=self.slits,
+                                                     scattlight=self.msscattlight,
                                                      calib_dir=self.calib_dir, setup=setup,
                                                      calib_id=calib_id)
         # Save the result
@@ -414,7 +418,8 @@ def get_align(self):
                                                  self.par['alignframe'], raw_files,
                                                  bias=self.msbias, bpm=self.msbpm,
                                                  dark=self.msdark, calib_dir=self.calib_dir,
-                                                 setup=setup, calib_id=calib_id)
+                                                 setup=setup, calib_id=calib_id,
+                                                 scattlight=self.msscattlight)
 
         # Instantiate
         # TODO: From JFH: Do we need the bpm here?  Check that this was in the previous code.
@@ -537,6 +542,95 @@ def get_bpm(self, frame=None):
         # Return
         return self.msbpm
 
+    def get_scattlight(self):
+        """
+        Load or generate the scattered light model.
+
+        Returns:
+            :class:`~pypeit.scattlight.ScatteredLight`: The processed calibration image including the model.
+        """
+        # Initialise to None
+        self.msscattlight = None
+
+        # Check for existing data
+        if not self._chk_objs(['msbpm', 'slits']):
+            msgs.warn('Must have the bpm and the slits defined to make a scattered light image!  '
+                      'Skipping and may crash down the line')
+            return self.msscattlight
+
+        # Check internals
+        self._chk_set(['det', 'calib_ID', 'par'])
+
+        # Prep
+        frame = {'type': 'scattlight', 'class': scattlight.ScatteredLight}
+        raw_scattlight_files, cal_file, calib_key, setup, calib_id, detname = \
+            self.find_calibrations(frame['type'], frame['class'])
+        scatt_idx = self.fitstbl.find_frames(frame['type'], calib_ID=self.calib_ID, index=True)
+
+        if len(raw_scattlight_files) == 0 and cal_file is None:
+            msgs.warn(f'No raw {frame["type"]} frames found and unable to identify a relevant '
+                      'processed calibration frame.  Continuing...')
+            return self.msscattlight
+
+        # If a processed calibration frame exists and we want to reuse it, do
+        # so:
+        if cal_file.exists() and self.reuse_calibs:
+            self.msscattlight = frame['class'].from_file(cal_file)
+            return self.msscattlight
+
+        # Scattered light model does not exist or we're not reusing it.
+        # Need to build everything from scratch.  Start with the trace image.
+        msgs.info('Creating edge tracing calibration frame using files: ')
+        for f in raw_scattlight_files:
+            msgs.prindent(f'{Path(f).name}')
+
+        # Reset the BPM
+        self.get_bpm(frame=raw_scattlight_files[0])
+
+        binning = self.fitstbl[scatt_idx[0]]['binning']
+        scattlightImage = buildimage.buildimage_fromlist(self.spectrograph, self.det,
+                                                         self.par['scattlightframe'], raw_scattlight_files,
+                                                         bias=self.msbias, bpm=self.msbpm,
+                                                         dark=self.msdark, calib_dir=self.calib_dir,
+                                                         setup=setup, calib_id=calib_id)
+
+        spatbin = parse.parse_binning(binning)[1]
+        pad = self.par['scattlight']['pad'] // spatbin
+        offslitmask = self.slits.slit_img(pad=pad, initial=True, flexure=None)==-1
+
+        model, modelpar, success = self.spectrograph.scattered_light(scattlightImage.image, offslitmask,
+                                                                     self.fitstbl[scatt_idx[0]])
+
+        if not success:
+            # Something went awry
+            msgs.warn('Scattered light modelling failed.  Continuing, but likely to fail soon...')
+            self.success = False
+            return self.msscattlight
+
+        # Now generate the DataModel
+        self.msscattlight = scattlight.ScatteredLight(PYP_SPEC=self.spectrograph.name,
+                                                      pypeline=self.spectrograph.pypeline,
+                                                      detname=scattlightImage.detector.name,
+                                                      nspec=scattlightImage.shape[0], nspat=scattlightImage.shape[1],
+                                                      binning=scattlightImage.detector.binning,
+                                                      pad=self.par['scattlight']['pad'],
+                                                      scattlight_raw=scattlightImage.image,
+                                                      scattlight_model=model,
+                                                      scattlight_param=modelpar)
+
+        # TODO :: Should we go back and recalculate the slit edges once the scattered light is known?
+
+        # Show the result if requested
+        if self.show:
+            self.msscattlight.show()
+
+        if self.msscattlight is not None:
+            self.msscattlight.set_paths(self.calib_dir, setup, calib_id, detname)
+            # Save the master scattered light model
+            self.msscattlight.to_file()
+
+        return self.msscattlight
+
     def get_flats(self):
         """
         Load or generate the flat-field calibration images.
@@ -622,7 +716,8 @@ def get_flats(self):
             pixel_flat = buildimage.buildimage_fromlist(self.spectrograph, self.det,
                                                         self.par['pixelflatframe'],
                                                         raw_pixel_files, dark=self.msdark,
-                                                        bias=self.msbias, bpm=self.msbpm)
+                                                        bias=self.msbias, bpm=self.msbpm,
+                                                        scattlight=self.msscattlight)
             if len(raw_lampoff_files) > 0:
                 # Reset the BPM
                 self.get_bpm(frame=raw_lampoff_files[0])
@@ -633,7 +728,7 @@ def get_flats(self):
                                                               self.par['lampoffflatsframe'],
                                                               raw_lampoff_files,
                                                               dark=self.msdark, bias=self.msbias,
-                                                              bpm=self.msbpm)
+                                                              bpm=self.msbpm, scattlight=self.msscattlight)
                 pixel_flat = pixel_flat.sub(lampoff_flat)
 
             # Initialise the pixel flat
@@ -656,7 +751,7 @@ def get_flats(self):
                 msgs.prindent(f'{Path(f).name}')
             illum_flat = buildimage.buildimage_fromlist(self.spectrograph, self.det,
                                                         self.par['illumflatframe'], raw_illum_files,
-                                                        dark=self.msdark, bias=self.msbias,
+                                                        dark=self.msdark, bias=self.msbias, scattlight=self.msscattlight,
                                                         flatimages=self.flatimages, bpm=self.msbpm)
             if len(raw_lampoff_files) > 0:
                 msgs.info('Subtracting lamp off flats using files: ')
@@ -667,7 +762,9 @@ def get_flats(self):
                                                                   self.par['lampoffflatsframe'],
                                                                   raw_lampoff_files,
                                                                   dark=self.msdark,
-                                                                  bias=self.msbias, bpm=self.msbpm)
+                                                                  bias=self.msbias,
+                                                                  scattlight=self.msscattlight,
+                                                                  bpm=self.msbpm)
                 illum_flat = illum_flat.sub(lampoff_flat)
 
             # Initialise the pixel flat
@@ -776,6 +873,7 @@ def get_slits(self):
         traceImage = buildimage.buildimage_fromlist(self.spectrograph, self.det,
                                                     self.par['traceframe'], raw_trace_files,
                                                     bias=self.msbias, bpm=self.msbpm,
+                                                    scattlight=self.msscattlight,
                                                     dark=self.msdark, calib_dir=self.calib_dir,
                                                     setup=setup, calib_id=calib_id)
         if len(raw_lampoff_files) > 0:
@@ -789,7 +887,8 @@ def get_slits(self):
             lampoff_flat = buildimage.buildimage_fromlist(self.spectrograph, self.det,
                                                           self.par['lampoffflatsframe'],
                                                           raw_lampoff_files, dark=self.msdark,
-                                                          bias=self.msbias, bpm=self.msbpm)
+                                                          bias=self.msbias, scattlight=self.msscattlight,
+                                                          bpm=self.msbpm)
             traceImage = traceImage.sub(lampoff_flat)
 
         edges = edgetrace.EdgeTraceSet(traceImage, self.spectrograph, self.par['slitedges'],
@@ -1254,7 +1353,7 @@ def default_steps():
         # Order matters!  And the name must match a viable "get_{step}" method
         # in Calibrations.
         # TODO: Does the bpm need to be done after the dark?
-        return ['bias', 'dark', 'bpm', 'slits', 'arc', 'tiltimg', 'wv_calib', 'tilts', 'flats']
+        return ['bias', 'dark', 'bpm', 'slits', 'arc', 'tiltimg', 'wv_calib', 'tilts', 'scattlight', 'flats']
 
 
 class IFUCalibrations(Calibrations):
@@ -1272,7 +1371,7 @@ def default_steps():
         """
         # Order matters!
         return ['bias', 'dark', 'bpm', 'arc', 'tiltimg', 'slits', 'wv_calib', 'tilts', 'align',
-                'flats']
+                'scattlight', 'flats']
 
 
 def check_for_calibs(par, fitstbl, raise_error=True, cut_cfg=None):

pypeit/core/framematch.py

@@ -29,9 +29,10 @@ def __init__(self):
                      ('pinhole', 'Pinhole observation used for tracing slit centers'),
                    ('pixelflat', 'Flat-field exposure used for pixel-to-pixel response'),
                    ('illumflat', 'Flat-field exposure used for illumination flat'),
-                   ('lampoffflats', 'Flat-field exposure with lamps off used to remove '
-                                   'persistence from lamp on flat exposures and/or thermal emission '
-                                   'from the telescope and dome'),
+                ('lampoffflats', 'Flat-field exposure with lamps off used to remove '
+                                 'persistence from lamp on flat exposures and/or thermal emission '
+                                 'from the telescope and dome'),
+                  ('scattlight', 'Frame (usually similar to pixelflat) used to determine scattered light model'),
                      ('science', 'On-sky observation of a primary target'),
                     ('standard', 'On-sky observation of a flux calibrator'),
                        ('trace', 'High-count exposure used to trace slit positions'),

pypeit/images/buildimage.py

@@ -77,6 +77,18 @@ class TiltImage(pypeitimage.PypeItCalibrationImage):
     calib_type = 'Tiltimg'
 
 
+class ScatteredLightImage(pypeitimage.PypeItCalibrationImage):
+    """
+    Simple DataContainer for the Scattered Light Image
+    """
+    # version is inherited from PypeItImage
+
+    # I/O
+    output_to_disk = ('SCATTLIGHT_IMAGE', 'SCATTLIGHT_FULLMASK', 'SCATTLIGHT_DETECTOR')
+    hdu_prefix = 'SCATTLIGHT_'
+    calib_type = 'ScattLight'
+
+
 class TraceImage(pypeitimage.PypeItCalibrationImage):
     """
     Simple DataContainer for the Trace Image
@@ -133,6 +145,7 @@ def construct_file_name(cls, calib_key, calib_dir=None, basename=None):
     arc=ArcImage,
     tilt=TiltImage,
     trace=TraceImage,
+    scattlight=ScatteredLightImage,
     align=AlignImage)
 """
 The list of classes that :func:`buildimage_fromlist` should use to decorate the
@@ -144,7 +157,7 @@ def construct_file_name(cls, calib_key, calib_dir=None, basename=None):
 
 
 def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm=None, dark=None,
-                        flatimages=None, maxiters=5, ignore_saturation=True, slits=None,
+                        scattlight=None, flatimages=None, maxiters=5, ignore_saturation=True, slits=None,
                         mosaic=None, calib_dir=None, setup=None, calib_id=None):
     """
     Perform basic image processing on a list of images and combine the results.
@@ -179,6 +192,8 @@ def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm=
         dark (:class:`~pypeit.images.buildimage.DarkImage`, optional):
             Dark-current image; passed directly to
             :func:`~pypeit.images.rawimage.RawImage.process` for all images.
+        scattlight (:class:`~pypeit.scattlight.ScatteredLight`, optional):
+            Scattered light model to be used to determine scattered light.
         flatimages (:class:`~pypeit.flatfield.FlatImages`, optional):
             Flat-field images for flat fielding; passed directly to
             :func:`~pypeit.images.rawimage.RawImage.process` for all images.
@@ -234,7 +249,7 @@ def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm=
 
     # Do it
     combineImage = combineimage.CombineImage(spectrograph, det, frame_par['process'], file_list)
-    pypeitImage = combineImage.run(bias=bias, bpm=bpm, dark=dark, flatimages=flatimages,
+    pypeitImage = combineImage.run(bias=bias, bpm=bpm, dark=dark, flatimages=flatimages, scattlight=scattlight,
                                    sigma_clip=frame_par['process']['clip'],
                                    sigrej=frame_par['process']['comb_sigrej'],
                                    maxiters=maxiters, ignore_saturation=ignore_saturation,

pypeit/images/combineimage.py

@@ -62,9 +62,8 @@ def __init__(self, spectrograph, det, par, files):
         if self.nfiles == 0:
             msgs.error('CombineImage requires a list of files to instantiate')
 
-    def run(self, bias=None, flatimages=None, ignore_saturation=False, sigma_clip=True,
-            bpm=None, sigrej=None, maxiters=5, slits=None, dark=None, combine_method='mean',
-            mosaic=False):
+    def run(self, bias=None, scattlight=None, flatimages=None, ignore_saturation=False, sigma_clip=True,
+            bpm=None, sigrej=None, maxiters=5, slits=None, dark=None, combine_method='mean', mosaic=False):
         r"""
         Process and combine all images.
 
@@ -137,6 +136,8 @@ def run(self, bias=None, flatimages=None, ignore_saturation=False, sigma_clip=Tr
             bias (:class:`~pypeit.images.buildimage.BiasImage`, optional):
                 Bias image for bias subtraction; passed directly to
                 :func:`~pypeit.images.rawimage.RawImage.process` for all images.
+            scattlight (:class:`~pypeit.scattlight.ScatteredLight`, optional):
+                Scattered light model to be used to determine scattered light.
             flatimages (:class:`~pypeit.flatfield.FlatImages`, optional):
                 Flat-field images for flat fielding; passed directly to
                 :func:`~pypeit.images.rawimage.RawImage.process` for all images.
@@ -194,7 +195,7 @@ def run(self, bias=None, flatimages=None, ignore_saturation=False, sigma_clip=Tr
             # Load raw image
             rawImage = rawimage.RawImage(ifile, self.spectrograph, self.det)
             # Process
-            pypeitImage = rawImage.process(self.par, bias=bias, bpm=bpm, dark=dark,
+            pypeitImage = rawImage.process(self.par, scattlight=scattlight, bias=bias, bpm=bpm, dark=dark,
                                            flatimages=flatimages, slits=slits, mosaic=mosaic)
 
             if self.nfiles == 1: