update csst/msc cosmic ray model

csst_cpic_sim/camera.py

+import os
 import math
 import numpy as np
 import scipy.ndimage as nd
@@ -7,6 +8,7 @@ from .config import config, S
 from .utils import region_replace, random_seed_select
 from .io import log
 from .optics import filter_throughput
+from .msc_cr import produceCR_Map
 
 cpism_refdata = config['cpism_refdata']
 MAG_SYSTEM = config['mag_system']
@@ -76,16 +78,19 @@ class CosmicRayFrameMaker(object):
         The pitch of the camera pixel in um.
     cr_rate : float
         The cosmic ray rate per second per cm2.
-
     """
 
-    def __init__(self) -> None:
+    def __init__(self, msc_model=False) -> None:
         self.tmp_size = [7, 101]
         self.freq_power = -0.9
         self.trail_std = 0.1
         self.depth = 10  # um
         self.pitch = 13  # um
         self.cr_rate = 1  # particle per s per cm2 from Miles et al. 2021
+        data_folder = os.path.dirname(__file__)
+        cr_path = os.path.join(data_folder, 'data/cr_path.txt')
+        self.attachedSize = np.loadtxt(cr_path)
+        self.msc_model = msc_model
 
     def make_CR(self, length, sigma, seed=-1):
         """
@@ -249,6 +254,15 @@ class CosmicRayFrameMaker(object):
         seed = random_seed_select(seed=seed)
         log.debug(f"cr position seed: {seed}")
 
+        if self.msc_model:
+            cr_map, cr_event_num = produceCR_Map(
+                xLen=shape[1], yLen=shape[0],
+                exTime=expt,
+                cr_pixelRatio=0.003*expt/600.,  # 5e-6/pixel , 1.3e-6 for cpic package
+                gain=1,  # gain = 1, means the unit of output image is electron
+                attachedSizes=self.attachedSizes)
+            return cr_map
+
         for i in range(len(cr_array)):
             cr = cr_array[i]
             x = np.random.rand() * sz[1]

csst_cpic_sim/config.py

@@ -92,10 +92,11 @@ config['pysyn_refdata'] = f'{cpism_refdata}/starmodel/grp/redcat/trds'
 config['catalog_folder'] = f'{cpism_refdata}/demo_catalog'
 config['csst_format'] = True
 config['nsample'] = 1
+config['msc_cr_model'] = False
 
 update_able_keys = [
     'apm_file', 'actuator_file', 'aberration',
-    'log_dir', 'log_level', 'catalog_folder',
+    'log_dir', 'log_level', 'catalog_folder', 'msc_cr_model'
     'nsample', 'csst_format', 'output', 'check_fits_header'
 ]
 

csst_cpic_sim/data/cr_path.txt

+# Graph from wfc-cr-attach, page 1
+0.00000	0.004684
+0.5031	0.004684
+0.5283	0.01873
+1.509	0.01873
+1.534	0.09327
+2.490	0.09327
+2.515	0.1034
+3.496	0.1034
+3.522	0.2440
+4.503	0.2440
+4.528	0.1107
+5.509	0.1107
+5.534	0.1013
+6.490	0.1013
+6.515	0.06090
+7.496	0.06090
+7.522	0.05834
+8.503	0.05834
+8.528	0.03875
+9.509	0.03875
+9.534	0.03066
+10.49	0.03066
+10.51	0.01788
+11.47	0.01788
+11.49	0.01831
+13.50	0.01831
+13.53	0.01235
+13.53	0.01235
+14.49	0.01235
+14.51	0.01064
+15.49	0.01064
+15.52	0.008091
+16.50	0.008091
+16.52	0.004684
+17.48	0.004684
+17.53	0.003833
+18.49	0.003833
+18.51	0.005536
+19.47	0.005536
+19.52	0.004684
+20.00	0.004684
\ No newline at end of file

csst_cpic_sim/main.py

@@ -342,6 +342,11 @@ def observation_simulation_from_config(obs_file, config_file):
     csst_format = default_config['csst_format']
     nsample = default_config['nsample']
 
+    if default_config['msc_cr_model']:
+        crmaker = CosmicRayFrameMaker(msc_model=True)
+    else:
+        crmaker = CosmicRayFrameMaker()
+
     obs_file = os.path.abspath(obs_file)
     file_list = []
 
@@ -417,6 +422,7 @@ def observation_simulation_from_config(obs_file, config_file):
                     output=output,
                     nsample=nsample,
                     csst_format=csst_format,
+                    crmaker=crmaker,
                     dataset=dataset,
                     prograss_bar=True)
             except Exception as e:

csst_cpic_sim/msc_cr.py

+import numpy as np
+import math
+
+# ---------- For Cosmic-Ray Simulation ------------
+# ---------- Zhang Xin ---------------
+# ---------- From MSC/CSST simulation package-------------
+# ---------- The initial code can be found at
+#        https://csst-tb.bao.ac.cn/code/liudezi/csst_msc_sim/-/blob/release_v2.0/ObservationSim/Instrument/Chip/Effects.py
+
+
+def getYValue(collection, x):
+    index = 0
+    if (collection.shape[1] == 2):
+        while (x > collection[index, 0] and index < collection.shape[0]):
+            index = index + 1
+        if (index == collection.shape[0] or index == 0):
+            return 0
+
+        deltX = collection[index, 0] - collection[index-1, 0]
+        deltY = collection[index, 1] - collection[index-1, 1]
+
+        if deltX == 0:
+            return (collection[index, 1] + collection[index-1, 1])/2.0
+        else:
+            a = deltY/deltX
+            return a * (x - collection[index-1, 0]) + collection[index-1, 1]
+    return 0
+
+
+def selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size):
+    normalRay = 0.90
+    nnormalRay = 1-normalRay
+    max_nrayLen = 100
+    pixelNum = int(xLen * yLen * cr_pixelRatio * normalRay)
+    pixelNum_n = int(xLen * yLen * cr_pixelRatio * nnormalRay)
+    CRPixelNum = 0
+
+    maxValue = max(attachedSizes[:, 1])
+    maxValue += 0.1
+
+    cr_event_num = 0
+    CRs = np.zeros(pixelNum)
+    while (CRPixelNum < pixelNum):
+        x = CR_max_size * np.random.random()
+        y = maxValue * np.random.random()
+        if (y <= getYValue(attachedSizes, x)):
+            CRs[cr_event_num] = np.ceil(x)
+            cr_event_num = cr_event_num + 1
+            CRPixelNum = CRPixelNum + round(x)
+
+    while (CRPixelNum < pixelNum + pixelNum_n):
+        nx = np.random.random()*(max_nrayLen-CR_max_size)+CR_max_size
+        CRs[cr_event_num] = np.ceil(nx)
+        cr_event_num = cr_event_num + 1
+        CRPixelNum = CRPixelNum + np.ceil(nx)
+
+    return CRs[0:cr_event_num]
+
+
+def defineEnergyForCR(cr_event_size, seed=12345):
+    import random
+    sigma = 0.6 / 2.355
+    mean = 3.3
+    random.seed(seed)
+    energys = np.zeros(cr_event_size)
+    for i in np.arange(cr_event_size):
+        energy_index = random.normalvariate(mean, sigma)
+        energys[i] = pow(10, energy_index)
+    return energys
+
+
+def convCR(CRmap=None, addPSF=None, sp_n=4):
+    sh = CRmap.shape
+
+    # sp_n = 4
+    subCRmap = np.zeros(np.array(sh)*sp_n)
+    pix_v0 = 1/(sp_n*sp_n)
+    for i in np.arange(sh[0]):
+        i_st = sp_n*i
+        for j in np.arange(sh[1]):
+            if CRmap[i, j] == 0:
+                continue
+            j_st = sp_n*j
+            pix_v1 = CRmap[i, j]*pix_v0
+            for m in np.arange(sp_n):
+                for n in np.arange(sp_n):
+                    subCRmap[i_st+m, j_st + n] = pix_v1
+
+    m_size = addPSF.shape[0]
+
+    subCRmap_n = np.zeros(np.array(subCRmap.shape) + m_size - 1)
+
+    for i in np.arange(subCRmap.shape[0]):
+        for j in np.arange(subCRmap.shape[1]):
+            if subCRmap[i, j] > 0:
+                convPix = addPSF*subCRmap[i, j]
+                subCRmap_n[i:i+m_size, j:j+m_size] += convPix
+
+    CRmap_n = np.zeros((np.array(subCRmap_n.shape)/sp_n).astype(np.int32))
+    sh_n = CRmap_n.shape
+
+    for i in np.arange(sh_n[0]):
+        i_st = sp_n*i
+        for j in np.arange(sh_n[1]):
+            p_v = 0
+            j_st = sp_n*j
+            for m in np.arange(sp_n):
+                for n in np.arange(sp_n):
+                    p_v += subCRmap_n[i_st+m, j_st + n]
+
+            CRmap_n[i, j] = p_v
+
+    return CRmap_n
+
+
+def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=20210317):
+    # Return: an 2-D numpy array
+    # attachedSizes = np.loadtxt('./wfc-cr-attachpixel.dat');
+    np.random.seed(seed)
+    CR_max_size = 20.0
+    cr_size = selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size)
+
+    cr_event_size = cr_size.shape[0]
+    cr_energys = defineEnergyForCR(cr_event_size, seed=seed)
+
+    CRmap = np.zeros([yLen, xLen])
+
+    # produce conv kernel
+    from astropy.modeling.models import Gaussian2D
+    o_size = 4
+    sp_n = 8
+
+    m_size = o_size*sp_n+1
+    m_cen = o_size*sp_n/2
+    sigma_psf = 0.2*sp_n
+    addPSF_ = Gaussian2D(1, m_cen, m_cen, sigma_psf, sigma_psf)
+    yp, xp = np.mgrid[0:m_size, 0:m_size]
+    addPSF = addPSF_(xp, yp)
+    convKernel = addPSF/addPSF.sum()
+
+    # ---------------------------------
+    for i in np.arange(cr_event_size):
+        xPos = round((xLen - 1) * np.random.random())
+        yPos = round((yLen - 1) * np.random.random())
+        cr_lens = int(cr_size[i])
+        if cr_lens == 0:
+            continue
+        pix_energy = cr_energys[i]/gain/cr_lens
+        pos_angle = 1/2*math.pi*np.random.random()
+
+        crMatrix = np.zeros([cr_lens+1, cr_lens + 1])
+
+        for j in np.arange(cr_lens):
+            x_n = int(np.cos(pos_angle)*j - np.sin(pos_angle)*0)
+            if x_n < 0:
+                x_n = x_n + cr_lens+1
+            y_n = int(np.sin(pos_angle)*j + np.cos(pos_angle)*0)
+            if x_n < 0 or x_n > cr_lens or y_n < 0 or y_n > cr_lens:
+                continue
+            crMatrix[y_n, x_n] = pix_energy
+
+        crMatrix_n = convCR(crMatrix, convKernel, sp_n)
+
+        xpix = np.arange(crMatrix_n.shape[0]) + int(xPos)
+        ypix = np.arange(crMatrix_n.shape[1]) + int(yPos)
+
+        sh = CRmap.shape
+        okx = (xpix >= 0) & (xpix < sh[1])
+        oky = (ypix >= 0) & (ypix < sh[0])
+
+        sly = slice(ypix[oky].min(), ypix[oky].max()+1)
+        slx = slice(xpix[okx].min(), xpix[okx].max()+1)
+        CRmap[sly, slx] += crMatrix_n[oky, :][:, okx]
+    return CRmap.astype(np.int32), cr_event_size

example/run_config.yaml

@@ -41,3 +41,4 @@ nsample: 1
 check_fits_header: true
 output: ./example_output/
 actuator_file: ${cpism_refdata}/optics/dark_zone_0011_volt.fits
+msc_cr_model: true

setup.py

@@ -28,7 +28,7 @@ setuptools.setup(
     package_dir={'csst_cpic_sim': 'csst_cpic_sim'},
     include_package_data=True,
     package_data={
-        'csst_cpic_sim': ['data/*.yaml', 'data/*.toml']
+        'csst_cpic_sim': ['data/*.yaml', 'data/*.toml', 'data/*.txt']
     },
 
     install_requires=requirements,

tests/test_main.py

@@ -54,8 +54,9 @@ class TestMain(unittest.TestCase):
             text = fid.read()
 
         text = text.replace('<output>', output)
+        config_text = text.replace('<msc_model>', "True")
         with open(output+'/run_config_new.yaml', 'w') as fid:
-            fid.write(text)
+            fid.write(config_text)
 
         observation_simulation_from_config(
             os.path.join(cases, '05_sci.yaml'),
@@ -66,6 +67,7 @@ class TestMain(unittest.TestCase):
         self.assertEqual(len(file), 4)
 
         clear_output()
+
         observation_simulation_from_config(
             os.path.join(cases, '05_sci.yaml'),
             None
@@ -74,6 +76,20 @@ class TestMain(unittest.TestCase):
         self.assertEqual(len(file), 1)
         self.assertEqual(file[0][:5], 'SCI')
 
+        clear_output()
+
+        config_text = text.replace('<msc_model>', "False")
+        with open(output+'/run_config_new.yaml', 'w') as fid:
+            fid.write(config_text)
+
+        observation_simulation_from_config(
+            os.path.join(cases, '05_sci.yaml'),
+            os.path.join(output, 'run_config_new.yaml')
+        )
+
+        file = os.listdir(output)
+        self.assertEqual(len(file), 4)
+
     def test_full_run_cmdline(self):
         clear_output()
         main(['run', os.path.join(cases, '05_sci.yaml')])

tests/testcases/run_config.yaml

@@ -40,3 +40,4 @@ csst_format: true
 nsample: 5
 check_fits_header: true
 output: <output>
+msc_cr_model: <msc_model>