BootCamp_Diffuse_Fit_Exercise.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"BootCamp_Diffuse_Fit_Exercise.ipynb","provenance":[],"collapsed_sections":[]},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.7"}},"cells":[{"cell_type":"markdown","metadata":{"colab_type":"text","id":"W2Nv6vLYV4zS"},"source":["# Diffuse Fit Workshop #\n","\n","Manuel Silva \n","\n","June 19, 2020\n","\n","\\\\\n","\n","*Many thanks to Austin for providing the dataset and Ibrahim for helping put this exercise together!*\n"]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"pMsPo7iBgu6W"},"source":["# Prepare the Environment #\n","* You need all these packages to play around with the data\n","* [numpy](https://www.numpy.org/) and [scipy](https://www.scipy.org/) are for managing and analyzing the data \n","* [matplotlib](https://matplotlib.org/) is good for making plots and visualizing data\n","* [json](https://pythonbasics.org/json/) is a data structure for storing and sharing data\n","* Start by importing them!\n","\n","\\\\\n","*if you are not too familiar with these packages, make sure to follow the links and look at some examples"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"kgvf5gISg1c6","colab":{}},"source":["import numpy as np\n","import scipy as sp\n","import scipy.stats\n","import matplotlib\n","%matplotlib inline\n","import matplotlib.pyplot as plt\n","import json"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"Eqp005RP4dZi"},"source":["# Download the dataset #\n","* Start by downloading the dataset to your machine\n","  * https://icecube.wisc.edu/~msilva/BootCamp_2019/toy_data.tar.gz\n","* If you aren't sure how to extract tarbells, copy and paste the following line:\n","  * tar -xzvf toy_data.tar.gz\n","* Make sure \"hese_toy_data.json\" and \"hese_toy_simulation.json\" can now be accessed on your machine"]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"NOrnyOChwWWf"},"source":["# Load the Simulation (eg expected rates) #\n","* The dataset was saved as a json file\n","* In the real world, you might use \".i3\", \".npy\", \".pkl\" datasets\n","  * Follow the example below, for loading the simulation dataset\n","  \n","* It is always a good idea to \"print\" every few lines and make sure the dataset is properly loaded into your python session\n"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"iQLrjYQbx09i","colab":{}},"source":["json_contents = json.load(open(path+'hese_toy_simulation.json', \"r\"))\n","simulation_mapping = json_contents[\"mapping\"]\n","simulation_events = np.array(json_contents[\"events\"])\n","del json_contents"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"colab_type":"code","id":"vpI4-zghYsVl","colab":{}},"source":["sim_map = sorted(simulation_mapping.items(), key=lambda x: x[1])\n","print (sim_map)"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"colab_type":"code","id":"K6vkE8NYnK89","colab":{}},"source":["print (simulation_events[:2])\n","[[(k, e[i]) for k, i in sim_map] for e in simulation_events[:2]]"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"itCJNUKhmXkJ"},"source":["# Load the Data (eg the observed rates) #\n","* Follow the same procedure to load the data\n","  * remember to change the variable names so you don't accidentally overwrite the simulation"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"yV4OUST0nbaE","colab":{}},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"colab_type":"code","id":"Bam_iHMpnhxV","colab":{}},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"colab_type":"code","id":"l97PkyW-nnnU","colab":{}},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"nEs27hZa1gRL"},"source":["# Computing the Expectation Value (eg $\\lambda$) #\n","\n","* Now let's run IceCube for 1 year, how many neutrinos do you expect to see?\n","  * what happens if we ask for number of neutrinos for energies > 100 TeV?\n","  * what about zeniths > 90?\n","* The best way to compute an expectation value is to use Monte Carlo!\n","* Start by generating neutrinos and \"inject\" them near IceCube\n","  * You usually inject neutrinos following chosen energy, zenith, azimuthal distributions, and over a certain area of IceCube (should we inject directly inside IceCube, further away, deepcore only?)\n","\n","\\\\\n","* Since you choose what events we generate (maybe you want high energy neutrinos vs low energy neutrinos), we need to normalize the simulated events\n","  * Take your ensemble of simulated events and compute a \"gen-weight\" for every event, refer to this \"gen-weight\" as a generated flux:\n","  \n","  $\\Phi_{generated} = N_{generated} [GeV^{-1} * sr^{-1} * cm^{-2}]$\n","  \n","* Now choose a theoretical model, either atmospheric or astrophysical flux model, lets proceed with astrophysical for now\n","  * Take 2 year MESE flux as an example\n","  \n","  $\\Phi_{model} = \\Phi_{MESE} = 2.06 *10^{-18}*(\\frac{E_{\\nu}}{100TeV})^{-2.46} [s^{-1} * GeV^{-1} * sr^{-1} * cm^{-2}]$\n","  \n","\\\\\n","* You can now combine the generated flux and theoretical fluxes to compute event rates\n","  * Each event will now have it's own expectation rate in units $s^{-1}$\n","  * This is literally the expectation (weight) that IceCube \"sees\" an event exactly like this every second\n","\n","  $w = \\frac{\\Phi_{model}}{\\Phi_{generated}}$\n","\n","* $w$ is the expectation value assuming a perfect detector with perfect simulation for that particular event with particular energy/zenith/etc..."]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"M9urtGtQ1m9G"},"source":["# Apply Realistic Weights to the Simulation! #\n","* After taking into account your event selection, detector effects, etc.... we now \"bin\" the simulation into bins such as energy/zenith/topology, take the sum of all the simulated event within a bin and denote as weight \"i\"\n","* If we choose a bin from 100 GeV-1 TeV, zenith < 90, cascade topology, the term $w_{i}$ would tell us how many neutrinos we expect in IceCube with these particular properties per second.\n","\n","  $w_{i} = \\Sigma_{bin-i}\\frac{\\Phi_{model}}{\\Phi_{generated}}$\n","  \n","* Since the HESE dataset measures the astrophysical and atmospheric neutrino fluxes, we need to take both expectations into account\n","\n","  $w_{i,atmos} = \\Sigma_{bin-i}\\frac{\\Phi_{atmos}}{\\Phi_{gen,atmos}}$,\n","\n","  $w_{i,astro} = \\Sigma_{bin-i}\\frac{\\Phi_{astro}}{\\Phi_{gen,astro}}$\n","\n","  $w_{i,total} = w_{i,atmos} + w_{i,astro}$\n","\n","\\\\\n","*  At this point, we are technically done computing our expecation value. \n","*  Recall tht we need to introduce variations to the simulation (our $\\Phi_{0}$ scan yesterday)\n","  \n","\n","# Let's write some code now... #\n","\n","* For the atmospheric component, we will add some overall normalization term ($N_{atmos}$)\n","\n","  $w'_{i,atmos}=N_{atmos}*w_{i,atmos}$\n","\n","* For the astrophysical component, we will add some overall normalization term ($N_{astro}$) along with an additional energy spectral index term ($\\gamma$)\n","  * scale your simulation by the following: $1.0⋅\\frac{E_{\\nu}}{100TeV}^{−2.5}$ \n","  * why did we do this?\n","\n","  $w'_{i,astro}=N_{astro}* \\frac{E_{\\nu}}{100TeV}^{−\\gamma}/ \\frac{E_{\\nu}}{100TeV}^{−2.5} * w_{i,astro}$\n","  \n","\\\\\n","* look at the simulation you loaded earlier (eg \"print simulation_events\")\n","#We gave you the weights!#\n","\n","* you want to takes these weights and rescales them to different $N_{atmos}$, $N_{astro}$, $\\gamma$\n","\n","* write a python function now \n","\n","def weight_event(event, atmos_norm, astro_norm, astro_gamma):\n","\n","    fill in your code here"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"COjJPleBC5wb","colab":{}},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"3vO_194te0J5"},"source":["# Bin the Data! #\n","* We can take all the HESE events and run a fit using the maximum likelihood technique you learned yesterday!\n","  * But remember, there are different detector effects for different energy ranges, zeniths, topologies, so we need to be a little more rigorous here!\n","*  Let's use energy, zenith, and topology as our observables and create some bins \n","\n","*  Start by creating bins for different energy ranges\n","  *  The HESE dataset is primarily for high energy, E = [60 TeV, 10 PeV]\n","  *  Do you think log space or linear space is preferable?\n","*  Now create bins for different zenith ranges\n","  *  HESE is all sky, $\\theta = [0, \\pi]$\n","  *  Do you think we should bin in linear space, log space\n","  *  What about binning it in $\\theta$ vs cos($\\theta$)\n","*  HESE can theoretically measure three different topolgies, but for today lets stick to tracks vs cascades \n","  * Create two bins for topology\n","\n","\\\\\n","\n","* Easiest way to create bins is by using \"[numpy.logspace](https://docs.scipy.org/doc/numpy/reference/generated/numpy.logspace.html)\" and \"[numpy.linspace](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linspace.html)\"\n","* When you are done, create a histogram in matplotlib using energy bins only\n"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"xuSTvs3VZSrF","colab":{}},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"colab_type":"code","id":"v9QOYSRaYC7c","colab":{}},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"QrnLN1qfFSfi"},"source":["# Bin the Data - Advanced #\n","* You can use the bins you defined above and apply it to both simulation and data\n","  * this means that there will now be a set of observed events and expected events, where each bin has its own set of observation and expectations\n","* What you should now do is for each energy, zenith, topologies bin defined above and scan over all possible $N_{atmos}$, $N_{astro}$, $\\gamma$ \n","  * As you can see, this has potential to get very memory heavy, very fast, so just use this code to bin the data and simulation"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"XKdet-gmLyx7","colab":{}},"source":["def make_bin_masks(energies, zeniths, topologies,\n","                   energy_bins=energy_bins, zenith_bins=zenith_bins, topology_bins=topology_bins):\n","    \n","    assert(len(energies) == len(zeniths))\n","    assert(len(energies) == len(topologies))\n","    \n","    n_energy_bins = len(energy_bins) - 1\n","    n_zenith_bins = len(zenith_bins) - 1\n","    n_topology_bins = len(topology_bins) - 1\n","    \n","    energy_mapping = np.digitize(energies, bins=energy_bins) - 1\n","    zenith_mapping = np.digitize(zeniths, bins=zenith_bins) - 1\n","    topology_mapping = np.digitize(topologies, bins=topology_bins) - 1\n","    bin_masks = []\n","    for i in range(n_topology_bins):\n","        for j in range(n_zenith_bins):\n","            for k in range(n_energy_bins):\n","                mask = topology_mapping == i\n","                mask = np.logical_and(mask, zenith_mapping == j)\n","                mask = np.logical_and(mask, energy_mapping == k)\n","                bin_masks.append(mask)\n","    return bin_masks\n","  \n","bin_masks = make_bin_masks(simulation_events[...,2], simulation_events[...,6], simulation_events[...,5])\n","data_masks = make_bin_masks(data_events[...,0], data_events[...,2], data_events[...,1])"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"s_bUmUcCNhT_"},"source":["# Pre-Fit Cross-Checks #\n","* At this point, you need to make sure that your set of observations and expectations make reasonable sense\n","* Write a function that takes the events, the bin masks, and our weights\n","  * Use the following parameters so you can double-check: atmos_norm=1.0, astro_norm=1.0, astro_gamma=2.5\n","\n","def get_expectation(events, masks, weighting):\n","\n","    fill in your code here\n","    \n","    return np.array(\"array of expectation value\")"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"IaEofOSrIBci","colab":{}},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"colab_type":"code","id":"BnLikmOkICNf","colab":{}},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"SPhTu3T0B6RP"},"source":["# Likelihood Evaluation #\n","*  Now we define our likelihood function, which will be a poissonian likelihood (like the previous example), the only difference is that the expectation is a bit more complicated, and we have 3 free parameters instead of one.\n","\n","\n","\n","$\\Large \\mathcal{L}(\\overrightarrow{\\theta | d})=\\prod_{i} \\frac{\\left(\\lambda_{i}\\right)^{k_{i}} e^{-\\lambda_{i}}}{k_i!}$\n","\n","\n","\n","$\\Large \\log \\mathcal{L}(\\theta \\overrightarrow{ |} d)=\\sum_{i} k_{i} \\cdot \\log \\left(\\lambda_{i}\\right)-\\lambda_{i}-\\log\\Gamma\\left[k_{i}+1\\right]$\n","\n","\\\\\n","\n","#Write a function to compute the loglikelihood#\n","\n","\\\\\n","\n","def llh(data, simulation_events, atmos_norm, astro_norm, astro_gamma):\n","\n","    fill in your code here\n","    \n","    return loglikelihood"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"OoS0xx8eCEUp","colab":{}},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"n27Dfk5FCQXt"},"source":["# Test using the data #\n","* Test the function using the data and some assumed values for your free parameters before you move on\n","  * Use the following parameters so you can double-check with us: atmos_norm=1.0, astro_norm=1.0, astro_gamma=2.5"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"sYhekmP-CKgV","colab":{}},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"IDwCAGD0CHo0"},"source":["# Minimize the llh function#\n","\n","Earlier today, we used a minimizer for 1 free parameter. Now use the same function, but modify it to fit this exercise."]},{"cell_type":"code","metadata":{"colab_type":"code","id":"F7dS9j6PCGHa","colab":{}},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"D4tw5wxYDLt1"},"source":["# Check your Final Results #\n","\n","* \"results\" contains the output of the fit\n","* \"results.x\" contains an array of the best-fit values"]},{"cell_type":"code","metadata":{"colab_type":"code","id":"zirwO2x7DRsC","colab":{}},"source":["print(results)\n","print(results.x)"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"colab_type":"code","id":"DD5kLBtSDXwH","colab":{}},"source":["print ('Atmospheric Normalization = ', results.x[0])\n","print ('Astrophysical Normalization = ', results.x[1])\n","print ('Astrophysical Gamma = ', results.x[2])"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"colab_type":"text","id":"X7pFeUQTcWo5"},"source":["# Congratulations! #\n","\n","You have officially completed a 3-dimensional diffuse fit!\n","Now go have a cup of coffee and relax...\n","\n","![alt text](https://amp.thisisinsider.com/images/5ae23026bd967120008b46ac-750-563.jpg)"]}]}