Single Photon Counting with FastCCD and CSXTOOLS

In [1]:

    

import numpy as np
import tifffile

%matplotlib inline
from matplotlib import pyplot as plt
plt.style.use('sbwnotebook')

from csxtools.fastccd import phocount
from csxtools.utils import make_panel_plot
from csxtools.ipynb import image_stack_to_movie, show_image_stack
import timeit


    

In [2]:

    

def get_photons(data):
    results = []
    histograms = []

    for n in range(1, 9):
        c = phocount.photon_count(data, (75,250), (0, 500), n, True)
        h = np.histogram(c[0], bins=100, range=(50, 350))
        results.append(c[0].ravel())
        histograms.append(h[0]) 
    
    hist_x = h[1][:-1]
    results = np.array(results)
    histograms = np.array(histograms)
    
    return (results, hist_x, histograms)


def plot_results(results, hist_x, histograms):
    color=iter(plt.cm.rainbow(np.linspace(0,1,9)))

    plt.subplot(121)
    for x,h in enumerate(histograms):
        plt.plot(hist_x, h, 'o-', c=next(color), label='{}'.format(x+1))
    plt.legend()

    m = np.nanmean(results, axis=1)
    s = np.nanstd(results, axis=1)
    
    ax1 = plt.subplot(122)
    ax1.plot(np.arange(1, 9), m,'bo')
    ax1.set_xlabel('ESUM Value')
    ax1.set_ylabel('Mean Integrated ADU')

    ax2 = ax1.twinx()
    ax2.plot(np.arange(1, 9), s,'rs')
    ax2.set_ylabel('Standard Deviation')


    

In [3]:

    

data = []
for x in range(1,661):
    tf = tifffile.TiffFile('/GPFS/xf23id/xf23id1/LBNL_Test_Data/20140410/Fe55-Exp100us-HC5837/bnl2_20140410_normalizedTiffs50count/image{}.tiff'.format(x))
    data.append(tf.pages[0].asarray())
data = np.array(data)
# Only take half the CCD
data = data[:,240:,:]
print(data.shape)


    

    




(660, 238, 960)

In [4]:

    

plt.figure()
plt.imshow(data[2], vmin=10, vmax=300, interpolation='none')
plt.xlabel('Columns')
plt.ylabel('Rows')


    

    
Out[4]:





<matplotlib.text.Text at 0x7f3fe4d93208>






    

In [5]:

    

t = timeit.timeit('phocount.photon_count(data, (150,250), (0, 300), 3)', globals=globals(), number=25)
print("Photon counting array of {} takes {:.3} seconds".format(data.shape, t / 25))


    

    




Photon counting array of (660, 238, 960) takes 0.335 seconds

In [6]:

    

h = np.histogram(data.ravel(), bins=300, range=(0, 300))
fig, ax = plt.subplots()
ax.semilogy(h[1][:-1], h[0], 'r-')


    

    
Out[6]:





[<matplotlib.lines.Line2D at 0x7f3fe78805f8>]






    

In [7]:

    

results, hist_x, histograms = get_photons(data)


    

In [8]:

    

plot_results(results, hist_x, histograms)


    

In [9]:

    

c = phocount.photon_count(data.astype(np.float32), (75,250), (170, 230), 3, True)
inti = c[0]
#inti[:,:,0::10] = np.nan
#inti[:,:,9::10] = np.nan
ff = np.nanmean(inti, axis=(0,1))
ffstd = np.nanstd(inti, axis=(0,1))
ffn = np.nansum(inti, axis=(0,1))


    

    




/GPFS/xf23id/users/swilkins/miniconda3/envs/analysis/lib/python3.5/site-packages/numpy/lib/nanfunctions.py:675: RuntimeWarning: Mean of empty slice
  warnings.warn("Mean of empty slice", RuntimeWarning)
/GPFS/xf23id/users/swilkins/miniconda3/envs/analysis/lib/python3.5/site-packages/numpy/lib/nanfunctions.py:1136: RuntimeWarning: Degrees of freedom <= 0 for slice.
  warnings.warn("Degrees of freedom <= 0 for slice.", RuntimeWarning)

In [10]:

    

fig = plt.figure(figsize=(18,18))
ax = make_panel_plot(n=3, fig=fig)

ax[2].plot(ff, 'r*-')
ax[1].plot(ffn, 'ks-')
ax[0].plot(ffstd, 'b+-')

ax[2].set_ylabel(r'$\leftangle I \rightangle$ [ADU]')
ax[2].set_title('Results by column ')
ax[1].set_ylabel('n [photons]')
ax[0].set_ylabel('$\sigma (I)$')
ax[0].set_xlabel('Column')

for a in ax:
    a.set_xlim((0,960))


    

In [11]:

    

import csxtools
csxtools.__version__


    

    
Out[11]:





'0.1.7+2.g67cddca'

In [ ]: