ObservationSim.py 22.2 KB
Newer Older
1
import os
Fang Yuedong's avatar
Fang Yuedong committed
2
3
4
5
6
import numpy as np
import mpi4py.MPI as MPI
import galsim
import logging
import psutil
Fang Yuedong's avatar
Fang Yuedong committed
7
import gc
8
9
10
from astropy.io import fits
from datetime import datetime

Fang Yuedong's avatar
Fang Yuedong committed
11
12
import traceback

13
14
15
from ObservationSim.Config import config_dir, ChipOutput
from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip
16
from ObservationSim.Instrument.Chip import Effects
17
from ObservationSim.Straylight import calculateSkyMap_split_g
18
19
from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp
from ObservationSim._util import get_shear_field, makeSubDir_PointingList
20
from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
Fang Yuedong's avatar
Fang Yuedong committed
21
22

class Observation(object):
23
    def __init__(self, config, Catalog, work_dir=None, data_dir=None):
24
        self.path_dict = config_dir(config=config, work_dir=work_dir, data_dir=data_dir)
Fang Yuedong's avatar
Fang Yuedong committed
25
        self.config = config
Fang Yuedong's avatar
Fang Yuedong committed
26
        self.tel = Telescope()
27
        self.focal_plane = FocalPlane(survey_type=self.config["obs_setting"]["survey_type"]) 
Fang Yuedong's avatar
Fang Yuedong committed
28
        self.filter_param = FilterParam() 
Fang Yuedong's avatar
Fang Yuedong committed
29
30
        self.chip_list = []
        self.filter_list = []
31
        self.all_filter = []
32
        self.Catalog = Catalog
Fang Yuedong's avatar
Fang Yuedong committed
33
34

        # Construct chips & filters:
Fang Yuedong's avatar
Fang Yuedong committed
35
        for i in range(self.focal_plane.nchips):
Fang Yuedong's avatar
Fang Yuedong committed
36
37
38
            chipID = i + 1

            # Make Chip & Filter lists
39
40
            chip = Chip(
                chipID=chipID, 
Fang Yuedong's avatar
Fang Yuedong committed
41
                config=self.config)
Fang Yuedong's avatar
Fang Yuedong committed
42
            filter_id, filter_type = chip.getChipFilter()
43
44
            filt = Filter(filter_id=filter_id, 
                filter_type=filter_type, 
45
                filter_param=self.filter_param)
46
47
48
49
            if not self.focal_plane.isIgnored(chipID=chipID):
                self.chip_list.append(chip)
                self.filter_list.append(filt)
            self.all_filter.append(filt)
Fang Yuedong's avatar
Fang Yuedong committed
50

51
    def run_one_chip(self, chip, filt, pointing, chip_output, wcs_fp=None, psf_model=None, shear_cat_file=None, cat_dir=None, sed_dir=None):
Fang Yuedong's avatar
Fang Yuedong committed
52

Fang Yuedong's avatar
Fang Yuedong committed
53
54
55
56
57
58
59
60
61
        chip_output.Log_info(':::::::::::::::::::Current Pointing Information::::::::::::::::::')
        chip_output.Log_info("RA: %f, DEC; %f" % (pointing.ra, pointing.dec))
        chip_output.Log_info("Time: %s" % datetime.utcfromtimestamp(pointing.timestamp).isoformat())
        chip_output.Log_info("Exposure time: %f" % pointing.exp_time)
        chip_output.Log_info("Satellite Position (x, y, z): (%f, %f, %f)" % (pointing.sat_x, pointing.sat_y, pointing.sat_z))
        chip_output.Log_info("Satellite Velocity (x, y, z): (%f, %f, %f)" % (pointing.sat_vx, pointing.sat_vy, pointing.sat_vz))
        chip_output.Log_info("Position Angle: %f" % pointing.img_pa.deg)
        chip_output.Log_info('Chip : %d' % chip.chipID)
        chip_output.Log_info(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::')
Fang Yuedong's avatar
Fang Yuedong committed
62

Fang Yuedong's avatar
Fang Yuedong committed
63
        if self.config["psf_setting"]["psf_model"] == "Gauss":
Fang Yuedong's avatar
Fang Yuedong committed
64
            psf_model = PSFGauss(chip=chip, psfRa=self.config["psf_setting"]["psf_rcont"])
Fang Yuedong's avatar
Fang Yuedong committed
65
        elif self.config["psf_setting"]["psf_model"] == "Interp":
Fang Yuedong's avatar
Fang Yuedong committed
66
            psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.path_dict["psf_dir"])
Fang Yuedong's avatar
Fang Yuedong committed
67
        else:
Fang Yuedong's avatar
Fang Yuedong committed
68
            chip_output.Log_error("unrecognized PSF model type!!", flush=True)
Fang Yuedong's avatar
Fang Yuedong committed
69

Fang Yuedong's avatar
Fang Yuedong committed
70
        # Figure out shear fields
Fang Yuedong's avatar
Fang Yuedong committed
71
        if shear_cat_file is not None:
72
            self.g1_field, self.g2_field, self.nshear = get_shear_field(config=self.config, shear_cat_file=shear_cat_file)
Fang Yuedong's avatar
Fang Yuedong committed
73

74
75
        # Apply astrometric simulation for pointing
        if self.config["obs_setting"]["enable_astrometric_model"]:
Fang Yuedong's avatar
Fang Yuedong committed
76
            dt = datetime.utcfromtimestamp(pointing.timestamp)
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
            date_str = dt.date().isoformat()
            time_str = dt.time().isoformat()
            ra_cen, dec_cen = on_orbit_obs_position(
                input_ra_list=[pointing.ra],
                input_dec_list=[pointing.dec],
                input_pmra_list=[0.],
                input_pmdec_list=[0.],
                input_rv_list=[0.],
                input_parallax_list=[1e-9],
                input_nstars=1,
                input_x=pointing.sat_x,
                input_y=pointing.sat_y,
                input_z=pointing.sat_z,
                input_vx=pointing.sat_vx,
                input_vy=pointing.sat_vy,
                input_vz=pointing.sat_vz,
Fang Yuedong's avatar
Fang Yuedong committed
93
                input_epoch="J2000",
94
95
96
97
98
99
100
101
                input_date_str=date_str,
                input_time_str=time_str
            )
            ra_cen, dec_cen = ra_cen[0], dec_cen[0]
        else:
            ra_cen = pointing.ra
            dec_cen = pointing.dec

Fang Yuedong's avatar
Fang Yuedong committed
102
103
        # Get WCS for the focal plane
        if wcs_fp == None:
104
            wcs_fp = self.focal_plane.getTanWCS(ra_cen, dec_cen, pointing.img_pa, chip.pix_scale)
Fang Yuedong's avatar
Fang Yuedong committed
105
106
107
108
109

        # Create chip Image
        chip.img = galsim.ImageF(chip.npix_x, chip.npix_y)
        chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
        chip.img.wcs = wcs_fp
110

Zhang Xin's avatar
Zhang Xin committed
111
112
        if self.config["obs_setting"]["enable_straylight_model"]:
            filt.setFilterStrayLightPixel(jtime = pointing.jdt, sat_pos = np.array([pointing.sat_x, pointing.sat_y, pointing.sat_z]), pointing_radec = np.array([pointing.ra,pointing.dec]), sun_pos = np.array([pointing.sun_x,pointing.sun_y,pointing.sun_z]))
113

Zhang Xin's avatar
Zhang Xin committed
114
        print("========================sky pix========================")
115
116
        print(filt.sky_background)

Fang Yuedong's avatar
Fang Yuedong committed
117
118
119
        if chip.survey_type == "photometric":
            sky_map = None
        elif chip.survey_type == "spectroscopic":
xin's avatar
xin committed
120
            # chip.loadSLSFLATCUBE(flat_fn='flat_cube.fits')
121
122
            flat_normal = np.ones_like(chip.img.array)
            if self.config["ins_effects"]["flat_fielding"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
123
                chip_output.Log_info("SLS flat preprocess,CHIP %d : Creating and applying Flat-Fielding"%chip.chipID)
124
                msg = str(chip.img.bounds)
Fang Yuedong's avatar
Fang Yuedong committed
125
                chip_output.Log_info(msg)
126
127
128
129
130
                flat_img = Effects.MakeFlatSmooth(
                    chip.img.bounds,
                    int(self.config["random_seeds"]["seed_flat"]))
                flat_normal = flat_normal * flat_img.array / np.mean(flat_img.array)
            if self.config["ins_effects"]["shutter_effect"] == True:
Fang Yuedong's avatar
Fang Yuedong committed
131
                chip_output.Log_info("SLS flat preprocess,CHIP %d : Apply shutter effect"%chip.chipID)
132
133
134
135
                shuttimg = Effects.ShutterEffectArr(chip.img, t_shutter=1.3, dist_bearing=735,
                                                    dt=1E-3)  # shutter effect normalized image for this chip
                flat_normal = flat_normal*shuttimg
                flat_normal = np.array(flat_normal,dtype='float32')
Fang Yuedong's avatar
Fang Yuedong committed
136
            sky_map = calculateSkyMap_split_g(
137
138
139
140
141
142
                skyMap=flat_normal,
                blueLimit=filt.blue_limit,
                redLimit=filt.red_limit,
                conf=chip.sls_conf,
                pixelSize=chip.pix_scale,
                isAlongY=0,
143
144
                flat_cube=chip.flat_cube, zoldial_spec = filt.zodical_spec)
            sky_map = sky_map+filt.sky_background
145
            del flat_normal
Fang Yuedong's avatar
Fang Yuedong committed
146

147
        if pointing.pointing_type == 'MS':
Fang Yuedong's avatar
Fang Yuedong committed
148
            # Load catalogues and templates
Fang Yuedong's avatar
Fang Yuedong committed
149
            self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, chip_output=chip_output, filt=filt)
150
            chip_output.create_output_file()
Fang Yuedong's avatar
Fang Yuedong committed
151
152
            self.nobj = len(self.cat.objs)

Fang Yuedong's avatar
Fang Yuedong committed
153
154
155
156
157
158
159
160
161
            for ifilt in range(len(self.all_filter)):
                temp_filter = self.all_filter[ifilt]
                # Update the limiting magnitude using exposure time in pointing
                temp_filter.update_limit_saturation_mags(exptime=pointing.exp_time, chip=chip)

                # Select cutting band filter for saturation/limiting magnitude
                if temp_filter.filter_type.lower() == self.config["obs_setting"]["cut_in_band"].lower():
                    cut_filter = temp_filter

Fang Yuedong's avatar
Fang Yuedong committed
162
163
164
165
166
            if self.config["ins_effects"]["field_dist"] == True:
                self.fd_model = FieldDistortion(chip=chip)
            else:
                self.fd_model = None

Fang Yuedong's avatar
Fang Yuedong committed
167
168
169
170
            # Loop over objects
            missed_obj = 0
            bright_obj = 0
            dim_obj = 0
Fang Yuedong's avatar
Fang Yuedong committed
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
 
            h_ext = generateExtensionHeader(
                chip=chip,
                xlen=chip.npix_x, 
                ylen=chip.npix_y, 
                ra=pointing.ra, 
                dec=pointing.dec, 
                pa=pointing.img_pa.deg, 
                gain=chip.gain, 
                readout=chip.read_noise, 
                dark=chip.dark_noise, 
                saturation=90000, 
                pixel_scale=chip.pix_scale, 
                pixel_size=chip.pix_size,
                xcen=chip.x_cen,
                ycen=chip.y_cen,
187
188
189
190
                extName='SCI',
                timestamp = pointing.timestamp,
                exptime = pointing.exp_time,
                readoutTime = 40.)
Fang Yuedong's avatar
Fang Yuedong committed
191

Fang Yuedong's avatar
Fang Yuedong committed
192
            for j in range(self.nobj):
193
194
                
                # (DEBUG)
195
                # if j >= 10:
Xin Zhang's avatar
Xin Zhang committed
196
                #     break
197

Fang Yuedong's avatar
Fang Yuedong committed
198
199
                obj = self.cat.objs[j]

200
                # load and convert SED; also caculate object's magnitude in all CSST bands
201
202
203
                try:
                    sed_data = self.cat.load_sed(obj)
                    norm_filt = self.cat.load_norm_filt(obj)
Fang Yuedong's avatar
Fang Yuedong committed
204
                    obj.sed, obj.param["mag_%s"%filt.filter_type.lower()], obj.param["flux_%s"%filt.filter_type.lower()] = self.cat.convert_sed(
205
206
207
208
209
                        mag=obj.param["mag_use_normal"],
                        sed=sed_data,
                        target_filt=filt, 
                        norm_filt=norm_filt,
                    )
Fang Yuedong's avatar
Fang Yuedong committed
210
                    _, obj.param["mag_%s"%cut_filter.filter_type.lower()], obj.param["flux_%s"%cut_filter.filter_type.lower()] = self.cat.convert_sed(
Fang Yuedong's avatar
Fang Yuedong committed
211
212
213
214
215
                        mag=obj.param["mag_use_normal"],
                        sed=sed_data,
                        target_filt=cut_filter, 
                        norm_filt=norm_filt,
                    )
216
                except Exception as e:
Fang Yuedong's avatar
Fang Yuedong committed
217
                    traceback.print_exc()
Fang Yuedong's avatar
Fang Yuedong committed
218
                    chip_output.Log_error(e)
219
                    continue
220
                
Fang Yuedong's avatar
Fang Yuedong committed
221
222
                # [TODO] Testing
                # chip_output.Log_info("mag_%s = %.3f"%(filt.filter_type.lower(), obj.param["mag_%s"%filt.filter_type.lower()]))
Fang Yuedong's avatar
Fang Yuedong committed
223
224

                # Exclude very bright/dim objects (for now)
225
226
227
                if cut_filter.is_too_bright(
                    mag=obj.param["mag_%s"%self.config["obs_setting"]["cut_in_band"].lower()],
                    margin=self.config["obs_setting"]["mag_sat_margin"]):
Fang Yuedong's avatar
Fang Yuedong committed
228
229
230
231
                    chip_output.Log_info("obj %s too birght!! mag_%s = %.3f"%(obj.id, cut_filter.filter_type, obj.param["mag_%s"%self.config["obs_setting"]["cut_in_band"].lower()]))
                    bright_obj += 1
                    obj.unload_SED()
                    continue
232
233
234
                if filt.is_too_dim(
                    mag=obj.getMagFilter(filt),
                    margin=self.config["obs_setting"]["mag_lim_margin"]):
Fang Yuedong's avatar
Fang Yuedong committed
235
                    chip_output.Log_info("obj %s too dim!! mag_%s = %.3f"%(obj.id, filt.filter_type, obj.getMagFilter(filt)))
Fang Yuedong's avatar
Fang Yuedong committed
236
                    dim_obj += 1
Fang Yuedong's avatar
Fang Yuedong committed
237
238
239
                    obj.unload_SED()
                    continue

Fang Yuedong's avatar
Fang Yuedong committed
240
                # Get corresponding shear values
Fang Yuedong's avatar
Fang Yuedong committed
241
                if self.config["shear_setting"]["shear_type"] == "constant":
Fang Yuedong's avatar
Fang Yuedong committed
242
                    if obj.type == 'star':
243
                        obj.g1, obj.g2 = 0., 0.
Fang Yuedong's avatar
Fang Yuedong committed
244
                    else:
245
                        obj.g1, obj.g2 = self.g1_field, self.g2_field
Fang Yuedong's avatar
Fang Yuedong committed
246
247
                elif self.config["shear_setting"]["shear_type"] == "extra":
                    try:
Fang Yuedong's avatar
Fang Yuedong committed
248
                        # [TODO]: every object with individual shear from input catalog(s)
249
                        obj.g1, obj.g2 = self.g1_field[j], self.g2_field[j]
Fang Yuedong's avatar
Fang Yuedong committed
250
                    except:
Fang Yuedong's avatar
Fang Yuedong committed
251
                        chip_output.Log_error("failed to load external shear.")
252
253
254
                elif self.config["shear_setting"]["shear_type"] == "catalog":
                    pass
                else:
Fang Yuedong's avatar
Fang Yuedong committed
255
                    chip_output.Log_error("Unknown shear input")
256
                    raise ValueError("Unknown shear input")
Fang Yuedong's avatar
Fang Yuedong committed
257
258
259
260
261
262
263

                # Get position of object on the focal plane
                pos_img, offset, local_wcs, real_wcs, fd_shear = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False, img_header=h_ext)

                # [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane
                # Otherwise they will be considered missed objects
                # if pos_img.x == -1 or pos_img.y == -1 or (not chip.isContainObj(x_image=pos_img.x, y_image=pos_img.y, margin=0.)):
Fang Yuedong's avatar
Fang Yuedong committed
264
                if pos_img.x == -1 or pos_img.y == -1:
Fang Yuedong's avatar
Fang Yuedong committed
265
266
                    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))
                    chip_output.Log_error("Objected missed: %s"%(obj.id))
Fang Yuedong's avatar
Fang Yuedong committed
267
268
269
                    missed_obj += 1
                    obj.unload_SED()
                    continue
Fang Yuedong's avatar
Fang Yuedong committed
270

Fang Yuedong's avatar
Fang Yuedong committed
271
272
                # Draw object & update output catalog
                try:
Fang Yuedong's avatar
Fang Yuedong committed
273
                    if self.config["run_option"]["out_cat_only"]:
Fang Yuedong's avatar
Fang Yuedong committed
274
                        isUpdated = True
Fang Yuedong's avatar
Fang Yuedong committed
275
                        obj.real_pos = obj.getRealPos(chip.img, global_x=obj.posImg.x, global_y=obj.posImg.y, img_real_wcs=obj.real_wcs)
276
                        pos_shear = 0.
277
                    elif chip.survey_type == "photometric" and not self.config["run_option"]["out_cat_only"]:
Fang Yuedong's avatar
Fang Yuedong committed
278
279
280
281
282
283
284
                        isUpdated, pos_shear = obj.drawObj_multiband(
                            tel=self.tel,
                            pos_img=pos_img, 
                            psf_model=psf_model, 
                            bandpass_list=filt.bandpass_sub_list, 
                            filt=filt, 
                            chip=chip, 
285
286
                            g1=obj.g1, 
                            g2=obj.g2, 
Fang Yuedong's avatar
Fang Yuedong committed
287
288
289
                            exptime=pointing.exp_time,
                            fd_shear=fd_shear)

Fang Yuedong's avatar
Fang Yuedong committed
290
                    elif chip.survey_type == "spectroscopic" and not self.config["run_option"]["out_cat_only"]:
Fang Yuedong's avatar
Fang Yuedong committed
291
292
293
294
295
296
297
                        isUpdated, pos_shear = obj.drawObj_slitless(
                            tel=self.tel, 
                            pos_img=pos_img, 
                            psf_model=psf_model, 
                            bandpass_list=filt.bandpass_sub_list, 
                            filt=filt, 
                            chip=chip, 
298
299
                            g1=obj.g1, 
                            g2=obj.g2, 
300
                            exptime=pointing.exp_time,
Fang Yuedong's avatar
Fang Yuedong committed
301
302
303
304
                            normFilter=norm_filt,
                            fd_shear=fd_shear)

                    if isUpdated == 1:
Fang Yuedong's avatar
Fang Yuedong committed
305
                        # TODO: add up stats
306
                        chip_output.cat_add_obj(obj, pos_img, pos_shear)
Fang Yuedong's avatar
Fang Yuedong committed
307
                        pass
Fang Yuedong's avatar
Fang Yuedong committed
308
309
310
                    elif isUpdated == 0:
                        missed_obj += 1
                        chip_output.Log_error("Objected missed: %s"%(obj.id))
Fang Yuedong's avatar
Fang Yuedong committed
311
                    else:
Fang Yuedong's avatar
Fang Yuedong committed
312
                        chip_output.Log_error("Draw error, object omitted: %s"%(obj.id))
Fang Yuedong's avatar
Fang Yuedong committed
313
314
                        continue
                except Exception as e:
Fang Yuedong's avatar
Fang Yuedong committed
315
                    traceback.print_exc()
Fang Yuedong's avatar
Fang Yuedong committed
316
317
318
319
320
321
322
                    chip_output.Log_error(e)
                    
                # # [C6 TEST]
                # chip_output.Log_info("check running:1: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ))
                # chip_output.Log_info('draw object %s'%obj.id)
                # chip_output.Log_info('mag = %.3f'%obj.param['mag_use_normal'])

Fang Yuedong's avatar
Fang Yuedong committed
323
324
325
                # Unload SED:
                obj.unload_SED()
                del obj
Fang Yuedong's avatar
Fang Yuedong committed
326
                gc.collect()
Fang Yuedong's avatar
Fang Yuedong committed
327

Fang Yuedong's avatar
Fang Yuedong committed
328
329
            del psf_model
            del self.cat
Fang Yuedong's avatar
Fang Yuedong committed
330
            gc.collect()
Fang Yuedong's avatar
Fang Yuedong committed
331

Fang Yuedong's avatar
Fang Yuedong committed
332
        chip_output.Log_info("check running:1: pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ))
Fang Yuedong's avatar
Fang Yuedong committed
333
334
335

        # Detector Effects
        # ===========================================================
Fang Yuedong's avatar
Fang Yuedong committed
336
        # whether to output zero, dark, flat calibration images.
337
338
339
340
341
342

        if not self.config["run_option"]["out_cat_only"]:
            chip.img = chip.addEffects(
                config=self.config, 
                img=chip.img, 
                chip_output=chip_output, 
343
                filt=filt,
344
345
346
                ra_cen=pointing.ra, 
                dec_cen=pointing.dec,
                img_rot=pointing.img_pa,
Fang Yuedong's avatar
Fang Yuedong committed
347
                exptime=pointing.exp_time,
348
349
350
351
352
353
354
                pointing_ID=pointing.id,
                timestamp_obs=pointing.timestamp,
                pointing_type=pointing.pointing_type,
                sky_map=sky_map, tel = self.tel,
                logger=chip_output.logger)
            
            if pointing.pointing_type == 'MS':
Fang Yuedong's avatar
Fang Yuedong committed
355
                datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
356
357
358
359
360
361
362
363
                date_obs = datetime_obs.strftime("%y%m%d")
                time_obs = datetime_obs.strftime("%H%M%S")
                h_prim = generatePrimaryHeader(
                    xlen=chip.npix_x, 
                    ylen=chip.npix_y, 
                    pointNum = str(pointing.id),
                    ra=pointing.ra, 
                    dec=pointing.dec, 
Fang Yuedong's avatar
Fang Yuedong committed
364
                    pixel_scale=chip.pix_scale,
365
366
367
368
369
                    date=date_obs,
                    time_obs=time_obs,
                    exptime=pointing.exp_time,
                    im_type='SCI',
                    sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z],
Fang Yuedong's avatar
Fang Yuedong committed
370
371
                    sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz],
                    chip_name=str(chip.chipID).rjust(2, '0'))
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
                h_ext = generateExtensionHeader(
                    chip=chip,
                    xlen=chip.npix_x,
                    ylen=chip.npix_y,
                    ra=pointing.ra,
                    dec=pointing.dec,
                    pa=pointing.img_pa.deg,
                    gain=chip.gain,
                    readout=chip.read_noise,
                    dark=chip.dark_noise,
                    saturation=90000,
                    pixel_scale=chip.pix_scale,
                    pixel_size=chip.pix_size,
                    xcen=chip.x_cen,
                    ycen=chip.y_cen,
                    extName='SCI',
                    timestamp=pointing.timestamp,
                    exptime=pointing.exp_time,
                    readoutTime=40.)
Fang Yuedong's avatar
Fang Yuedong committed
391
                
392
393
                chip.img = galsim.Image(chip.img.array, dtype=np.uint16)
                hdu1 = fits.PrimaryHDU(header=h_prim)
Fang Yuedong's avatar
Fang Yuedong committed
394
                hdu1.add_checksum()
395
396
                hdu1.header.comments['CHECKSUM'] = 'HDU checksum'
                hdu1.header.comments['DATASUM'] = 'data unit checksum'
397
                hdu2 = fits.ImageHDU(chip.img.array, header=h_ext)
Fang Yuedong's avatar
Fang Yuedong committed
398
                hdu2.add_checksum()
399
400
401
                hdu2.header.comments['XTENSION'] = 'extension type'
                hdu2.header.comments['CHECKSUM'] = 'HDU checksum'
                hdu2.header.comments['DATASUM'] = 'data unit checksum'
402
403
404
                hdu1 = fits.HDUList([hdu1, hdu2])
                fname = os.path.join(chip_output.subdir, h_prim['FILENAME'] + '.fits')
                hdu1.writeto(fname, output_verify='ignore', overwrite=True)
Fang Yuedong's avatar
Fang Yuedong committed
405
406
407
408

                chip_output.Log_info("# objects that are too bright %d out of %d"%(bright_obj, self.nobj))
                chip_output.Log_info("# objects that are too dim %d out of %d"%(dim_obj, self.nobj))
                chip_output.Log_info("# objects that are missed %d out of %d"%(missed_obj, self.nobj))
Fang Yuedong's avatar
Fang Yuedong committed
409
410
        del chip.img

Fang Yuedong's avatar
Fang Yuedong committed
411
        chip_output.Log_info("check running:2: pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ))
Fang Yuedong's avatar
Fang Yuedong committed
412

413
    def runExposure_MPI_PointingList(self, pointing_list, shear_cat_file=None, chips=None, use_mpi=False):
Fang Yuedong's avatar
Fang Yuedong committed
414
415
416
417
        if use_mpi:
            comm = MPI.COMM_WORLD
            ind_thread = comm.Get_rank()
            num_thread = comm.Get_size()
Fang Yuedong's avatar
Fang Yuedong committed
418

Fang Yuedong's avatar
Fang Yuedong committed
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
        if chips is None:
            nchips_per_fp = len(self.chip_list)
            run_chips = self.chip_list
            run_filts = self.filter_list
        else:
            # Only run a particular set of chips
            run_chips = []
            run_filts = []
            nchips_per_fp = len(chips)
            for ichip in range(len(self.chip_list)):
                chip = self.chip_list[ichip]
                filt = self.filter_list[ichip]
                if chip.chipID in chips:
                    run_chips.append(chip)
                    run_filts.append(filt)
Fang Yuedong's avatar
Fang Yuedong committed
434

435
        for ipoint in range(len(pointing_list)):
Fang Yuedong's avatar
Fang Yuedong committed
436
437
            for ichip in range(nchips_per_fp):
                i = ipoint*nchips_per_fp + ichip
438
439
                pointing = pointing_list[ipoint]
                pointing_ID = pointing.id
Fang Yuedong's avatar
Fang Yuedong committed
440
441
442
                if use_mpi:
                    if i % num_thread != ind_thread:
                        continue
Fang Yuedong's avatar
Fang Yuedong committed
443
444
445
446
447

                pid = os.getpid()

                sub_img_dir, prefix = makeSubDir_PointingList(path_dict=self.path_dict, config=self.config, pointing_ID=pointing_ID)

Fang Yuedong's avatar
Fang Yuedong committed
448
449
                chip = run_chips[ichip]
                filt = run_filts[ichip]
Fang Yuedong's avatar
Fang Yuedong committed
450
                # chip_output.Log_info("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid))
Fang Yuedong's avatar
Fang Yuedong committed
451
452
453
454
455
                chip_output = ChipOutput(
                    config=self.config, 
                    focal_plane=self.focal_plane, 
                    chip=chip, 
                    filt=filt,  
456
457
                    exptime=pointing.exp_time,
                    pointing_type=pointing.pointing_type,
Fang Yuedong's avatar
Fang Yuedong committed
458
459
460
                    pointing_ID=pointing_ID,  
                    subdir=sub_img_dir,
                    prefix=prefix)
Fang Yuedong's avatar
Fang Yuedong committed
461
                chip_output.Log_info("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid))
462
                self.run_one_chip(
Fang Yuedong's avatar
Fang Yuedong committed
463
464
465
                    chip=chip, 
                    filt=filt, 
                    chip_output=chip_output, 
Fang Yuedong's avatar
Fang Yuedong committed
466
467
                    pointing=pointing)
                chip_output.Log_info("finished running chip#%d..."%(chip.chipID))
468
469
                for handler in chip_output.logger.handlers[:]:
                    chip_output.logger.removeHandler(handler)
Fang Yuedong's avatar
Fang Yuedong committed
470
                gc.collect()