Monte Carlo Monday

This is Day 1 of my outreach week challenge : what are Monte Carlo samples and why do we need them? Learn how to use simulated data to hunt for new particles.

Gambling in physics?

Let’s get the obvious pun out of the way and focus on the science. No, it’s not about that Monte Carlo…

Where would you rather be?..

Rather than gambling, physicists are often interested interested in working out probabilities, and relying on the Law of large numbers to see them realised (or not) in nature. Take for instance the example of tossing a coin : if the coin is “fair”, then the probability of getting heads is 50%. So is the probability of getting tails, by the way. What does that tell us about the outcome of a single coin toss? Not much…

The underlying idea of classical mechanics is that of determinism : given a set of physical laws (the dynamics of a coin flipping through the air as it moves up and down), and an exact knowledge of the initial conditions (how you flipped the coin), the outcome of that coin toss is in theory computable. The only problem is that it is way too complicated! Too many variables, known with too little precision, come into play; the deterministic character of the physical system is lost behind a veil of uncertainties and errors.

But, were you to repeat that coin toss in roughly similar conditions (again, the experiment is clearly not perfect), these errors and unkowns would average out as you repeat your tosses. That’s how the Law of large numbers saves the day : with enough statistics, you approach the true probabilistic distribution. That is, although you might very well get three tails in a row at first, over a hundred or a thousand tosses you should observe about half of them to be heads and half of them to be tails (if not, you’d conclude there’s something wrong with your probabilites, i.e. with your modelling of the coin).

The case for randomness

This is a very simple example, yet one that has far-reaching consequences. Take the case of particle physics. In the so-called Standard Model of particle physics we have a vast number of free parameters, complicated mathematics and quantum oddities that render the concept of one-off observations completely irrelevant (in fact, quantum mechanics does away with the notion of determinism, but that’s another story entirely). We still need to be able to work out predictions of our model and compare them with data - coming up with sophisticated approaches of simulating reality is therefore crucial, and that’s where Monte Carlo methods come in.

These techniques rely on two core principles : parameterized ignorance, and randomness. The latter is quite straightforward and simply consists in generating random numbers, something we now how to do quite well in computer science; by checking, via the first principle, whether or not they meet some condition, and by repeating the process many times, we can approximate the probability that is the prediction of our model. Indeed, it is often the case that computing these probabilities by hand (i.e. figuring out the 50-50 heads/tails distribution in a coin toss) is either really hard or downright impossible - you have to thank the concept of infinite series for that…

So what’s this “parameterized ignorance” about? It’s simply recognizing that, although we can’t predict a single outcome of a physical process, we have some idea of what might be going on. Take the example of the Higgs boson (to which I’ll come back in a bit) : we don’t observe it directly, rather we detect the lighter particles it decays to. Problem 1 : it doesn’t decay to just one set of particles. The first probability we have to work with is what the final state we observe is. Problem 2 : we don’t know in advance what the properties of these particles are going to be (in particular their energies and velocities). Although we know that taken all together, they must give some total amount of energy, say, we have no idea how this quota is divided amongst the daughter particles. Problem 3 : all this can be influenced by how specifically the Higgs boson decayed, and that added uncertainty (or probability) messes up with all the above.

You might think it’s a downright mess - but sprinkle a little Monte Carlo fairy dust over it all, and things start falling in place nicely. Example.

Rice pi(e)

Draw a one-by-one square on the floor. Inside that square, draw a perfect circle that is tangent to it. Now throw grains of rice at random over the floor, and get this : by counting how many grains fall within the circle, you can estimate $\pi$.

How? Well, there’s a one-to-one correspondence between the surface area of the square (or the circle) and how many (same sized) grains of rice can fit in. That’s what the whole concept of area is about. We know, quite precisely, that the area of the square is $A_{\square}=1^2=1$ and that of the circle is $A_{\bigcirc}=\pi\times \frac{1}{2}^2=\pi/4$. How much is $\pi$ is our parameterized ignorance here. Although we know what the formula is, we’re still unsure (or pretending to be) about how exactly it is realized in nature.

So throw your rice, count how many land on the square, how many land on the circle, do some arithmetic and you’ll get an approximation of $\pi\simeq 3.1415$ :

Here only done on a single quadrant. "After placing 30,000 random points, the estimate for π is within 0.07% of the actual value." (Wikimedia)

Why not try it at home in Python?

from __future__ import division
import random

def MCpi(n):
    return 4 / n * sum(1 for k in range(n) if random.random()**2 + random.random()**2 <= 1)
for npoints in [10,20,30,40,50,100,1000,10000]:                                                                                            
    print "Sampling ",npoints," points : pi ~ ",MCpi(npoints) 
Sampling  10  points : pi ~  2.8
Sampling  20  points : pi ~  3.4
Sampling  30  points : pi ~  3.3333
Sampling  40  points : pi ~  2.9
Sampling  50  points : pi ~  2.88
Sampling  100  points : pi ~  3.12
Sampling  1000  points : pi ~  3.076
Sampling  10000  points : pi ~  3.1564

Crucially, you see that the method converges to the true value of $\pi$ for large enough sampling :

Convergence to pi

MC at the LHC

Although these techniques can be applied to a wide range of problems, from mathematics and physics to biology and finance, I’ll focus here on its applications in high energy particle physics at the Large Hadron Collider, which is the machine I’m working on.

From theory to experiment…

As I’ve already mentioned, a lot of the calculations involved in making a specific prediction are incredibly hard (usually because they involve summing series with infinitely many terms!), and the more unknowns you have (undiscovered yet theorised particles, very rare phenomena, unknown interferences with the experimental environment, etc.) the harder. Turning to numerical solutions is therefore the way to go. Another wide use of Monte Carlo techniques is for detector simulation : emulating the performance of the actual apparatus used to record the data, in order to estimate its sensitivity, efficiency, losses, etc. In a world where the tiniest quantum effects could be signs of major new physics, getting a firm handle on both theoretical and experimental uncertainties is of the upmost importance.

It is also quite convenient to study “individual” (so to speak) processes. I mentioned earlier the case of the Higgs boson decaying to some final state - that would be considered one sample. You’d input as parameters everything you know about the Higgs boson, those final state particles, how they interact, how the detector behaves with them, everything else that might come into play, then start generating unique events, based on random distributions. You seed some randomness into the process (how energetic the Higgs was, how fast the daughter particles are going to fly off from each other), then let the machine compute, following the rules you input, and give you an answer : “the decay products are two electrons and two positrons, with the following properties - blurts out technical data”. That’s coin toss 1. Do that, say, a million times, and you have an idea of the behaviour of that particular process (Higgs decaying to two electrons and two positrons) over a million different events.

That would be one sample, corresponding maybe to the Standard Model situation - if you’re investigating new physics theories, you might want to adjust your parameters and simulate events again. If you’re sticking to the Standard Model, you might want to simulate other processes that end up with this particular final state, or other final states that can originate from the Higgs. And so on and so forth. At the end of the day, you end up with a massive collection of Monte Carlo samples, describing all the relevant physics that correspond to your analysis, in an LHC-ready way (remember that this is simulated data as seen by the detector). Time for the analysis itself…

… and from data to discovery

How was the Higgs boson discovered at the LHC in 2012? Well, at that point we had Monte Carlo simulations for relevant processes without the Higgs. In particular, final states with four leptons (electrons, or their heavier cousins the muons) were of interest, as that’s where the Higgs was thought to show the most clearly if it existed. Once we were satisfied with how well our Higgs-less samples fitted the Higgs-less cases - that is, that they represented accurately physical processes up to the possible existence of this new particle - we started looking at the data.

The GIF below shows the collected data (black dots with error bars), and the estimated background (red columns) corresponding to four leptons that did not originate from a Higgs. The bottom pad then subtracts the data from the prediction, and you can clearly see something fishy happening around 125 GeV… There are more events with four leptons because we “forgot” to include a source of these final states : the Higgs boson! This is then confirmed by zooming in and comparing with a full Standard Model Monte Carlo (blue, including a 125 GeV Higgs particle) - and yes, it does match, we have a discovery!

Building up the Higgs

comments powered by Disqus