Experimental data analysis on foil open area

Brian Larsen, ISR-1

Data provided by Phil Fernandes, ISR-1 2016-9-14

The setup is a foil in its holder mounted to a foil holder meant to bock incident ions. The foil has a ~0.6mm hole in it to provide a baseline. The goal is to use the relative intensity of the witness hole to determine the intensity of holes in the foil.

A quick summary:

Foil is placed 0.66” from front of MCP surface
Beam is rastered to cover full foil and “witness” aperture
Beam is 1.0 keV Ar+, slightly underfocused
Accumulate data for set period of time (either 60s or 180s, identified in spreadsheet)
Total_cts is the # of counts through the foil and the witness aperture
Witness_cts is the # of counts in the witness aperture only
Foil_cts = total_cts – witness_cts
Open area OA = (foil_cts/witness_cts) * (witness_area/foil_area)



In [1]:

    
import itertools
from pprint import pprint
from operator import getitem

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import numpy as np
import spacepy.plot as spp
import pymc as mc
import tqdm

from MCA_file_viewer_v001 import GetMCAfile









    



/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/matplotlib/__init__.py:913: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))






    



---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-1-6e1255deb835> in <module>()
      7 import numpy as np
      8 import spacepy.plot as spp
----> 9 import pymc as mc
     10 import tqdm
     11 

ModuleNotFoundError: No module named 'pymc'



In [ ]:

    
def plot_box(x, y, c='r', lw=0.6, ax=None):
    if ax is None:
        plt.plot((xind[0], xind[0]), (yind[0], yind[1]), lw=lw, c=c)
        plt.plot((xind[1], xind[1]), (yind[0], yind[1]), lw=lw, c=c)
        plt.plot((xind[0], xind[1]), (yind[0], yind[0]), lw=lw, c=c)
        plt.plot((xind[0], xind[1]), (yind[1], yind[1]), lw=lw, c=c)
    else:
        ax.plot((xind[0], xind[0]), (yind[0], yind[1]), lw=lw, c=c)
        ax.plot((xind[1], xind[1]), (yind[0], yind[1]), lw=lw, c=c)
        ax.plot((xind[0], xind[1]), (yind[0], yind[0]), lw=lw, c=c)
        ax.plot((xind[0], xind[1]), (yind[1], yind[1]), lw=lw, c=c)



In [ ]:

    
ZZ, XX, YY = GetMCAfile('16090203.mca')
# It is believed as of 2016-09-19 that the MCA records 2 counts for each count. 
#    This means all data are even and all the data can be divided by 2 to give the
#    right number of counts. Per emails Larsen-Fernandes 2016-09-17
#    These data are integres and care muct be taken to assure that /2 does not
#    lead to number that are not representable in float
ZZ = ZZ.astype(float)
ZZ /= 2
XX = XX.astype(np.uint16) # as they all should be integers anyway



In [ ]:

    
xind = (986, 1003)
yind = (492, 506)

fig = plt.figure(figsize=(20,8))
ax1 = fig.add_subplot(131)
ax2 = fig.add_subplot(132)
ax3 = fig.add_subplot(133)

pc = ax1.pcolormesh(XX, YY, ZZ, norm=LogNorm())
plt.colorbar(pc, ax=ax1)
plot_box(xind, yind, ax=ax1)

ax2.hist(ZZ.flatten(), 20)
ax2.set_yscale('log')

ax3.hist(ZZ.flatten(), 20, normed=True)
ax3.set_yscale('log')

Do some calculations to try and match Phil's analysis

Phil's data:

File name Witness cts Total cts Foil cts Open area

16090203 658 4570 3912 0.00102



In [ ]:

    
total_cnts = ZZ.sum()
print('Total counts:{0}  -- Phil got {1} -- remember /2'.format(total_cnts, 4570/2)) # remember we did a /2



In [ ]:

    
# Is the whitness hole at x=1000, y=500?
XX.shape, YY.shape, ZZ.shape



In [ ]:

    
print(ZZ[yind[0]:yind[1], xind[0]:xind[1]])
plt.figure()
plt.pcolormesh(XX[xind[0]:xind[1]], YY[yind[0]:yind[1]], ZZ[yind[0]:yind[1], xind[0]:xind[1]] , norm=LogNorm())
plt.colorbar()

witness_counts = ZZ[yind[0]:yind[1], xind[0]:xind[1]].sum()

print('Witness counts: {0}, Phil got {1}/2={2}'.format(witness_counts, 658, 658/2))
wit_pixels = 46
print('There {0} pixels in the witness peak'.format(wit_pixels))

total_counts = ZZ.sum()
print("There are a total of {0} counts".format(total_counts))

Can we get a noise estimate?

1) Try all pixels with a value where a neighbor does not. This assumes that real holes are large enough to have a point spread function and therefore cannot be in a single pixel.



In [ ]:

    
def neighbor_inds(x, y, xlim=(0,1023), ylim=(0,1023), center=False, mask=False):
    """
    given an x and y index return the 8 neighbor indices
    
    if center also return the center index
    if mask return a boolean mask over the whole 2d array
    """
    xi = np.clip([x + v for v in [-1, 0, 1]], xlim[0], xlim[1])
    yi = np.clip([y + v for v in [-1, 0, 1]], ylim[0], ylim[1])
    ans = [(i, j) for i, j in itertools.product(xi, yi)]
    if not center:
        ans.remove((x,y))
    if mask:
        out = np.zeros((np.diff(xlim)+1, np.diff(ylim)+1), dtype=np.bool)
        for c in ans:
            out[c] = True
    else:
        out = ans
    return np.asarray(out)

print(neighbor_inds(2,2))
print(neighbor_inds(2,2, mask=True))
print(ZZ[neighbor_inds(500, 992, mask=True)])



In [ ]:

    
def get_alone_pixels(dat):
    """
    loop over all the data and store the value of all lone pixels
    """
    ans = []
    for index, x in tqdm.tqdm_notebook(np.ndenumerate(dat)):
        if (np.sum([ZZ[i, j] for i, j in neighbor_inds(index[0], index[1])]) == 0) and x != 0:
            ans.append((index, x))
    return ans
# print((neighbor_inds(5, 4)))
alone = get_alone_pixels(ZZ)
pprint(alone)
# ZZ[neighbor_inds(5, 4)[0]].shape
# print((neighbor_inds(5, 4))[0])
# print(ZZ[(neighbor_inds(5, 4))[0]].shape)
# ZZ[4,3]



In [ ]:

    
ZZ[(965, 485)]



In [ ]:

    
print(neighbor_inds(4,3)[0])
print(ZZ[neighbor_inds(4,3)[0]])
print(ZZ[3,2])

ni = neighbor_inds(4,3)[0]
print(ZZ[ni[0], ni[1]])



In [ ]:

    
(ZZ % 2).any() # not all even any longer

Noise estimates

Not we assume that all lone counts are noise that can be considered random and uniform over the MCP. This then provides a number of counts per MCA pixel that we can use.



In [ ]:

    
n_noise = np.sum([v[1] for v in alone])
n_pixels = 1024*1024
noise_pixel = n_noise/n_pixels
print("There were a total of {0} random counts over {1} pixels, {2} cts/pixel".format(n_noise, n_pixels, noise_pixel))

Maybe we should consider just part of the MCP, lets get the min,max X and min,max Y where there are counts and just use that area. This will increase the cts/pixel.



In [ ]:

    
minx_tmp = ZZ.sum(axis=0)
minx_tmp.shape
print(minx_tmp)

miny_tmp = ZZ.sum(axis=1)
miny_tmp.shape
print(miny_tmp)

Looks to go all the way to all sides in X-Y.

Work to total open area calculations

Now we can model the total open area of the foil given the noise estimate per pixel and the pixels that are a part of the witness sample and the total area.

We model the observed background as Poisson with center at the real background:

$obsnbkg \sim Pois(nbkg)$

We model the observed witness sample, $obswit$, as Poisson with center of background per pixel times number of pixels in peak plus the number of real counts:

$obswit \sim Pois(nbkg/C + witc)$, $C = \frac{A_w}{A_t}$

This then leaves the number of counts in open areas of the system (excluding witness) as a Poisson with center of background per pixel times number of pixels in the system (less witness) plus the real number of counts.

$obsopen \sim Pois(nbkg/D + realc)$, $D=\frac{A_t - A_w}{A_t}$

Then then the open area is given by the ratio number of counts, $realc$, over an unknown area, $A_o$, as related to witness counts, $witc$, to the witness area, $A_w$, which is assumed perfect as as 0.6mm hole.

$\frac{A_o}{realc}=\frac{A_w}{witc} => A_o = \frac{A_w}{witc}realc $



In [ ]:

    
Aw = np.pi*(0.2/2)**2 # mm**2
Af = 182.75 # mm**2  this is the area of the foil
W_F_ratio = Aw/Af

print(Aw, Af, W_F_ratio)

C = wit_pixels/n_pixels
D = (n_pixels-wit_pixels)/n_pixels
print('C', C, 'D', D)


nbkg = mc.Uniform('nbkg', 1, n_noise*5) # just 1 to some large number
obsnbkg = mc.Poisson('obsnbkg', nbkg, observed=True, value=n_noise)

witc = mc.Uniform('witc', 0, witness_counts*5) # just 0 to some large number
obswit = mc.Poisson('obswit', nbkg*C + witc, observed=True, value=witness_counts)

realc = mc.Uniform('realc', 0, total_counts*5) # just 0 to some large number
obsopen = mc.Poisson('obsopen', nbkg*D + realc, observed=True, value=total_counts-witness_counts)

@mc.deterministic(plot=True)
def open_area(realc=realc, witc=witc):
    return realc*Aw/witc/Af

model = mc.MCMC([nbkg, obsnbkg, witc, obswit, realc, obsopen, open_area])



In [ ]:

    
model.sample(200000, burn=100, thin=30, burn_till_tuned=True)
mc.Matplot.plot(model)

# 1000, burn=100, thin=30    0.000985 +/- 0.000058
# 10000, burn=100, thin=30   0.000982 +/- 0.000061
# 100000, burn=100, thin=30  0.000984 +/- 0.000059
# 200000, burn=100, thin=30  0.000986 +/- 0.000059
# 1000000, burn=100, thin=30 0.000985 +/- 0.000059



In [ ]:

    
print("Foil 1 \n")

witc_mean = np.mean(witc.trace()[...])
witc_std = np.std(witc.trace()[...])

print("Found witness counts of {0} turn into {1} +/- {2}  ({3:.2f}%)\n".format(witness_counts, witc_mean, witc_std, witc_std/witc_mean*100))

realc_mean = np.mean(realc.trace()[...])
realc_std = np.std(realc.trace()[...])

print("Found non-witness counts of {0} turn into {1} +/- {2}  ({3:.2f}%)\n".format(total_counts-witness_counts, realc_mean, realc_std, realc_std/realc_mean*100))

nbkg_mean = np.mean(nbkg.trace()[...])
nbkg_std = np.std(nbkg.trace()[...])

print("Found noise counts of {0} turn into {1} +/- {2}  ({3:.2f}%)\n".format(0, nbkg_mean, nbkg_std, nbkg_std/nbkg_mean*100))

OA_median = np.median(open_area.trace()[...])
OA_mean = np.mean(open_area.trace()[...])
OA_std = np.std(open_area.trace()[...])
print("The open area fraction is {0:.6f} +/- {1:.6f}   ({2:.2f}%) at the 1 stddev level from 1 measurement\n".format(OA_mean, OA_std,OA_std/OA_mean*100 ))
print("Phil got {0} for 1 measurement\n".format(0.00139))
print("The ratio Brian/Phil is: {0:.6f} or {1:.6f}".format(OA_mean/0.00139, 0.00139/OA_mean))

Run again allowing some uncertainity on witness and foil areas



In [ ]:

    
_Aw = np.pi*(0.2/2)**2 # mm**2
_Af = 182.75 # mm**2  this is the area of the foil

Aw = mc.Normal('Aw', _Aw, (_Aw*0.2)**-2) # 20%
Af = mc.Normal('Af', _Af, (_Af*0.1)**-2)  # 10%

print(_Aw, _Af)

C = wit_pixels/n_pixels
D = (n_pixels-wit_pixels)/n_pixels
print('C', C, 'D', D)


nbkg = mc.Uniform('nbkg', 1, n_noise*5) # just 1 to some large number
obsnbkg = mc.Poisson('obsnbkg', nbkg, observed=True, value=n_noise)

witc = mc.Uniform('witc', 0, witness_counts*5) # just 0 to some large number
obswit = mc.Poisson('obswit', nbkg*C + witc, observed=True, value=witness_counts)

realc = mc.Uniform('realc', 0, total_counts*5) # just 0 to some large number
obsopen = mc.Poisson('obsopen', nbkg*D + realc, observed=True, value=total_counts-witness_counts)

@mc.deterministic(plot=True)
def open_area(realc=realc, witc=witc, Aw=Aw, Af=Af):
    return realc*Aw/witc/Af

model = mc.MCMC([nbkg, obsnbkg, witc, obswit, realc, obsopen, open_area, Af, Aw])



In [ ]:

    
model.sample(200000, burn=100, thin=30, burn_till_tuned=True)



In [ ]:

    
mc.Matplot.plot(nbkg)
mc.Matplot.plot(witc)
mc.Matplot.plot(realc)
# mc.Matplot.plot(open_area)
mc.Matplot.plot(Aw)

_ = spp.plt.hist(open_area.trace(), 20)



In [ ]:

    
print("Foil 1 \n")

witc_mean = np.mean(witc.trace()[...])
witc_std = np.std(witc.trace()[...])

print("Found witness counts of {0} turn into {1} +/- {2}  ({3:.2f}%)\n".format(witness_counts, witc_mean, witc_std, witc_std/witc_mean*100))

realc_mean = np.mean(realc.trace()[...])
realc_std = np.std(realc.trace()[...])

print("Found non-witness counts of {0} turn into {1} +/- {2}  ({3:.2f}%)\n".format(total_counts-witness_counts, realc_mean, realc_std, realc_std/realc_mean*100))

nbkg_mean = np.mean(nbkg.trace()[...])
nbkg_std = np.std(nbkg.trace()[...])

print("Found noise counts of {0} turn into {1} +/- {2}  ({3:.2f}%)\n".format(0, nbkg_mean, nbkg_std, nbkg_std/nbkg_mean*100))

OA_median = np.median(open_area.trace()[...])
OA_mean = np.mean(open_area.trace()[...])
OA_std = np.std(open_area.trace()[...])
print("The open area fraction is {0:.6f} +/- {1:.6f}   ({2:.2f}%) at the 1 stddev level from 1 measurement\n".format(OA_mean, OA_std,OA_std/OA_mean*100 ))
print("Phil got {0} for 1 measurement\n".format(0.00139))
print("The ratio Brian/Phil is: {0:.6f} or {1:.6f}".format(OA_mean/0.00139, 0.00139/OA_mean))



In [ ]:

    
mc.Matplot.plot(Aw)



In [ ]: