-
Notifications
You must be signed in to change notification settings - Fork 0
/
example_script_pc.py
88 lines (70 loc) · 2.64 KB
/
example_script_pc.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
# -*- coding: utf-8 -*-
"""Example photon counting script."""
import os
from pathlib import Path
import warnings
import numpy as np
import matplotlib.pyplot as plt
from emccd_detect.emccd_detect import EMCCDDetect
from PhotonCount.corr_photon_count import get_count_rate
here = Path(os.path.abspath(os.path.dirname(__file__)))
def imagesc(data, title=None, vmin=None, vmax=None, cmap='viridis',
aspect='equal', colorbar=True):
"""Plot a scaled colormap."""
fig, ax = plt.subplots()
im = ax.imshow(data, vmin=vmin, vmax=vmax, cmap=cmap, aspect=aspect)
if title:
ax.set_title(title)
if colorbar:
fig.colorbar(im, ax=ax)
return fig, ax
if __name__ == '__main__':
# Specify relevant detector properties
emccd = EMCCDDetect(
em_gain=5000.,
full_well_image=60000., # e-
full_well_serial=100000., # e-
dark_current=3e-5, # e-/pix/s
cic=1.3e-3, # e-/pix/frame
read_noise=100., # e-/pix/frame
bias=10000., # e-
qe=0.8,
cr_rate=0., # hits/cm^2/s
pixel_pitch=13e-6,
eperdn=7.,
nbits=14,
numel_gain_register=604
)
fluxmap = np.load(Path(here, 'fluxmap.npy'))
# Simulate frames
# Set frametime to get a good output, like 0.1 phot/pix or less
frametime = 10 # s
frame_e_list = []
frame_e_dark_list = []
nframes = 100
for i in range(nframes):
# Simulate bright
frame_dn = emccd.sim_sub_frame(fluxmap, frametime)
# Simulate dark
frame_dn_dark = emccd.sim_sub_frame(np.zeros_like(fluxmap), frametime)
# Convert from dn to e- and bias subtract
frame_e = frame_dn * emccd.eperdn - emccd.bias
frame_e_dark = frame_dn_dark * emccd.eperdn - emccd.bias
frame_e_list.append(frame_e)
frame_e_dark_list.append(frame_e_dark)
frame_e_cube = np.stack(frame_e_list)
# Photon count, co-add, and correct for photometric error
thresh = 500. # see warnings below
if emccd.read_noise <=0:
warnings.warn('read noise should be greater than 0 for effective '
'photon counting')
if thresh < 4*emccd.read_noise:
warnings.warn('thresh should be at least 4 or 5 times read_noise for '
'accurate photon counting')
mean_rate = get_count_rate(frame_e_cube, thresh, emccd.em_gain)
# Plot images
imagesc(fluxmap, 'input flux map')
imagesc(np.max(frame_e_cube, axis=0), 'max frame (pixel-by-pixel)')
imagesc(frame_e_cube[0], 'random frame')
imagesc(mean_rate, 'get_count_rate')
plt.show()