Tools from the Couch: MCMC and all that


Presented virtually June 25 2020

Handwritten notes

Some references:
Roberto Trotta: Bayesian Methods in Cosmology
My own (very) slightly longer notes from 2019

Resources:
emcee: an ensemble sampler (you need this for the tutorial!)
pymultinest: a python implementation of the classic nested sampling algorithm
pyGTC: an easy-to-use program for making Giant Triangle Confusograms
corner: another GTC maker
getdist: yet another GTC-maker
Pippi: a more sophisticated program for parsing and plotting MCMC chains

Final exercise:

You checked the arXiv this morning and it appears that the POCO experiment has recorded an unexplained excess number of events in their detector!

POCO is a dark matter direct detection experiment located deep underground in SNAILAB, in Sadbury, Ontario.  It consists of  80 kg of 1,10-Diiododeca-1,3,5,7,9-pentayne (that is,  C10I2). In their exuberance (and possibly due to the deadly fumes), POCO published their data with a world-class chemistry analysis, but forgot to perform a dark matter analysis! Fortunately, you have some code from a previous project that allows you to calculate a projected event rate (including their very well-modeled background), which provides the expected number of events per energy bin, specifically for this experiment! All you need to do now is write a likelihood function (you can use a Poisson likelihood), and run an MCMC analysis. Your job is to find the likely DM mass and cross section to explain these events (which are definitely not due to the nearby potassium mine), and plot the 2d posterior distribution. Once your paper is on hep-ph tomorrow, that phone call from the King of Sweden is pretty much in the bag.

You can set a flat prior on the dark matter mass between 0 and 100 GeV. For the dark matter nucleon cross section, you know it should be somewhere between 10-35 and 10-45 cm2 . A log-flat prior should be sufficient. Finally, the rate code allows for two extra inputs (efficiency and background normalization), you can set those to 1.

You’ll need the rate code here , it uses WIMpy, which you should download and install.

The data points you’ll need (don’t worry about the energy bins, the code takes care of that):

[0,214,360,123,89,84,74,93,115,138,132,143,149,157,146,147,136,126,112,76]