Exercise 2: Muon stopping power¶¶

The stopping power of a particle usually refers to the energy loss rate \(dE/dx\) when it passes through matter. When charged particles travel through our LArTPC detector, they interact with the argon and lose energy.

Note

Currently this tutorial only computes dQ/dx distribution. Need refresh to make a 2D plot of dQ/dx vs residual range.

MIP and Bragg peak¶¶

Minimally ionizing particles (MIP) are charged particles which lose energy when passing through matter at a rate close to minimal. Particles such as muons often have energy losses close to the MIP level and are treated in practice as MIP. The only exception is when the muon comes to a stop and experiences a Bragg peak.

Fig. 12 Example of muon Bragg peak. The muon is travelling from bottom left to top right. The color scale represents energy deposition. Red means more energy deposited. The sudden increase in deposited (lost) energy towards the end of the muon trajectory is the Bragg peak. From MicroBooNE (arxiv: 1704.02927)¶

I. Motivation¶

We know that the energy loss rate of a MIP in argon is ~2 MeV/cm. Hence our goal is to carefully isolate the MIP-like sections of muons (i.e. leaving out the ends of the track), and compute the (reconstructed) energy loss along these trajectories \(dQ/dx\). This can inform the detector calibration effort, for example, since we can compare the peak of the \(dQ/dx\) histogram with the theoretical expected value (although there are many detector effects that make this not straightforward). We can also study the spatial uniformity of the detector by looking at MIP \(dQ/dx\) values in different regions of the detector, etc. If we plot the dQ/dx change along the so-called “residual range” (i.e. distance to the end of the muon trajectory), we get a characteristic plot (due to the Bragg peak). In this tutorial we will focus on reproducing this plot.

Fig. 13 Example of what we expect for a muon dQ/dx versus residual range 2d histogram. The sudden increase in deposited (lost) energy towards the end of the muon trajectory (= low residual range) is the Bragg peak. From MicroBooNE (arxiv: 2010.02390)¶

II. Setup¶

Again, we start by setting our working environment. Some necessary boilerplate code:

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 pandas as pd
import yaml

from plotly.offline import init_notebook_mode, iplot
import plotly.graph_objs as go
init_notebook_mode(connected=False)

c. MLRECO specific imports for model loading and configuration setup¶

The imports below load some auxiliary functions and classes required to run the full chain in interactive mode.

import warnings
warnings.filterwarnings('ignore')
from mlreco.main_funcs import process_config, prepare

Let’s load the configuration file inference.cfg to setup the full chain architecture and weights. This file uses the keyword DATA_DIR to symbolize the path to validation set and checkpoint file. We need to replace it with the actual location defined previously.

cfg = yaml.load(open('%s/inference.cfg' % DATA_DIR, 'r').read().replace('DATA_DIR', DATA_DIR),Loader=yaml.Loader)

To keep things simple we will ask the chain to stay quiet and we can load a bigger dataset than the 100-events default one:

cfg['model']['modules']['chain']['verbose'] = False
# cfg['iotool']['dataset']['data_keys'] = ['/sdf/group/neutrino/ldomine/mpvmpr_012022/test.root']

So far cfg was simply a dictionary loaded from a YAML configuration file. It needs to be consumed by a helper function process_config to be ready for the full chain usage:

process_config(cfg, verbose=False)

Config processed at: Linux tur012 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.¶

This next cell loads the dataset and the model to the notebook environment. hs stands for “handlers”. It contains various useful access points, such as the I/O iterator hs.data_io_iter to directly access the dataset or the trainer instance hs.trainer which will enable us to actually run the network.

# prepare function configures necessary "handlers"
hs = prepare(cfg)
# Optionally, you can specifiy a list of images by listing the image entry ID numbers:
# hs = prepare(cfg, event_list=[0, 1, 2])
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
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)

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

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.

As usual, the model is now ready to be used (check for successful weight loading). Let’s do one forward iteration to retrieve a handful of events.

data, result = hs.trainer.forward(dataset)

III. Setup the analysis tools¶

The particles in a given image are classified into one of five particle types: photon \(\gamma\), electron \(e\), muon \(\mu\), pion \(\pi\), and proton \(p\). Obtaining particles from the full chain is quite simple: we initialize the FullChainEvaluator for this batch of events and examine the particle composition through the get_particles (for true particles, get_true_particles) method.

from analysis.classes.ui import FullChainEvaluator

# Only run this cell once!
predictor = FullChainEvaluator(data, result, cfg, deghosting=True)

entry = 7    # Batch ID for current sample

IV. Visualize true and predicted particles¶

Let’s first import plotting functions from lartpc_mlreco3d for easier visualization:

from mlreco.visualization import scatter_points, scatter_cubes, plotly_layout3d
from mlreco.visualization.plotly_layouts import white_layout, trace_particles, trace_interactions, dualplot

The evaluator will handle ghost masking internally, so all you have to do is retrieve the predicted and true semantic segmentation labels for visualization and analysis:

pred = predictor.get_predicted_label(entry, 'segment')
truth = predictor.get_true_label(entry, 'segment')
points = predictor.data_blob['input_data'][entry]
# Check if dimensions agree
assert (pred.shape[0] == truth.shape[0])

Let’s plot the voxel-level true and predicted semantic labels side-by-side:

trace1, trace2 = [], []
edep = points[:, 5]

trace1 += scatter_points(points,
                        markersize=1,
                        color=truth,
                        cmin=0, cmax=5, colorscale='rainbow')

trace2 += scatter_points(points,
                        markersize=1,
                        color=pred,
                        cmin=0, cmax=5, colorscale='rainbow')

fig = dualplot(trace1, trace2, titles=['True semantic labels (true no-ghost mask)', 'Predicted semantic labels (predicted no-ghost mask)' ])

iplot(fig)

lartpc_mlreco3d Tutorials