init

BKGscript.py

+from concurrent import futures
+from glob import glob
+from astropy.io import fits
+import os
+import numpy as np
+from time import time
+from concurrent.futures import ThreadPoolExecutor
+
+A, B = 2, 2
+
+
+def split(data: np.ndarray, i: int) -> np.ndarray:
+    h, w = data.shape
+    x, y = h//A, w//B
+    n, m = i//B, i % B
+    data = data[n*x:(n+1)*x,
+                m*y:(m+1)*y].copy()
+    return data
+
+
+def concatenate(data_list):
+    data_u = np.concatenate(data_list[:B], axis=1)
+    data_d = np.concatenate(data_list[B:], axis=1)
+    data = np.concatenate([data_u, data_d], axis=0)
+    return data
+
+
+def array_combine(ndarray, mode="mean") -> np.ndarray:
+    """ Function to combine 3-D data array
+
+    Parameters
+    ----------
+    ndarray: array, input data cube (3D)
+    model: mean, median, sum, mean_clip, median_clip, default is mean
+    """
+    if mode == "median":
+        array = np.median(ndarray, axis=0)
+    elif mode == "median_clip":
+        ndarray = np.sort(ndarray, axis=0)[1:-1]
+        array = np.median(ndarray, axis=0)
+    elif mode == "sum":
+        array = np.sum(ndarray, axis=0)
+    elif mode == "mean":
+        array = np.mean(ndarray, axis=0)
+    elif mode == "mean_clip":
+        ndarray = np.sort(ndarray, axis=0)[1:-1]
+        array = np.mean(ndarray, axis=0)
+    return array
+
+
+def load_bias(path: str, i: int) -> np.ndarray:
+    with fits.open(path) as hdul:
+        du = hdul[1].data
+    du = du.astype(int)
+    du = split(du, i)
+    return du
+
+
+def load_dark(path: str, i: int, bias: np.ndarray) -> np.ndarray:
+    with fits.open(path) as hdul:
+        du = hdul[1].data
+        hu = hdul[0].header
+    du = du.astype(int)
+    du = du - bias
+    du = du / hu["EXPTIME"]
+    du = split(du, i)
+    return du
+
+
+def load_flat(path: str, i: int, bias: np.ndarray, dark: np.ndarray) -> np.ndarray:
+    with fits.open(path) as hdul:
+        du = hdul[1].data
+        hu = hdul[0].header
+    du = du.astype(int)
+    du = du - bias - dark * hu["EXPTIME"]
+    du = du / hu["EXPTIME"]
+    du = du / np.median(du)
+    du = split(du, i)
+    return du
+
+
+def save_fits(func, mode: str, path_list, save_path, *args) -> np.ndarray:
+    ch_list = []
+    with ThreadPoolExecutor() as tpe:
+        for i in range(A*B):
+            futures = [tpe.submit(func, path, i, *args) for path in path_list]
+            du_list = [future.result() for future in futures]
+            du = array_combine(du_list, mode)
+            ch_list.append(du)
+    du = concatenate(ch_list)
+    du_output = os.path.basename(path_list[0])
+    du_output = du_output.replace("raw", "combine")
+    du_output = os.path.join(save_path, du_output)
+    du = du.astype(np.float32)
+    du_fits = fits.HDUList([fits.PrimaryHDU(data=du)])
+    du_fits.writeto(du_output, overwrite=True)
+    return du
+
+
+def main(input_path : str, save_path : str, number_list : list):
+    for number in number_list:
+        print(number, end=' ')
+        bias_path = glob(os.path.join(input_path, "MSC*/*CLB*_" + number + '_*'))
+        bias = save_fits(load_bias, "median", bias_path, save_path)
+        print('bias finish', end=' ')
+        dark_path = glob(os.path.join(input_path, "MSC*/*CLD*_" + number + '_*'))
+        dark = save_fits(load_dark, "median", dark_path, save_path, bias)
+        print('dark finish', end=' ')
+        flat_path = glob(os.path.join(input_path, "MSC*/*CLF*_" + number + '_*'))
+        flat = save_fits(load_flat, "median", flat_path, save_path, bias, dark)
+        print('flat finish')
+
+
+if __name__ == "__main__":
+    input_path = "/data/test20211012/ref/"
+    save_path = "/data/test20211012/ref/"
+    number_list = ['06', '07', '08', '09', '11', '12', '13', '14',
+                   '15', '16', '17', '18', '19', '20', '22', '23', '24', '25']
+    main(input_path, save_path, number_list)
+

Dockerfile

+FROM hub.cstcloud.cn/csst/python:3.8
+
+RUN pip install --no-deps ccdproc && pip install deepCR
+
+RUN apt-get update && apt-get install -y vim
+
+RUN mkdir -p /home/csstpipeline/code/csst-msc-instrument-effect/csst_msc_instrument_effect
+# COPY run.sh /app/bin/
+COPY instrument_effect_wrapper.py setup.py requirements.txt \ 
+    /home/csstpipeline/code/csst-msc-instrument-effect/
+
+COPY csst_msc_instrument_effect/* \ 
+    /home/csstpipeline/code/csst-msc-instrument-effect/csst_msc_instrument_effect/
+
+COPY MSC_crmask.ini /home/csstpipeline/csst-msc-instrument-effect-master/
+
+RUN cd /home/csstpipeline/code/csst-msc-instrument-effect/ \
+    && python setup.py install
+
+COPY cleanup.sh /app/bin/
+COPY csst_msc_iec.sh /app/bin/run.sh

MSC_crmask.ini

+[global]
+;which method used for cosmic ray detection
+;model=deepCR
+
+update_flag=True
+;whether save CR cleaned image
+save_flag=True
+;flag image suffix
+flag_suffix=flg
+;CR masked image suffix
+data_suffix=crclean
+;whether append CR masked image to the output FITS file
+append_flag=False
+;mask image suffix
+mask_suffix=crmask
+;whether use GPU
+gpu_flag=False
+;whether fill appropriate value for CR contaminated pixels
+fill_flag=True
+;filling method
+fill_method=meanmask
+;torch thread
+torch_thread=1
+
+[deepCR]
+threshold=0.5
+inpaint=True
+binary=True
+patch=256
+segment=True
+parallel=True
+;Cosmic ray clean training model
+clean_model=/home/csstpipeline/csst-msc-instrument-effect-main/CSST_2021-12-30_CCD23_epoch20.pth
+;Inpaint model
+inpaint_model=ACS-WFC-F606W-2-32
+n_jobs=4
+;train parameters
+;ignore pixel, e.g., bad pixel, saturation, npy or fits
+;ignore=ignore.fits
+;sky=sky.fits
+;aug_sky_min=0
+;aug_sky_max=0
+;name=xxxx
+hidden=50
+gpu=False
+epoch=50
+batch=16
+lr=0.005
+auto_lr_decay=True
+lr_decay_patience=4
+lr_decay_factor=0.1
+save_after=1
+plot_every=10
+use_tqdm=False
+use_tqdm_notebook=False
+plot_roc=False
+directory=./data/
+
+[lacosmic]
+sigclip=4.5
+sigfrac=0.3
+objlim=5.0
+gain=1.0
+readnoise=6.5
+satlevel=65535.0
+pssl=0.0
+niter=4
+sepmed=True
+cleantype=meanmask
+fsmode=median
+psfmodel=gauss
+psffwhm=2.5
+psfsize=7
+psfbeta=4.765
+gain_apply=True

Makefile

+IMAGE=hub.cstcloud.cn/csst/detector-effect-correction
+
+build:
+	docker build --network=host -t $(IMAGE) .
+
+push:
+	docker push $(IMAGE)
+
+clean:
+	docker rmi $(IMAGE)
+
+rmi:
+	for H in h0 h1 h2 ; do echo $$H ; ssh $$H docker rmi $(IMAGE) ; done

README.md

+# csst_msc_instrument_effect 仪器效应改正
+
+用于改正CSST的0级数据中的仪器效应问题，如本底、暗流、平场，并标记出坏像元、热像元、暖像元、溢出像元、宇宙线。也会运算出权重weight。
+
+## 在命令行运行
+
+```bash
+sh csst_msc_iec.sh [id]
+```
+
+## 使用python运行
+```python
+from csst_msc_instrument_effect.msc_iec import InstrumentEffectCorrection
+
+
+iec = InstrumentEffectCorrection(
+    data_path=data_filename,    # 科学图fits文件路径
+    bias_path=bias_filename,    # 本底fits文件路径(可以是路径列表)
+    dark_path=dark_filename,    # 暗场fits文件路径(可以是路径列表)
+    flat_path=flat_filename,    # 平场fits文件路径(可以是路径列表)
+    cray_path=cray_filename,    # 模拟宇宙线文件路径(因为模拟数据的问题, 暂时添加, 之后会去除)
+    output_path=output_path,    # 输出路径
+    config_path=config_path,    # crmask配置文件路径, 指的就是MSC_crmask.ini
+)
+iec.run()
+```
+生成的文件会有img, flag, weight, header四个
+
+## 使用BKGscript.py
+
+运行BKGscript.py 更改main函数内input_path与save_path与number_list三个变量
+input_path为参考源文件路径
+save_path为输出路径
+number_list为需要进行合并的参考文件编号

cleanup.sh

+#!/bin/bash
+
+# 当前0级数据处理完成，告知下阶段继续处理
+touch /next-stage/$1

csst_msc_iec.sh

+#!/bin/bash
+
+cd /
+python /home/csstpipeline/csst-msc-instrument-effect-main/instrument_effect_wrapper.py $1

csst_msc_instrument_effect/__init__.py

+# init file
\ No newline at end of file

csst_msc_instrument_effect/msc_iec.py

+import os
+from typing import List
+from astropy.io.fits import header
+from datetime import datetime
+
+import numpy as np
+from astropy.io import fits
+from deepCR import deepCR
+import torch
+
+from csst_msc_instrument_effect.msc_crmask import CRMask
+
+
+class InstrumentEffectCorrection:
+    def __init__(
+        self,
+        data_input,
+        bias_input,
+        dark_input,
+        flat_input,
+        output_path,
+        config_path,
+        cray_input=None,
+    ) -> None:
+        self.data_input = data_input
+        self.bias_input = bias_input
+        self.dark_input = dark_input
+        self.flat_input = flat_input
+        self.cray_input = cray_input
+        self.output = output_path
+        self.config_path = config_path
+        # RDNOISE
+        self.RDNOISE = "RDNOISE1"
+        # GAIN
+        self.GAIN = "GAIN1"
+        # EXPTIME
+        self.EXPTIME = "EXPTIME"
+        # output
+        if not os.path.exists(self.output):
+            os.mkdir(self.output)
+
+    def set_data(self):
+        # data
+        with fits.open(self.data_input) as hdul:
+            self.primary_header = hdul[0].header
+            self.data_raw = hdul[1].data.astype(int)
+            self.image_header = hdul[1].header
+
+    def set_bias(self, mode):
+        # bias
+        def func(path):
+            du = fits.getdata(path)
+            return du.astype(int)
+
+        if isinstance(self.bias_input, list) and len(self.bias_input) != 1:
+            bias = np.array([*map(func, self.bias_input)])
+            self.bias, var = array_combine(bias, mode)
+        elif isinstance(self.bias_input, list) and len(self.bias_input) == 1:
+            self.bias = func(self.bias_input[0])
+        elif isinstance(self.bias_input, np.ndarray) and len(self.bias_input.shape) == 2:
+            self.bias = self.bias_input
+        else:
+            self.bias = func(self.bias_input)
+
+    def set_dark(self, mode):
+        # dark
+        def func(path):
+            with fits.open(path) as hdul:
+                du = hdul[1].data
+                hu = hdul[0].header
+            du = du.astype(int)
+            du = du - self.bias
+            du = du - self.cray / self.image_header[self.GAIN]  # tmp
+            return du / hu[self.EXPTIME]
+
+        if isinstance(self.dark_input, list) and len(self.dark_input) != 1:
+            dark = np.array([*map(func, self.dark_input)])
+            self.dark, var = array_combine(dark, mode)
+        elif isinstance(self.dark_input, list) and len(self.dark_input) == 1:
+            self.dark = func(self.dark_input[0])
+        elif isinstance(self.dark_input, np.ndarray) and len(self.dark_input.shape) == 2:
+            self.dark = self.dark_input
+        else:
+            self.dark = func(self.dark_input)
+
+    def set_flat(self, mode):
+        # flat
+        def func(path):
+            with fits.open(path) as hdul:
+                du = hdul[1].data
+                hu = hdul[0].header
+            du = du.astype(int)
+            du = du - self.bias
+            du = du / hu[self.EXPTIME]
+            du = du - self.dark
+            return du / np.median(du)
+
+        if isinstance(self.flat_input, list) and len(self.flat_input) != 1:
+            flat = np.array([*map(func, self.flat_input)])
+            self.flat, var = array_combine(flat, mode)
+        elif isinstance(self.flat_input, list) and len(self.flat_input) == 1:
+            self.flat = func(self.flat_input[0])
+        elif isinstance(self.flat_input, np.ndarray) and len(self.flat_input.shape) == 2:
+            self.flat = self.flat_input
+        else:
+            self.flat = func(self.flat_input)
+
+    def set_cray(self):
+        # cray
+        if not self.cray_input == None:
+            self.cray = fits.getdata(self.cray_input).astype(int)
+        else:
+            self.cray = 0
+
+    def fix_data(self,):
+        self.data_fix0 = np.divide(
+            self.data_raw - self.bias - self.dark * self.primary_header[self.EXPTIME],
+            self.flat,
+            out=np.zeros_like(self.data_raw, float),
+            where=(self.flat != 0),
+        )
+
+    def set_flag(self,):
+        flag = np.zeros_like(self.data_raw, dtype=np.uint16)
+        # 00000001:   坏像元
+        # 因为探测器本身原因造成不可用于有效科学研究的像元, 像元响应小于0.5中值 大于1.5中值
+        med = np.median(self.flat)
+        flg = (self.flat < 0.5 * med) | (1.5 * med < self.flat)
+        flag = flag | (flg * 1)
+        # 00000010:   热像元
+        # 因为探测器本身原因造成的影响科学研究结果的像元. 像元在150秒积分时间内, 暗流计数大于探测器平均读出噪声的平方.
+        dark = self.dark.copy() * self.primary_header[self.EXPTIME]
+        dark[dark < 0] = 0
+        flg = 1 * self.image_header[self.RDNOISE] ** 2 <= dark  # 不确定是否包含 暂定包含
+        flag = flag | (flg * 2)
+        # 00000100:   暖像元
+        # 因为探测器本身原因造成的影响科学研究结果的像元. 像元的150秒积分时间内, 暗流计数大于0.5被读出噪声的平方, 但小于1被读出噪声的平方.
+        flg = (0.5 * self.image_header[self.RDNOISE] ** 2 < dark) & (
+            dark < 1 * self.image_header[self.RDNOISE] ** 2
+        )
+        flag = flag | (flg * 4)
+        # 00001000:   饱和溢出像元
+        # 饱和像元及流量溢出污染的像元.
+        flg = self.data_raw == 65535
+        flag = flag | (flg * 8)
+        # 00010000:   宇宙线像元
+        # 宇宙线污染的像元
+        crobj = CRMask(self.data_fix0, config_path=self.config_path)
+        flag_fits, data_fits = crobj.cr_mask()
+        flag_fits = flag_fits[1].data.astype(np.uint16)
+        flag = flag | (flag_fits * 16)
+        # 00100000:   卫星或者人造移动天体轨迹污染的像元.
+        pass
+        # 01000000:   鬼像污染的像元.
+        pass
+        # 10000000:   散射光污染的像元(包括dragon breath)
+        pass
+        self.flag = flag
+        # self.data_fix1 = data_fits[1].data
+
+    def set_weight(self):
+        data = self.data_fix0.copy()
+        data[self.data_fix0 < 0] = 0
+        weight_raw = 1. / (self.image_header[self.GAIN] * data + self.image_header[self.RDNOISE] ** 2)
+        bias_weight = np.std(self.bias)
+        weight = 1. / (1. / weight_raw + 1. / bias_weight) * (self.primary_header[self.EXPTIME] ** 2)
+        weight[self.flag > 0] = 0
+        self.weight = weight
+
+    def save(self):
+        data_filename = self.data_input
+        data_basename0 = os.path.basename(data_filename).replace("L0", "L1")
+        data_basename = os.path.basename(data_basename0).replace(".fits", "_img.fits")
+        data_output = os.path.join(self.output, data_basename)
+        data_header = self.image_header.copy()
+        data_header["EXTNAME"] = "img"
+        data_header[self.GAIN] = (
+            data_header[self.GAIN] * self.primary_header[self.EXPTIME]
+        )
+        data_header["BUNIT"] = "e/s"
+        data_fits = fits.HDUList(
+            [
+                fits.PrimaryHDU(header=self.primary_header),
+                fits.ImageHDU(
+                    data=self.data_fix0.astype(np.float32)
+                    / self.primary_header[self.EXPTIME],
+                    header=data_header,
+                ),
+            ]
+        )
+        data_fits.writeto(data_output, overwrite=True)
+        self.data_output = data_output
+
+        flag_output = data_output.replace("img", "flg")
+        flag_header = self.image_header.copy()
+        flag_header["EXTNAME"] = "flg"
+        flag_fits = fits.HDUList(
+            [
+                fits.PrimaryHDU(header=self.primary_header),
+                fits.ImageHDU(data=self.flag.astype(np.uint16), header=flag_header),
+            ]
+        )
+        flag_fits.writeto(flag_output, overwrite=True)
+        self.flag_output = flag_output
+
+        weight_output = data_output.replace("img", "wht")
+        weight_header = self.image_header.copy()
+        weight_header["EXTNAME"] = "wht"
+        weight_fits = fits.HDUList(
+            [
+                fits.PrimaryHDU(header=self.primary_header),
+                fits.ImageHDU(
+                    data=self.weight.astype(np.float32), header=weight_header
+                ),
+            ]
+        )
+        weight_fits.writeto(weight_output, overwrite=True)
+        self.weight_output = weight_output
+
+        header_name = data_basename.replace(".fits", ".head")
+        header_output = os.path.join(self.output, header_name)
+        hu = self.image_header
+        hu["INSTRU_S"] = (0, "instrument effect status")
+        hu["INST_TOL"] = (1234, "instrument effect operation time")
+        hu["INSTRU_V"] = ("0.1", "instrument effect version")
+        hu["INSTRU_P"] = ("test.conf", "instrument effect config file")
+        header_fits = fits.HDUList([fits.PrimaryHDU(header=hu)])
+        header_fits.writeto(header_output, overwrite=True)
+        self.header_output = header_output
+
+        # bias_output = data_output.replace("img", "bias")
+        # bias_fits = fits.HDUList([fits.PrimaryHDU(data=self.bias)])
+        # bias_fits.writeto(bias_output, overwrite=True)
+        # dark_output = data_output.replace("img", "dark")
+        # dark_fits = fits.HDUList([fits.PrimaryHDU(data=self.dark)])
+        # dark_fits.writeto(dark_output, overwrite=True)
+        # flat_output = data_output.replace("img", "flat")
+        # flat_fits = fits.HDUList([fits.PrimaryHDU(data=self.flat)])
+        # flat_fits.writeto(flat_output, overwrite=True)
+
+    def run(self):
+        self.set_data()
+        self.set_cray()
+        self.set_bias(mode="median_clip")
+        self.set_dark(mode="mean_clip")
+        self.set_flat(mode="mean_clip")
+        self.fix_data()
+        self.set_flag()
+        self.set_weight()
+        self.save()
+
+
+def array_combine(ndarray, mode="mean"):
+    """ Function to combine 3-D data array
+
+    Parameters
+    ----------
+    ndarray: array, input data cube (3D)
+    model: mean, median, sum, mean_clip, median_clip, default is mean
+    """
+    var = np.std(ndarray, axis=0)
+    ndarray_sort = np.sort(ndarray, axis=0)  ### sort the flux for each pixel
+    ndarray_cut = ndarray_sort[1:-1]  ### cut the min and max flux for each pixel
+
+    if mode == "median":
+        array = np.median(ndarray_sort, axis=0)
+    elif mode == "median_clip":
+        array = np.median(ndarray_cut, axis=0)
+    elif mode == "sum":
+        array = np.sum(ndarray_sort, axis=0)
+    elif mode == "mean":
+        array = np.mean(ndarray_sort, axis=0)
+    elif mode == "mean_clip":
+        array = np.mean(ndarray_cut, axis=0)
+    return array, var
+

demo.py

+from csst_msc_instrument_effect.msc_iec import InstrumentEffectCorrection
+
+
+
+iec = InstrumentEffectCorrection(
+    data_input="../data/test/MSC_MS_210525160000_100000022_24_raw.fits",#"../data/test/MSC_MS_210525123000_100000001_07_raw.fits",#"/data/L0/MSC_MS_210525121500_100000001_08_raw.fits",
+    bias_input="/data/BKG/MSC_CLB_210525120000_100000000_24_raw.fits",
+    dark_input="/data/BKG/MSC_CLD_210525121000_100000000_24_raw.fits",
+    flat_input="/data/BKG/MSC_CLF_210525120500_100000000_24_raw.fits",
+    cray_input="/data/BKG/MSC_CRD_210525121000_100000000_24_raw.fits",
+    output_path="/data/test/",
+    config_path="MSC_crmask.ini",
+)
+iec.run()

instrument_effect_wrapper.py

+import sys
+from datetime import datetime
+
+from csst_dfs_api.facility.level0 import Level0DataApi
+from csst_dfs_api.facility.level0prc import Level0PrcApi
+from csst_dfs_api.facility.calmerge import CalMergeApi
+from csst_msc_instrument_effect.msc_iec import InstrumentEffectCorrection
+
+l0Api = Level0DataApi()
+prcApi = Level0PrcApi()
+calApi = CalMergeApi()
+
+# -----------------------------------------------------
+# get data
+l0_id = sys.argv[1]
+
+rec = l0Api.get(level0_id=l0_id)
+print(rec)
+print("Find Level0 Data :" + rec.data.file_path)
+
+bias_rec = calApi.get_latest_by_l0(level0_id=l0_id, ref_type="bias")
+dark_rec = calApi.get_latest_by_l0(level0_id=l0_id, ref_type="dark")
+flat_rec = calApi.get_latest_by_l0(level0_id=l0_id, ref_type="flat")
+cray_rec = calApi.get_latest_by_l0(level0_id=l0_id, ref_type="cray")  # temp
+
+print("Find Level0 Bias :" + bias_rec.data.file_path)
+print("Find Level0 Dark :" + dark_rec.data.file_path)
+print("Find Level0 Flat :" + flat_rec.data.file_path)
+# ------------------------------------------------------
+# to-do data processing here
+iec = InstrumentEffectCorrection(
+    data_input=rec.data.file_path,
+    bias_input=bias_rec.data.file_path,
+    dark_input=dark_rec.data.file_path,
+    flat_input=flat_rec.data.file_path,
+    output_path="/home/csstpipeline/data/L05_test/",
+    config_path="/home/csstpipeline/csst-msc-instrument-effect-master/MSC_crmask.ini",
+    cray_input=cray_rec.data.file_path,  # temp
+)
+iec.run()
+
+params_file_path = iec.config_path
+result_file_path = ', '.join([
+    iec.data_output,
+    iec.flag_output,
+    iec.weight_output,
+    iec.header_output,
+])
+dt = datetime.now()
+# ------------------------------------------------------
+# write result
+rec = prcApi.write(
+    level0_id=l0_id,
+    pipeline_id="P1",
+    prc_module="MSC-IE",  # csst-msc-instrument-effect
+    params_file_path=params_file_path,
+    prc_status=0,
+    prc_time=dt.strftime('%Y-%m-%d %H:%M:%S'),
+    result_file_path=result_file_path)
+print('Write Result:', rec)

requirements.txt

+astropy
+ccdproc
+deepCR
\ No newline at end of file

run.py

+import os
+from glob import glob
+import argparse
+from astropy.io.fits import getdata, getheader
+from time import time
+
+from csst_msc_instrument_effect.msc_iec import InstrumentEffectCorrection
+
+
+def run(data_path, output, ref, conf):
+    header = getheader(data_path)
+    number = header['DETECTOR'][-2:]
+    # number = os.path.basename(data_path).split('_')[4]  # 原来版本
+    number_list = ['06', '07', '08', '09', '11', '12', '13', '14', '15',
+                   '16', '17', '18', '19', '20', '22', '23', '24', '25']
+    if number not in number_list:
+        print(number, 'not in')
+        return
+    bias = getdata(glob(os.path.join(ref, "*CLB*_" + number + '_*'))[0])
+    dark = getdata(glob(os.path.join(ref, "*CLD*_" + number + '_*'))[0])
+    flat = getdata(glob(os.path.join(ref, "*CLF*_" + number + '_*'))[0])
+    iec = InstrumentEffectCorrection(
+        data_input=data_path,
+        bias_input=bias,
+        dark_input=dark,
+        flat_input=flat,
+        output_path=output,
+        config_path=conf
+    )
+    iec.run()
+    print(number, 'done')
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        '-i', '--input', help='fits file path or folder path', required=True)
+    parser.add_argument(
+        '-o', '--output', help='output path', default='/data/L05/')
+    parser.add_argument(
+        '-r', '--ref', help='ref path', default='/data/ref/')
+    parser.add_argument(
+        '-c', '--conf', help='config path', default="/home/csstpipeline/csst-msc-instrument-effect-main/MSC_crmask.ini")
+    args = parser.parse_args()
+    if os.path.isdir(args.input):
+        for path in glob(os.path.join(args.input.strip(), '*_MS_SCI_*.fits')):
+        #for path in glob(os.path.join(args.input.strip(), '*_MS_*.fits')):
+            run(path, args.output, args.ref, args.conf)
+    else:
+        print('run test')
+        begin = time()
+        run(args.input.strip(), args.output.strip(), args.ref.strip(), args.conf.strip())
+        print(time() - begin)
+
+
+if __name__ == '__main__':
+    # python run.py -i "/data/test20211012/30s/MSC_MS_210525220000_100000020_13_raw.fits" -o "/data/test20211012/output/150s/"
+    # python run.py -i "/data/test20211012/150s/" -o '/data/test20211012/output/150s/'
+    # python run.py -i "/data/cali_20211012/L0/150s/" -o '/data/cali_20211012/output/'
+    # python run.py -i '/data/test20211012/30s/' -o '/data/test20211012/output/30s/'
+    # python run.py -i "/data/sim_data/20220413/MSC_0000100/CSST_MSC_MS_CRS_20270810081950_20270810082220_100000100_13_L0_1.fits" -o "/home/csstpipeline/test/"
+    main()

run.sh

+# 运行: sh /home/csstpipeline/csst-msc-instrument-effect-main/run.sh [输入文件路径] [输出文件路径] [背景文件路径] [配置文件路径]
+# 结果默认保存在/data/cali_20211012/output/
+# 背景文件路径默认在/data/test20211012/ref/
+# 配置文件路径默认在/home/csstpipeline/csst-msc-instrument-effect-main/MSC_crmask.ini
+input=${1}
+output=${2:-'/data/cali_20211012/output/'}
+ref=${3:-'/data/test20211012/ref/'}
+conf=${4:-"/home/csstpipeline/csst-msc-instrument-effect-main/MSC_crmask.ini"}
+echo ${input}
+echo ${output}
+echo ${ref}
+echo ${conf}
+python /home/csstpipeline/csst-msc-instrument-effect-main/run.py -i ${input} -o ${output} -r ${ref} -c ${conf}

run_18.py

+import os
+os.environ['KMP_DUPLICATE_LIB_OK'] = 'TRUE'
+from glob import glob
+from astropy.io.fits import getdata, getheader
+from concurrent.futures import ThreadPoolExecutor
+from time import time
+
+from csst_msc_instrument_effect.msc_iec import InstrumentEffectCorrection
+
+
+def run(data_path, output, ref, conf):
+    header = getheader(data_path)
+    number = header['DETECTOR'][-2:]
+    number_list = ['06', '07', '08', '09', '11', '12', '13', '14', '15',
+                   '16', '17', '18', '19', '20', '22', '23', '24', '25']
+    if number not in number_list:
+        print(number, 'not in')
+        return
+    bias = getdata(glob(os.path.join(ref, "*CLB*_" + number + '_*'))[0])
+    dark = getdata(glob(os.path.join(ref, "*CLD*_" + number + '_*'))[0])
+    flat = getdata(glob(os.path.join(ref, "*CLF*_" + number + '_*'))[0])
+    iec = InstrumentEffectCorrection(
+        data_input=data_path,
+        bias_input=bias,
+        dark_input=dark,
+        flat_input=flat,
+        output_path=output,
+        config_path=conf
+    )
+    iec.run()
+    print(number, 'done')
+
+
+if __name__ == '__main__':
+    input_path = 'D:/Desktop/data/'
+    output_path = 'D:/Desktop/data/output/'
+    ref_path = r'D:/Desktop/data/ref/'
+    conf_path = 'D:/Desktop/test/csst-msc-instrument-effect/CSST_2021-12-30_CCD23_epoch20.pth'
+    with ThreadPoolExecutor(18) as tpe:
+        for path in glob(os.path.join(input_path, '*_MS_SCI_*.fits')):
+            tpe.submit(run, path, output_path, ref_path, conf_path)

setup.py

+from setuptools import setup
+
+
+setup(
+    name='csst_msc_instrument_effect',
+    author='xiezhou',
+    version='0.1',
+    packages=['csst_msc_instrument_effect'],
+    entry_points={
+        'console_scripts': [
+            'csst-msc-iec=csst_msc_instrument_effect.msc_iec:main'
+        ]
+    }
+)