Machine learning in COMET experiment (part II)

In previous post I've written a short explanation of COMET - Japanese experiment in particle physics. This second post will be devoted to machine learning approach we developed for tracking (currently this is uses the information which is probably unavailable in online, but I still hope we will be able to get some online information and thus apply part of the approach as high-level trigger).

Local wire features

Anyway, the data we are going to use in classification is:

• energy deposit for each wire
• time of deposit for each wire

Let's look at their distributions:

Energy deposited on each 'wire' by signal and background tracks

Time elapsed starting from the moment of triggering till the time when particles are detected.

How do we obtain these values?
COMET uses straw chambers (which we call wires for simplicity), these are long tubes, which fill detector. They are filled with gas, which is ionized by moving charged particles (mostly, electrons).

After ionization, electrons and ions are moving to opposite directions, so we probably can estimate the moment of ionization by drift time (which is of course much greater compared to time of particle flight in detector). Drift time also gives approximate information about the distance between center of wire and track.

Straw tubes used in many detector systems

Finally, each wire has one more characteristic, namely, the distance to target:

Signal hits are more frequent on the inner layers by construction of experiment.

The time is counted from the moment of when trigger worked (it's another subdetector system). The event is 'recorded' starting from that moment.

The basic algorithm that was proposed (baseline algorithm) is using the cut on energy deposition. As you can see, there is really significant difference: energy deposited by background hits is higher. The reason is that background particles are moving slower, so they ionize more particles.

ROC AUC when we use the only feature - energy, is around 0.95, which seems to be very high. Nevertheless, it's not enough, since we have around 4400 wires, of which around 1000 gets activated (this number is called occupancy) within each measurement (event), while the signal track usually contains around 80 points.

So the event is represented using 4400 pairs (energy, time), of which most are zeros.

So the noise which passes through such basic filter is still very significant and there is large room for improvements.

First, let's combine all the information we have about the single wire (distance from center, time and energy deposition), let's call them wire features:

ROC curves (using physical notation here), we see that usage of wire features made us able to twice decrease the background efficiency. 

Neighbors features

Let me remind you, how the whole picture of hits in COMET detector looks like (we use here projection on the plane, orthogonal to beam line):

Blue dots designate background hits, most of which are tracks of different charged particles, say protons; but some of them are simply noise.
Signal hits, which we are looking for,  are red points forming an arc of circle. It's approximate radius is known ahead.

Simple but useful observation one can see in this plot: signal hits nearly always come together.
This implies that we can try using features collected from neighbours of point.

This drives to dramatic improvement of  classification: we almost get rid of random noise tracks, still there are some coupled misidentified tracks. AUC is about 0.995, fpr (background efficiency) decreased by 4-5 times. In principle it suffices to use only information from left and right wires from the same layer.

Ok, I believe that was simple. Is there still something we missed and that could be improved?
For sure, this is the whole shape of the track. It's time to try using this information.

Hough transform and circle shape

Hough transform was initially developed to detect lines and circles. Being quite trivial, it is one of the most effective algorithms in high energy physics.

After using GBDT trained on wire features and features of its neighbours, we are getting quite clean picture with very few false positives. All we want is to detect and remove possible isolated 'islands' with misidentified background hits.

This is done by very approximate reconstruction of track centers. Since we know approximate radius of track center, we can use the Hough transform with fixed radius. It looks like:

Visualization of Hough transform for circles with fixed radius. We are trying to reconstruct the center of track, going through red points. Assuming that we know the radius of fitted track, all possible centers are laying on the circle with center in red point.

We discretize the space of possible track centers and for each point we reconstruct how likely it is the center of track. It is done using sparse matrices and some normalization + regularization (because otherwise tracks with few or many points will have very low/high probabilities).

Once we computed Hough transform, we leave only those centers, where hough transform is high, applying some nonlinear transformation there and applying inverse hough + some filtering. This way we obtain for each wire the probability that it belongs to some signal track.

Finally, we collect all the information we get for each wire:

• local features (energy deposition, timing, layer_id)
• features collected from neighbors
• result of inverse hough transform

And train GBDT on these features to obtain final classifier. It's ROC AUC is 0.9993 (100 times less probability of misordering)

Final classifier is red and it is very close to ideal one. The ROC AUC is about 0.9993

When we are comparing ROC curves at the threshold of interest (with very high signal sensitivity), things are bit worse, but still very impressing: 

At stable benchmark, the background yield decreased by factor of 34. Original ROC curve is not seen on the plot, since it is much lower.

Visualization of all steps

Initial picture of hits in CyDET. Red are signal tracks, blue are background ones.

After we apply initial GBDT (which uses wire features + neighbors), we have some bck hits (to the right, for instance). See the 'island' of misidentified background to the right.
By greed dots we denote the possible centers of tracks. The bigger the point, the greater value of Hough transform image.

Now we apply some non-linear transformations to leave only centers with very high probability. Then applying inverse hough transform and apply second GBDT, which incorporates also information from inverse Hough transform.

Conclusion

I've described how simple machine learning techniques coupled with well-known algorithms can produce a very good result, superior to many complex approaches.
These results are to be checked on better simulation data, when it will be ready. Also we have some ideas about possibilities of classification for online triggers, which use same ideas.

Links

1. Repository with algorithms and results of tracking.
   Most of plots are taken from it, thanks to Ewen Gillies.
2. Reproducible experiment platform used in experiments, gradient boosting was used from scikit-learn
3. Review of straw chambers. Straw chambers are used within many different experiments due to their high resolution and cheapness.

Machine learning used in tracking of COMET (japanese particle physics experiment), part I

Recently I worked together with intern from Imperial College over the tracking system of COMET for two months, and here I'm going to briefly sum our (impressing) results. To begin with, let's explain what this experiment about.

From physics courses we know that there are conservation laws, most famous are conservation of  energy and momenta, but there are other. For instance, conservation of different leptonic family numbers.

For instance, electronic number is number of electrons + number of electronic neutrinos. For some time these numbers (electronic number, muonic number, tauonic number) were considered to be preserved, but in the late 1960's it was proved, that neutrino oscillations change the leptonic number of system (neutrino can become electronic from muonic, for instance). This observation is called lepton flavour violation. However, in the Standard Model overall leptonic number (which is sum of three named family numbers) is conserved.

Since that time physicists are searching for other processes which change the leptonic number and include charged leptons (electrons, muons, tauons), which could be a key for new physics, however no success in this direction and the results are currently mostly upper limits we can prove for some processes. Such possible processes have the name of CLFV (charged lepton flavour violation).

Details of COMET experiment

COMET is small experiment (which is only in plans at this moment) built to detect one specific decay (muon to electron conversion on nucleus), which is suppressed in the standard model:
$$ \mu^{-} + N \to e^{-} + N.$$

The frequency of this process is extremely small in SM: it's probability is about $10^{-52}$, which in particular means that we are sure this will not happen during the experiment.

COMET is sending lots of muons (with fixed energy) on aluminum, and makes the muons to interact with it. The process we are looking for is
$$ \mu^{-} + Al \to e^{-} + Al,$$
while this isn't the only way we can obtain electron from muon
$$ \mu^{-} \to e^{-} + \nu_{\mu} + \overline{\nu}_e.$$

Pay attention,  that second process is 'normal' in the sense it doesn't violate the leptonic conservation laws.

The only way to determine whether we had decay-in-orbit (DIO) or not is to measure the energy of resulting electron (since neutrinos are very elusive particles and we are unable to detect it). Fortunately, the distributions of energies of electron is quite different for these processes and background electrons usually have lower energy.

Distribution of energies for signal process and main background (Decay-in-orbit) in COMET

Some other illustrations of COMET experiment

The scheme of planned experiment COMET. Detector (with aluminum target) will be in the end of this pipeline.
This is what planned during phase II, but for the first time experiment will be simpler

COMET during phase I will be much simpler.

Detector of COMET (CyDET). The detector is filled with sensitive wires (white).
When muon hits the aluminum target (in the center), the electron produced travels over helix trajectories of larger radius in magnetic field and hits wires. The radius of trajectory depends on the transverse momentum of electron.

Data collected from COMET

The main information (apart one from trigger) is the data taken from wires, which are registering when charged particle flies nearby (by adsorbing the particles after ionization) and also it measures the time when this happened.

Another characteristic, which helps in distinguishing electrons, is the radius of trajectory, which is computed (in constant magnetic field) by formula
$$r = \dfrac{p_T}{eB}, $$
so for electron we know the distribution of gyroradius from distribution of energies. Gyroradius is how we actually will reconstruct energy of electron.

Typical signal event in COMET, particle is helixing in magnetic field, leaving energy in drift cambers.

In orthogonal projection of CyDET, the collected information from wires will look like this:

Typical picture of signal event in experiment. Red dots are corresponding to signal hits and form a circle of needed radius.
Also one can see here many background hits, which come from other particles flying through detector.

But the shape of circle with some radius is not the only useful information. During signal hits usually less energy is disposed (which may sound a bit contradictory to what I wrote earlier, but the faster particle flies, the less energy it disposes).

The energy deposited is very strong feature, though leaves significant amount of bck hits, which can spoil the picture, when we have high occupancy.

Apart from this, we can use the 'hit time' (measured from the moment when trigger detected particle).

The time after trigger tends to be smaller, while background hits are distributed uniformly across the time.
This is actually to be checked with better monte-carlo simulations (which are quite poor at the moment).

The main goal of our part is classification of signal vs background hits.

In the next post I'll write about how we applied machine learning to this problem and achieved promising results.

PS. There is second part of this post, devoted to machine learning.

Links

1. Official COMET site.
2. Detailed presentation about CLFV, COMET (and e2mu).
3. Раритеты микромира, an article in russian about rare processes (CLFV is one of examples)  
4. COMET phase-I proposal
5. Repository with algorithms and results of tracking.

Automatic reweighting with gradient boosting

In high energy physics there is frequently used trick: reweighting of distribution. The most frequent reason to reweight data is to introduce some kind of constistency between real data and simulated one. The typical picture of how it looks like you can find below:

The aim is to find such weights for MC, that makes original distribution look like target distribution. I've found a way to use here gradient boosting over regression trees, though it is quite different from usual GBRT in loss, update rule and splitting criterion.

Update rule of reweighting GB
Here is the definition of weights:
$$w = \begin{cases}w, \text{event from target distribution} \\ e^{pred} w, \text{event from original distribution} \end{cases}, $$
so as you see, we are looking for multiplier $e^{pred}$, which will reweight the original distribution.

Pred is raw prediction of gradient boosting.

Splitting criterion of reweighting GB
The most obvious way to select splitting is to maximize binned $\chi^2$ statistics, so we are looking for the tree, which splits the space on the most 'informative' cells, where the difference is significant.

$$ \chi^2 = \sum_{\text{bins}} \dfrac{ (w_{target} - w_{original} )^2}{w_{target} + w_{original} } $$

Here I selected a bit more symmetric version, though this was not necessary.

Computing optimal value in the leaf
since we are going to remove difference in distributions, the optimal value is obviously:
$$ \text{leaf_value} = \log{\dfrac{w_{target}}{ w_{original} }} $$