<jupyter_start><jupyter_text>Investigating Detectors=====**Version 0.1**Understanding the behavior of the CCDs in a camera requires digging deep into calibration exposures. That is where you can uncover effects that might not be noticible in on-sky exposures, but may subtly contaminate the data if left uncorrected. It is also how camera engineering teams optimize and debug the performance of the camera when it's still in the lab. We're going to look at two test exposures taken with one of the Rubin Observatory CCDs. They're both biases; each image has a zero second exposure time and the detector was not illuminated.Please download a tarball of the images for this notebook: [investigating_detectors.tar.gz](https://northwestern.box.com/s/nfa3owg2wpyknca415gl017ny6phx6w6). As a reminder, you can unpack these files via `tar -zxvf investigating_detectors.tar.gz`* * *By <NAME> (University of Washington)<jupyter_code>from astropy.io import fits

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

matplotlib.rcParams['figure.dpi'] = 120<jupyter_output><empty_output><jupyter_text>Photon Transfer Curve====1) Simulated Images-----The "Photon Transfer Curve" is the name given to the relationship between the signal level and the noise level in a sensor. We're going to do a few experiments to show how it works in principle, and then we'll look at some real images and make some diagnostic measurements.First we need a model of the noise in CCD image. I'm going to give this to you so we all start out on the same page<jupyter_code>def simulated_image(signal_level, read_noise, gain):
    """
    Return a 1-D simulated "image" with the noise properties of
    a CCD sensor. The image is always 1000 pixels long.
    
    signal_level is the mean number of electrons in each pixel.
    read_noise is the noise of the readout amplifier, in electrons.
    gain is the number of electrons per ADU.
    
    """
    return (1/gain) * (read_noise*np.random.randn(1000) +  np.random.poisson(signal_level, size=1000))<jupyter_output><empty_output><jupyter_text>Before diving in to programming, take a careful look at the components in the simulated image. What are the two noise sources, and why do they have that functional form? We're going to be looking a lot at the image "gain"; does it make sense how that is applied?Let's make some simulations. What we want to do is loop over a set of input levels light levels, from zero to "full well" capacity (on order of 10,000 electrons). For each simulated image, we want to measure the mean signal level (because that's what we see as users of a CCD) and the standard deviation of that image. Save those in two lists, but at the end convert those back to numpy arrays to make downstream usage easier.For right now, set the read noise to 5, and the gain to 0.8.<jupyter_code># Question
noise_level_list = []
measured_signal_level_list = []
input_signal_levels = []

for i in range(10000):
    input_signal_levels.append(i)
    noise = simulated_image(i, 5, 0.8)
    noise_level_list.append(np.std(noise))
    measured_signal_level_list.append(np.mean(noise + i))
noise_levels = np.array(noise_level_list)
measured_signal_levels = np.array(measured_signal_level_list)<jupyter_output><empty_output><jupyter_text>Plot the noise vs. the measured signal level, on a log-log plot.<jupyter_code># Question
plt.plot(measured_signal_levels, noise_levels)
# Complete
plt.xscale('log')
plt.yscale('log')
plt.ylabel("RMS")
plt.xlabel("Measured Signal Level")<jupyter_output><empty_output><jupyter_text>What is the behavior you see? What are the two different noise regimes?Fit a straight line to the "bright" portion of the data (high signal levels) and print the resulting coefficients. Remember that you're looking at a log-log plot, and so you want to fit the logs of the variables. You can add this to the plot in the cell above.Why does the line have that value of the slope? Now we're going to plot something slightly different. Plot the *variance* this time, and on a linear plot instead of log-log (again vs. measured signal level). Fit a straight line to the data (in linear space) and print the coefficients. Also print the *reciprocal* of the slope.Where did this slope come from?<jupyter_code># Question

plt.plot(measured_signal_levels, noise_levels**2)
    

plt.ylabel("Variance")
plt.xlabel("Measured Signal Level")<jupyter_output><empty_output><jupyter_text>The slope here is related to the gain (either proportionally or inversely, depending on how one chooses to define gain). This can be summarized as$$ \frac{1}{\textrm{gain}} =  \langle \frac{\textrm{Variance}}{\textrm{Mean Signal Level}} \rangle $$It's a clever and useful trick, or at least it seems like a trick, because the standard deviation plot wasn't affected by the gain at all. Go back and try varying the gain and re-run the plots, and you'll see what does and doesn't change.One way to think of it is that the measured signal level is affected linearly by the gain, but the variance is affected by the square of the gain. Dividing these two gives you a linear relation back, but when dividing the square root of the variance, the gain cancels out. 2) Looking at real bias frames------Remember that a bias frame is an image that is exposed for zero seconds; it's just immediately read-out without being exposed to light. You might think that is a pretty boring image, particularly if you're at the telescope and getting ready for a night of observing. But to a camera engineer, bias frames hold lots of information about how the camera is operating.We're going to look at images from one example LSST sensor; this was taken on a test stand and not the actual camera, so don't take it as representative of real camera performance.Our first step is, as usual, to look at the image and make sure it seems reasonable.<jupyter_code>bias1_file = fits.open("00258334360-S10-det003.fits")
bias1_data = bias1_file[1].data
# Question
plt.imshow(bias1_data)
plt.colorbar()<jupyter_output><empty_output><jupyter_text>Notice that when we plotted `bias1_file[1].data`, the image we get is 2048 by 576 pixels. Because LSST sensors have 16 separate amplifiers, the data from each one of them is put in a different "header data unit" (HDU) in the FITS file. You can get to them by substituting n in `bias1_file[n]`, where n is the amplifier number. 3) Looking for structure------The bias looks mostly like Gaussian noise, but if you look carefully some parts of the image look like they have some "structure". Let's make a few plots: try plotting the mean of the data *along columns* in one plot, and along *rows* in another.Start with just a single amplifier, but if you like you can learn more by plotting each amplifier as a different line. Hint: the amplifiers each have different mean levels that you probably want to subtract off.<jupyter_code>len(bias1_file)
for hdu in bias1_file[1:16]:
    plt.plot(np.mean(hdu.data, axis=0) - np.mean(hdu.data))
for hdu in bias1_file[1:16]:
    plt.plot(np.mean(hdu.data, axis=1) - np.mean(hdu.data))
plt.plot(np.mean(bias1_data, axis=1))<jupyter_output><empty_output><jupyter_text>These "simple" bias frames turn out to have a lot of structure in them, particularly at the start of columns. This isn't something we can dive much further into, because it's really an electronics problem (that was known about at the time). It's also worth noting that it's fractionally a small effect. We will have to make sure our subsequent analyses are not affected by the issue though. 4) Measuring the noise------Bias images usually have some repeatable structure to them, so a useful trick is to use the *difference* of two bias frames taken close in time. Let's measure the standard deviation for the differences between the biases, doing so separately for each amplifier. This isn't the final read noise value yet, because it's still in ADU and not in electrons. We will store the results in a dictionary for later use. We load the second image:<jupyter_code>bias2_file = fits.open("00258334672-S10-det003.fits")
# Question
measured_stddevs = {}

bias1_decs = bias1_file[1:16]
for i, hdu in enumerate(bias2_file[1:16]):
    hdu_difference = hdu.data - bias1_decs[i].data
    
    stddev =  np.std(hdu_difference)

    measured_stddevs[i] = stddev
measured_stddevs<jupyter_output><empty_output><jupyter_text>5) Measuring the gain---We have just one more step before we can report the read noise. We need to measure the gains so we can convert the noise in ADU into electrons. To do that, we're going to use the trick we saw at the start of this notebook. We need to add two things though: we want to use pairs of images, to cancel out any fixed spatial patterns, and we need images with significant counts in them so that we're not just measuring read noise. The formula we want to implement is thus:$$ \frac{1}{\textrm{gain}} =  \langle \frac{(I_1 - I_2)^2}{I_1 + I_2} \rangle $$where $I_1$ and $I_2$ are the pixel values from each image, and the $\langle$ $\rangle$ brackets denote taking the mean of this ratio over all pixels.We have some flat field images from those same sensors that we can use:<jupyter_code>flat1_file = fits.open("00258342968-S10-det003.fits")
flat2_file = fits.open("00258343136-S10-det003.fits")<jupyter_output><empty_output><jupyter_text>Since each amplifier can have a slightly different gain, we want to run this per-HDU and output a table of values. Since we're looping over the HDUs, we can also print the finished read noise values at the same time. Note that those have a factor of $\sqrt{2}$ because we took the difference of two bias frames, so the noise is greater than a single image.<jupyter_code># Question
import sys
for hdu in range(1,16):
    flat1_data = flat1_file[hdu].data
    flat2_data = flat2_file[hdu].data
    
    debiased_flat1 = flat1_data - bias1_file[hdu].data
    debiased_flat2 = flat2_data - bias1_file[hdu].data
    squared_noise = (debiased_flat1 - debiased_flat2)**2
    summed_intensity = debiased_flat1 + debiased_flat2
    # Some pixels with low counts are likely artifacts and can skew the measurement.
    # It helps to only keep pixels that have significant flux; you can experiment with this cutoff
    ok_values = (summed_intensity > 150000)# complete
    # Remember that as we defined gain above, the formula returns 1/gain.
    reciprocal_gain = np.mean(squared_noise[ok_values]/summed_intensity[ok_values])

    print(hdu, reciprocal_gain, 1/reciprocal_gain, measured_stddevs[hdu-1]/reciprocal_gain/np.sqrt(2))<jupyter_output>1 1.3988061749758016 0.7148953285234814 4.78107872980674
2 1.3910039804044478 0.7189052037861461 4.813561304524701
3 1.381810649623231 0.7236881553001949 4.789324376693804
4 1.3576246393258946 0.7365806210592443 4.901565269877944
5 1.3726965726319666 0.7284931134363004 6.305726481503794
6 1.3619033468978834 0.7342664971620639 5.119340105536518
7 1.3557778825803555 0.7375839456067621 5.039591292713493
8 1.3589476948805257 0.7358634947961825 5.12379636481906
9 1.4007110978570168 0.7139230934415564 5.572846959208616
10 1.3935363505960843 0.7175987907113084 5.46685053679957
11 1.3829211313725187 0.7231070357624242 4.965580436323176
12 1.3897076872652003 0.7195757850112319 5.313490201803283
13 1.3839613148861132 0.7225635494604056 5.0566653112161575
14 1.3831262707896557 0.7229997875964564 5.579632396410896
15 1.376421609723124 0.7265215780803937 5.333729730786731
{'fraction_non_alphanumeric': 0.06897155361050328, 'fraction_numerical': 0.08437636761487965, 'mean_word_length': 4.880597014925373, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 42, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '66528', 'license': 'MIT', 'n_tokens_mistral': 4097, 'n_tokens_neox': 3358, 'n_words': 1486, 'path': 'Session11/Day1/InvestigatingDetectors-Copy1.ipynb', 'repo_name': 'yupinghuang/LSSTC-DSFP-Sessions', 'size': '11427'}