1 """
2 Evolve the CCD, generating fake dark current map using the
3 smoothed-broken-powerlaw spectral model fitting (from fit_evol.py).
4
5 First use the dates of actual dark current cals (for direct comparison), then
6 go once per year.
7 """
8
9 import dark_models
10 from Chandra.Time import DateTime
11 import matplotlib.ticker as ticker
12 import os
13 import re
14 import ParseTable
15 import numpy as np
16
18 """Plot 1999223 actual, degradation model, 2008273 actual and 2008273 predicted"""
19 figure(1, figsize=(7,5))
20 clf()
21 plot_darkcal('aca_dark_cal/1999223/imd.fits', -10, '-g', label='1999 SSD closed', step=True)
22 plot_darkcal('aca_dark_cal/2008273/imd.fits', -19, '-b', label='2008 actual', step=True)
23 plot_darkcal('worst/from1999223/2008273.fits', -19, '.b', label='2008 predicted')
24 plot_darkcal('worst/from1999223/2019182.fits', -18, '.r', label='2019 predicted')
25
26 pars = dark_models.nompars(dyear=1.0)
27 darkmodel = dark_models.smooth_twice_broken_pow(pars, dark_models.xall)
28 plot(dark_models.xall, darkmodel, label='Degradation model')
29
30 xlim(1, 1e4)
31 ylim(0.5, 3e5)
32 xlabel('Dark current (e-/sec)')
33 ylabel('Number per bin')
34 title('Dark current distributions')
35 legend(loc=3)
36 savefig('dark_overlays.png')
37
38 -def plot_darkcal(darkfile, Tccd, plotsymcol='-', label=None, step=None):
39 """
40 Plot a real dark cal, scaled from temp C{Tccd} to -19
41
42 @param darkfile: FITS file containing a dark current map image.
43 @param Tccd: CCD temperature (degC)
44 @param plotsymcol: Matplotlib symbol/color specifierr
45 @param label: Matplotlib plot label
46 @param step: Matplotlib plot step specifier
47
48 @return: C{(dark, x, y)}
49 - C{dark}: Dark current map (scaled to temperature C{Tccd})
50 - C{x, y}: Histogram values plotted
51 """
52 hdus = pyfits.open(darkfile)
53
54 dark = hdus[0].data.flatten() / dark_models.temp_scalefac(Tccd)
55 x, y = plot_dark(dark, label, plotsymcol, step=step)
56 return dark, x, y
57
58 -def plot_dark(dark, label=None, plotsymcol='.', step=False):
59 """Plot dark cal data as a histogram. 'dark' is a flat array."""
60 dark = dark + np.random.uniform(-0.5, 0.5, len(dark))
61 y, x = histogram(dark, dark_models.xbins, new=True)
62 if step:
63 x, y = make_steps(x, y)
64 else:
65 x = (x[1:] + x[:-1])/2
66 y = np.array(y, dtype=float)
67 y[y < 1e-5] = 0.1
68 loglog(x, y, plotsymcol, label=label)
69 return x, y
70
72 """For x as bin edges (n+1 vals) and y as bin vals (n vals), make the corresponding
73 arrays that can be plotted as a line."""
74 lx = x[:-1]
75 rx = x[1:]
76 newx = np.array([lx, rx]).flatten(2)
77 newy = np.array([y, y]).flatten(2)
78 return newx, newy
79
81 mxdates = (DateTime(x).mxDateTime for x in alldarkcals['date'])
82 yrs = [x.year + x.day_of_year / yeardays for x in mxdates]
83 fracs = alldarkcals[thresh] / npix
84
85 figure(1, figsize=(6, 4.5))
86 ax = subplot(1, 1, 1)
87 plot(yrs, fracs, '.')
88 formatter = ticker.FormatStrFormatter('%d')
89 ax.xaxis.set_major_formatter(formatter)
90 plot(warmpix['year'], np.array(warmpix[thresh]) / npix)
91 xlim(2000, 2020)
92 ylim(-0.005, 0.251)
93 xlabel('Year')
94 ylabel('Fraction')
95 title('Warm pixel fraction (Nominal and Worst)')
96 draw()
97
98
99 npix = 1024.0**2
100 secs2k = DateTime('2000:001:00:00:00').secs
101 sec2year = 1 / (86400. * 365.25)
102 yeardays = 365.25
103 alldarkcals = ParseTable.parse_table('darkcal_stats.dat')
104
106 """Make dark current maps by degrading the CCD using the best-fit model
107 from fit_evol.py"""
108
109 date0 = DateTime('1999-05-23T18:00:00')
110
111 if initial == 'pristine':
112
113
114 dark = dark_models.pristine_ccd()
115 darkcals = alldarkcals
116
117 elif re.match('from|zero', initial):
118 darkmap, year, doy = re.match(r'(\S+)(\d{4})(\d{3})', initial).groups()
119 yeardoy = year + ':' + doy
120 date0 = DateTime(yeardoy)
121 darkcals = alldarkcals[alldarkcals['date'] > yeardoy]
122
123 if darkmap == 'zero':
124 dark = dark_models.zero_dark_ccd()
125 else:
126 ok = alldarkcals['date'] == yeardoy
127 tccd = alldarkcals[ok]['tccd'][0]
128 scalefac = dark_models.temp_scalefac(tccd)
129 print 'Scaling dark cal', yeardoy, 'by', '%.2f' % scalefac
130 hdus = pyfits.open(os.path.join('aca_dark_cal', year+doy, 'imd.fits'))
131 dark = hdus[0].data.flatten() / scalefac
132
133 warmpix = dict(year=[], n100=[], n200=[], n2000=[])
134
135 print 'Calculating for', case, 'case starting from', initial, 'CCD'
136
137
138 dates = [x['date'] for x in darkcals] + [str(yr) + ':182' for yr in range(2009, 2020)]
139
140 datelast = date0
141 for datestr in dates:
142 date = DateTime(datestr)
143 datemx = date.mxDateTime
144 print date.date
145
146
147 dyear = (datemx - datelast.mxDateTime).days / yeardays
148 dark_models.degrade_ccd(dark, dyear)
149
150
151 try:
152 ok = darkcals['date'] == datestr
153 tccd = darkcals[ok]['tccd'][0]
154 except IndexError:
155 if datemx.year > 2014:
156 if case == 'nominal':
157 tccd = -18
158 else:
159 tccd = -15
160 else:
161 tccd = -19
162
163
164 scalefac = dark_models.temp_scalefac(tccd)
165 dark_tccd = dark * scalefac
166 outdir = os.path.join(case, initial)
167 if not os.path.exists(outdir):
168 os.makedirs(outdir)
169 outfile = os.path.join(outdir, '%04d%03d.fits' % (datemx.year, datemx.day_of_year))
170
171 if os.path.exists(outfile):
172 os.unlink(outfile)
173 hdu = pyfits.PrimaryHDU(np.float32(dark_tccd.reshape(1024,1024)))
174 hdu.writeto(outfile)
175
176
177 warmpix['year'].append(datemx.year + datemx.day_of_year / yeardays)
178 warmpix['n100'].append(len(where(dark_tccd > 100)[0]))
179 warmpix['n200'].append(len(where(dark_tccd > 200)[0]))
180 warmpix['n2000'].append(len(where(dark_tccd > 2000)[0]))
181
182 datelast = date
183
184
185
186
187
188
189 return warmpix
190