Switch to use romanisim

pyproject.toml

@@ -23,6 +23,7 @@ dependencies = [
   "skycatalogs",
   "packaging",
   "setuptools",
+  "romanisim",
 ]
 
 [project.urls]

roman_imsim/bandpass.py

@@ -1,4 +1,5 @@
-import galsim.roman as roman
+import romanisim.models as models
+
 from galsim.config import BandpassBuilder, GetAllParams, RegisterBandpassType
 
 
@@ -22,10 +23,16 @@ def buildBandpass(self, config, base, logger):
             the constructed Bandpass object.
         """
         req = {"name": str}
-        kwargs, safe = GetAllParams(config, base, req=req)
+        opt = {"SCA": int}
+        kwargs, safe = GetAllParams(config, base, req=req, opt=opt)
 
         name = kwargs["name"]
-        bandpass = roman.getBandpasses(red_limit=2000)[name]
+        if "SCA" in kwargs:
+            bandpass = models.bandpass.getBandpasses(red_limit=2000, sca=kwargs["SCA"], bandnames=[name])[
+                name
+            ]
+        else:
+            bandpass = models.bandpass.getBandpasses(red_limit=2000, bandnames=[name])[name]
 
         return bandpass, safe
 

roman_imsim/config/default.yaml

@@ -32,7 +32,8 @@ image:
         mjd: { type: ObSeqData, field: mjd }
 
     bandpass:
-        type: RomanBandpass
+        type: RomanBandpassTrimmed
+        SCA: '@image.SCA'
         name: { type: ObSeqData, field: filter }
 
     # When you want to have multiple images generate the same random galaxies, then
@@ -65,7 +66,7 @@ image:
     #     last: 18
     #     repeat: 3  # repeat each SCA num 3 times before moving on, for the 3 filters.
     #
-    SCA: 10
+    SCA: 4
     mjd: { type: ObSeqData, field: mjd }
     filter: { type: ObSeqData, field: filter }
     exptime: { type: ObSeqData, field: exptime }
@@ -79,10 +80,32 @@ image:
     # is even more pronounced in this case.
     use_fft_bright: True
 
-    # Override this in your config to enable noise, e.g. with type: RomanNoise
     noise:
         # These are all by default turned on, but you can turn any of them off if desired:
-        type: NoNoise
+        type: RomanNoise
+        mjd: { type: ObSeqData, field: mjd }
+        stray_light: True
+        thermal_background: True
+        reciprocity_failure: True
+        dark_current:
+            turn_on: True
+            use_crds: True
+        nonlinearity:
+            turn_on: True
+            use_crds: True
+        ipc:
+            turn_on: True
+            use_crds: True
+        read_noise:
+            turn_on: True
+            use_crds: True
+        gain:
+            turn_on: True
+            use_crds: True
+        sky_subtract: False
+    # noise:
+    #     # These are all by default turned on, but you can turn any of them off if desired:
+    #     type: NoNoise
 
 stamp:
     type: Roman_stamp
@@ -105,24 +128,24 @@ gal:
 
 input:
     obseq_data:
-        dir: /hpc/group/cosmology/OpenUniverse2024/
-        # file_name: Roman_WAS_obseq_11_1_23.fits
-        visit: 1
+        # dir: /hpc/group/cosmology/OpenUniverse2024/
+        file_name: /home/yuedong/Work/SimHackJune2025/Roman_WAS_obseq_11_1_23.fits
+        visit: 12909
         SCA: '@image.SCA'
     RomanPSFInterpolator:
         kind: corners
         n_waves: 5
     sky_catalog:
-        file_name: config/skyCatalog.yaml
+        file_name: /home/yuedong/Work/SimHackJune2025/skyCatalog/skyCatalog.yaml
         edge_pix: 512
         mjd: { type: ObSeqData, field: mjd }
-        exptime: { type: ObSeqData, field: exptime }          
+        exptime: { type: ObSeqData, field: exptime }
         obj_types: ['diffsky_galaxy', 'snana', 'star', 'gaia_star']
 
 output:
 
     nfiles: 1
-    dir: output
+    dir: /home/yuedong/Work/SimHackJune2025/output/RomanWAS_new/images/truth
     file_name:
         type: FormattedStr
         format: "Roman_LSST_true_%s_%s_%s.fits.gz"
@@ -132,7 +155,7 @@ output:
             - '@image.SCA'
 
     truth:
-        dir: output
+        dir: output/RomanWAS_new/truth
         file_name:
             type: FormattedStr
             format: "Roman_LSST_index_%s_%s_%s.txt"

roman_imsim/detector_physics.py

@@ -172,7 +172,15 @@ 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)
+        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):
 
@@ -309,7 +317,13 @@ def add_effects_scafile(self, im, wt, pointing):
             self.sky_mean,
             sky_noise,
         )
-        return im, self.sky[self.sky.bounds & im.bounds] - self.sky_mean, dq, self.sky_mean, sky_noise
+        return (
+            im,
+            self.sky[self.sky.bounds & im.bounds] - self.sky_mean,
+            dq,
+            self.sky_mean,
+            sky_noise,
+        )
 
     def add_effects_galsim(self, im, wt, pointing):
         """
@@ -439,7 +453,13 @@ def add_effects_simple(self, im, wt, pointing):
             self.sky_mean,
             sky_noise,
         )
-        return im, self.sky[self.sky.bounds & im.bounds] - self.sky_mean, dq, self.sky_mean, sky_noise
+        return (
+            im,
+            self.sky[self.sky.bounds & im.bounds] - self.sky_mean,
+            dq,
+            self.sky_mean,
+            sky_noise,
+        )
 
     def set_diff(self, im=None):
         if self.params["save_diff"]:
@@ -509,7 +529,16 @@ def bfe(self, im):
                 a = (a + np.fliplr(a) + np.flipud(a) + np.flip(a)) / 4.0
 
                 r = 0.5 * (3.25 / 4.25) ** (1.5) / 1.5  # source-boundary projection
-                B = (a[2, 2], a[3, 2], a[2, 3], a[3, 3], a[4, 2], a[2, 4], a[3, 4], a[4, 4])
+                B = (
+                    a[2, 2],
+                    a[3, 2],
+                    a[2, 3],
+                    a[3, 3],
+                    a[4, 2],
+                    a[2, 4],
+                    a[3, 4],
+                    a[4, 4],
+                )
 
                 A = np.array(
                     [
@@ -584,7 +613,6 @@ def bfe(self, im):
 
         for gj in range(num_grids):
             for gi in range(num_grids):
-
                 a_components_pad = np.zeros(
                     (
                         4,
@@ -595,7 +623,13 @@ def bfe(self, im):
                     )
                 )  # (4,5,5,sub_grid,sub_grid)
                 a_components_pad = np.zeros(
-                    (4, 2 * nbfe + 1, 2 * nbfe + 1, bin_size * n_sub + 2 * nbfe, bin_size * m_sub + 2 * nbfe)
+                    (
+                        4,
+                        2 * nbfe + 1,
+                        2 * nbfe + 1,
+                        bin_size * n_sub + 2 * nbfe,
+                        bin_size * m_sub + 2 * nbfe,
+                    )
                 )  # (4,5,5,sub_grid,sub_grid)
 
                 for comp in range(4):
@@ -858,7 +892,6 @@ def dark_current(self, im):
             self.im_dark = im - self.im_dark
 
         else:
-
             dark_current_ = self.dark_current_.clip(0)
 
             # opt for numpy random geneator instead for speed
@@ -989,7 +1022,6 @@ def add_persistence(self, im, pointing):
                 im.array[:, :] += galsim.roman.roman_detectors.fermi_linear(x.array, dt) * roman.exptime
 
         else:
-
             # setup parameters for persistence
             Q01 = self.df["PERSIST"].read_header()["Q01"]
             Q02 = self.df["PERSIST"].read_header()["Q02"]
@@ -1023,7 +1055,7 @@ def add_persistence(self, im, pointing):
 
                 # Do linear interpolation
                 a = np.zeros(x.shape)
-                a += ((x < Q01)) * x / Q01
+                a += (x < Q01) * x / Q01
                 a += ((x >= Q01) & (x < Q02)) * (Q02 - x) / (Q02 - Q01)
                 im.array[:, :] += a * self.df["PERSIST"][0, :, :][0] * fac_dt
 
@@ -1049,7 +1081,7 @@ def add_persistence(self, im, pointing):
 
                 a = np.zeros(x.shape)
                 a += ((x >= Q05) & (x < Q06)) * (x - Q05) / (Q06 - Q05)
-                a += ((x >= Q06)) * (x / Q06) ** alpha  # avoid fractional power of negative values
+                a += (x >= Q06) * (x / Q06) ** alpha  # avoid fractional power of negative values
                 im.array[:, :] += a * self.df["PERSIST"][5, :, :][0] * fac_dt
 
         return im
@@ -1132,7 +1164,6 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel):
 
                     for dy in range(-1, 2):
                         for dx in range(-1, 2):
-
                             array_out[
                                 gj * grid_size : (gj + 1) * grid_size,
                                 gi * grid_size : (gi + 1) * grid_size,
@@ -1258,7 +1289,6 @@ def finalize_sky_im(self, im, pointing):
             im = self.e_to_ADU(im)
             im.quantize()
         else:
-
             bound_pad = galsim.BoundsI(xmin=1, ymin=1, xmax=4096, ymax=4096)
             im_pad = galsim.Image(bound_pad)
             im_pad.array[4:-4, 4:-4] = im.array[:, :]