Icecube Research

This page contains research notes and documentation on work for Icecube in the Spring Semester of 2016. A presentation of nice pictures and whatnot can be found here.

1. h1 = plt.hist2d(cos(data["zenith"]),data["azimuth"],bins=100) //now we can do some 2d histogram plotting

   1. This is what we see from this. We see the same distribution as above horizontally (that's the zenith angle) but looking vertically we get a much more of a homogeneous distribution. That means that there is little difference how many neutrinos the detector sees depending on what direction tangent to the surface of the south pole you look. That is, except for a few spikes in the down-going neutrino side. Let's take a closer look.

1. h1 = plt.hist(data["azimuth"],bins=100)//let's look at the azimuthal histogram

   1. This is interesting. We see 6 very distinct spikes in the number of events at even spacing in azimuth angle. This only happens to the down-going data, as we saw above. I guess that this is because the detector is a hexagonal grid and that may put some directional asymmetry in the data. I wonder if they should have cut this out of the final data. But they know what they're doing.

Week Five-Six Notes (3/10/16-3/24/16)

1. Continuing from that histogram of the azimuth, the reason why we got those 6 bumps is because tracks that go right down a row of DOMs will get better sensitivity than otherwise. This is why we use right ascension instead of the azimuth because azimuth is based on the orientation of the detector but right ascension stays constant relative to the stars and icecube rotates in it. So when we plot right ascension (via "h1 = plt.hist(data["ra"],bins=100)") and get:

• 1. and that looks much more uniform. Then we can plot that in 2d and get a more true histogram (by "h1 = plt.hist2d(cos(data["zenith"]),data["ra"],bins=100)") and get this:

1. Now that looks more true to life. The next step is to try to plot this in a mollweide projection. We can use healpy (sudo apt-get install python-healpy) which does a good job of indexing the sky so that each bin is equal sky.

2. You know, this code is getting a bit complicated. We should start making it so we can do all this in one step. I hear shell scripts coming on! To make one, we type "touch scriptname.sh" and then use "emacs scriptname.sh".

   1. Fun emacs shortcuts:

      1. C-<cha> means press ctrl and whatever character at the same time. M-<cha> means press alt and whatever character

      2. C-x C-c ends the session

      3. C-g quits a command

      4. C-v is page down, M-v is page up, C-l (L) centers the cursor

      5. C-p (previous) moves up one line, C-n (next) moves down one line, C-b (back) moves back one space, C-f (forward) moves forward one space. These replace arrow keys

      6. M-f and M-b move whole words at a time forwards or backwards

      7. C-a and C-e move to the beginning and end of a line. M-a and M-e move to the beginning and end of a sentence

   2. Here I wrote a shell script which calls a python program to do all these things we just talked about with one line. I type in "./2dhist_radec.sh" and it runs. Admittedly this is not super elegant. It also has the problem that I can't get the histogram to stay displayed. It appears then disappears immediately, so I worked around it by saving a png of the plot then opening it with eog ("sudo apt-get install eog"). Oh also you'll need to change permissions of the file with "chmod 777 <scriptname.sh>"

1. Of course, this still has the problem that we are plotting onto a rectangle. But the Earth is a circle. We like circles. Let's plot more circl-ey. Luckily healpy is a great library for plotting mollweide projections of data. I used Erik's code here to do this (we can color-code the code in emacs by typing "M-x python-mode"). That code looks like:

   1.  

• 1. and creates:

1. Now we can consider the p-values of each pixel. That is, how likely is each pixel to be a point source, given the background. We use the data itself to try and guess how many events we expect for a given pixel in a given declination, then we compare the expected value to what we actually see in each pixel and calculate that p-value assuming that the probability follows a poisson distribution.

   1. Here is the code that does that:

• 1. The line to calculate n[i] is admittedly confusing because we have to normalize our histogram to account for the number of pixels in a bin of declination. The math behind this horrid line is here:

• 1. And we get this graph of log(p-value) vs log(frequency):

• 1. And this map of pvalues (notice that there are no noticeable patterns)

Week Seven-Eight Notes (3/25/16-4/6/16)

1. Week seven was a bit light on work because Erik (our mentor) left for a conference. We completed the work we had to do mostly in the eighth week.

2. We learned a better way to run a python script. Once you have started python ("ipython -pylab") we can run our script ("%run -i <scriptname>.py") and it runs from inside python, then whatever variables (or plots) were made inside the script are now in the work space. So if pval.py makes a histogram called h1, we can look at it by running "h1.draw()" from inside python.

3. One other edit I made to the pvalue code is that I also looked at the nearest neighbors to a given pixel and added them in to the consideration of whether there was a point source or not there. This is to accommodate for the fact that the uncertainty in degree of a given event is on the order of the width of a pixel, so we want to account for that in the pvalue test. This required two additions to the code:

• 1. A couple effects arose from this. First, there are a few more very low p-value events as can be seen in the first graph below, and second, near the equator and at a couple other declinations, there are noticeable strips of lower pvalues. I think that this is because there are really big changes between bins in declination at those lines, and so when we look at several bin that overlap across those shifts in declination, the poisson distribution sees a dramatic divergence from what it expects. This is seen in the second graph

1. The next thing we have to do is to run a trials corrected p-value test on the data. To do this, we need to see what the distribution of p-values is in a sample of data that we know there are no point sources in. The easiest way to check is to scramble the actual data in right ascension so we know there are no real point sources and see what the p-value distribution looks like.

   1. The edits to the code make it become:

• 1. And makes a pval distribution given no point sources that looks like this after scrambling in right ascension ten times and adding all the pvalues into one histogram:

• 1. And scrambling 1000 times gives this distribution:

Week Nine Notes (4/9/16-4/13/16)

1. So on the agenda: we need to fix the expectation calculations when looking at nearest-neighbor pixels, we also need to change the histogram of the distribution of p-values to be a probability distribution. This means that we need to normalize the histogram so that it has an integrated area of one.

2. To fix the nearest neighbors calculations, we need to use the expectation of each pixel and sum them together instead of multiplying the expected value for the central pixel by 9 (because there are 8 neighbors and the pixel itself). We can do this by implementing the following code:

• 1. Which will produce a new map that no longer has the statistical artifacts from rapid changes in expectation between declinations. This can clearly be seen in the skymaps below. The left hand map shows the problem and the right hand one shows what the revised code does.

• 1. And the pvalue distribution looks like this

1. Ok, now we can implement that in the pvalue trials correction code and then normalize the resulting histogram. This can be accomplished by adding the condition "normed = True" to our plt.hist line that creates the histogram. Our distribution of pvalues from data that was scrambled in right ascension, after running 1000 trials, looks like this:

• 1. One interesting note about this is that the slope is much more linear than it looked before we were accounting for nearest neighbors, so I suspect that one ramification of the nearest neighbors calculation is to smooth out the pvalues in some sense.

Week Ten-Eleven Notes (4/14/16-4/27/16)

1. So the first thing we're going to try this week is to adjust how we are getting an expected value for our poisson distribution in each declination band. The old way was to make a histogram with 100 bins in declination, try to normalize that with some messy trig to get the rough number of events per pixel area in a band. Then when you wanted to find an expectation for a given pixel later, you had to round its declination to the nearest bin. This whole situation has lots of over complexity. What we are going to try to do now is to get the expectation by finding all the pixels in a given declination by healpy, add the events in all those pixels up, then dividing by the number of pixels in that band. This will give us a more accurate expectation, and will make finding the right expectation later more simple. The changes to the code to do this are here:

• 1. Which gives me this skymap of p-values:

• 1. And the distribution of p-values looks like this:

• 1. So this new technique gives more low p-values than normal. However this brings us to a rather big problem that we've just realized we have to deal with...

1. These p-values are way too low for what we expect to see. Given that there are papers saying there are no point sources, and we know there are no point sources in this data, a p-value of 10^-13 should only happen once in ten trillion pixels. That map only has 12288 pixels, so we shouldn't see many points above around 10^-4. What this means is that somehow we are underestimating our expected background on a significant number of points. Our expectation could just be two or three times too low and we could have this problem. Where is the error coming from?

   1. We tried adding the expectation up for each pixel and comparing that to the total number of events in the data set. This should be the same if our expectations are correct. We get 121089.3 from the expectation and 122691 from the data. That difference can be explained away by events we don't account for in the polar regions. So no help there.

   2. Next we can look at what distribution of values we get for the number of events in a given declination and then see if it looks very poisson-ey to us.

      1. We make these plots and compare them to the poisson distribution we are using to describe the event distribution in that declination and, as you can see below (these are just 3 of the 124 declination distributions we made), the fits are very believable:

• 1. Now that that's solved we need to do one more big thing. We need to start injecting signal into our scrambled datasets to see how many neutrinos we need to see from a point source before it becomes statistically signifigant. To do this, we will need a distribution of the highest p-values we get from a bunch of scrambles to do trials corrections, and we need code that gives us the p-value of an injection of x neutrinos and y declinations.

     1. We plotted the lowest p-value from 1000 scrambles and get this distribution:

• • 1. And we can see here that we don't see any p-values above 10^-11 from random fluctuations in the background.

Week Twelve-Thirteen Notes (4/28/16-5/12/16)

• 1. So now that we know we're calculating the p-values correctly, we need to figure out how many neutrinos we'd need to see from a hypothetical point source above the background before we see it as the most significant point in the sky 50% of the time (that will be our discovery potential). We will only look at 7 or so declinations for simplicity, and because the background doesn't change too incredibly drastically in declination. So basically we're going to try to make some code to scramble all our data in right ascension, add a few extra events to one pixel, and see if it's the highest p-value on the map and has high enough statistics to be distinguishable reliably. Then we repeat this a whole bunch of times for 7 declination with different amounts of added events.

  2. So how much is significant? Well after 1000 background scrambles to get a trials adjustment distribution (see above) we find that the highest p-value from random background was 1.713e-11. This means the right tail .1% (or 3 sigma) cuttoff is about 1.7e-11, so if our injection point is both the highest point on the map and if it is above that 1.7e-11, then we consider that a success.

  3. Next we wrote some code ("injection_pval.py") that adds extra events over and over to scrambled maps and creates a distribution of "what percent of injections are above the threshold" indicated by color plotted against the number of injections on the x axis and the declination on the y axis.

• 1. Awesome! It is interesting to note that this actually closely resembles the distribution of events in declination we got at the very beginning, but this makes sense as our test statistic is based on a poisson distribution, who's standard deviation is related to the number of events in a bin. Funny how things work out. Now that we have that, we need to convert "number of events we need to see in icecube" to neutrino flux since that is a more universal value that anyone can use. We got conversion factors from Erik (which vary in declination and depend on the cross sectional area of the detector versus the energy distribution of neutrinos and the cutting criteria for events). Then we used those values to create a discovery potential for IceCube at the 3 sigma level (we'd do 5 sigma, but that requires orders of magnitude more computer time).

• • 1. Which is very similar in shape to IceCube's own results after 7 years and much more sophisticated analysis (albeit they are about 100 times more sensitive and 2 sigma above us, but they have 7 times more data and much more complicated algorithms).

• 1. For just three months of work by 3 undergraduates, I'd say we did pretty good!

  2. Attached at the bottom are our final code, a zip with all our final plots, and a pdf of our poster and (eventually) our paper. Unfortunately the data we used is not publicly available. All publicly available data can be found here: http://icecube.wisc.edu/science.

References

1. This is the "about us" page on the Icecube public webpage. It has general information about who they are and what they do. http://icecube.wisc.edu/about/overview

2. This is a link to the main Icecube page for the UMD part of the collaboration. http://www.icecube.umd.edu/

3. These are some cool facts about Icecube. http://icecube.wisc.edu/about/facts

4. This is a page with research highlights of the experiment. http://icecube.wisc.edu/science/highlights

5. This has some great information on neutrinos: https://icecube.wisc.edu/outreach/neutrinos