Chip.py 48.4 KB
Newer Older
Fang Yuedong's avatar
Fang Yuedong committed
1
2
3
import galsim
import os
import numpy as np
Fang Yuedong's avatar
Fang Yuedong committed
4
5
import pickle
import json
6
import ObservationSim.Instrument._util as _util
Fang Yuedong's avatar
Fang Yuedong committed
7
8
from astropy.table import Table
from numpy.random import Generator, PCG64
Fang Yuedong's avatar
Fang Yuedong committed
9
10
from astropy.io import fits
from datetime import datetime
Fang Yuedong's avatar
Fang Yuedong committed
11

12
13
14
from ObservationSim.Instrument.Chip import Effects as effects
from ObservationSim.Instrument.FocalPlane import FocalPlane
from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
Fang Yuedong's avatar
Fang Yuedong committed
15
from ObservationSim.Instrument._util import rotate_conterclockwise
16
from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
17

Wei Chengliang's avatar
Wei Chengliang committed
18
19
from ObservationSim.Instrument.Chip.libCTI.CTI_modeling import CTI_sim

Fang Yuedong's avatar
Fang Yuedong committed
20
21
22
23
24
25
try:
    import importlib.resources as pkg_resources
except ImportError:
    # Try backported to PY<37 'importlib_resources'
    import importlib_resources as pkg_resources

Fang Yuedong's avatar
Fang Yuedong committed
26
class Chip(FocalPlane):
27
    def __init__(self, chipID, ccdEffCurve_dir=None, CRdata_dir=None, sls_dir=None, config=None, treering_func=None, logger=None):
Fang Yuedong's avatar
Fang Yuedong committed
28
29
        # Get focal plane (instance of paraent class) info
        super().__init__()
30
31
32
        self.nsecy = 2
        self.nsecx = 8
        self.gain_channel = np.ones(self.nsecy* self.nsecx)
Fang Yuedong's avatar
Fang Yuedong committed
33
        self._set_attributes_from_config(config)
Fang Yuedong's avatar
Fang Yuedong committed
34

35
36
        self.logger = logger

Fang Yuedong's avatar
Fang Yuedong committed
37
38
        # A chip ID must be assigned
        self.chipID = int(chipID)
Fang Yuedong's avatar
Fang Yuedong committed
39
        self.chip_name = str(chipID).rjust(2, '0')
Fang Yuedong's avatar
Fang Yuedong committed
40
41
42
43
44

        # Get corresponding filter info
        self.filter_id, self.filter_type = self.getChipFilter()
        self.survey_type = self._getSurveyType()

Fang Yuedong's avatar
Fang Yuedong committed
45
46
47
        if self.filter_type != "FGS":
            self._getChipRowCol()

48
        # Set the relavent specs for detectors
Fang Yuedong's avatar
Fang Yuedong committed
49
50
51
52
53
54
55
56
57
58
        try:
            with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath("chip_definition.json") as chip_definition:
                with open(chip_definition, "r") as f:
                    chip_dict = json.load(f)[str(self.chipID)]
        except AttributeError:
            with pkg_resources.path('ObservationSim.Instrument.data.ccd', "chip_definition.json") as chip_definition:
                with open(chip_definition, "r") as f:
                    chip_dict = json.load(f)[str(self.chipID)]
        for key in chip_dict:
            setattr(self, key, chip_dict[key])
59
        
Fang Yuedong's avatar
Fang Yuedong committed
60
        if self.filter_type == "FGS":
61
            if ("field_dist" in config) and (config["ins_effects"]["field_dist"]) == False:
Fang Yuedong's avatar
Fang Yuedong committed
62
63
64
65
66
67
68
69
70
71
72
73
74
                self.fdModel = None
            else:
                fgs_name = self.chip_name[0:4]
                try:
                    with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_pr4_%s.pickle"%(fgs_name.lower())) as field_distortion:
                        with open(field_distortion, "rb") as f:
                            self.fdModel = pickle.load(f)
                except AttributeError:
                    with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_pr4_%s.pickle"%(fgs_name.lower())) as field_distortion:
                        with open(field_distortion, "rb") as f:
                            self.fdModel = pickle.load(f)
        else:
            # Get the corresponding field distortion model
75
            if ("field_dist" in config) and (config["ins_effects"]["field_dist"] == False):
Fang Yuedong's avatar
Fang Yuedong committed
76
77
78
                self.fdModel = None
            else:
                try:
79
                    with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModel_v2.0.pickle") as field_distortion:
Fang Yuedong's avatar
Fang Yuedong committed
80
81
82
83
84
85
86
                        with open(field_distortion, "rb") as f:
                            self.fdModel = pickle.load(f)
                except AttributeError:
                    with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
                        with open(field_distortion, "rb") as f:
                            self.fdModel = pickle.load(f)

Fang Yuedong's avatar
Fang Yuedong committed
87
88
        # Get boundary (in pix)
        self.bound = self.getChipLim()
89
        
Fang Yuedong's avatar
Fang Yuedong committed
90
91
        self.ccdEffCurve_dir = ccdEffCurve_dir
        self.CRdata_dir = CRdata_dir
92
93
        
        slsconfs = chip_utils.getChipSLSConf(chipID=self.chipID)
Fang Yuedong's avatar
Fang Yuedong committed
94
        if np.size(slsconfs) == 1:
95
96
97
98
99
100
            try:
                with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs) as conf_path:
                    self.sls_conf = str(conf_path)
            except AttributeError:
                with pkg_resources.path('ObservationSim.Instrument.data.sls_conf', slsconfs) as conf_path:
                    self.sls_conf = str(conf_path)
Fang Yuedong's avatar
Fang Yuedong committed
101
        else:
Fang Yuedong's avatar
Fang Yuedong committed
102
103
            # self.sls_conf = [os.path.join(self.sls_dir, slsconfs[0]), os.path.join(self.sls_dir, slsconfs[1])]
            self.sls_conf = []
104
105
106
107
108
109
110
111
112
113
114
115
            try:
                with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs[0]) as conf_path:
                    self.sls_conf.append(str(conf_path))
            except AttributeError:
                with pkg_resources.path('ObservationSim.Instrument.data.sls_conf', slsconfs[0]) as conf_path:
                    self.sls_conf.append(str(conf_path))
            try:
                with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs[1]) as conf_path:
                    self.sls_conf.append(str(conf_path))
            except AttributeError:
                with pkg_resources.path('ObservationSim.Instrument.data.sls_conf', slsconfs[1]) as conf_path:
                    self.sls_conf.append(str(conf_path))
Fang Yuedong's avatar
Fang Yuedong committed
116
117
118
119
        
        self.effCurve = self._getChipEffCurve(self.filter_type)
        self._getCRdata()

Fang Yuedong's avatar
Fang Yuedong committed
120
        # Define the sensor model
121
        if "bright_fatter" in config["ins_effects"] and config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric":
Fang Yuedong's avatar
Fang Yuedong committed
122
            self.sensor = galsim.SiliconSensor(strength=self.df_strength, treering_func=treering_func)
Fang Yuedong's avatar
Fang Yuedong committed
123
124
125
        else:
            self.sensor = galsim.Sensor()

126
127
        self.flat_cube = None # for spectroscopic flat field cube simulation

Fang Yuedong's avatar
Fang Yuedong committed
128
129
130
131
132
133
134
135
136
137
138
    def _set_attributes_from_config(self, config):
        # Default setting
        self.read_noise = 5.0   # e/pix
        self.dark_noise = 0.02  # e/pix/s
        self.rotate_angle = 0.
        self.overscan   = 1000
        # Override default values
        for key in ["gain", "bias_level, dark_exptime", "flat_exptime", "readout_time", "full_well", "read_noise", "dark_noise", "overscan"]:
            if key in config["ins_effects"]:
                setattr(self, key, config["ins_effects"][key])

Fang Yuedong's avatar
Fang Yuedong committed
139
140
141
142
143
144
145
146
147
    def _getChipRowCol(self):
        self.rowID, self.colID = self.getChipRowCol(self.chipID)

    def getChipRowCol(self, chipID):
        rowID = ((chipID - 1) % 5) + 1
        colID = 6 - ((chipID - 1) // 5)
        return rowID, colID

    def _getSurveyType(self):
148
        if self.filter_type in _util.SPEC_FILTERS:
Fang Yuedong's avatar
Fang Yuedong committed
149
            return "spectroscopic"
150
        elif self.filter_type in _util.PHOT_FILTERS:
Fang Yuedong's avatar
Fang Yuedong committed
151
            return "photometric"
Fang Yuedong's avatar
Fang Yuedong committed
152
153
        # elif self.filter_type in ["FGS"]:
        #     return "FGS"
Fang Yuedong's avatar
Fang Yuedong committed
154
155
156

    def _getChipEffCurve(self, filter_type):
        # CCD efficiency curves
157
        if filter_type in ['NUV', 'u', 'GU']: filename = 'UV0.txt'
Fang Yuedong's avatar
Fang Yuedong committed
158
        if filter_type in ['g', 'r', 'GV', 'FGS']: filename = 'Astro_MB.txt' # TODO, need to switch to the right efficiency curvey for FGS CMOS
Fang Yuedong's avatar
Fang Yuedong committed
159
        if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt'
160
161
162
163
164
165
        try:
            with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath(filename) as ccd_path:
                table = Table.read(ccd_path, format='ascii')
        except AttributeError:
            with pkg_resources.path('ObservationSim.Instrument.data.ccd', filename) as ccd_path:
                table = Table.read(ccd_path, format='ascii')
166
        throughput = galsim.LookupTable(x=table['col1'], f=table['col2'], interpolant='linear')
Fang Yuedong's avatar
Fang Yuedong committed
167
168
169
170
        bandpass = galsim.Bandpass(throughput, wave_type='nm')
        return bandpass

    def _getCRdata(self):
171
172
173
174
175
176
        try:
            with pkg_resources.files('ObservationSim.Instrument.data').joinpath("wfc-cr-attachpixel.dat") as cr_path:
                self.attachedSizes = np.loadtxt(cr_path)
        except AttributeError:
            with pkg_resources.path('ObservationSim.Instrument.data', "wfc-cr-attachpixel.dat") as cr_path:
                self.attachedSizes = np.loadtxt(cr_path)
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
    
    def loadSLSFLATCUBE(self, flat_fn='flat_cube.fits'):
        try:
            with pkg_resources.files('ObservationSim.Instrument.data').joinpath(flat_fn) as data_path:
                flat_fits = fits.open(data_path, ignore_missing_simple=True)
        except AttributeError:
            with pkg_resources.path('ObservationSim.Instrument.data', flat_fn) as data_path:
                flat_fits = fits.open(data_path, ignore_missing_simple=True)

        fl = len(flat_fits)
        fl_sh = flat_fits[0].data.shape
        assert fl == 4, 'FLAT Field Cube is Not 4 layess!!!!!!!'
        self.flat_cube = np.zeros([fl, fl_sh[0], fl_sh[1]])
        for i in np.arange(0, fl, 1):
            self.flat_cube[i] = flat_fits[i].data
Fang Yuedong's avatar
Fang Yuedong committed
192

193
    def getChipFilter(self, chipID=None):
Fang Yuedong's avatar
Fang Yuedong committed
194
195
        """Return the filter index and type for a given chip #(chipID)
        """
196
        filter_type_list = _util.ALL_FILTERS
Fang Yuedong's avatar
Fang Yuedong committed
197
198
199
200
        if chipID == None:
            chipID = self.chipID

        # updated configurations
Fang Yuedong's avatar
Fang Yuedong committed
201
202
        if chipID>42 or chipID<1: raise ValueError("!!! Chip ID: [1,42]")
        if chipID in [6, 15, 16, 25]:  filter_type = "y"
Fang Yuedong's avatar
Fang Yuedong committed
203
        if chipID in [11, 20]:         filter_type = "z"
Fang Yuedong's avatar
Fang Yuedong committed
204
        if chipID in [7, 24]:          filter_type = "i"
Fang Yuedong's avatar
Fang Yuedong committed
205
        if chipID in [14, 17]:         filter_type = "u"
Fang Yuedong's avatar
Fang Yuedong committed
206
        if chipID in [9, 22]:          filter_type = "r"
207
        if chipID in [12, 13, 18, 19]: filter_type = "NUV"
Fang Yuedong's avatar
Fang Yuedong committed
208
209
210
211
212
        if chipID in [8, 23]:          filter_type = "g"
        if chipID in [1, 10, 21, 30]:  filter_type = "GI"
        if chipID in [2, 5, 26, 29]:   filter_type = "GV"
        if chipID in [3, 4, 27, 28]:   filter_type = "GU"
        if chipID in range(31, 43):    filter_type = 'FGS'
Fang Yuedong's avatar
Fang Yuedong committed
213
214
215
216
217
218
        filter_id = filter_type_list.index(filter_type)

        return filter_id, filter_type

    def getChipLim(self, chipID=None):
        """Calculate the edges in pixel for a given CCD chip on the focal plane
Fang Yuedong's avatar
Fang Yuedong committed
219
        NOTE: There are 5*4 CCD chips in the focus plane for photometric / spectroscopic observation.
Fang Yuedong's avatar
Fang Yuedong committed
220
221
222
223
224
225
        Parameters:
            chipID:         int
                            the index of the chip
        Returns:
            A galsim BoundsD object
        """
226
227
228
229
230
231
232
233
        xmin, xmax, ymin, ymax = 1e10, -1e10, 1e10, -1e10
        xcen = self.x_cen / self.pix_size
        ycen = self.y_cen / self.pix_size
        sign_x = [-1., 1., -1., 1.]
        sign_y = [-1., -1., 1., 1.]
        for i in range(4):
            x = xcen + sign_x[i] * self.npix_x / 2.
            y = ycen + sign_y[i] * self.npix_y / 2.
234
            x, y = _util.rotate_conterclockwise(x0=xcen, y0=ycen, x=x, y=y, angle=self.rotate_angle)
235
236
237
            xmin, xmax = min(xmin, x), max(xmax, x)
            ymin, ymax = min(ymin, y), max(ymax, y)
        return galsim.BoundsD(xmin, xmax, ymin, ymax)
Fang Yuedong's avatar
Fang Yuedong committed
238
239
240
241



    def getSkyCoverage(self, wcs):
Fang Yuedong's avatar
Fang Yuedong committed
242
        # print("In getSkyCoverage: xmin = %.3f, xmax = %.3f, ymin = %.3f, ymax = %.3f"%(self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax))
Fang Yuedong's avatar
Fang Yuedong committed
243
244
245
246
247
248
249
250
251
252
        return super().getSkyCoverage(wcs, self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax)


    def getSkyCoverageEnlarged(self, wcs, margin=0.5):
        """The enlarged sky coverage of the chip
        """
        margin /= 60.0
        bound = self.getSkyCoverage(wcs)
        return galsim.BoundsD(bound.xmin - margin, bound.xmax + margin, bound.ymin - margin, bound.ymax + margin)

Fang Yuedong's avatar
Fang Yuedong committed
253
    def isContainObj(self, ra_obj=None, dec_obj=None, x_image=None, y_image=None, wcs=None, margin=1):
Fang Yuedong's avatar
Fang Yuedong committed
254
        # magin in number of pix
Fang Yuedong's avatar
Fang Yuedong committed
255
256
257
258
259
260
261
262
        if (ra_obj is not None) and (dec_obj is not None):
            if wcs is None:
                wcs = self.img.wcs
            pos_obj = wcs.toImage(galsim.CelestialCoord(ra=ra_obj*galsim.degrees, dec=dec_obj*galsim.degrees))
            x_image, y_image = pos_obj.x, pos_obj.y
        elif (x_image is None) or (y_image is None):
            raise ValueError("Either (ra_obj, dec_obj) or (x_image, y_image) should be given")

Fang Yuedong's avatar
Fang Yuedong committed
263
264
        xmin, xmax = self.bound.xmin - margin, self.bound.xmax + margin
        ymin, ymax = self.bound.ymin - margin, self.bound.ymax + margin
Fang Yuedong's avatar
Fang Yuedong committed
265
        if (x_image - xmin) * (x_image - xmax) > 0.0 or (y_image - ymin) * (y_image - ymax) > 0.0:
Fang Yuedong's avatar
Fang Yuedong committed
266
267
268
269
270
271
272
            return False
        return True

    def getChipNoise(self, exptime=150.0):
        noise = self.dark_noise * exptime + self.read_noise**2
        return noise

Zhang Xin's avatar
Zhang Xin committed
273
    def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='SCI', sky_map=None, post_flash_map=None, tel=None, logger=None):
274
        # Set random seeds
Fang Yuedong's avatar
Fang Yuedong committed
275
276
277
278
279
280
        SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"])
        SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"])
        SeedRnNonuni = int(config["random_seeds"]["seed_rnNonUniform"])
        SeedBadColumns = int(config["random_seeds"]["seed_badcolumns"])
        SeedDefective = int(config["random_seeds"]["seed_defective"])
        SeedCosmicRay = int(config["random_seeds"]["seed_CR"])
Fang Yuedong's avatar
Fang Yuedong committed
281
        fullwell = int(self.full_well)
Fang Yuedong's avatar
Fang Yuedong committed
282
        if config["ins_effects"]["add_hotpixels"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
283
284
285
            BoolHotPix = True
        else:
            BoolHotPix = False
Fang Yuedong's avatar
Fang Yuedong committed
286
        if config["ins_effects"]["add_deadpixels"]== True:
Fang Yuedong's avatar
Fang Yuedong committed
287
288
289
            BoolDeadPix = True
        else:
            BoolDeadPix = False
290
        self.logger = logger
Fang Yuedong's avatar
Fang Yuedong committed
291

292
        # Get Poisson noise generator
293
294
        rng_poisson, poisson_noise = chip_utils.get_poisson(
            seed=int(config["random_seeds"]["seed_poisson"]) + pointing_ID*30 + self.chipID, sky_level=0.)
295

296
        # Add sky background
Fang Yuedong's avatar
Fang Yuedong committed
297
        if config["ins_effects"]["add_back"] == True:
298
299
            img, sky_map = chip_utils.add_sky_background(img=img, filt=filt, exptime=exptime, sky_map=sky_map, tel=tel)
            del sky_map
300

Fang Yuedong's avatar
Fang Yuedong committed
301
        # Apply flat-field large scale structure for one chip
Fang Yuedong's avatar
Fang Yuedong committed
302
        if config["ins_effects"]["flat_fielding"] == True:
303
304
305
            chip_utils.log_info(msg="  Creating and applying Flat-Fielding", logger=self.logger)
            chip_utils.log_info(msg=str(img.bounds), logger=self.logger)
            flat_img, flat_normal = chip_utils.get_flat(img=img, seed=int(config["random_seeds"]["seed_flat"]))
306
307
            if self.survey_type == "photometric":
                img *= flat_normal
Fang Yuedong's avatar
Fang Yuedong committed
308
            del flat_normal
Fang Yuedong's avatar
Fang Yuedong committed
309
            if config["output_setting"]["flat_output"] == False:
Fang Yuedong's avatar
Fang Yuedong committed
310
311
                del flat_img

Zhang Xin's avatar
Zhang Xin committed
312
313
314
        if post_flash_map is not None:
            img = img + post_flash_map

Fang Yuedong's avatar
Fang Yuedong committed
315
        # Apply Shutter-effect for one chip
Fang Yuedong's avatar
Fang Yuedong committed
316
        if config["ins_effects"]["shutter_effect"] == True:
317
            chip_utils.log_info(msg="  Apply shutter effect", logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
318
            shuttimg = effects.ShutterEffectArr(img, t_shutter=1.3, dist_bearing=735, dt=1E-3)    # shutter effect normalized image for this chip
319
320
            if self.survey_type == "photometric":
                img *= shuttimg
Fang Yuedong's avatar
Fang Yuedong committed
321
            if config["output_setting"]["shutter_output"] == True:    # output 16-bit shutter effect image with pixel value <=65535
Fang Yuedong's avatar
Fang Yuedong committed
322
323
324
325
                shutt_gsimg = galsim.ImageUS(shuttimg*6E4)
                shutt_gsimg.write("%s/ShutterEffect_%s_1.fits" % (chip_output.subdir, self.chipID))
                del shutt_gsimg
            del shuttimg
326
327
328
329
330
        # # Add Poisson noise to the resulting images
        # # (NOTE): this can only applied to the slitless image
        # # since it dose not use photon shooting to draw stamps
        # if self.survey_type == "spectroscopic":
        #     img.addNoise(poisson_noise)
Fang Yuedong's avatar
Fang Yuedong committed
331
332

        # Add cosmic-rays
333
        if config["ins_effects"]["cosmic_ray"] == True and pointing_type=='SCI':
334
335
336
337
338
            chip_utils.log_info(msg="  Adding Cosmic-Ray", logger=self.logger)
            img, crmap_gsimg, cr_event_num = chip_utils.add_cosmic_rays(img=img, chip=self, exptime=exptime, 
                                                    seed=SeedCosmicRay+pointing_ID*30+self.chipID)
            chip_utils.outputCal(
                chip=self,
Fang Yuedong's avatar
Fang Yuedong committed
339
340
341
342
343
344
345
                img=crmap_gsimg,
                ra_cen=ra_cen,
                dec_cen=dec_cen,
                img_rot=img_rot,
                im_type='CRS',
                pointing_ID=pointing_ID,
                output_dir=chip_output.subdir,
346
                exptime=exptime,
347
348
                project_cycle=config["project_cycle"],
                run_counter=config["run_counter"],
349
                timestamp=timestamp_obs)
Fang Yuedong's avatar
Fang Yuedong committed
350
351
            del crmap_gsimg

352
        # Apply PRNU effect and output PRNU flat file:
Fang Yuedong's avatar
Fang Yuedong committed
353
        if config["ins_effects"]["prnu_effect"] == True:
354
355
356
            chip_utils.log_info(msg="  Applying PRNU effect", logger=self.logger)
            img, prnu_img = chip_utils.add_PRNU(img=img, chip=self, 
                                seed=int(config["random_seeds"]["seed_prnu"]+self.chipID))
Fang Yuedong's avatar
Fang Yuedong committed
357
            if config["output_setting"]["prnu_output"] == True:
358
                prnu_img.write("%s/FlatImg_PRNU_%s.fits" % (chip_output.subdir,self.chipID))
Fang Yuedong's avatar
Fang Yuedong committed
359
            if config["output_setting"]["flat_output"] == False:
360
361
                del prnu_img

362
363
364
365
366
367
        # # Add dark current
        # if config["ins_effects"]["add_dark"] == True:
        #     dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(exptime+0.5*self.readout_time)))
        #     img.addNoise(dark_noise)

        # Add dark current & Poisson noise
368
        InputDark = False
Fang Yuedong's avatar
Fang Yuedong committed
369
        if config["ins_effects"]["add_dark"] == True:
370
371
372
373
            if InputDark:
                img = chip_utils.add_inputdark(img=img, chip=self, exptime=exptime)
            else:
                img, _ = chip_utils.add_poisson(img=img, chip=self, exptime=exptime, poisson_noise=poisson_noise)
374
375
376
377
378
379
        else:
            img, _ = chip_utils.add_poisson(img=img, chip=self, exptime=exptime, poisson_noise=poisson_noise, dark_noise=0.)

        # Add diffusion & brighter-fatter effects
        if config["ins_effects"]["bright_fatter"] == True:
            img = chip_utils.add_brighter_fatter(img=img)
380
381
382
383
384
385
386

        # Add Hot Pixels or/and Dead Pixels
        rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
        badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
        img = effects.DefectivePixels(img, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)

        # Apply Bad lines 
Fang Yuedong's avatar
Fang Yuedong committed
387
        if config["ins_effects"]["add_badcolumns"] == True:
388
            img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
389
390

        # Apply Nonlinearity on the chip image
Fang Yuedong's avatar
Fang Yuedong committed
391
        if config["ins_effects"]["non_linear"] == True:
392
            chip_utils.log_info(msg="  Applying Non-Linearity on the chip image", logger=self.logger)
393
394
395
            img = effects.NonLinearity(GSImage=img, beta1=5.e-7, beta2=0)

        # Apply CCD Saturation & Blooming
Fang Yuedong's avatar
Fang Yuedong committed
396
        if config["ins_effects"]["saturbloom"] == True:
397
            chip_utils.log_info(msg="  Applying CCD Saturation & Blooming", logger=self.logger)
398
399
400
            img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell)

        # Apply CTE Effect
Wei Chengliang's avatar
Wei Chengliang committed
401
402
403
404
405
406
407
408
409
        ###if config["ins_effects"]["cte_trail"] == True:
        ###    chip_utils.log_info(msg="  Apply CTE Effect", logger=self.logger)
        ###    img = effects.CTE_Effect(GSImage=img, threshold=27)

        pre1 = self.prescan_x  #27
        over1= self.overscan_x #71
        pre2 = self.prescan_y  #0 #4
        over2= self.overscan_y #84 #80

Fang Yuedong's avatar
Fang Yuedong committed
410
        if config["ins_effects"]["cte_trail"] == True:
411
            chip_utils.log_info(msg="  Apply CTE Effect", logger=self.logger)
Wei Chengliang's avatar
Wei Chengliang committed
412
413
414
415
416
417
418
419
420
421
422
423
424
            ###img = effects.CTE_Effect(GSImage=img, threshold=27)
            ###CTI_modeling
            ### 2*8 -> 1*16 img-layout
            img = chip_utils.formatOutput(GSImage=img)
            self.nsecy = 1
            self.nsecx = 16

            img_arr = img.array
            ny, nx = img_arr.shape
            dx = int(nx/self.nsecx)
            dy = int(ny/self.nsecy)
            newimg = galsim.Image(nx, int(ny+over2), init_value=0)
            for ichannel in range(16):
Wei Chengliang's avatar
Wei Chengliang committed
425
                print('\n***add CTI effects: pointing-{:} chip-{:} channel-{:}***'.format(pointing_ID, self.chipID, ichannel+1))
Wei Chengliang's avatar
Wei Chengliang committed
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
                #nx,ny,noverscan,nsp,nmax = 4608,4616,84,3,10
                noverscan,nsp,nmax = over2,3,10
                beta,w,c = 0.478,84700,0
                t = np.array([0.74,7.7,37],dtype=np.float32)
                rho_trap = np.array([0.6,1.6,1.4],dtype=np.float32)
                trap_seeds = np.array([0,1000,10000],dtype=np.int32) + ichannel + self.chipID*16
                release_seed = 50 + ichannel + pointing_ID*30  + self.chipID*16
                newimg.array[:, 0+ichannel*dx:dx+ichannel*dx] = CTI_sim(img_arr[:, 0+ichannel*dx:dx+ichannel*dx],dx,dy,noverscan,nsp,nmax,beta,w,c,t,rho_trap,trap_seeds,release_seed)
            newimg.wcs = img.wcs
            del img
            img = newimg

            ### 1*16 -> 2*8 img-layout
            img = chip_utils.formatRevert(GSImage=img)
            self.nsecy = 2
            self.nsecx = 8
442
443
444

        ### prescan & overscan
        if config["ins_effects"]["add_prescan"] == True:
Wei Chengliang's avatar
Wei Chengliang committed
445
446
447
448
449
450
            chip_utils.log_info(msg="  Apply pre/over-scan", logger=self.logger)
            if config["ins_effects"]["cte_trail"] == False:
                img = chip_utils.AddPreScan(GSImage=img, pre1=pre1, pre2=pre2, over1=over1, over2=over2)
            if config["ins_effects"]["cte_trail"] == True:
                img = chip_utils.AddPreScan(GSImage=img, pre1=pre1, pre2=pre2, over1=over1, over2=0)

451
452
        ### 1*16 output
        if config["ins_effects"]["format_output"] == True:
Wei Chengliang's avatar
Wei Chengliang committed
453
            chip_utils.log_info(msg="  Apply 1*16 format", logger=self.logger)
454
455
456
457
458
            img = chip_utils.formatOutput(GSImage=img)
            self.nsecy = 1
            self.nsecx = 16
                    

459
460
        # Add Bias level
        if config["ins_effects"]["add_bias"] == True:
461
            chip_utils.log_info(msg="  Adding Bias level and 16-channel non-uniformity", logger=self.logger)
462
463
464
            if config["ins_effects"]["bias_16channel"] == True:
                img = effects.AddBiasNonUniform16(img, 
                    bias_level=float(self.bias_level), 
465
                    nsecy = self.nsecy, nsecx=self.nsecx, 
466
467
468
469
                    seed=SeedBiasNonuni+self.chipID,
                    logger=self.logger)
            elif config["ins_effects"]["bias_16channel"] == False:
                img += self.bias_level
470
471

        # Add Read-out Noise
Fang Yuedong's avatar
Fang Yuedong committed
472
473
        if config["ins_effects"]["add_readout"] == True:
            seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID
474
475
476
477
478
            rng_readout = galsim.BaseDeviate(seed)
            readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
            img.addNoise(readout_noise)

        # Apply Gain & Quantization
479
        chip_utils.log_info(msg="  Applying Gain (and 16 channel non-uniformity) & Quantization", logger=self.logger)
480
        if config["ins_effects"]["gain_16channel"] == True:
481
            img, self.gain_channel = effects.ApplyGainNonUniform16(
482
                img, gain=self.gain, 
483
                nsecy = self.nsecy, nsecx=self.nsecx, 
484
485
486
487
488
                seed=SeedGainNonuni+self.chipID,
                logger=self.logger)
        elif config["ins_effects"]["gain_16channel"] == False:
            img /= self.gain
            
489
490
491
492
493
494
495
        img.array[img.array > 65535] = 65535
        img.replaceNegative(replace_value=0)
        img.quantize()

        ######################################################################################
        # Output images for calibration pointing
        ######################################################################################
Fang Yuedong's avatar
Fang Yuedong committed
496
        # Bias output
497
        if config["ins_effects"]["add_bias"] == True and config["output_setting"]["bias_output"] == True and pointing_type=='CAL':
498
499
500
501
            if self.logger is not None:
                self.logger.info("  Output N frame Bias files")
            else:
                print("  Output N frame Bias files", flush=True)
Fang Yuedong's avatar
Fang Yuedong committed
502
            NBias = int(config["output_setting"]["NBias"])
Fang Yuedong's avatar
Fang Yuedong committed
503
504
505
            for i in range(NBias):
                BiasCombImg, BiasTag = effects.MakeBiasNcomb(
                    self.npix_x, self.npix_y, 
Fang Yuedong's avatar
Fang Yuedong committed
506
                    bias_level=float(self.bias_level), 
Fang Yuedong's avatar
Fang Yuedong committed
507
                    ncombine=1, read_noise=self.read_noise, gain=1,
508
509
                    seed=SeedBiasNonuni+self.chipID,
                    logger=self.logger)
510
                # Readout noise for Biases is not generated with random seeds. So readout noise for bias images can't be reproduced.
Fang Yuedong's avatar
Fang Yuedong committed
511
512
513
                if config["ins_effects"]["cosmic_ray"] == True:
                    if config["ins_effects"]["cray_differ"] == True:
                        cr_map, cr_event_num = effects.produceCR_Map(
Fang Yuedong's avatar
Fang Yuedong committed
514
515
                            xLen=self.npix_x, yLen=self.npix_y, 
                            exTime=0.01, 
Fang Yuedong's avatar
Fang Yuedong committed
516
                            cr_pixelRatio=0.003*(0.01+0.5*self.readout_time)/150., 
Fang Yuedong's avatar
Fang Yuedong committed
517
518
519
520
521
522
523
                            gain=self.gain, 
                            attachedSizes=self.attachedSizes,
                            seed=SeedCosmicRay+pointing_ID*30+self.chipID+1)
                            # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
                    BiasCombImg += cr_map
                    del cr_map

Fang Yuedong's avatar
Fang Yuedong committed
524
                # Non-Linearity for Bias
Fang Yuedong's avatar
Fang Yuedong committed
525
                if config["ins_effects"]["non_linear"] == True:
526
527
528
529
                    if self.logger is not None:
                        self.logger.info("  Applying Non-Linearity on the Bias image")
                    else:
                        print("  Applying Non-Linearity on the Bias image", flush=True)
530
531
                    BiasCombImg = effects.NonLinearity(GSImage=BiasCombImg, beta1=5.e-7, beta2=0)

Wei Chengliang's avatar
Wei Chengliang committed
532
533
534
535
536
537
538
539
540
541
542
543
544
545
                ###########################START
                ### prescan & overscan
                if config["ins_effects"]["add_prescan"] == True:
                    chip_utils.log_info(msg="  Apply pre/over-scan", logger=self.logger)
                    BiasCombImg = chip_utils.AddPreScan(GSImage=BiasCombImg, pre1=pre1, pre2=pre2, over1=over1, over2=over2)

                ### 1*16 output
                if config["ins_effects"]["format_output"] == True:
                    chip_utils.log_info(msg="  Apply 1*16 format", logger=self.logger)
                    BiasCombImg = chip_utils.formatOutput(GSImage=BiasCombImg)
                    self.nsecy = 1
                    self.nsecx = 16
                ###########################END

546
                # Apply Bad lines 
Fang Yuedong's avatar
Fang Yuedong committed
547
                if config["ins_effects"]["add_badcolumns"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
548
                    BiasCombImg = effects.BadColumns(BiasCombImg-float(self.bias_level)+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(self.bias_level)-5
Fang Yuedong's avatar
Fang Yuedong committed
549

550
                BiasCombImg, self.gain_channel = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain,
Wei Chengliang's avatar
Wei Chengliang committed
551
                    nsecy = self.nsecy, nsecx=self.nsecx, 
552
553
                    seed=SeedGainNonuni+self.chipID,
                    logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
554
555
                # BiasCombImg = effects.AddOverscan(
                #     BiasCombImg, 
Fang Yuedong's avatar
Fang Yuedong committed
556
                #     overscan=float(config["ins_effects"]["bias_level"])-2, gain=self.gain, 
Fang Yuedong's avatar
Fang Yuedong committed
557
558
559
560
                #     widthl=27, widthr=27, widtht=8, widthb=8)
                BiasCombImg.replaceNegative(replace_value=0)
                BiasCombImg.quantize()
                BiasCombImg = galsim.ImageUS(BiasCombImg)
561
                timestamp_obs += 10 * 60
562
563
                chip_utils.outputCal(
                    chip=self,
Fang Yuedong's avatar
Fang Yuedong committed
564
565
566
567
                    img=BiasCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
Xin Zhang's avatar
Xin Zhang committed
568
                    im_type='BIAS',
Fang Yuedong's avatar
Fang Yuedong committed
569
570
                    pointing_ID=pointing_ID,
                    output_dir=chip_output.subdir,
571
                    exptime=0.0,
572
573
                    project_cycle=config["project_cycle"],
                    run_counter=config["run_counter"],
574
                    timestamp=timestamp_obs)
Fang Yuedong's avatar
Fang Yuedong committed
575
576
577
            del BiasCombImg

        # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file
578
        if config["ins_effects"]["flat_fielding"] == True and config["output_setting"]["flat_output"] == True and pointing_type=='CAL':
579
580
581
582
            if self.logger is not None:
                self.logger.info("  Output N frame Flat-Field files")
            else:
                print("  Output N frame Flat-Field files", flush=True)
Fang Yuedong's avatar
Fang Yuedong committed
583
            NFlat = int(config["output_setting"]["NFlat"])
Fang Yuedong's avatar
Fang Yuedong committed
584
            if config["ins_effects"]["add_bias"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
585
586
                biaslevel = self.bias_level
                overscan = biaslevel-2
587
            elif config["ins_effects"]["add_bias"] == False:
Fang Yuedong's avatar
Fang Yuedong committed
588
589
                biaslevel = 0
                overscan = 0
590
            darklevel = self.dark_noise*(self.flat_exptime+0.5*self.readout_time)
Fang Yuedong's avatar
Fang Yuedong committed
591
592
593
594
595
596
597
598
            for i in range(NFlat):
                FlatSingle = flat_img * prnu_img + darklevel
                FlatCombImg, FlatTag = effects.MakeFlatNcomb(
                    flat_single_image=FlatSingle, 
                    ncombine=1, 
                    read_noise=self.read_noise,
                    gain=1, 
                    overscan=overscan, 
599
                    biaslevel=0,
600
601
                    seed_bias=SeedDefective+self.chipID,
                    logger=self.logger
Fang Yuedong's avatar
Fang Yuedong committed
602
                    )
Fang Yuedong's avatar
Fang Yuedong committed
603
604
605
                if config["ins_effects"]["cosmic_ray"] == True:
                    if config["ins_effects"]["cray_differ"] == True:
                        cr_map, cr_event_num = effects.produceCR_Map(
Fang Yuedong's avatar
Fang Yuedong committed
606
                            xLen=self.npix_x, yLen=self.npix_y, 
607
                            exTime=self.flat_exptime+0.5*self.readout_time, 
Fang Yuedong's avatar
Fang Yuedong committed
608
                            cr_pixelRatio=0.003*(self.flat_exptime+0.5*self.readout_time)/150., 
Fang Yuedong's avatar
Fang Yuedong committed
609
610
611
612
                            gain=self.gain, 
                            attachedSizes=self.attachedSizes,
                            seed=SeedCosmicRay+pointing_ID*30+self.chipID+3)
                            # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
Fang Yuedong's avatar
Fang Yuedong committed
613
                    FlatCombImg += cr_map
Fang Yuedong's avatar
Fang Yuedong committed
614
                    del cr_map
Fang Yuedong's avatar
Fang Yuedong committed
615

Fang Yuedong's avatar
Fang Yuedong committed
616
                if config["ins_effects"]["non_linear"] == True:
617
618
619
620
                    if self.logger is not None:
                        self.logger.info("  Applying Non-Linearity on the Flat image")
                    else:
                        print("  Applying Non-Linearity on the Flat image", flush=True)
621
                    FlatCombImg = effects.NonLinearity(GSImage=FlatCombImg, beta1=5.e-7, beta2=0)
Fang Yuedong's avatar
Fang Yuedong committed
622

Wei Chengliang's avatar
Wei Chengliang committed
623
624
625
626
627
628
629
                ###if config["ins_effects"]["cte_trail"] == True:
                ###    FlatCombImg = effects.CTE_Effect(GSImage=FlatCombImg, threshold=3)
                ###########################START
                pre1 = self.prescan_x  #27
                over1= self.overscan_x #71
                pre2 = self.prescan_y  #0 #4
                over2= self.overscan_y #84 #80
Fang Yuedong's avatar
Fang Yuedong committed
630
                if config["ins_effects"]["cte_trail"] == True:
Wei Chengliang's avatar
Wei Chengliang committed
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
                    chip_utils.log_info(msg="  Apply CTE Effect", logger=self.logger)
                    ###img = effects.CTE_Effect(GSImage=img, threshold=27)
                    ###CTI_modeling
                    ### 2*8 -> 1*16 img-layout
                    FlatCombImg = chip_utils.formatOutput(GSImage=FlatCombImg)
                    self.nsecy = 1
                    self.nsecx = 16

                    img_arr = FlatCombImg.array
                    ny, nx = img_arr.shape
                    dx = int(nx/self.nsecx)
                    dy = int(ny/self.nsecy)
                    newimg = galsim.Image(nx, int(ny+over2), init_value=0)
                    for ichannel in range(16):
                        print('\n***add CTI effects: pointing-{:} chip-{:} channel-{:}***'.format(pointing_ID, self.chipID, ichannel+1))
                        #nx,ny,noverscan,nsp,nmax = 4608,4616,84,3,10
                        noverscan,nsp,nmax = over2,3,10
                        beta,w,c = 0.478,84700,0
                        t = np.array([0.74,7.7,37],dtype=np.float32)
                        rho_trap = np.array([0.6,1.6,1.4],dtype=np.float32)
                        trap_seeds = np.array([0,1000,10000],dtype=np.int32) + ichannel + self.chipID*16
                        release_seed = 50 + ichannel + pointing_ID*30  + self.chipID*16
                        newimg.array[:, 0+ichannel*dx:dx+ichannel*dx] = CTI_sim(img_arr[:, 0+ichannel*dx:dx+ichannel*dx],dx,dy,noverscan,nsp,nmax,beta,w,c,t,rho_trap,trap_seeds,release_seed)
                    newimg.wcs = FlatCombImg.wcs
                    del FlatCombImg
                    FlatCombImg = newimg

                    ### 1*16 -> 2*8 img-layout
Wei Chengliang's avatar
Wei Chengliang committed
659
                    FlatCombImg = chip_utils.formatRevert(GSImage=FlatCombImg)
Wei Chengliang's avatar
Wei Chengliang committed
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
                    self.nsecy = 2
                    self.nsecx = 8

                ### prescan & overscan
                if config["ins_effects"]["add_prescan"] == True:
                    chip_utils.log_info(msg="  Apply pre/over-scan", logger=self.logger)
                    if config["ins_effects"]["cte_trail"] == False:
                        FlatCombImg = chip_utils.AddPreScan(GSImage=FlatCombImg, pre1=pre1, pre2=pre2, over1=over1, over2=over2)
                    if config["ins_effects"]["cte_trail"] == True:
                        FlatCombImg = chip_utils.AddPreScan(GSImage=FlatCombImg, pre1=pre1, pre2=pre2, over1=over1, over2=0)

                ### 1*16 output
                if config["ins_effects"]["format_output"] == True:
                    chip_utils.log_info(msg="  Apply 1*16 format", logger=self.logger)
                    FlatCombImg = chip_utils.formatOutput(GSImage=FlatCombImg)
                    self.nsecy = 1
                    self.nsecx = 16
                ###########################END
Fang Yuedong's avatar
Fang Yuedong committed
678

679
680
681
682
683
                # Add Hot Pixels or/and Dead Pixels
                rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
                badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
                FlatCombImg = effects.DefectivePixels(FlatCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)

Fang Yuedong's avatar
Fang Yuedong committed
684
                # Apply Bad lines 
Fang Yuedong's avatar
Fang Yuedong committed
685
                if config["ins_effects"]["add_badcolumns"] == True:
686
                    FlatCombImg = effects.BadColumns(FlatCombImg, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
687

688
                # Add Bias level
Fang Yuedong's avatar
Fang Yuedong committed
689
                if config["ins_effects"]["add_bias"] == True:
690
691
692
693
                    if self.logger is not None:
                        self.logger.info("  Adding Bias level and 16-channel non-uniformity")
                    else:
                        print("  Adding Bias level and 16-channel non-uniformity")
Fang Yuedong's avatar
Fang Yuedong committed
694
                    # img += float(config["ins_effects"]["bias_level"])
695
696
                    FlatCombImg = effects.AddBiasNonUniform16(FlatCombImg, 
                        bias_level=biaslevel, 
Wei Chengliang's avatar
Wei Chengliang committed
697
                        nsecy = self.nsecy, nsecx=self.nsecx, 
698
699
                        seed=SeedBiasNonuni+self.chipID,
                        logger=self.logger)
700
701
                
                # Add Read-out Noise
Fang Yuedong's avatar
Fang Yuedong committed
702
                if config["ins_effects"]["add_readout"] == True:
703
                    seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID + 3
704
705
706
                    rng_readout = galsim.BaseDeviate(seed)
                    readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
                    FlatCombImg.addNoise(readout_noise)
Fang Yuedong's avatar
Fang Yuedong committed
707

708
                FlatCombImg, self.gain_channel = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain,
Wei Chengliang's avatar
Wei Chengliang committed
709
                    nsecy = self.nsecy, nsecx=self.nsecx, 
710
711
                    seed=SeedGainNonuni+self.chipID,
                    logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
712
713
714
715
                # FlatCombImg = effects.AddOverscan(FlatCombImg, overscan=overscan, gain=self.gain, widthl=27, widthr=27, widtht=8, widthb=8)
                FlatCombImg.replaceNegative(replace_value=0)
                FlatCombImg.quantize()
                FlatCombImg = galsim.ImageUS(FlatCombImg)
716
                timestamp_obs += 10 * 60
717
718
                chip_utils.outputCal(
                    chip=self,
Fang Yuedong's avatar
Fang Yuedong committed
719
720
721
722
                    img=FlatCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
Xin Zhang's avatar
Xin Zhang committed
723
                    im_type='FLAT',
Fang Yuedong's avatar
Fang Yuedong committed
724
725
                    pointing_ID=pointing_ID,
                    output_dir=chip_output.subdir,
726
                    exptime=self.flat_exptime,
727
728
                    project_cycle=config["project_cycle"],
                    run_counter=config["run_counter"],
729
                    timestamp=timestamp_obs)
Fang Yuedong's avatar
Fang Yuedong committed
730

Fang Yuedong's avatar
Fang Yuedong committed
731
732
733
734
735
736
737
            del FlatCombImg, FlatSingle, prnu_img
            # flat_img.replaceNegative(replace_value=0)
            # flat_img.quantize()
            # galsim.ImageUS(flat_img).write("%s/FlatImg_Vignette_%s.fits" % (chip_output.subdir, self.chipID))
            del flat_img

        # Export Dark current images
738
        if config["ins_effects"]["add_dark"] == True and config["output_setting"]["dark_output"] == True and pointing_type=='CAL':
739
740
741
742
            if self.logger is not None:
                self.logger.info("  Output N frame Dark Current files")
            else:
                print("  Output N frame Dark Current files", flush=True)
Fang Yuedong's avatar
Fang Yuedong committed
743
            NDark = int(config["output_setting"]["NDark"])
Fang Yuedong's avatar
Fang Yuedong committed
744
            if config["ins_effects"]["add_bias"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
745
746
                biaslevel = self.bias_level
                overscan = biaslevel-2
747
            elif config["ins_effects"]["add_bias"] == False:
Fang Yuedong's avatar
Fang Yuedong committed
748
749
750
751
752
                biaslevel = 0
                overscan = 0
            for i in range(NDark):
                DarkCombImg, DarkTag = effects.MakeDarkNcomb(
                    self.npix_x, self.npix_y, 
753
                    overscan=overscan, bias_level=0, darkpsec=0.02, exptime=self.dark_exptime+0.5*self.readout_time,
Fang Yuedong's avatar
Fang Yuedong committed
754
                    ncombine=1, read_noise=self.read_noise, 
755
756
                    gain=1, seed_bias=SeedBiasNonuni+self.chipID,
                    logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
757
758
759
                if config["ins_effects"]["cosmic_ray"] == True:
                    if config["ins_effects"]["cray_differ"] == True:
                        cr_map, cr_event_num = effects.produceCR_Map(
Fang Yuedong's avatar
Fang Yuedong committed
760
                            xLen=self.npix_x, yLen=self.npix_y, 
761
                            exTime=self.dark_exptime+0.5*self.readout_time, 
Fang Yuedong's avatar
Fang Yuedong committed
762
                            cr_pixelRatio=0.003*(self.dark_exptime+0.5*self.readout_time)/150., 
Fang Yuedong's avatar
Fang Yuedong committed
763
764
765
766
                            gain=self.gain, 
                            attachedSizes=self.attachedSizes,
                            seed=SeedCosmicRay+pointing_ID*30+self.chipID+2)
                            # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
Fang Yuedong's avatar
Fang Yuedong committed
767
                    DarkCombImg += cr_map
Fang Yuedong's avatar
Fang Yuedong committed
768
769
770
771
                    cr_map[cr_map > 65535] = 65535
                    cr_map[cr_map < 0] = 0
                    crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
                    del cr_map
Wei Chengliang's avatar
Wei Chengliang committed
772
773
774
775
776
777
778
779
780
781
782
783
784
                    ###########################START
                    ### prescan & overscan
                    if config["ins_effects"]["add_prescan"] == True:
                        chip_utils.log_info(msg="  Apply pre/over-scan", logger=self.logger)
                        crmap_gsimg = chip_utils.AddPreScan(GSImage=crmap_gsimg, pre1=pre1, pre2=pre2, over1=over1, over2=over2)

                    ### 1*16 output
                    if config["ins_effects"]["format_output"] == True:
                        chip_utils.log_info(msg="  Apply 1*16 format", logger=self.logger)
                        crmap_gsimg = chip_utils.formatOutput(GSImage=crmap_gsimg)
                        self.nsecy = 1
                        self.nsecx = 16
                    ###########################END
785
786
                    chip_utils.outputCal(
                        chip=self,
Fang Yuedong's avatar
Fang Yuedong committed
787
788
789
790
791
792
793
                        img=crmap_gsimg,
                        ra_cen=ra_cen,
                        dec_cen=dec_cen,
                        img_rot=img_rot,
                        im_type='CRD',
                        pointing_ID=pointing_ID,
                        output_dir=chip_output.subdir,
794
                        exptime=self.dark_exptime,
795
796
                        project_cycle=config["project_cycle"],
                        run_counter=config["run_counter"],
797
                        timestamp=timestamp_obs)
Fang Yuedong's avatar
Fang Yuedong committed
798
                    del crmap_gsimg
Fang Yuedong's avatar
Fang Yuedong committed
799
800

                # Non-Linearity for Dark
Fang Yuedong's avatar
Fang Yuedong committed
801
                if config["ins_effects"]["non_linear"] == True:
802
803
804
805
                    if self.logger is not None:
                        self.logger.info("  Applying Non-Linearity on the Dark image")
                    else:
                        print("  Applying Non-Linearity on the Dark image", flush=True)
806
                    DarkCombImg = effects.NonLinearity(GSImage=DarkCombImg, beta1=5.e-7, beta2=0)
Fang Yuedong's avatar
Fang Yuedong committed
807

Wei Chengliang's avatar
Wei Chengliang committed
808
809
810
811
812
813
814
                ###if config["ins_effects"]["cte_trail"] == True:
                ###    DarkCombImg = effects.CTE_Effect(GSImage=DarkCombImg, threshold=3)
                ###########################START
                pre1 = self.prescan_x  #27
                over1= self.overscan_x #71
                pre2 = self.prescan_y  #0 #4
                over2= self.overscan_y #84 #80
Fang Yuedong's avatar
Fang Yuedong committed
815
                if config["ins_effects"]["cte_trail"] == True:
Wei Chengliang's avatar
Wei Chengliang committed
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
                    chip_utils.log_info(msg="  Apply CTE Effect", logger=self.logger)
                    ###img = effects.CTE_Effect(GSImage=img, threshold=27)
                    ###CTI_modeling
                    ### 2*8 -> 1*16 img-layout
                    DarkCombImg = chip_utils.formatOutput(GSImage=DarkCombImg)
                    self.nsecy = 1
                    self.nsecx = 16

                    img_arr = DarkCombImg.array
                    ny, nx = img_arr.shape
                    dx = int(nx/self.nsecx)
                    dy = int(ny/self.nsecy)
                    newimg = galsim.Image(nx, int(ny+over2), init_value=0)
                    for ichannel in range(16):
                        print('\n***add CTI effects: pointing-{:} chip-{:} channel-{:}***'.format(pointing_ID, self.chipID, ichannel+1))
                        #nx,ny,noverscan,nsp,nmax = 4608,4616,84,3,10
                        noverscan,nsp,nmax = over2,3,10
                        beta,w,c = 0.478,84700,0
                        t = np.array([0.74,7.7,37],dtype=np.float32)
                        rho_trap = np.array([0.6,1.6,1.4],dtype=np.float32)
                        trap_seeds = np.array([0,1000,10000],dtype=np.int32) + ichannel + self.chipID*16
                        release_seed = 50 + ichannel + pointing_ID*30  + self.chipID*16
                        newimg.array[:, 0+ichannel*dx:dx+ichannel*dx] = CTI_sim(img_arr[:, 0+ichannel*dx:dx+ichannel*dx],dx,dy,noverscan,nsp,nmax,beta,w,c,t,rho_trap,trap_seeds,release_seed)
                    newimg.wcs = DarkCombImg.wcs
                    del DarkCombImg
                    DarkCombImg = newimg

                    ### 1*16 -> 2*8 img-layout
                    DarkCombImg = chip_utils.formatRevert(GSImage=DarkCombImg)
                    self.nsecy = 2
                    self.nsecx = 8

                ### prescan & overscan
                if config["ins_effects"]["add_prescan"] == True:
                    chip_utils.log_info(msg="  Apply pre/over-scan", logger=self.logger)
                    if config["ins_effects"]["cte_trail"] == False:
                        DarkCombImg = chip_utils.AddPreScan(GSImage=DarkCombImg, pre1=pre1, pre2=pre2, over1=over1, over2=over2)
                    if config["ins_effects"]["cte_trail"] == True:
                        DarkCombImg = chip_utils.AddPreScan(GSImage=DarkCombImg, pre1=pre1, pre2=pre2, over1=over1, over2=0)

                ### 1*16 output
                if config["ins_effects"]["format_output"] == True:
                    chip_utils.log_info(msg="  Apply 1*16 format", logger=self.logger)
                    DarkCombImg = chip_utils.formatOutput(GSImage=DarkCombImg)
                    self.nsecy = 1
                    self.nsecx = 16
                ###########################END
Fang Yuedong's avatar
Fang Yuedong committed
863

864
865
866
867
868
                # Add Hot Pixels or/and Dead Pixels
                rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
                badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
                DarkCombImg = effects.DefectivePixels(DarkCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)

Fang Yuedong's avatar
Fang Yuedong committed
869
                # Apply Bad lines 
Fang Yuedong's avatar
Fang Yuedong committed
870
                if config["ins_effects"]["add_badcolumns"] == True:
871
                    DarkCombImg = effects.BadColumns(DarkCombImg, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
872

873
                # Add Bias level
Fang Yuedong's avatar
Fang Yuedong committed
874
                if config["ins_effects"]["add_bias"] == True:
875
876
877
878
                    if self.logger is not None:
                        self.logger.info("  Adding Bias level and 16-channel non-uniformity")
                    else:
                        print("  Adding Bias level and 16-channel non-uniformity")
Fang Yuedong's avatar
Fang Yuedong committed
879
                    # img += float(config["ins_effects"]["bias_level"])
880
881
                    DarkCombImg = effects.AddBiasNonUniform16(DarkCombImg, 
                        bias_level=biaslevel, 
Wei Chengliang's avatar
Wei Chengliang committed
882
                        nsecy = self.nsecy, nsecx=self.nsecx, 
883
884
                        seed=SeedBiasNonuni+self.chipID,
                        logger=self.logger)
885
886

                # Add Read-out Noise
Fang Yuedong's avatar
Fang Yuedong committed
887
                if config["ins_effects"]["add_readout"] == True:
888
                    seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID + 2
889
890
891
                    rng_readout = galsim.BaseDeviate(seed)
                    readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
                    DarkCombImg.addNoise(readout_noise)
Fang Yuedong's avatar
Fang Yuedong committed
892

893
                DarkCombImg, self.gain_channel = effects.ApplyGainNonUniform16(
Fang Yuedong's avatar
Fang Yuedong committed
894
                    DarkCombImg, gain=self.gain, 
Wei Chengliang's avatar
Wei Chengliang committed
895
                    nsecy = self.nsecy, nsecx=self.nsecx, 
896
897
                    seed=SeedGainNonuni+self.chipID,
                    logger=self.logger)
Fang Yuedong's avatar
Fang Yuedong committed
898
899
900
901
902
903
904
                # DarkCombImg = effects.AddOverscan(
                #     DarkCombImg, 
                #     overscan=overscan, gain=self.gain, 
                #     widthl=27, widthr=27, widtht=8, widthb=8)
                DarkCombImg.replaceNegative(replace_value=0)
                DarkCombImg.quantize()
                DarkCombImg = galsim.ImageUS(DarkCombImg)
905
                timestamp_obs += 10 * 60
906
                chip_utils.outputCal(
907
                    chip=self,
Fang Yuedong's avatar
Fang Yuedong committed
908
909
910
911
                    img=DarkCombImg,
                    ra_cen=ra_cen,
                    dec_cen=dec_cen,
                    img_rot=img_rot,
Xin Zhang's avatar
Xin Zhang committed
912
                    im_type='DARK',
Fang Yuedong's avatar
Fang Yuedong committed
913
914
                    pointing_ID=pointing_ID,
                    output_dir=chip_output.subdir,
915
                    exptime=self.dark_exptime,
916
917
                    project_cycle=config["project_cycle"],
                    run_counter=config["run_counter"],
918
                    timestamp = timestamp_obs)
Fang Yuedong's avatar
Fang Yuedong committed
919
920
921
922
            del DarkCombImg
        # img = galsim.ImageUS(img)

        # # 16 output channel, with each a single image file
Fang Yuedong's avatar
Fang Yuedong committed
923
        # if config["ins_effects"]["readout16"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
924
925
926
927
928
929
930
931
932
933
934
935
        #     print("  16 Output Channel simulation")
        #     for coli in [0, 1]:
        #         for rowi in range(8):
        #             sub_img = effects.readout16(
        #                 GSImage=img, 
        #                 rowi=rowi, 
        #                 coli=coli, 
        #                 overscan_value=self.overscan)
        #             rowcoltag = str(rowi) + str(coli)
        #             img_name_root = chip_output.img_name.split(".")[0]
        #             sub_img.write("%s/%s_%s.fits" % (chip_output.subdir, img_name_root, rowcoltag))
        #     del sub_img
Zhang Xin's avatar
Zhang Xin committed
936
        return img
937