Incorporate det into sca

.gitignore

@@ -0,0 +1,3 @@
+/build/*
+/dist/*
+/roman_imsim.egg-info/*

config/was.yaml

@@ -88,16 +88,42 @@ image:
     # read_noise: False
     # sky_subtract: False
 
-    ignore_noise: False
-    stray_light: True
-    thermal_background: True
-    reciprocity_failure: True
-    dark_current: True
-    nonlinearity: True
-    ipc: True
-    read_noise: True
+    add_effects:
+        Background:
+            model: simple_model
+            thermal_background: True
+            stray_light: True
+        QuantumEfficiency: 
+            model: lab_model
+        BrighterFatter: 
+            model: lab_model
+        Persistence:
+            model: lab_model
+        # ReciprocityFailure:
+        #     model: simple_model
+        DarkCurrent: 
+            model: lab_model
+        Saturation:
+            model: lab_model
+        Nonlinearity: 
+            model: lab_model
+        IPC:
+            model: lab_model
+        DeadPixel:
+            model: lab_model
+        VTPE:
+            model: lab_model
+        ReadNoise: 
+            model: lab_model
+        Gain:
+            model: lab_model
+        Bias:
+            model: lab_model
     sky_subtract: False
 
+    sca_filepath: /hpc/group/cosmology/phy-lsst/roman_sensors
+    save_diff: False
+
     # nobjects: 500
 
 stamp:
@@ -124,17 +150,14 @@ gal:
 
 input:
     obseq_data:
-        # file_name: /lus/grand/projects/RomanDESC/final_roman_runfiles/was/Roman_WAS_obseq_11_1_23.fits
         file_name: /hpc/group/cosmology/OpenUniverse2024/RomanWAS/Roman_WAS_obseq_11_1_23.fits
-        visit: 1
+        visit: 4906
         SCA: '@image.SCA'
     roman_psf:
         SCA: '@image.SCA'
         n_waves: 5
     sky_catalog:
-        # file_name: /lus/grand/projects/RomanDESC/roman_rubin_cats_v1.1.2_faint/skyCatalog.yaml 
-        # file_name: /hpc/group/cosmology/OpenUniverse2024/roman_rubin_cats_v1.1.2_faint/skyCatalog.yaml 
-        file_name: /hpc/home/yf194/Work/projects/roman_imsim/config/skyCatalog.yaml
+        file_name: /hpc/group/cosmology/OpenUniverse2024/roman_rubin_cats_v1.1.2_faint/skyCatalog.yaml
         edge_pix: 512
         mjd: { type: ObSeqData, field: mjd }
         exptime: { type: ObSeqData, field: exptime }

roman_imsim/__init__.py

@@ -16,3 +16,6 @@
 from .skycat import *
 from .stamp import *
 from .wcs import *
+from .skycat import *
+from .photonOps import *
+from .bandpass import *

roman_imsim/detector_physics.py

@@ -83,15 +83,6 @@ def __init__(self, params, visit, SCA):
         )[self.sca]
         self.radec = self.WCS.toWorld(galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2)))
 
-        self.WCS = roman.getWCS(
-            world_pos=galsim.CelestialCoord(ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"]),
-            PA=obseq_data.ob["pa"],
-            date=self.date,
-            SCAs=self.sca,
-            PA_is_FPA=True,
-        )[self.sca]
-        self.radec = self.WCS.toWorld(galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2)))
-
 
 class modify_image(object):
     """
@@ -121,8 +112,6 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_
         if sca_filepath is not None:
             self.df = fio.FITS(sca_filepath + "/" + sca_number_to_file[self.pointing.sca])
             print("------- Using SCA files --------")
-            self.df = fio.FITS(sca_filepath + "/" + sca_number_to_file[self.pointing.sca])
-            print("------- Using SCA files --------")
         else:
             self.df = None
             print("------- Using simple detector model --------")
@@ -139,16 +128,6 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_
             "truth", self.get_path_name(use_galsim=use_galsim)
         )
 
-        b = galsim.BoundsI(xmin=1, ymin=1, xmax=roman.n_pix, ymax=roman.n_pix)
-        old_filename = os.path.join(self.params["output"]["dir"], imfilename)
-        if not os.path.exists(
-            self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim))
-        ):
-            os.mkdir(self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim)))
-        new_filename = os.path.join(self.params["output"]["dir"], imfilename).replace(
-            "truth", self.get_path_name(use_galsim=use_galsim)
-        )
-
         b = galsim.BoundsI(xmin=1, ymin=1, xmax=roman.n_pix, ymax=roman.n_pix)
         im = fio.FITS(old_filename)[-1].read()
         im = galsim.Image(im, bounds=b, wcs=self.pointing.WCS)
@@ -160,7 +139,6 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_
                 force_cvz = True
         self.setup_sky(im, self.pointing, rng, visit * sca, force_cvz=force_cvz)
 
-        img, err, dq, sky_mean, sky_noise = self.add_effects(im, None, self.pointing, use_galsim=use_galsim)
         img, err, dq, sky_mean, sky_noise = self.add_effects(im, None, self.pointing, use_galsim=use_galsim)
 
         write_fits(
@@ -172,7 +150,6 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_
             self.pointing.sca,
             sky_mean=sky_mean,
         )
-        write_fits(old_filename, new_filename, img, sky_noise, dq, self.pointing.sca, sky_mean=sky_mean)
 
     def get_path_name(self, use_galsim=False):
 
@@ -562,10 +539,6 @@ def bfe(self, im):
         # The img is clipped by the saturation level here to cap the brighter fatter effect
         # and avoid unphysical behavior
 
-        array_pad = self.saturate(im.copy()).array[4:-4, 4:-4]  # img of interest 4088x4088
-        array_pad = np.pad(
-            array_pad, [(4 + nbfe, 4 + nbfe), (4 + nbfe, 4 + nbfe)], mode="symmetric"
-        )  # 4100x4100 array
         array_pad = self.saturate(im.copy()).array[4:-4, 4:-4]  # img of interest 4088x4088
         array_pad = np.pad(
             array_pad, [(4 + nbfe, 4 + nbfe), (4 + nbfe, 4 + nbfe)], mode="symmetric"
@@ -697,8 +670,6 @@ def get_eff_sky_bg(self, pointing, radec):
         radec               : World coordinate position of image
         """
 
-        sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date)
-        sky_level *= (1.0 + roman.stray_light_fraction) * roman.pixel_scale**2
         sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date)
         sky_level *= (1.0 + roman.stray_light_fraction) * roman.pixel_scale**2
 
@@ -739,9 +710,6 @@ def setup_sky(self, im, pointing, rng, rng_iter, force_cvz=False):
             self.dark_current_ = (
                 roman.dark_current * roman.exptime + self.df["DARK"][:, :].flatten() * roman.exptime
             )
-            self.dark_current_ = (
-                roman.dark_current * roman.exptime + self.df["DARK"][:, :].flatten() * roman.exptime
-            )
         if self.df is None:
             self.gain = roman.gain
         else:
@@ -755,8 +723,6 @@ def setup_sky(self, im, pointing, rng, rng_iter, force_cvz=False):
             radec = pointing.radec
         sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date)
         sky_level *= 1.0 + roman.stray_light_fraction
-        sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date)
-        sky_level *= 1.0 + roman.stray_light_fraction
         # Make a image of the sky that takes into account the spatially variable pixel scale. Note
         # that makeSkyImage() takes a bit of time. If you do not care about the variable pixel
         # scale, you could simply compute an approximate sky level in e-/pix by multiplying
@@ -864,8 +830,6 @@ def dark_current(self, im):
             # opt for numpy random geneator instead for speed
             self.im_dark = self.rng_np.poisson(dark_current_).reshape(im.array.shape).astype(im.dtype)
             im.array[:, :] += self.im_dark
-            self.im_dark = self.rng_np.poisson(dark_current_).reshape(im.array.shape).astype(im.dtype)
-            im.array[:, :] += self.im_dark
 
         # NOTE: Sky level and dark current might appear like a constant background that can be
         # simply subtracted. However, these contribute to the shot noise and matter for the
@@ -960,11 +924,6 @@ def add_persistence(self, im, pointing):
         p_list = np.array([get_pointing(self.params, i, pointing.sca) for i in dither_list_selected])
         dt_list = np.array([(pointing.date - p.date).total_seconds() for p in p_list])
         p_pers = p_list[np.where((dt_list > 0) & (dt_list < roman.exptime * 10))]
-        dither_list_selected = dither_sca_array[dither_sca_array[:, 1] == pointing.sca, 0]
-        dither_list_selected = dither_list_selected[np.abs(dither_list_selected - pointing.visit) < 10]
-        p_list = np.array([get_pointing(self.params, i, pointing.sca) for i in dither_list_selected])
-        dt_list = np.array([(pointing.date - p.date).total_seconds() for p in p_list])
-        p_pers = p_list[np.where((dt_list > 0) & (dt_list < roman.exptime * 10))]
 
         if self.df is None:
             # iterate over previous exposures
@@ -986,7 +945,6 @@ def add_persistence(self, im, pointing):
                 x = x.clip(0)  # remove negative stimulus
 
                 im.array[:, :] += galsim.roman.roman_detectors.fermi_linear(x.array, dt) * roman.exptime
-                im.array[:, :] += galsim.roman.roman_detectors.fermi_linear(x.array, dt) * roman.exptime
 
         else:
 
@@ -1106,8 +1064,6 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel):
             num_grids = 4  # num_grids <= 8
             grid_size = 4096 // num_grids
 
-            array_pad = im.array[4:-4, 4:-4]  # it's an array instead of img
-            array_pad = np.pad(array_pad, [(5, 5), (5, 5)], mode="symmetric")  # 4098x4098 array
             array_pad = im.array[4:-4, 4:-4]  # it's an array instead of img
             array_pad = np.pad(array_pad, [(5, 5), (5, 5)], mode="symmetric")  # 4098x4098 array
 
@@ -1127,8 +1083,6 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel):
                         for i in range(3):
                             tmp = (t.dot(K[j, i, :, :])).dot(t.T)  # grid_sizexgrid_size
                             K_pad[j, i, :, :] = np.pad(tmp, [(1, 1), (1, 1)], mode="symmetric")
-                            tmp = (t.dot(K[j, i, :, :])).dot(t.T)  # grid_sizexgrid_size
-                            K_pad[j, i, :, :] = np.pad(tmp, [(1, 1), (1, 1)], mode="symmetric")
 
                     for dy in range(-1, 2):
                         for dx in range(-1, 2):
@@ -1177,11 +1131,6 @@ def add_read_noise(self, im):
                 self.rng_np.normal(loc=0.0, scale=read_noise).reshape(im.array.shape).astype(im.dtype)
             )
             im.array[:, :] += self.im_read
-            read_noise = self.df["READ"][2, :, :].flatten()  # flattened 4096x4096 array
-            self.im_read = (
-                self.rng_np.normal(loc=0.0, scale=read_noise).reshape(im.array.shape).astype(im.dtype)
-            )
-            im.array[:, :] += self.im_read
 
             # noise_array = self.rng_np.normal(loc=0., scale=read_noise)
             # 4088x4088 img

roman_imsim/effects/__init__.py

@@ -0,0 +1,15 @@
+from .roman_effects import RomanEffects, setup_sky
+from .brighter_fatter import BrighterFatter
+from .nonlinearity import Nonlinearity
+from .background import Background
+from .quantum_efficiency import QuantumEfficiency
+from .persistence import Persistence
+from .reciprocity_failure import ReciprocityFailure
+from .dark_current import DarkCurrent
+from .saturation import Saturation
+from .ipc import IPC
+from .dead_pixel import DeadPixel
+from .vtpe import VTPE
+from .read_noise import ReadNoise
+from .gain import Gain
+from .bias import Bias

roman_imsim/effects/background.py

@@ -0,0 +1,54 @@
+import os
+import galsim
+from . import RomanEffects
+import galsim.roman as roman
+
+
+class Background(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+        self.thermal_background = (
+            self.params["thermal_background"] if "thermal_background" in self.params else False
+        )
+        self.stray_light = self.params["stray_light"] if "stray_light" in self.params else False
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        if self.save_diff:
+            orig = image.copy()
+
+        self.logger.warning("Simple model will be applied for background.")
+        pointing = self.pointing
+        # Build current specification sky level if sky level not given
+        if self.force_cvz:
+            radec = self.translate_cvz(pointing.radec)
+        else:
+            radec = pointing.radec
+        sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date)
+        self.logger.debug("Adding sky_level = %s", sky_level)
+
+        if self.stray_light:
+            self.logger.debug("Stray light fraction = %s", roman.stray_light_fraction)
+            sky_level *= 1.0 + roman.stray_light_fraction
+        # Create sky image
+        self.sky = galsim.Image(bounds=image.bounds, wcs=pointing.WCS)
+        pointing.WCS.makeSkyImage(self.sky, sky_level)
+        if self.thermal_background:
+            tb = roman.thermal_backgrounds[pointing.filter] * pointing.exptime
+            self.logger.debug("Adding thermal background: %s", tb)
+            self.sky += tb
+        self.sky.addNoise(self.noise)
+
+        # [TODO] Not entirely sure about this block, since the 'auto' option is meant to
+        # let the software choose which drawing method to use based on the total flux.
+        if self.base["image"]["draw_method"] not in ["phot", "auto"]:
+            image.addNoise(self.noise)
+
+        # Adding sky level to the image.
+        image += self.sky[self.sky.bounds & image.bounds]
+        if self.save_diff:
+            prev = image.copy()
+            diff = prev - orig
+            diff.write(os.path.join(self.diff_dir, "sky_a.fits"))
+        return image

roman_imsim/effects/bias.py

@@ -0,0 +1,27 @@
+import os
+import fitsio as fio
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class Bias(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        self.logger.warning("No bias will be applied.")
+        return image
+
+    def lab_model(self, image):
+        if self.sca_filepath is None:
+            self.logger.warning("No bias data provided; no bias will be applied.")
+            return image
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        self.logger.warning("Lab measured model will be applied for bias.")
+        bias = self.df["BIAS"][:, :]  # 4096x4096 img
+
+        image.array[:, :] += bias
+        return image