Exercise 3: electron vs gamma separation¶
1. Introduction¶
Fig. 14 Difference between an electron vs photon is at the start of the electromagnetic shower, where the photon has a gap. From Wouter Van De Pontseele, ICHEP 2020.¶
Electrons are visible in a LArTPC detector because of the electromagnetic showers that they trigger.
Photons, on the other hand, are neutral (no charge) and thus remain invisible to the LArTPC eyes until they convert into electrons (pair production) or Compton scatter. In both cases, the visible outcome will be an electromagnetic shower.
How can we differentiate the two, then? The answer is in the very beginning of the EM shower. For an electron, this shower will be topologically connected to the interaction vertex where the electron was produced. For a photon, there will be a gap (equal to the photon travel path) until the EM shower start (when the photon becomes indirectly visible through pair production or Compton scatter). That seems simple enough, right? Wrong, of course.
Energetic photons could interact at a distance short enough from the interaction vertex, that we would not be able to see the gap. Or, the hadronic activity might be invisible, because it includes neutral particles or because the particles are too low energy to be seen. In that case the interaction vertex might be hard to identify, and the notion of a gap goes away too. For such cases, fortunately, there is another way to tell electrons from gamma showers. Another major difference is in the energy loss rate at the start of the EM shower. An electron would leave ionization corresponding to a single ionizing particle, whereas a pair of electron + positron coming from a photon pair production would add up to two ionizing particle. Thus, we expect the dE/dx at the beginning of the shower to be roughly twice larger in the case of a gamma-induced shower compared to an electron-induced shower.
Fig. 15 Example from MicroBooNE. Left is the shower \(dE/dx\), right is the gap between the vertex and shower start. From Wouter Van De Pontseele, ICHEP 2020.¶
Why do we care? The difference becomes significant if, for example, you are looking for electron neutrinos. One of the key signatures you would be looking for are electrons.
In this exercise, we will focus on finding the start of EM showers and computing the reconstructed dQ/dx in these segments. Optionally, you could compare that to the result of using automatic PID as predicted by the chain.
2. Setup¶
a. Software and data directory¶
import os, sys
SOFTWARE_DIR = '%s/lartpc_mlreco3d' % os.environ.get('HOME')
DATA_DIR = os.environ.get('DATA_DIR')
# Set software directory
sys.path.append(SOFTWARE_DIR)
b. Numpy, Matplotlib, and Plotly for Visualization and data handling.¶
import numpy as np
import matplotlib.pyplot as plt
import seaborn
seaborn.set(rc={
'figure.figsize':(15, 10),
})
seaborn.set_context('talk')
import plotly
import plotly.graph_objs as go
from plotly.subplots import make_subplots
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=False)
c. MLRECO specific imports for model loading and configuration setup¶
from mlreco.main_funcs import process_config, prepare
import warnings, yaml
warnings.filterwarnings('ignore')
cfg = yaml.load(open('%s/inference.cfg' % DATA_DIR, 'r').read().replace('DATA_DIR', DATA_DIR),Loader=yaml.Loader)
process_config(cfg, verbose=False)
/usr/local/lib/python3.8/dist-packages/MinkowskiEngine/__init__.py:36: UserWarning:
The environment variable `OMP_NUM_THREADS` not set. MinkowskiEngine will automatically set `OMP_NUM_THREADS=16`. If you want to set `OMP_NUM_THREADS` manually, please export it on the command line before running a python script. e.g. `export OMP_NUM_THREADS=12; python your_program.py`. It is recommended to set it below 24.
Config processed at: Linux ampt017 3.10.0-1160.42.2.el7.x86_64 #1 SMP Tue Sep 7 14:49:57 UTC 2021 x86_64 x86_64 x86_64 GNU/Linux
$CUDA_VISIBLE_DEVICES="0"
d. Initialize and load weights to model using Trainer.¶
# prepare function configures necessary "handlers"
hs = prepare(cfg)
dataset = hs.data_io_iter
Welcome to JupyROOT 6.22/09
Loading file: /sdf/home/l/ldomine/lartpc_mlreco3d_tutorials/book/data/mpvmpr_062022_test_small.root
Loading tree sparse3d_reco
Warning in <TClass::Init>: no dictionary for class larcv::EventNeutrino is available
Warning in <TClass::Init>: no dictionary for class larcv::NeutrinoSet is available
Warning in <TClass::Init>: no dictionary for class larcv::Neutrino is available
Loading tree sparse3d_reco_chi2
Loading tree sparse3d_reco_hit_charge0
Loading tree sparse3d_reco_hit_charge1
Loading tree sparse3d_reco_hit_charge2
Loading tree sparse3d_reco_hit_key0
Loading tree sparse3d_reco_hit_key1
Loading tree sparse3d_reco_hit_key2
Loading tree sparse3d_pcluster_semantics_ghost
Loading tree cluster3d_pcluster
Loading tree particle_pcluster
Loading tree particle_mpv
Loading tree sparse3d_pcluster_semantics
Loading tree sparse3d_pcluster
Loading tree particle_corrected
Found 101 events in file(s)
Shower GNN: True
Track GNN: True
Particle GNN: False
Interaction GNN: True
Kinematics GNN: False
Cosmic GNN: False
Since one of the GNNs are turned on, process_fragments is turned ON.
Fragment processing is turned ON. When training CNN models from
scratch, we recommend turning fragment processing OFF as without
reliable segmentation and/or cnn clustering outputs this could take
prohibitively large training iterations.
Shower GNN: True
Track GNN: True
Particle GNN: False
Interaction GNN: True
Kinematics GNN: False
Cosmic GNN: False
Since one of the GNNs are turned on, process_fragments is turned ON.
Fragment processing is turned ON. When training CNN models from
scratch, we recommend turning fragment processing OFF as without
reliable segmentation and/or cnn clustering outputs this could take
prohibitively large training iterations.
Freezing 82 weights for a sub-module ppn
Freezing 141 weights for a sub-module uresnet_lonely
Freezing 141 weights for a sub-module uresnet_deghost
Freezing 146 weights for a sub-module graph_spice
Freezing 120 weights for a sub-module grappa_track
Freezing 120 weights for a sub-module grappa_shower
Restoring weights for from /sdf/home/l/ldomine/lartpc_mlreco3d_tutorials/book/data/weights_full_mpvmpr_062022.ckpt...
Done.
Let’s load one iteration worth of data into our notebook:
data, result = hs.trainer.forward(dataset)
Deghosting Accuracy: 0.9830
Segmentation Accuracy: 0.9900
PPN Accuracy: 0.8843
Clustering Accuracy: 0.2691
Clustering Edge Accuracy: 0.1252
Shower fragment clustering accuracy: 0.9581
Shower primary prediction accuracy: 0.9434
Track fragment clustering accuracy: 0.9937
Interaction grouping accuracy: 0.9763
Particle ID accuracy: 0.8409
Primary particle score accuracy: 0.9755
e. Setup Evaluator¶
from analysis.classes.ui import FullChainEvaluator
# Only run this cell once!
evaluator = FullChainEvaluator(data, result, cfg, deghosting=True)
print(evaluator)
FullChainEvaluator(num_images=10)
entry = 4 # Batch ID for current sample
print("Batch ID = ", evaluator.index[entry])
Batch ID = 4
3. Identifying Shower Primaries¶
Step 1: Get shower primary fragments¶
By using the only_primaries=True option, we can select out primary particles in this image. We will also load true_particles for comparison.
particles = evaluator.get_particles(entry, only_primaries=True)
true_particles = evaluator.get_true_particles(entry, only_primaries=True)
from pprint import pprint
pprint(particles)
[Particle( Image ID=0 | Particle ID=6 | Semantic_type: Track | PID: Proton | Primary: 1 | Score = 99.99% | Interaction ID: 6 | Size: 528 ),
Particle( Image ID=0 | Particle ID=7 | Semantic_type: Track | PID: Muon | Primary: 1 | Score = 99.95% | Interaction ID: 6 | Size: 2006 ),
Particle( Image ID=0 | Particle ID=8 | Semantic_type: Track | PID: Proton | Primary: 1 | Score = 99.98% | Interaction ID: 6 | Size: 440 )]
Alternatively, as you may have noticed, the primariness information is also stored in the Particle instance as an attribute with name is_primary. If you prefer to view the full image and then select out primaries manually:
particles = evaluator.get_particles(entry, only_primaries=False)
true_particles = evaluator.get_true_particles(entry, only_primaries=False)
from pprint import pprint
pprint(particles)
[Particle( Image ID=0 | Particle ID=0 | Semantic_type: Shower Fragment | PID: Photon | Primary: 0 | Score = 96.25% | Interaction ID: 0 | Size: 220 ),
Particle( Image ID=0 | Particle ID=4 | Semantic_type: Shower Fragment | PID: Photon | Primary: 0 | Score = 92.49% | Interaction ID: 4 | Size: 235 ),
Particle( Image ID=0 | Particle ID=5 | Semantic_type: Track | PID: Photon | Primary: 0 | Score = 80.81% | Interaction ID: 5 | Size: 1664 ),
Particle( Image ID=0 | Particle ID=6 | Semantic_type: Track | PID: Proton | Primary: 1 | Score = 99.99% | Interaction ID: 6 | Size: 528 ),
Particle( Image ID=0 | Particle ID=7 | Semantic_type: Track | PID: Muon | Primary: 1 | Score = 99.95% | Interaction ID: 6 | Size: 2006 ),
Particle( Image ID=0 | Particle ID=8 | Semantic_type: Track | PID: Proton | Primary: 1 | Score = 99.98% | Interaction ID: 6 | Size: 440 ),
Particle( Image ID=0 | Particle ID=9 | Semantic_type: Track | PID: Photon | Primary: 0 | Score = 74.19% | Interaction ID: 9 | Size: 596 ),
Particle( Image ID=0 | Particle ID=10 | Semantic_type: Track | PID: Photon | Primary: 0 | Score = 74.67% | Interaction ID: 10 | Size: 1245 ),
Particle( Image ID=0 | Particle ID=11 | Semantic_type: Track | PID: Photon | Primary: 0 | Score = 67.31% | Interaction ID: 11 | Size: 1250 ),
Particle( Image ID=0 | Particle ID=12 | Semantic_type: Track | PID: Pion | Primary: 0 | Score = 49.21% | Interaction ID: 12 | Size: 90 ),
Particle( Image ID=0 | Particle ID=13 | Semantic_type: Track | PID: Muon | Primary: 0 | Score = 78.96% | Interaction ID: 13 | Size: 3354 ),
Particle( Image ID=0 | Particle ID=16 | Semantic_type: Delta Ray | PID: Electron | Primary: 0 | Score = 64.23% | Interaction ID: 10 | Size: 28 ),
Particle( Image ID=0 | Particle ID=17 | Semantic_type: Delta Ray | PID: Electron | Primary: 0 | Score = 46.07% | Interaction ID: 13 | Size: 21 )]
Let’s quickly plot the particles and visualize which ones are predicted as primaries. Here is one way to do it with the trace_particles function:
from mlreco.visualization.plotly_layouts import white_layout, trace_particles, trace_interactions
traces = trace_particles(particles, color='is_primary', colorscale='rdylgn') # is_primary for coloring with respect to primary label
traces_true = trace_particles(true_particles, color='is_primary', colorscale='rdylgn')
fig = make_subplots(rows=1, cols=2,
specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],
horizontal_spacing=0.05, vertical_spacing=0.04)
fig.add_traces(traces, rows=[1] * len(traces), cols=[1] * len(traces))
fig.add_traces(traces_true, rows=[1] * len(traces_true), cols=[2] * len(traces_true))
fig.layout = white_layout()
fig.update_layout(showlegend=False,
legend=dict(xanchor="left"),
autosize=True,
height=600,
width=1500,
margin=dict(r=20, l=20, b=20, t=20))
iplot(fig)
The green voxels are predicted primary particles, while red indicates non-primary.
It is often easier to further break down the shower into different fragments and locate which of the shower fragments actually correspond to a predicted primary.
fragments = evaluator.get_fragments(entry)
traces = trace_particles(fragments, color='is_primary', colorscale='rdylgn') # is_primary for coloring with respect to primary label
traces_right = trace_particles(fragments, color='id', colorscale='rainbow') # This time, we'll plot the predicted particle
fig = make_subplots(rows=1, cols=2,
specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],
horizontal_spacing=0.05, vertical_spacing=0.04)
fig.add_traces(traces, rows=[1] * len(traces), cols=[1] * len(traces))
fig.add_traces(traces_right, rows=[1] * len(traces_right), cols=[2] * len(traces_right))
fig.layout = white_layout()
fig.update_layout(showlegend=False,
legend=dict(xanchor="left"),
autosize=True,
height=600,
width=1500,
margin=dict(r=20, l=20, b=20, t=20))
iplot(fig)
# TODO: Plot true fragment labels
Step 2: Identify the startpoint of the shower primary¶
During initialization of the Particle instance, PPN predictions are assigned to each particle if the distance between then is less than a predetermined threshold (attaching_threshold). PPN predictions that are matched to particles in this way are then stored in each Particle instance as attributes (ppn_candidates)
print("Minimum voxel distance required to assign ppn prediction to particle fragment = ", evaluator.attaching_threshold)
Minimum voxel distance required to assign ppn prediction to particle fragment = 2
fragments = evaluator.get_fragments(entry, only_primaries=False)
The first three columns are the \((x,y,z)\) coordinates of the PPN points. The fourth column is the PPN prediction score, and the last column indicates the predicted semantic type of the point.
We first visualize whether the predicted ppn candidates accurately locate the shower fragment start:
traces = trace_particles(fragments, color='id', size=1, scatter_ppn=True, highlight_primaries=True) # Set scatter_ppn=True for plotting PPN information
traces_true = trace_particles(true_particles, color='id', size=1)
fig = make_subplots(rows=1, cols=2,
specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],
horizontal_spacing=0.05, vertical_spacing=0.04)
fig.add_traces(traces, rows=[1] * len(traces), cols=[1] * len(traces))
fig.add_traces(traces_true, rows=[1] * len(traces_true), cols=[2] * len(traces_true))
fig.layout = white_layout()
fig.update_layout(showlegend=False,
legend=dict(xanchor="left"),
autosize=True,
height=600,
width=1500,
margin=dict(r=20, l=20, b=20, t=20))
iplot(fig)