Merge branch 'add_ghost' into 'develop' (c105b14b) · Commits · csst-sims / csst_msc_sim

config/obs_config_SCI.yaml

+10 −4

Original line number	Diff line number	Diff line
		@@ -9,13 +9,14 @@
		###############################################

		# Observation type
		obs_type: null
		obs_type: "WIDE"
		obs_type_code: "101"
		obs_id: "00000001" # this setting will only be used if pointing list file is not given

		# Define list of chips
		run_chips: [6,7,8,9,11,12,13,14,15,16,17,18,19,20,22,23,24,25] # Photometric chips
		# run_chips: [6,7,8,9,11,12,13,14,15,16,17,18,19,20,22,23,24,25] # Photometric chips
		#run_chips: [1,2,3,4,5,10,21,26,27,28,29,30] # Spectroscopic chips
		run_chips: [17]

		# Define observation sequence
		call_sequence:
		@@ -27,6 +28,11 @@ call_sequence:
		shutter_effect: YES
		flat_fielding: YES
		field_dist: YES

		ghost_SCI:
		mag_threshold: 16
		field_dist: YES

		# Accumulate fluxes from sky background
		sky_background:
		# [Optional]: exposure time of the pointing will be used as default.
		@@ -83,5 +89,5 @@ call_sequence:
		gain_16channel: YES
		# Output the final image
		quantization_and_output:
		format_output: YES
		format_output: NO
		...

observation_sim/mock_objects/Ghost.py

0 → 100644

+199 −0

Original line number	Diff line number	Diff line
		"""
		Ghost image calculation (single-class refactor)

		This module keeps the original behavior but wraps everything into ONE class: `Ghost`.

		All lengths are in meters.
		"""

		from __future__ import annotations

		from typing import Dict, Tuple

		import numpy as np
		import galsim


		class Ghost:
		"""
		Single-class ghost image model.

		Usage
		-----
		ghost = Ghost() # uses built-in defaults
		result = ghost.calculate(detector_id=13, x0=0.0, y0=0.0)

		The returned dict contains:
		- "ghosts": ndarray shape (3, 4) with rows [x, y, radius, relative_energy]
		- "angle_x_deg": float
		- "angle_y_deg": float
		- "wave_name": str (band name)
		"""

		# ----------- Default optical constants (meters) -----------
		DEFAULT_ZP = 4.593
		DEFAULT_RP = 0.16367
		DEFAULT_D1 = 0.005
		DEFAULT_D2 = 0.008
		DEFAULT_PIXEL = 1.0e-5 # kept for completeness; unused by formula

		def __init__(
		self,
		*,
		Zp: float = DEFAULT_ZP,
		Rp: float = DEFAULT_RP,
		D1: float = DEFAULT_D1,
		D2: float = DEFAULT_D2,
		pixel: float = DEFAULT_PIXEL,
		) -> None:
		self.Zp = float(Zp)
		self.Rp = float(Rp)
		self.D1 = float(D1)
		self.D2 = float(D2)
		self.pixel = float(pixel)

		# ---- Wave table (from original script) ----
		# Stored as: wave_num -> (wave_name, wave_length, nd, tf1, rf1, tf2, rf2, TD, RD)
		self._waves: Dict[int, Tuple[str, float, float, float, float, float, float, float, float]] = {
		0: ("nuv", 2.9e-7, 1.495, 0.825, 0.0093, 0.825, 0.0093, 0.555, 0.15),
		1: ("u", 3.6e-7, 1.473, 0.954, 0.0138, 0.954, 0.0138, 0.56, 0.15),
		2: ("g", 4.8e-7, 1.463, 0.98, 0.0197, 0.98, 0.0197, 0.8, 0.2),
		3: ("r", 6.2e-7, 1.458, 0.98, 0.015, 0.98, 0.015, 0.91, 0.09),
		4: ("i", 7.5e-7, 1.454, 0.985, 0.01, 0.985, 0.01, 0.889, 0.11),
		5: ("z", 9.1e-7, 1.452, 0.99, 0.0085, 0.99, 0.0085, 0.6, 0.1),
		6: ("y", 9.7e-7, 1.451, 0.99, 0.0035, 0.99, 0.0035, 0.36, 0.1),
		}

		# ---- Detector -> wave_num mapping (from original script) ----
		self._detector_to_wave_num: Dict[int, int] = {
		6: 6,
		7: 4,
		8: 2,
		9: 3,
		11: 5,
		12: 0,
		13: 0,
		14: 1,
		15: 6,
		16: 6,
		17: 1,
		18: 0,
		19: 0,
		20: 5,
		22: 3,
		23: 2,
		24: 4,
		25: 6,
		}

		# -----------------------------
		# Helpers
		# -----------------------------

		def _get_wave_tuple(self, detector_id: int) -> Tuple[str, float, float, float, float, float, float, float, float]:
		"""
		Return wave properties tuple for a detector_id.
		Raises KeyError with a friendly message if unknown.
		"""
		try:
		wave_num = self._detector_to_wave_num[detector_id]
		except KeyError as e:
		known = ", ".join(map(str, sorted(self._detector_to_wave_num.keys())))
		raise KeyError(f"Unknown detector_id={detector_id}. Known ids: {known}") from e

		return self._waves[wave_num]

		# -----------------------------
		# Public API
		# -----------------------------

		def calculate(self, detector_id: int, x0: float, y0: float) -> Dict[str, object]:
		"""
		Calculate ghost images for a given detector and focal-plane coordinate.

		Parameters
		----------
		detector_id : int
		Detector identifier.
		x0, y0 : float
		Absolute focal-plane coordinates of the image point [m].

		Returns
		-------
		dict with keys:
		ghosts: ndarray (3, 4) rows [ghost_x, ghost_y, ghost_radius, relative_energy]
		angle_x_deg: float
		angle_y_deg: float
		wave_name: str
		"""
		wave_name, _wave_length, nd, tf1, rf1, tf2, rf2, TD, RD = self._get_wave_tuple(detector_id)

		# Slope defined by pupil geometry
		u = self.Rp / self.Zp

		# Ghost radii (physical lengths)
		r1 = 2.0 * self.D1 * u / nd
		r2 = 2.0 * self.D2 * u
		r3 = 2.0 * (self.D1 / nd + self.D2) * u

		# Relative energy terms
		m1 = rf1 * rf2
		m2 = rf2 * RD
		m3 = rf1 * RD * (tf2 ** 2)

		# NOTE: Original script uses a fixed +0.313 x-offset; preserve behavior.
		denom = self.Zp - self.D1 * (1.0 - 1.0 / nd)
		angle_x = (float(x0) + 0.313) / denom
		angle_y = float(y0) / denom

		# Angles inside the filter
		af_x = angle_x / nd
		af_y = angle_y / nd

		# Lateral shifts on the focal plane
		xd1 = np.tan(af_x) * 2.0 * self.D1
		yd1 = np.tan(af_y) * 2.0 * self.D1
		xd2 = np.tan(angle_x) * 2.0 * self.D2
		yd2 = np.tan(angle_y) * 2.0 * self.D2

		# Combined shifts
		xd3 = xd1 + xd2
		yd3 = yd1 + yd2

		# Ghost positions
		x1, y1 = float(x0) + xd1, float(y0) + yd1
		x2, y2 = float(x0) + xd2, float(y0) + yd2
		x3, y3 = float(x0) + xd3, float(y0) + yd3

		ghosts = [
		[x1, y1, r1, m1],
		[x2, y2, r2, m2],
		[x3, y3, r3, m3],
		]

		return {
		"ghosts": ghosts,
		"angle_x_deg": float(angle_x * 180.0 / np.pi),
		"angle_y_deg": float(angle_y * 180.0 / np.pi),
		"wave_name": wave_name,
		}

		def draw_on_chip(self, chip, pos_pix, r_pix, flux=1.0):
		x_pix, y_pix = pos_pix.x, pos_pix.y
		profile = galsim.TopHat(radius=r_pix, flux=flux)
		chip.img.setOrigin(0, 0)
		profile.drawImage(
		image = chip.img,
		center = galsim.PositionD(x_pix, y_pix),
		add_to_image=True
		)
		chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)


		if __name__ == "__main__":
		g = Ghost()

		# Original sanity checks
		print(g.calculate(13, 0.0, 0.0))
		print(g.calculate(6, 0.19233, -0.24876))
		print(g.calculate(25, -0.19233, 0.24876))

observation_sim/mock_objects/init.py

+1 −0

Original line number	Diff line number	Diff line
		@@ -6,3 +6,4 @@ from .Star import Star
		from .Stamp import Stamp
		from .FlatLED import FlatLED
		from .ExtinctionMW import ExtinctionMW
		from .Ghost import Ghost

observation_sim/sim_steps/init.py

+2 −0

Original line number	Diff line number	Diff line
		@@ -17,10 +17,12 @@ class SimSteps:
		from .add_brighter_fatter_CTE import add_brighter_fatter, apply_CTE
		from .readout_output import add_prescan_overscan, add_readout_noise, apply_gain, quantization_and_output, add_crosstalk
		from .add_LED_flat import add_LED_Flat
		from .add_ghost import add_ghosts_SCI


		SIM_STEP_TYPES = {
		"scie_obs": "add_objects",
		"ghost_SCI": "add_ghosts_SCI",
		"sky_background": "add_sky_background",
		"cosmic_rays": "add_cosmic_rays",
		"PRNU_effect": "apply_PRNU",

observation_sim/sim_steps/add_ghost.py

0 → 100644

+85 −0

Original line number	Diff line number	Diff line
		import galsim

		from observation_sim.mock_objects import Ghost
		from observation_sim.psf import FieldDistortion

		def add_ghosts_SCI(self, chip, filt, tel, pointing, catalog, obs_param):
		# Load catalogues
		if catalog is None:
		self.chip_output.Log_error(
		"Catalog interface class must be specified for SCIE-OBS")
		raise ValueError(
		"Catalog interface class must be specified for SCIE-OBS")
		cat = catalog(config=self.overall_config, chip=chip,
		pointing=pointing, chip_output=self.chip_output, filt=filt)

		# Apply field distortion model
		if obs_param["field_dist"] is True:
		fd_model = FieldDistortion(chip=chip, img_rot=pointing.img_pa.deg)
		else:
		fd_model = None

		# Get chip WCS
		if not hasattr(self, 'h_ext'):
		_, _ = self.prepare_headers(chip=chip, pointing=pointing)

		chip_wcs = galsim.FitsWCS(header=self.h_ext)
		ghost_model = Ghost()

		# Loop over objects
		for j in range(len(cat.objs)):
		obj = cat.objs[j]

		try:
		sed_data = cat.load_sed(obj)
		norm_filt = cat.load_norm_filt(obj)

		obj.sed, obj.param["mag_%s" % filt.filter_type.lower()], obj.param["flux_%s" % filt.filter_type.lower()] = cat.convert_sed(
		mag=obj.param["mag_use_normal"],
		sed=sed_data,
		target_filt=filt,
		norm_filt=norm_filt,
		mu=obj.mu
		)
		except Exception as e:
		traceback.print_exc()
		self.chip_output.Log_error(e)
		continue

		# Select only bright stars
		if obj.type != 'star' or obj.getMagFilter(filt) >= obs_param["mag_threshold"]:
		continue

		# Get position of object on the focal plane
		pos_img, _, _, _, fd_shear = obj.getPosImg_Offset_WCS(
		img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext, ra_offset=self.ra_offset, dec_offset=self.dec_offset)

		# [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane
		if pos_img is None:
		self.chip_output.Log_info('obj_ra = %.6f, obj_dec = %.6f, obj_ra_orig = %.6f, obj_dec_orig = %.6f' % (
		obj.ra, obj.dec, obj.ra_orig, obj.dec_orig))
		self.chip_output.Log_error("Object missed: %s" % (obj.id))
		obj.unload_SED()
		continue

		# Get number of photons for total photons for the object
		nphotons_tot = obj.getElectronFluxFilt(filt, tel, pointing.exp_time)

		x_m = pos_img.x * chip.pix_size * 1e-3
		y_m = pos_img.y * chip.pix_size * 1e-3

		ghost_list = ghost_model.calculate(chip.chipID, x_m, y_m)["ghosts"]

		for ghost in ghost_list:
		x, y, r, ratio = ghost
		factor = (1e-3) * chip.pix_size
		x_pix = x / factor
		y_pix = y / factor
		r_pix = r / factor
		ghost_pos = obj.getRealPos(chip.img, global_x=x_pix, global_y=y_pix,
		img_real_wcs=obj.chip_wcs)
		# print(f"star at {obj.getRealPos(chip.img, global_x=pos_img.x, global_y=pos_img.y, img_real_wcs=obj.chip_wcs)}")
		# print("ghost info: ", ghost_pos.x, ghost_pos.y, r_pix, ratio * nphotons_tot)
		ghost_model.draw_on_chip(chip, ghost_pos, r_pix, flux=ratio * nphotons_tot)

		return chip, filt, tel, pointing
		No newline at end of file