LightMap Hit Pattern Examples¶
In this notebook, we go over a method for modeling a hit pattern — the photon collection efficiency as a function of SiPM, given a photon generation position inside a detector (in this case, a TPC).
All caveats regarding the total light treatment in this project apply equally here.
Jump to:
Imports¶
import numpy as np
import matplotlib.pyplot as plt
import histlite as hl
import LightMap
%matplotlib inline
# plotting aesthetics
plt.rc('figure', dpi=150, figsize=(4,3), facecolor='w')
plt.rc('savefig', dpi=150, facecolor='w')
plt.rc('lines', linewidth=1.5)
Loading an existing LightMap¶
We'll start by loading up a LightMap previously computed in the total-light tutorial. This is useful for two reasons:
- We can ensure that our hitpattern is defined for the same geometry.
- We can optionally use the total light map in conjunction with the hitpattern.
The models are stored under:
import os
home_dir = os.getenv('HOME')
model_dir = '{}/lightmap/total/tutorial'.format(home_dir)
And we recover the LightMap itself by, e.g.:
lm = LightMap.load_model(model_dir, 'LightMapSplineXYZ')
lm
The geometry of that lightmap is:
tpc = lm.tpc
tpc
HitPattern — low statistics¶
We'll work with the same MC dataset as in the total light tutorial:
mc_dir = '/p/lustre2/nexouser/optical_sample/g4'
filenames = [
'{}/PhotonBomb_gps_ActiveOptical_seed{}.nEXOevents.root'.format(mc_dir, i)
for i in range(1, 1001)]
High statistics MC is required to parameterize the hit pattern with reasonable accuracy. And now, we need not only the per-"event" summary information, but the details of where each collected photon landed. For this reason, the fitting is done in terms of filenames rather than data pre-loaded into RAM.
The HitPattern is initially constructed like so:
# assumes 384x120 SiPM array
hp1 = LightMap.hitpattern.HitPatternHist(tpc)
We'll start by fitting with a single file:
hp1.fit(filenames[:1])
Over 15 seconds per file is not the best, so for this tutorial we'll save immediately after every fit:
LightMap.save_model(model_dir, 'HitPatternHist__0001_files', hp1)
We reload the lightmap to make sure everything is working:
hp1a = LightMap.load_model(model_dir, 'HitPatternHist__0001_files')
A quick check that the underlying histograms match:
np.all(hp1.H.values == hp1a.H.values)
HitPattern — decent statistics¶
We won't be able to do much with one file. How long are we willing to wait? Let's say 60 mins. How many files will that get us?
60 * 60 / 16. # mins * sec-per-min / sec-per-file
Off we go, then, processing 200 files:
hp = LightMap.hitpattern.HitPatternHist(tpc)
hp.fit(filenames[:200])
Again, we save right away:
LightMap.save_model(model_dir, 'HitPatternHist__0200_files', hp)
#hp = LightMap.load_model(model_dir, 'HitPatternHist__0200_files')
HitPattern Visualization¶
Before we try anything else, we should make sure the hitpattern makes sense visually.
We'll define some helper variables and a function:
figsize = (1.5*3.9, 1.5*1.24)
figsize6 = (2*figsize[0], 3 * figsize[1])
def hp_plot(ax, h):
out = hl.plot2d(ax, h / h.values.max(), cmap='Greys_r', vmin=0)
ax.set_xlabel(r'SiPM $\theta$ bin')
ax.set_ylabel(r'SiPM $z$ bin')
return out
And now, here is an example of continuous evaluation in $r_\text{gen}$:
q = hp.centers_theta[17]
z = hp.centers_z[4]
rs = np.linspace(hp.centers_r[6], hp.centers_r[7], 6)
fig, axs = plt.subplots(3, 2, figsize=figsize6)
axs = np.ravel(axs)
title_fmt = (r'$(r_\mathsf{{gen}},\theta_\mathsf{{gen}},z_\mathsf{{gen}})'
r'\ =\ (\mathbf{{{:.0f}}},\ {:.0f}^\circ,\ {:.0f})$')
for (i,r) in enumerate(rs):
ax = axs[i]
h = hp(r, q, z, cyl=True)
hp_plot(ax, h)
ax.set_title(title_fmt.format(r, np.degrees(q), z))
plt.tight_layout()
and in $\theta_\text{gen}$:
r = hp.centers_r[6]
z = hp.centers_z[4]
qs = np.linspace(hp.centers_theta[16], hp.centers_theta[17], 6)
fig, axs = plt.subplots(3, 2, figsize=figsize6)
axs = np.ravel(axs)
title_fmt = (r'$(r_\mathsf{{gen}},\theta_\mathsf{{gen}},z_\mathsf{{gen}})'
r'\ =\ ({:.0f},\ \mathbf{{{:.0f}^\circ}},\ {:.0f})$')
for (i,q) in enumerate(qs):
ax = axs[i]
h = hp(r, q, z, cyl=True)
hp_plot(ax, h)
ax.set_title(title_fmt.format(r, np.degrees(q), z))
plt.tight_layout()
and finally, in $z_\text{gen}$:
r = hp.centers_r[6]
q = hp.centers_theta[17]
zs = np.linspace(hp.centers_z[4], hp.centers_z[5], 6)
fig, axs = plt.subplots(3, 2, figsize=figsize6)
axs = np.ravel(axs)
title_fmt = (r'$(r_\mathsf{{gen}},\theta_\mathsf{{gen}},z_\mathsf{{gen}})'
r'\ =\ ({:.0f},\ {:.0f}^\circ,\ \mathbf{{{:.0f}}})$')
for (i,z) in enumerate(zs):
ax = axs[i]
h = hp(r, q, z, cyl=True)
hp_plot(ax, h)
ax.set_title(title_fmt.format(r, np.degrees(q), z))
plt.tight_layout()
HitPattern value access¶
Ok, so what are these HitPatternHist objects? If called like a function, they behave similarly to LightMaps, but instead of accepting arrays of generation positions and returning corresponding floats, they accept one generation position at a time and return the entire hl.Hist of detection probabilities.
Note: the below plot revealed that there is either a new, or a pre-existing but previously undetected, bug in the SiPMID $\leftrightarrow$ (col,row) conversion; see the overly bright column on the right side. We're going to ignore it for now, but obviously this needs to be fixed — and in general, it would be good to get some clarity on how exactly LightMap.hitpattern.sipm_id_to_col_row() and its TODO sibling, LightMap.hitpattern.col_row_to_sipm_id() should be expressed.
fig, ax = plt.subplots(figsize=figsize)
hp_plot(ax, hp(142, -142, -400))
To be clear,
hp(142, -142, -400)
Either Cartesian (default) or cylindrical coordinates can be used:
hp(400, np.radians(45), -400) # silliness: why y in rad?
hp(400, np.radians(45), -400, cyl=True) # ok
hp(400/np.sqrt(2), 400/np.sqrt(2), -400) # ok
We can capture the total collection probability (summing over SiPMs):
hp(400, np.radians(45), -400, cyl=True).values.sum()
We can also specify a total lightmap from which to obtain the overall normalization:
hp(400, np.radians(45), -400, cyl=True, total=lm).values.sum()
Of course, this will match the collection efficiency reported by lm:
lm(400, np.radians(45), -400, cyl=True)
Random sampling¶
We can sample from the hitpattern:
col, row = hp.sample(30, 400, np.radians(-45), -400, cyl=True, total=lm, seed=1)
col, row
And we can convert the sampled values to a histogram:
hp.sample_to_image((col, row))
We can also obtain a randomly sampled histogram directly:
hp.sample_image(30, 400, np.radians(-45), -400, cyl=True, total=lm, seed=1)
The latter is useful if, for example, you want to go straight to plotting:
fig, ax = plt.subplots(figsize=figsize)
h = hp.sample_image(10000, 400, np.radians(-45), -400, cyl=True, total=lm, seed=2)
hp_plot(ax, h)
Finally, note that you can either specify a quantum efficiency, reducing the number of observed photons:
hp.sample(30, 400, np.radians(-45), -400, cyl=True, total=lm, seed=1, qe=.15)
Or you can sample a specific number of observed rather than initial photons:
hp.sample(30, 400, np.radians(-45), -400, cyl=True, total=lm, seed=1, all=True)
Towards maximum likelihood¶
This material probably deserves to be formalized as its own tool, as part of LightMap or elsewhere.
Let's arbitrarily suppose that we have an event that initially generates 60k photons, fairly close to the skin:
x, y, z = 500, 100, -400
im = hp.sample_image(60000, x, y, z, total=lm, qe=.15, seed=1)
im.values.sum()
The event data look like so:
fig, ax = plt.subplots(figsize=figsize)
hp_plot(ax, im)
Compare with the expectation for an event at this position:
fig, ax = plt.subplots(figsize=figsize)
hp_plot(ax, hp(x, y, z))
Let's write down a simple binned-Poissonian log likelihood for the position:
def minus_llh(im, x, y, z):
# expectation must be nonzero, and we're still working with low stats
# so, smooth the expectation
if not hp.tpc.mask(x, y, z):
return 1e5
expected = hp(x, y, z)#.gaussian_filter(2)
obs = im.values.ravel()
pred = im.values.sum() * expected.values.ravel()
# logl = -pred + obs*log(pred)
mask = (pred > 0)
return -np.sum(-pred[mask] + obs[mask] * np.log(pred[mask]))
vminus_llh = np.vectorize(minus_llh)
This looks like it should be "enough information" to enable a position reconstruction; let's try. First, we test the TS at the true generated position, and compare with spots in the middle and on the other side of the TPC:
minus_llh(im, x, y, z), minus_llh(im, 0, 0, z), minus_llh(im, -x, -y, z)
Next, let's scan the TS on a grid:
X, Y, Z, V = tpc.get_grid_cartesian()
X.size
%time grid_ts = vminus_llh(im, X, Y, Z)
The best fit position so far is:
ibest = np.argmin(grid_ts)
ibest
X[ibest], Y[ibest], Z[ibest]
We take this as a seed for a likelihood fit:
from scipy import optimize
%time optimize.fmin(lambda xyz: minus_llh(im, *xyz), x0=(X[ibest], Y[ibest], Z[ibest]))
Sooo... this basically works, but it's pretty slow.
What's next?¶
Resolving the hit pattern as a function of event generation position is a relatively high-dimensional (5D) problem. Here we have only a preliminary solution. Things that can be done with this code:
- Experiment with light-only reconstructions
- Experiment with light-charge matching
- If we can sort out the
sipm_id_to_col_row()andcol_row_to_sipm_id()conversions, then the above.samplecan be used to do per-SiPM "fastlight" MC.
Here are some open issues:
- Histlite was written more for convenience than for speed; many of its sanity checks are unnecessary and slow for performance sensitive code. Once the parameterization is built, evaluation for given generated positions would be faster with dedicated code (unlike for total light, I think C++ might help a lot here).
- This feels like an excellent use case for deep learning methods. In IceCube we have a "cascade generator" that can substitute for light propagation spline tables by effectively simulating (coarsely) each event hypothesis, rather than looking up light expectations. Maybe we can do the same for nEXO.