Spring 2021—IceCube Research

Ok, not so scientific, but it illustrates a few things. One is that there is an absurd amount of these particles traveling through a given volume or, more precisely, a given surface area at any given time. The other incredible part is how few interact as they flow through us and our surroundings. Professor Jordan Goodman was also at the presentation. He said that he calculated the number of neutrinos a person might expect to interact with their body to be 1 or 2 during their entire lifetime.

The experiment has to deal with a lot of background noise due to many different reasons arising from the fact that neutrinos are indirectly detected. The "noise" amounts to about 100,000 neutrinos of little interest per year, while the signal—the number of interactions the scientists are interested in—amounts to about 100 neutrino detections per year.

    Will asked:

Does the direction of a neutrino's origin or source affect the results detected at IceCube?

    2. When neutrinos travel through space, the probability that they will preserve their initial flavor upon interaction changes with the distance travelled. This phenomena is known as neutrino oscillation.

Here are two graphs depicting probability (the vertical axis) as a function of distance travelled per unit energy (the horizontal axis). Forget about the units of energy for now. You can think of the point on the very upper left corner of each graph as the initial interaction that created the neutrino. You'll notice that at that point the probability is 1.0, or 100%, for a particular line on each plot. On the left it's the blue line which represents the probability that a neutrino will retain its initial flavor over a distance—the red line represents the probability its flavor changes to another kind. On the right it's the black line which assumes an electron type neutrino was initially generated. The blue and red lines plot the probabilities that the electron neutrino switches to a muon type or tau type, respectively. You might also notice that the distances are very long for probability changes with the minimal chance of an electron retaining its flavor occurring around 17,000 km. One thing that stuck out to me is that the probability that an electron type neutrino becomes another type never quite drops to zero and that the probability it becomes a muon type is always higher than the chances it becomes a tau type.

    3. Scientists know that neutrinos have mass because they oscillate.

Initially, the standard model predicted the neutrino would not have any mass. The observation of neutrino oscillation proves that prediction wrong, though. The explanation has to do with quantum mechanics so things quickly get above my pay grade, but essentially a neutrino would not oscillate if it were massless. The experimental findings on oscillation were only recently observed starting in the early 2000s and with Nobel Prize winning work done by the Super-Kamiokande and Sudbury Neutrino Observatories awarded in 2015. The theory was first predicted in the 1950s by Bruno Pontecorvo.

    4. Energy conservation always applies! Even with neutrinos that originate as lower mass lepton type and later interact to produce higher mass leptons.

If an electron type neutrino generated in a star galaxies away from us traveled to IceCube and interacted to create a tau particle, how could that mass discrepancy be reconciled? Conservation of Energy! If that neutrino arrived at the South Pole and struck a water nucleus with enough energy to create a tau particle and the probability decided a tau particle would emerge from the reaction, the energy of the original neutrino would be transformed in a way that accounted for the tau's mass, speed, and whatever other variables that make up it's total energy. If that original electron neutrino didn't have enough energy for the process, there would be no interaction, regardless of the probability at that point. It would just continue traveling until it had another chance to interact in a way that allows for the conservation of energy.

————————————————————————————————————————————————————————————————————————————————————————————

Setting-up Docker and JupyterLab on an Ubuntu Server – 2/28/2021

Docker and JupyterLab will be needed regardless of the specific research I undertake this semester. I installed both of these on an Ubuntu server running Focal 20.04 (LTS) on my local network. 

Docker has a lot of documentation so their website is the place to start. Here's the installation guide link for an Ubuntu OS. It has all the specs and requirements needed for a successful installation. I had no issue following the instructions there. I used the repository method of installation and chose to stick with the default "stable" version.

Getting JupyterLab to run was as easy as copying and pasting the command Erik sent me: docker run -ti --rm -v /Users/blaufuss/jupyter-notebooks:/home/jovyan -p 8888:8888 icecube/icetray-notebook:latest  start.sh jupyter lab. This might require sudo depending on the user privileges and the highlighted part is specific to the filesystem on the server. So, the one here is specific to Erik's filesystem, mine is "home/bryan".

Accessing the JupyterLab requires ssh tunneling. OpenSSH is already installed on an Ubuntu server out of the box...or at least I'm pretty sure it is, so it was just a matter of logging on from my Mac laptop using the correct command:

    ssh -L 5454:localhost:8888 <user@address>

did the trick. To actually access JupyterLab on the server through my laptop, I just had to open a browser tab and enter

    localhost:5454

and copy-paste my token from the terminal connected to the server when prompted. The readout from the docker command in the terminal can be seen below with my token blocked out since it's basically a password.

————————————————————————————————————————————————————————————————————————————————————————————

Binned Analysis: In or Out

Step one to understanding point source searches was to grasp the binned analysis. A binned analysis is kind of what it sounds like. Datum points fall into a bin and are accounted for strictly within that bin. The bin will have some dimensions. For example, if you took the celestial sphere and sliced it into horizontal bands, you would end up with declination bands of a certain width. If you divided it into 25 equally spaced bands, each one would be 7.2 degrees tall—each of those bands is a bin. Now start examining points that correlate to reconstructed neutrino events coming from spots in the sky and you can start populating each of the declination bands with information. Voila, you have bins and information about neutrinos inside them. Here's a plot generated from python code which outputs the number of events per declination bin (vertical, horizontal axes). The sine of the declination was used instead of the declination angle itself.

This graph was made with a bit of python code Erik gave us and a library called matplotlib. It looks nice and all but when it comes to point sources we're dealing with the celestial sphere here—3D! So there should be a way to bin points in the sky. Here's an example.

The map was generated using a different block of code in the same Jupyter Notebook but works based on a different python library called healpy. 

Essentially the idea is that we can look at how many events happen in a bin and compare that to pseudo-randomly generated background. The background is actually made from real data scrambled across declination bands. Since moving a neutrino event to another RA value is essentially equivalent to changing its arrival time in the day, we could randomly assign a new RA value to every neutrino datum while preserving its declination angle. If we plotted all the new data, we would expect to see a randomly distributed background. Now compare the real data to the scrambled data and see if any bin has more events than what you would expect from background only. That's how a binned analysis works in a nutshell.

Below I have inserted the code, outputs, and personal notes describing what the code does in markdown text boxes in between code blocks and outputs. This example only scratches the surface, binning points into pixels on a skymap like the one you see above.

————————————————————————————————————————————————————————————————————————————————————————————

IMPORTANT! GitHub Code and Viewing Instructions

From here on out, the code we're using gets extremely long and this particular version of Google Sites can't handle the amount of HTML and rendering required to keep outputting similar embedded code blocks like the one above. Screenshotting all of it just isn't feasible either due to the length of the scripts and my technical issues. I'll include the snippets that I can muster but you'll have to follow the links provided here to my GitHub account to see all of what Erik provided us and the complete outputs of my code runs.

The main branch for my work with IceCube can be found here: https://github.com/BryanFRezende/IceCube-Point-Source-Analysis. All of the code and files are available from this homepage.

I will link specific file pages as I reference them in my descriptions below so that you can just click on them and see the specific file without having to navigate from the homepage.

IMPORTANT IF THE FILE IS NOT LOADING:

Sometimes Jupyter Notebooks don't load on GitHub. It's an issue with their servers and it's pretty sporadic. If reloading the page a couple of times doesn't solve the issue, it's not because the file isn't there, the servers are just encountering this well-documented issue. To view the file do the following:

• 1. Copy the url on the GitHub page you are trying to view.

  2. Navigate to the Jupyter nbviewer website: https://nbviewer.jupyter.org/

  3. Paste the url into the bar on that page and hit "Go!"

That should display the file. If not, let me know in the comments section below.

Speaking of GitHub, if you don't have an account yet, you need to make one, and it's worth the hour-long intensive to learn how to work with it locally on your terminal. You'll probably use it for the rest of your life if this class was even remotely near what you're interested in getting into as a career. Here's the link to the video I used to learn how to use git.

————————————————————————————————————————————————————————————————————————————————————————————

Artifacts That Tear Up the Sky...and My Data

There are a few problems with the binned analysis. For one it doesn't take into consideration the energy and angular uncertainties of each neutrino event. Another problem is that the contents of neighboring bins aren't accounted for in the contributions to the probability that a spot in the sky is seeing higher than background-only event rates. A third problem that becomes quite obvious as we worked more with the code Erik gave us is that the binned method introduces artifacts, or errors, in the data due to the way the sky is chopped up. You can choose different bin sizes to try to minimize the effects, but there will always be some artifacts. Look at this plot here to see what I mean:

This skymap shows the probabilities of events above or below expected background on a sort of log scale depicted by colors. Blue means lower probability and yellow means higher, but the focus is on the bands of artifacts introduced by the binning method. Here is the code that generated the same map as above but with only 25 declination bands instead of 100.

• Follow this link for the full code and outputs.

Near the bottom of the code on my GitHub, there are a couple of blocks I wrote to plot the most significant point generated from a simulation on a separate skymap. It also follows a similar pattern to the plot above because it's based on the flawed data that generated it.

I also tried writing a script that would produce skymaps with varying levels of declination bands to see which numbers produced the least amount of artifacts, but that became unnecessary once we decided to go for a full likelihood, unbinned analysis.

————————————————————————————————————————————————————————————————————————————————————————————

Unbinned, Full Likelihood Analysis

This is where things get interesting. Instead of proceeding with the bins and investing time into finding parameters that minimize artifacts and account for neighboring contributions, we just moved toward a full likelihood analysis, an unbinned method, if you will. This method is far superior since no artifacts are introduced and the weighted contributions from all neutrinos from a chosen dataset are taken into account to produce the most accurate analysis. The downsides?

    1. It takes way more computing power

    2. We have to unpack this thing:

and this...

Not the easiest task but we were promised it was doable so I trusted Erik.

Essentially, those formulas calculate a test statistic (TS) that boils the whole experiment down to a number. So I ask a question, "Is this particular point in the sky interesting?" I run the code which generates the TS based on the formulas above and my chosen dataset. But what does that value mean? Without context, nothing.

To make sense of the TS I have to develop a framework around the question I'm asking. We chose to go with a dedicated point source search, which asks the question I posed above—does a particular point in the sky have more neutrino events than the expected number of events from background only? Defining the question matters because the amount of work necessary to set up the study varies according to the type of search. An all sky search or a catalog search would be more work, though if I understood correctly, it would be similar.

Anyway, set-up for a dedicated search requires defining a sensitivity and discovery threshold. Think of a maximum and minimum bound, respectively, where the sensitivity is the maximum level where the point source is consistent with a "background only" hypothesis and the discovery is the minimum level required to declare a point source candidate has been identified. Candidate is in bold here because even if a spot in the sky returns a TS value above discovery, there are other factors that need to be accounted for before declaring a discovery to the scientific community.

This image was helpful for me to understand how the thresholds are defined graphically.

The vertical red line is the sensitivity TS threshold and the blue line is the discovery. But where do all the colored curves come from and what do the percentages mean? Allow me to explain: 

As I set up the study, I use my real data during a time window of my choice and a set of Monte Carlo data, to aid in simulation runs. Then code is run that generates background, and various different levels of simulated signal events at the point I'm interested in. I might run dozens of different levels of "injected signal", simulations that randomly produce a number of events above the expected background from the point in the sky I'm analyzing. If I run this simulation, say at an injected signal level of 5 events, 100,000 times, a nice histogram where the fraction of trials (on the vertical axis) for a given TS level (horizontal axis) can be plotted. Observe:

The two graphs are both the result of a set of simulations with the indicated level of signal injections on the upper right-hand corner of each plot. Each of the signal levels (N=5,10, and 50) was generated with 15,000 simulations each. The background was made with 75,000 runs (N=0). The graph on the left plots the fraction of trials returning a particular TS value, while the plot on the right shows the number of events seen for those same trials. Remember, this is in the context of just one spot in the sky. If you drew a vertical line through the peaks of the NS graph, you would see that the line corresponds to the number of injected signal events chosen for that simulation run. But there are deviations! Since each simulation is random, sometimes more signal events are seen and sometimes less. The area under each curve is 1, or in other words, 100% of the events from a chosen injection level are contained in the curve. Duh. But the interesting part comes when using this information to establish the thresholds. Look back at the plot with the red and blue curves. This is what we're looking for with the simulated injection levels—where 90% of the events from some injection level fall above the mean of the background only curve. That's sensitivity. The curve which has 50% of its events fall above the 5 sigma level (5 sigma being 3x10^-7 number of total events in that little portion of the tail) is the injection level that quantifies discovery.

• Look at the full code Erik gave us that generated these plots here.

Since there are so many simulations involved, time is a concern. I added in timestamp functionality to keep track of how long it would take to run a simulation of a given number of trials (a trial is a single simulation for a specified injection level). Since the computing time is pretty linear, I could estimate how long a larger simulation would take and optimize the run based on how much time I had.

————————————————————————————————————————————————————————————————————————————————————————————

PSO J030947.49+271757.31—The Oldest Known Blazar

We had to make a decision on which point in the sky to analyze. Meghan, being the wisest of the group, decided we should stick to the sensible option of analyzing the TXS source that Erik had showed us earlier in the semester. TXS 0506+056 is a blazar that IceCube identified as an interesting neutrino point source in 2017 (Wikipedia). It made sense since we about 3 weeks (being generous here) before our poster session. I of course am far too curious, so I googled the most recently discovered blazar and found the PSO J030... you see it in the title, I'll call it PSO for short (SIMBAD Database, Wikipedia). Dr. Silvia Belladitta and her team over at the University of Insubria discovered it using radio and x-ray telescope data (NVSS and Swift-XRT) in 2019 (Discovery Paper). This thing is like 30 billion light-years away...so I had to suggest we analyze this source as well. Luckily, the team said ok, but I essentially doubled our work. Worth it.

We met with Erik again and go our code running smoothly. Instead of using injected signal numbers we used injected flux normalizations (FN). The FN is just a time-integrated version of injected signal events— it's flux, an amount over time. A small conceptual difference, a big computing power difference. We considered using the HEPCMS cluster to run the code since using FN would slow things down so much, but we ended up deciding the amount of overhead to copy all the data and libraries we needed into the cluster wouldn't justify the return. So we specialized, Meghan spearheaded the TXS point source and I focused on PSO.

Things got very serious, very quick. A lot of background information on the components generating the simulation, TS values, and the statistics behind the analysis. We had to do some coding to make the script fully functional and store our data. I was able to code the following scripts:

    Code to save the simulation outputs as NumPy arrays: 

    Code to generate statistics:

• The full code can be found here.

————————————————————————————————————————————————————————————————————————————————————————————

The Final Results

Lots of work and multi-hour-long simulations later we had done it. Meghan finished the TXS analysis and I had finished compiling the PSO results.

TXS had a significance of 3.5 sigma. Right in line with what we thought it should.

And for PSO?

..a TS of 0. No flux above the expected background rate was detected at that spot in the sky. I guess I knew that going into it because Erik had said that general region in the sky wasn't of particular interest according to their all sky searches. It was still a good learning experience!

Here are the outputs of the results and the full code that generated them. 

 BACKGROUND TS values (Mean, 3 sigma): (0.22467422894638492, 6.498487314264295)    SIGNAL TS vales by flux norm (Mean, 10th percentile): flux_norm = 10^-10.85: (8.705993589288152, 0.0) flux_norm = 10^-10.8: (9.82466105741314, 0.040355408488844696) flux_norm = 10^-10.7: (13.660711126483214, 0.9643661636011038) flux_norm = 10^-10.495: (25.485703253869627, 6.017902801776261) flux_norm = 10^-10.485: (25.857644574847267, 7.288778054381336) flux_norm = 10^-10.48: (26.224145663460703, 6.60428706438666)

Each of the flux levels was generated using 1,200 simulated trials and the background was run using 13,000 trials. These numbers were chosen so that the 90% levels of the flux curves could be sufficiently quantified and so that the background 3 sigma level could be quantified. 3 sigma is something like 15 out of 10,000 so 13,000 was enough for background. We could have run more if we had more time, but this wouldn't have significantly changed the results.

Another thing I want to point out—we used slightly different standards to define the sensitivity and discovery thresholds than what I explained a few logs before this one. We were experiencing statistical noise which placed our sensitivity level above the discovery threshold. To combat the issue, we defined discovery as the FN where 90% of events were above 3 sigma of the background. Notice we chose 3 sigma instead of 5. To quantify 5 sigma we would have needed way more background trials... for TXS. Again, PSO had a TS of zero, so there would be no need to define a stricter discovery threshold, it's consistent with background only.

Since our analysis found TXS to be 3.5 sigma, there is now cause to be more rigorous and generate more background. I think Meghan will likely have more about this on her logbook, but another thing to be aware of here is that since we knew TXS was an interesting source, to begin with, we've introduced a bias into our analysis. By choosing PSO by virtue of it being a recently discovered blazar doesn't incur quite as much of a penalty since we don't have any information as to what it would be doing with its neutrino emissions. Since we knew TXS had significant activity, the 3.5 sigma result is a bit overstated. There are other analyses that IceCube does to account for biases like this one and get more conservative numbers on the significance of their findings. But that's a conversation to be had in another class.

Thanks for reading my logbook! If you have questions or corrections please comment at the bottom of the page. The last bit of information I have are my group's final deliverables you can find below.

Here is our poster:

A huge thanks to Dr. Erik Blaufuss for being so patient and all-around awesome! This semester's work was a major learning experience in manipulating big data and applying statistics to scientific questions.

Dr. Jabeen a massive thanks to you This class was amazing and I have a whole lot of work to show for it. Thanks for pushing us to become better presenters, coders, and researchers-in-training.

————————————————————————————————————————————————————————————————————————————————————————————

END