# Simulating Cosmic Ray Showers using Monte Carlo in Python 3

This is a final project using python to simulate cosmic ray showers. Upper division undergraduate physics knowledge required, specifically in special relativity and quantum mechanics. Python knowledge in Monte Carlo simulation (important), discrete derivatives, ODEs (Euler / Runge-Kutta), and NumPy and mathplotlib arguments required. Please add comments throughout the code. For the ipynb files, I use nteract.

The project entails building a Monte Carlo simulation in 2 dimensions that will reproduce the cosmic-ray air shower effect. The overall algorithm includes:

1. Keeping track of the following information for all particles: identity (proton, muon,...), current position, and current velocity.
2. At each timestamp dt, use Monte Carlo to test for random decay/interaction, and split into multiple particles if such an event happens. If it does happen, use sub-algorithm:
3. 2.a Choose a random direction for one of the decay products
4. 2.b Momentum vectors of all decay products determined using conservation of energy and momentum
5. 2.c Using known velocity of decaying particle, apply Lorentz boost to find momentum vectors of all decay products in lab frame.
6. 2.d Concert from momentum to velocity, and update list of particles.
7. Otherwise, advance the positions of all particles as (dx,dy) = (vx*dt,vy*dt) (This is all detailed in "cosmic_assign.ipynb)

I have added the files necessary. Look at "cosmic_assign.ipynb" first. This file gives complete description of concept, physics equations, the main loop, and required API. Read it carefully. All functions in API must be used. "Cosmic_API_tests_full.ipynb" is what it sounds like, it contains tests for the API functions. "Cosmic_presentation. ipynb" contains guidelines for the final presentation of results. "cosmic.py" is the module you will add code to to eventually call in the Monte Carlo sim (so far only contains "import NumPy as np"). The two files that will be turned in are "cosmic.py" and "cosmic_presentation.ipynb".

VERY IMPORTANT: I have set the due date for this project to be May 29th. However, I request a project checkpoint no later than May 19th. The complete project does not have to be completed at this date. Just a checkpoint to see the progress. Please let me know ASAP if any information is missing or if you need any files in a different format (I included each as a .ipynb and .pdf files)

## Get Help With a similar task to - Simulating Cosmic Ray Showers using Monte Carlo in Python 3

In [ ]: %matplotlib inline from cosmic import * import matplotlib.pyplot as plt 1. lorentz_beta_gamma Verify that lorentz_beta_gamma computes the correct values of and . Test input: (beta_x, beta_y) = (0.3, 0.4) Expected output: In [ ]: test_beta_xy = (0.30, 0.40) test_beta, test_gamma = lorentz_beta_gamma(*test_beta_xy) print(test_beta, test_gamma) assert np.abs(test_beta - 0.5) <= 1e-6 assert np.abs(test_gamma - 1.15470) <= 1e-6 2. boost_momentum Using the test beta/gamma from above, boost a known momentum to another frame of reference. Test input: (p_x, p_y) = (40, -40) MeV/c, for a pion (mass m = 139.6 MeV/c**2.) Same (beta_x, beta_y) = (0.3, 0.4) from above. Expected output: (calculated by hand) Energy: MeV In [ ]: test_p_xy = (40., -40.) test_p = np.sqrt(test_p_xy**2 + test_p_xy**2) test_m = 139.6 test_E = np.sqrt(test_p**2 + test_m**2 ) print(test_E) px_boost, py_boost = boost_momentum(test_m, test_p_xy, test_p_xy, test_beta_xy, test_beta_xy) print([px_boost, py_boost] ) assert np.abs(test_E - 150.63)/150.63 <= 1e-4 assert np.abs(px_boost + 12.92)/12.92 <= 1e-4 assert np.abs(py_boost + 110.56)/110.56 <= 1e-4 3. check_for_interaction Several test cases for interactions with very high/low probability, so the outcome should always be True/False. In [ ]: ## Unyt version h = 15 * 1000 / 3e8 # km --> light-seconds test_p1 = velocity_to_momentum(particle_mass('electron'), 0.3, 0.4) test_part_1 = ('electron', test_p1, test_p1, [], []) print("Electron interacts: (must be False)") print(check_for_interaction(test_part_1, 1e6, h)) test_p2 = velocity_to_momentum(particle_mass('muon'), 0.3, 0.4) test_part_2 = ('muon', test_p2, test_p2, [], []) print("Muon interacts in 1e-10 seconds: (False)") print(check_for_interaction(test_part_2, 1e-10, h)) print("Muon interacts in 100 seconds: (True)") print(check_for_interaction(test_part_2, 10, h)) test_p3 = velocity_to_momentum(particle_mass('proton'), 0.7, 0.4) test_part_3 = ('proton', test_p3, test_p3, [], []) print("Proton interacts in 1 second: (True)") print(check_for_interaction(test_part_3, 1 , h)) test_p4 = velocity_to_momentum(particle_mass('proton'), 0.3, 0.2) test_part_4 = ('proton', test_p4, test_p4, [], []) print("Slow proton doesn't interact (below threshold): (False)") print(check_for_interaction(test_part_4, 1, h )) 4. proton_downscatter Verify that proton downscattering returns decay-product tuples, and that total momentum is zero in the proton's rest frame. In [ ]: part1, part2, part3 = proton_downscatter() print(part1,'\n', part2, '\n', part3) # Should be a proton and two pions as decay product tuples # Check momentum conservation: should print 0 twice print(part1 + part2 + part3) print(part1 + part2 + part3) 5. Main loop test In [ ]: import unyt as u print(1*u.m) # Define light-second distance u.define_unit('ls', 1*u.s * u.physical_constants.c) (1*u.ls).to('km') In [ ]: # Proton shower px, py = velocity_to_momentum(particle_mass('proton'), 0.959999, -0.28) particles = [ ('proton', px, py, ( * u.m).to_value('ls'), ([20e3] * u.m).to_value('ls')), ] particles = run_cosmic_MC(particles, 1e-5, 50) pc = { 'proton': 'black', 'pion': 'red', 'muon': 'green', 'electron': 'blue', 'neutrino': 'gold', } #print(particles) for particle in particles: plt.plot(particle, particle, color=pc[particle]) In [ ]: # Muon shower px, py = velocity_to_momentum(particle_mass('muon'), 0.92, -0.24) init_particles = [ ('muon', px, py, np.array(), np.array([1e-4]) ) ] # ~ 30 km height in light-sec particles = run_cosmic_MC(init_particles, 2e-7, 100) for particle in particles: plt.plot(particle, particle, color=pc[particle]) β γ β = 0.5 γ = 1/ = 2/ ≈ 1.15471 − β2‾ ‾‾‾‾‾√ 3‾√ E = ≈ 150.63+m2 p2‾ ‾‾‾‾‾‾‾√ ( , ) = (−12.92, −110.56)p′x p′y

import numpy as np

Final Proj/cosmic_API_tests_full.ipynb { "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "from cosmic import *\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. lorentz_beta_gamma\n", "\n", "Verify that lorentz_beta_gamma computes the correct values of $\\beta$ and $\\gamma$.\n", "\n", "__Test input:__ (beta_x, beta_y) = (0.3, 0.4)\n", "\n", "__Expected output:__\n", "\n", "$$\n", "\\beta = 0.5 \\\\\n", "\\gamma = 1/\\sqrt{1-\\beta^2} = 2/\\sqrt{3} \\approx 1.1547\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "test_beta_xy = (0.30, 0.40)\n", "\n", "test_beta, test_gamma = lorentz_beta_gamma(*test_beta_xy)\n", "print(test_beta, test_gamma)\n", "\n", "assert np.abs(test_beta - 0.5) <= 1e-6\n", "assert np.abs(test_gamma - 1.15470) <= 1e-6" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. boost_momentum\n", "\n", "Using the test beta/gamma from above, boost a known momentum to another frame of reference.\n", "\n", "__Test input:__ (p_x, p_y) = (40, -40) MeV/c, for a pion (mass m = 139.6 MeV/c**2.) Same (beta_x, beta_y) = (0.3, 0.4) from above.\n", "\n", "__Expected output:__ (calculated by hand)\n", "\n", "Energy: $E = \\sqrt{m^2 + p^2} \\approx 150.63$ MeV\n", "\n", "$(p_x', p_y') = (-12.92, -110.56)$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "test_p_xy = (40., -40.)\n", "test_p = np.sqrt(test_p_xy**2 + test_p_xy**2)\n", "test_m = 139.6 \n", "\n", "test_E = np.sqrt(test_p**2 + test_m**2 )\n", "print(test_E)\n", "\n", "px_boost, py_boost = boost_momentum(test_m, test_p_xy, test_p_xy, test_beta_xy, test_beta_xy)\n", "print([px_boost, py_boost] )\n", "\n", "assert np.abs(test_E - 150.63)/150.63 <= 1e-4\n", "assert np.abs(px_boost + 12.92)/12.92 <= 1e-4\n", "assert np.abs(py_boost + 110.56)/110.56 <= 1e-4\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. check_for_interaction\n", "\n", "Several test cases for interactions with very high/low probability, so the outcome should always be True/False." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Unyt version\n", "\n", "h = 15 * 1000 / 3e8 # km --> light-seconds\n", "\n", "test_p1 = velocity_to_momentum(particle_mass('electron'), 0.3, 0.4)\n", "test_part_1 = ('electron', test_p1, test_p1, [], [])\n", "print(\"Electron interacts: (must be False)\")\n", "print(check_for_interaction(test_part_1, 1e6, h))\n", "\n", "test_p2 = velocity_to_momentum(particle_mass('muon'), 0.3, 0.4)\n", "test_part_2 = ('muon', test_p2, test_p2, [], [])\n", "print(\"Muon interacts in 1e-10 seconds: (False)\")\n", "print(check_for_interaction(test_part_2, 1e-10, h))\n", "\n", "print(\"Muon interacts in 100 seconds: (True)\")\n", "print(check_for_interaction(test_part_2, 10, h))\n", "\n", "test_p3 = velocity_to_momentum(particle_mass('proton'), 0.7, 0.4)\n", "test_part_3 = ('proton', test_p3, test_p3, [], [])\n", "print(\"Proton interacts in 1 second: (True)\")\n", "print(check_for_interaction(test_part_3, 1 , h))\n", "\n", "test_p4 = velocity_to_momentum(particle_mass('proton'), 0.3, 0.2)\n", "test_part_4 = ('proton', test_p4, test_p4, [], [])\n", "print(\"Slow proton doesn't interact (below threshold): (False)\")\n", "print(check_for_interaction(test_part_4, 1, h ))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. proton_downscatter\n", "\n", "Verify that proton downscattering returns decay-product tuples, and that total momentum is zero in the proton's rest frame." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "part1, part2, part3 = proton_downscatter()\n", "print(part1,'\\n', part2, '\\n', part3)\n", "# Should be a proton and two pions as decay product tuples\n", "\n", "# Check momentum conservation: should print 0 twice\n", "print(part1 + part2 + part3)\n", "print(part1 + part2 + part3)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Main loop test" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import unyt as u\n", "\n", "print(1*u.m)\n", "\n", "# Define light-second distance\n", "u.define_unit('ls', 1*u.s * u.physical_constants.c)\n", "(1*u.ls).to('km')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Proton shower\n", "\n", "px, py = velocity_to_momentum(particle_mass('proton'), 0.959999, -0.28)\n", "\n", "particles = [\n", " ('proton', px, py, ( * u.m).to_value('ls'), ([20e3] * u.m).to_value('ls')),\n", "]\n", "\n", "particles = run_cosmic_MC(particles, 1e-5, 50)\n", "\n", "pc = {\n", " 'proton': 'black',\n", " 'pion': 'red',\n", " 'muon': 'green',\n", " 'electron': 'blue',\n", " 'neutrino': 'gold',\n", "}\n", "\n", "#print(particles)\n", "\n", "\n", "for particle in particles:\n", " plt.plot(particle, particle, color=pc[particle])\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Muon shower\n", "px, py = velocity_to_momentum(particle_mass('muon'), 0.92, -0.24)\n", "\n", "init_particles = [ ('muon', px, py, np.array(), np.array([1e-4]) ) ] # ~ 30 km height in light-sec\n", "particles = run_cosmic_MC(init_particles, 2e-7, 100)\n", "\n", "for particle in particles:\n", " plt.plot(particle, particle, color=pc[particle])" ] } ], "metadata": { "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.8" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 } __MACOSX/Final Proj/._cosmic_API_tests_full.ipynb Final Proj/cosmic_assign.ipynb { "cells": [ { "cell_type": "markdown", "source": [ "# PHYS 2600 - Final Project 1: cosmic\n", "\n", "## Monte Carlo Simulation of Cosmic Ray Showers \n" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "1eac413f25017b99f6479e25bc88bd8c", "grade": false, "grade_id": "cell-142a0448af844469", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "This is one possible _final project_; __you must choose one (and only one) project to complete and hand in.__ \n", "\n", "The deadline for this project is __12:00 midnight (Boulder time) at the end of Thursday, 30 April.__ Since final grades must be submitted to the university shortly after the end of the final exam period, to allow for timely grading _no late submissions will be accepted._\n", "\n", "You are permitted to discuss the final projects with other students in the class, including both the physics and computing aspects, but __your submitted project must be entirely your own__; no direct copying of Python code is permitted, and you must write up the presentation in your own words.\n", "\n", "When you submit your project, you must provide __two (2) files__:\n", "\n", "* cosmic.py: a Python module (see lecture 22) which implements the functions listed in the Application Progamming Interface (API) below.\n", "* cosmic_presentation.ipynb: a Jupyter notebook which uses the code from your Python module to answer the physics questions. Use of Markdown text, MathJax equations, and plots are all encouraged to make your presentation clear and compelling!\n", "\n", "The rubric for grading will be split 50/50 between the code (in your .py module and Jupyter notebook) and the response to the physics questions (provided in cosmic_presentation.ipynb):\n", "* Basic functionality of code (passes provided tests, follows API specification): __30%__\n", "* __Six (6)__ additional API tests _you write_ (tests are provided and are correct and useful): __10%__\n", "* Documentation of code (docstrings and comments make the code clear): __10%__\n", "* Correct and detailed answers to physics questions: __40%__\n", "* Quality of presentation (clear and readable, appropriate use of Markdown, equations, and plots): __10%__\n", "\n", "A _bonus of up to 10%_ is available for the extra-credit \"challenge\" physics question at the end. (These challenge questions are meant to be really hard! But partial credit will be awarded if you make some progress on them.)" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "548c5a59204c1bd5cf17cafa7eb3c9bc", "grade": false, "grade_id": "cell-1abec7f2b0d59fa6", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "# Overview\n", "\n", "Every second, the Earth is bombarded by __cosmic rays__: high-energy particles flying at us from the cosmos. These naturally-accelerated particles have been a treasure trove for research, leading to discoveries of new elementary particles such as the positron and the muon in the early 20th century. The collisions of a cosmic proton with Earth's atmosphere lead to a spectacular multi-pronged trail of new particles called a __cosmic ray shower__.\n", "\n", "For this project, you will build a Monte Carlo simulation (in two dimensions to make things simpler) that will reproduce the cosmic-ray air shower effect. The set of interactions and particles you study will be limited to five: the proton, electron, pion, muon, and neutrino. You'll need to account for special relativity, which is a crucial ingredient since the cosmic ray particles are moving very fast!\n", "\n", "In addition to making nice simulated pictures of a spectacular physical phenomenon, you'll be able to directly investigate how the flux of cosmic rays in the upper atmosphere translates into detection of muons in ground-based observatories." ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "d19418b93c3d733f735ae89b67783fce", "grade": false, "grade_id": "cell-d0bddfdb0dd97095", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "# Physics background and important equations\n", "\n", "For this project, we will adopt a _simplified model_ of particle physics, which omits a lot of complicated details but includes enough physics to get the right qualitative behavior of cosmic ray showers.\n", "\n", "We will consider _only_ the following set of particles, and we will ignore their electric charges completely:\n", "\n", "| name | symbol | mass |\n", "|------|---------|-----------|\n", "| proton | $p$ | 938.3 MeV/c^2 |\n", "| electron | $e$ | 0.5110 MeV/c^2 |\n", "| pion | $\\pi$ | 139.6 MeV/c^2 |\n", "| muon | $\\mu$ | 105.7 MeV/c^2 |\n", "| neutrino | $\\nu$ | $\\sim$ 0 |\n", "\n", "(\"MeV\" stands for \"mega-electron-volt\", a unit of energy; electron-volt based units are common in particle physics.) The proton and electron are familiar; the muon is a heavier (unstable) cousin of the electron, with the same charge. Neutrinos are ghostly particles with nearly zero mass and very weak interactions. Finally, the pion is a bound state of the strong nuclear force, like the proton.\n", "\n", "We're missing some details from the [real Standard Model of particle physics](http://www.particleadventure.org/standard-model.html), including the distinction between different types of neutrinos, the anti-particles for each of the above particles, and neutrons. But this is enough to get reasonable cosmic ray physics!\n", "\n", "We will consider __only the following interactions__ in our model:\n", "\n", "1. Proton downscattering: $p \\rightarrow p + \\pi + \\pi$\n", "2. Pion decay: $\\pi \\rightarrow \\mu + \\nu$\n", "3. Muon decay: $\\mu \\rightarrow e + \\nu$\n", "\n", "The first process is due to collisions of the cosmic ray proton with the nuclei inside of air molecules in the Earth's atmosphere; this inelastic collision produces a pair of pions (due to charge conservation: the final state is really $p + \\pi^+ + \\pi^-$.) The other two processes are decays of unstable particles, and would happen even in a vacuum.\n", "\n", "Strictly speaking, there should be _two_ neutrinos in the decay of the muon - but we'll simplify by just considering one. The pion can also decay as $\\pi \\rightarrow e + \\nu$, but this is sufficiently rare that we can just ignore it.\n", "\n", "We can break our understanding of all of these processes into three smaller questions\"\n", "\n", "* What is the probability of a process happening, per unit time?\n", "* What does the process look like in the _rest frame_ of the initial particle?\n", "* Given the process in the rest frame, what does it look like in the _lab frame_ where we observe it?\n", "\n", "### Probability of decay/downscattering\n", "\n", "As we saw in detail when dealing with Monte Carlo simulations, for a decay process with lifetime $\\tau$, the probability that we observe a decay in time interval $dt$ is\n", "\n", "$$\n", "p(\\tau, dt) = 1 - e^{-dt/\\tau}.\n", "$$\n", "\n", "In our (laboratory) frame of reference, the three lifetimes we will need are:\n", "\n", "$$\n", "\\tau_p = \\exp \\left( \\frac{h}{7\\ {\\rm km}} \\right) \\times \\frac{1}{\\beta} \\times (1.14 \\times 10^{-5}\\ \\rm{s}) \\\\\n", "\\tau_{\\pi} = \\gamma \\times (2.603 \\times 10^{-8}\\ \\rm{s}) \\\\\n", "\\tau_{\\mu} = \\gamma \\times (2.197 \\times 10^{-6}\\ \\rm{s}) \\\\\n", "$$\n", "\n", "where $\\beta$ and $\\gamma$ are the usual relativistic factors,\n", "\n", "$$\n", "\\mathbf{\\beta} = \\frac{\\mathbf{v}}{c}, \\\\\n", "\\gamma = \\frac{1}{\\sqrt{1-v^2/c^2}}.\n", "$$\n", "\n", "Note that $\\mathbf{\\beta} = (\\beta_x, \\beta_y)$ is a vector since we're working in two dimensions. The variable $h$ is the height above the Earth's surface; the proton downscattering rate depends strongly on the density of the air, and thus on the height.\n", "\n", "The first process, $p \\rightarrow p + \\pi + \\pi$, arises from scattering off of air molecules; as a result, the time between scatters is inversely proportional to the speed of the proton, $\\beta$. Our model assumes a simple model for the density for air consistent (see appendix C). The other two processes are decays; in the rest frame of a $\\pi$ or $\\mu$ the lifetime is constant, but when they move at high speeds their apparent lifetime is extended due to __time dilation__.\n", "\n", "One more complication: the proton scattering requires the proton to exceed a minimum \"__threshold energy__\" for the interaction to take place at all. (A lot of energy is used up to convert to the mass of the two pions that are created.) Assuming the nucleus being recoiled from is very heavy, the threshold for this interaction is\n", "\n", "$$\n", "E_p \\geq (m_p + 2m_\\pi) c^2 \\approx 1218\\ \\rm{MeV}.\n", "$$\n", "\n", "__For a proton with less energy than the threshold, the downscattering _cannot occur_;__ in our model we will just assume that such a proton doesn't interact anymore.\n", "\n", "\n", "### Finding the decay/scattering products in the rest frame\n", "\n", "Let's start with decay (processes 2 and 3). In the rest frame of the mother particle, the system before and after decay looks like this:\n", "\n", "<img src=\"https://physicscourses.colorado.edu/phys2600/phys2600_fa18/img/proj1_decay.png\" width=400px />\n", "\n", "Conservation of momentum works for __any value__ of the angle $\\theta$; quantum mechanics predicts the probability of any particular $\\theta$. We will assume the decay is __isotropic__, which means that all $\\theta$s are equally likely. Once $\\theta$ is chosen, conservation of momentum tells us the rest - a quick calculation (see appendix) gives\n", "\n", "$$\n", "|\\mathbf{p}| = \\frac{c(M_X^2 - M_Y^2)}{2M_X}\n", "$$\n", "\n", "where $M_X$ is the mass of the mother particle, and $M_Y$ is the mass of the heavy daughter particle. (The other daughter particle is always a neutrino, so we take $M_Z = 0$.)\n", "\n", "Now the other process. This time, there are _three_ particles in the final state:\n", "\n", "<img src=\"https://physicscourses.colorado.edu/phys2600/phys2600_fa18/img/proj1_downscatter.png\" width=400px />\n", "\n", "\n", "In this case, it turns out that on top of the angles, the _momentum of each pion_ is random and determined by quantum effects. (It has to be, since momentum conservation alone doesn't tell us what the pion momentum should be!) A simple approximate model for the pion momentum distribution is\n", "\n", "$$\n", "p(E_\\pi) = b e^{-b(E_\\pi - m_\\pi c^2)}\n", "$$\n", "\n", "where the coefficient $b = 1/(500\\ \\rm{MeV})$, and $E_\\pi = \\sqrt{|\\mathbf{p}_\\pi|^2 c^2 + m_\\pi^2 c^4} \\geq m_\\pi c^2$. We can draw from this distribution by drawing $y$ from the uniform distribution over $[0,1]$, and then computing\n", "\n", "$$\n", "E_\\pi = m_\\pi c^2 - \\frac{1}{b} \\log (1-y),\n", "$$\n", "\n", "and then the pion momentum is\n", "\n", "$$\n", "p_\\pi = \\sqrt{E_\\pi^2 - m_\\pi^2 c^4} / c\n", "$$\n", "\n", "Once the two pion momenta are known, then momentum conservation tells us that\n", "\n", "$$\n", "\\mathbf{p}_p = -\\mathbf{p}_{\\pi,1} - \\mathbf{p}_{\\pi,2}.\n", "$$\n", "\n", "If you check the energy of the proton, you'll find that conservation of energy is violated! This is because we're ignoring the heavy nucleus that the proton is recoiling off of - so we don't expect energy conservation in the subsystem we're looking at.\n", "\n", "\n", "### Boosting back to the lab frame\n", "\n", "Once we have momentum vectors in the rest frame of an interaction, we have to apply a __Lorentz transformation__ (a \"boost\") to go back to the lab frame. In two dimensions, the lab components of the momentum vector are given by\n", "\n", "$$\n", "p_{x,\\rm{lab}} = p_x - \\gamma \\frac{E}{c} \\beta_x + (\\gamma - 1) \\beta_x \\frac{\\beta_x p_x + \\beta_y p_y}{\\beta^2} \\\\\n", "p_{y,\\rm{lab}} = p_y - \\gamma \\frac{E}{c} \\beta_y + (\\gamma - 1) \\beta_y \\frac{\\beta_x p_x + \\beta_y p_y}{\\beta^2}\n", "$$\n", "\n", "now defining two _directional_ velocity factors $(\\beta_x, \\beta_y) = (v_x / c, v_y / c)$, so that $\\beta = \\sqrt{\\beta_x^2 + \\beta_y^2}$. The right-hand side requires $p_x, p_y$, and relativistic energy $E$ in the rest frame of the decaying particle, and then $\\beta_x, \\beta_y, \\gamma$ are taken from the lab-frame velocity of the decaying particle.\n", "\n", "Finally, to get the speed from the momentum, starting from the relativistic formula for a massive particle\n", "\n", "$$\n", "\\mathbf{p} = \\gamma m \\mathbf{v} = \\frac{m\\mathbf{v}}{\\sqrt{1 - v^2/c^2}},\n", "$$\n", "\n", "we can show (see appendix A) that\n", "\n", "$$\n", "\\mathbf{v} = \\frac{\\mathbf{p}/m}{\\sqrt{1 + (p/mc)^2}}\n", "$$\n", "\n", "for a massive particle. (If we have a _massless_ particle like a neutrino, then its speed is always just $c$, pointing in the same direction as $\\mathbf{p}$. Your code should treat this as a special case!)" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "a6705cfbb6a1b95edb70264553743b2e", "grade": false, "grade_id": "cell-f201306b13816de8", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "# Computational Strategy and Algorithms\n", "\n", "There are a lot of parts to keep track of in a cosmic ray shower! The good news is that once we implement each small piece of physics above, we can use the magic of Monte Carlo simulation to chain them all together easily! Aside from the Monte Carlo, all we need are some basic kinematics (force-free motion) to keep track of particle positions. Here is the overall algorithm we want to use:\n", "\n", "1. Keep track of the following info for all particles: __identity__ (proton, muon, ...), current __position__, and current __velocity__.\n", "2. At each timestep $dt$, we use Monte Carlo to test for random decay/interaction, and split into multiple particles if such an event happens. (For the proton scattering only, we check to see if $E_p$ is above the threshold first.)\n", "3. Otherwise, we advance the positions of all particles as $(\\Delta x, \\Delta y) = (v_x dt, v_y dt)$.\n", "\n", "All of the complications happen in that second step. Whenever one of the three interaction processes from above happens, we find the outcome with the following sub-algorithm:\n", "\n", "1. Choose a random direction for one of the decay products (for the proton, choose two random directions and two random momenta for each of the two pions.)\n", "2. The momentum vectors of all decay products are then determined using conservation of energy and momentum.\n", "3. Using the known velocity of the decaying particle, apply a Lorentz boost to find the momentum vectors of all decay products in the lab frame.\n", "4. Convert from momentum to velocity, and update the list of particles.\n" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "80dea3f022aa4e8b16ca65600d8c9d07", "grade": false, "grade_id": "cell-f2ba414b33b68bf4", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "# Application Programming Interface (API) specification\n", "\n", "__You must implement all of the below functions, according to the specifications given.__ \n", "\n", "Two commonly-used tuple structures in our API are the __particle__ and the __decay product__:\n", "\n", "\n", "my_particle = (name, p_x, p_y, traj_x, traj_y)\n", "my_product = (name, p_x, p_y)\n", "\n", "\n", "Here name is a string, one of: 'proton', 'pion', 'muon', 'electron', or 'neutrino'. traj_x and traj_y are arrays of distances, corresponding to the history of a particle's position. Finally, p_x and p_y are the momentum components of a particle in some reference frame. \n", "\n", "(Why not just use \"particle\" tuples for everything? A design decision to avoid potential confusion. Product tuples will be created _outside_ the lab frame, and they won't have any position assigned to them until we do some work.)\n", "\n", "You can use __any set of units for this project__, and things will work _as long as you're self-consistent!_ I recommend \"particle physics units\", in which the speed of light is set to $c = 1$: this puts energy, momentum, and mass all in common units of energy. My code uses MeV for energy and s for time, which automatically requires that all distances are given in units of light-seconds (since $c$ is 1.)\n", "\n", "### Special relativity:\n", "* lorentz_beta_gamma(beta_x, beta_y): Computes and returns the dimensionless Lorentz factors (beta, gamma) from the components of velocity.\n", "* boost_momentum(m, px, py, beta_x, beta_y): Computes and returns (px_lab, py_lab), the components of momentum in the new frame of reference after a boost by $(\\beta_x, \\beta_y)$.\n", "* momentum_to_velocity(m, px, py): Returns the dimensionless velocity (beta_x, beta_y) corresponding to the given momentum and mass.\n", "\n", "### Particle interactions:\n", "* check_for_interaction(particle, dt, h): Checks to see if a particle will interact within time dt at height h above the Earth's surface. Returns True if a decay or downscatter occurs, and False otherwise. (As a convention, we assume that $h=0$ __corresponds to the Earth's surface.__)\n", "* particle_decay(name): Simulates the decay of a pion or muon in its rest frame. Returns a tuple (product_Y, product_Z) containing _decay product_ tuples for each decay product.\n", "* proton_downscatter(): Simulates the downscattering of a proton in its rest frame. Returns a tuple (product_pi_1, product_pi_2, product_p) containing _decay product_ tuples for each downscattering product.\n", "\n", "### Other functions:\n", "* boost_product_to_lab(product, beta_x, beta_y): Given a decay product tuple (name, m, px, py) (see above) and the two components of velocity (beta_x, beta_y) of the parent particle, returns the tuple (name, beta_x_lab, beta_y_lab) giving the velocity in the lab frame, i.e. after boosting the momentum by $(-\\beta_x, -\\beta_y)$ (_note the minus signs since we're boosting __to the lab frame__.)_\n", "* particle_mass(name): Given a particle name, return the corresponding mass in whatever units you have chosen.\n", "* velocity_to_momentum(m, beta_x, beta_y): Returns the momentum vector components (p_x, p_y) corresponding to the given dimensionless velocity and mass. Not used in the main loop, but useful for initializing particles with a given speed. (__Note:__ this function should give an error if $m$ is zero!)\n", "* run_API_tests(): A custom function that should use assertions to test the other functions implemented in the API. If all tests are passed, the function should simply return the None object. You should implement __at least six (6) tests__ inside this function on your own; additional tests will be provided after the checkpoint is due.\n", "\n", "## Main loop\n", "\n", "The below code implements the \"main loop\", running the Monte Carlo given an initial state, timestep, and number of timesteps to evolve. Once you have implemented the API above, add this code to your cosmic.py module, then call it to run the Monte Carlo simulation!" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "8d9488eb4eb7161dde3aabd939572d6a", "grade": false, "grade_id": "cell-f60d14d26a183edc", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "code", "source": [ "def run_cosmic_MC(particles, dt, Nsteps):\n", " \"\"\"\n", " Main loop for cosmic ray Monte Carlo project.\n", " \n", " Arguments:\n", " =====\n", " * particles: a list of any length, containing \"particle\" tuples.\n", " A particle tuple has the form:\n", " (name, p_x, p_y, traj_x, traj_y)\n", " where p_x and p_y are momentum vector components, and\n", " traj_x and traj_y are NumPy arrays of distance.\n", " \n", " * dt: time between Monte Carlo steps, in seconds.\n", " * Nsteps: Number of steps to run the Monte Carlo before returning.\n", " \n", " Returns:\n", " =====\n", " A list of particle tuples.\n", " \n", " Example usage:\n", " =====\n", " \n", " >> init_p_x, init_p_y = velocity_to_momentum(particle_mass('muon'), 0.85, -0.24) # beta ~ 0.8, relativistic\n", " >> init_particles = [ ('muon', init_p_x, init_p_y, np.array(), np.array([1e-4]) ) ] # ~ 30 km height in light-sec\n", " >> particles = run_cosmic_MC(init_particles, 1e-5, 100)\n", " \n", " The 'particles' variable after running should contain three particle tuples: the\n", " initial muon, an electron, and a neutrino. (Even though the muon decays,\n", " we keep it in the final particle list for its trajectory.)\n", " \n", " \"\"\"\n", " \n", " stopped_particles = []\n", " for step_i in range(Nsteps):\n", " updated_particles = []\n", " \n", " for particle in particles:\n", " # Unpack particle tuple\n", " (name, p_x, p_y, traj_x, traj_y) = particle\n", " (beta_x, beta_y) = momentum_to_velocity(particle_mass(name), p_x, p_y)\n", "\n", " # Check for interaction\n", " does_interact = check_for_interaction(particle, dt, traj_y[-1])\n", " \n", " if does_interact: \n", " if name == 'proton':\n", " decay_products = proton_downscatter()\n", " else:\n", " stopped_particles.append(particle)\n", " decay_products = particle_decay(name)\n", " \n", " # Transform products back to lab frame\n", " for product in decay_products:\n", " (product_name, product_p_x, product_p_y) = boost_product_to_lab(product, beta_x, beta_y)\n", " \n", " # If this was a proton scatter, then the \"new\" proton is\n", " # the same as the original, so keep track of its trajectory!\n", " if name == 'proton' and product_name == 'proton':\n", " product_traj_x = traj_x\n", " product_traj_y = traj_y\n", " else:\n", " product_traj_x = np.array([traj_x[-1]])\n", " product_traj_y = np.array([traj_y[-1]])\n", " \n", " # Make new particle tuple and append\n", " product_particle = (product_name, product_p_x, product_p_y, \n", " product_traj_x, product_traj_y)\n", " updated_particles.append( product_particle )\n", " else:\n", " # Doesn't interact, so compute motion\n", " traj_x = np.append(traj_x, traj_x[-1] + beta_x * dt)\n", " traj_y = np.append(traj_y, traj_y[-1] + beta_y * dt)\n", "\n", " updated_particles.append( (name, p_x, p_y, traj_x, traj_y) )\n", " \n", " \n", " # Run next timestep\n", " particles = updated_particles\n", " \n", " # Add stopped particles back to list and return\n", " particles.extend(stopped_particles)\n", " return particles" ], "outputs": [], "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "24f437310677c7b08413a08e049ab569", "grade": false, "grade_id": "cell-aaac1fc6dc7bf4fb", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "---\n", "---\n", "---\n", "\n", "The below appendices are included for your information, but you don't need to understand them in gory detail for the project - just use the resulting formulas as given above. (You may need some of them for the challenge question, though.)\n", "\n", "## Appendix A: Kinematics for decay/scattering\n", "\n", "Here I show the full derivations for the equations governing decay and scattering in our simplified model. Let's start with the simpler case of $X \\rightarrow Y + \\nu$ decay for processes 2 and 3 above. With only two particles in the final state, momentum conservation requires them to be back-to-back. Once the angle $\\theta$ is fixed (by randomly choosing it), the only quantity left to determine is the magnitude $p$ of the momentum vectors.\n", "\n", "For a relativistic system, the expression for energy is $E = \\sqrt{p^2c^2 + m^2c^4}$. The mass of the $\\nu$ is zero, so applying conservation of energy to the initial and final states, we find\n", "\n", "$$\n", "E_X = E_Y + E_Z \\\\\n", "m_X c^2 = \\sqrt{p^2 c^2 + m_Y^2 c^4} + pc\n", "$$\n", "\n", "which we can solve to find\n", "\n", "$$\n", "p = \\frac{c(m_X^2 - m_Y^2)}{2m_X}.\n", "$$\n", "\n", "\n", "Now, proton downscattering. This is actually nice and simple: since we are taking both the pion angles _and_ momenta to be determined randomly,\n", "\n", "$$\n", "p_{\\pi,i} = (p_\\pi \\cos \\theta_i, p_\\pi \\sin \\theta_i),\n", "$$\n", "\n", "then the proton momentum is immediately given by\n", "\n", "$$\n", "p_p = -p_{\\pi,1} - p_{\\pi,2} = (-p_\\pi [\\cos \\theta_1 + \\cos \\theta_2], -p_\\pi [\\sin \\theta_1 + \\sin \\theta_2]).\n", "$$\n", "\n", "Conservation of energy will be violated unless we consider the full 4-object process $p + X \\rightarrow p + X + \\pi + \\pi$ - but that's fine. The one assumption we need to add is that $X$ (a nucleus from the Earth's atmosphere) is _not_ traveling relativistically, so our incoming proton must have enough energy __in the lab frame__ to produce a pair of pions. This minimum is called the __threshold energy__. \n", "\n", "For $p \\rightarrow p + \\pi + \\pi$, the minimum possible energy occurs if all three products are at rest, which yields\n", "\n", "$$\n", "E_p \\geq (m_p + 2m_\\pi) c^2 \\approx 1218\\ \\rm{MeV}.\n", "$$\n", "\n", "(If we work this out carefully including the air nucleus we're scattering off of, we find the same result in the limit that the nucleus is very heavy.)" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "4d64c9989f68f897ac00ee550de5e5fa", "grade": false, "grade_id": "cell-ab84ffc94108b1c1", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "## Appendix B: Lorentz boosts in two spatial dimensions \n", "\n", "Most intro textbooks that discuss special relativity just show the effects of a Lorentz boost only along one of the coordinate axes. For this project, we need the more complicated expression for a coordinate change in an arbitrary direction. The mixing between coordinates can be described by a __boost matrix__ of the form\n", "\n", "$$\n", "\\left(\\begin{array}{c}ct' \\\\ x' \\\\ y'\\end{array} \\right) = \n", "\\left( \\begin{array}{ccc} \n", " \\gamma & -\\gamma \\beta_x & -\\gamma \\beta_y \\\\ \n", " -\\gamma \\beta_x & 1 + \\frac{\\gamma-1}{|\\beta|^2} \\beta_x^2 & \\frac{\\gamma-1}{|\\beta|^2} \\beta_x \\beta_y \\\\\n", " -\\gamma \\beta_y & \\frac{\\gamma-1}{|\\beta|^2} \\beta_x \\beta_y & 1 + \\frac{\\gamma-1}{|\\beta|^2} \\beta_y^2\n", "\\end{array} \\right) \\cdot \\left( \\begin{array}{c} ct \\\\ x \\\\ y \\end{array} \\right)\n", "$$" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "cccf5b628febaff0c8884b194ab330b7", "grade": false, "grade_id": "cell-2451b9ea5a1c09b6", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "Textbook version: for boost by speed $\\beta_x = v_x/c$ along the x-axis, we have\n", "\n", "$$\n", "t' = \\gamma(t - \\beta_x x) \\\\\n", "x' = \\gamma(x - \\beta_x ct) \\\\\n", "y' = y\n", "$$\n", "\n", "In general, for a boost by $\\boldsymbol{\\beta} = (\\beta_x, \\beta_y)$ in a general direction, we can write a vector version of these equations:\n", "\n", "$$\n", "t' = \\gamma(t - \\boldsymbol{\\beta} \\cdot \\mathbf{x}) \\\\\n", "\\mathbf{x}' = \\mathbf{x} - \\gamma ct \\boldsymbol{\\beta} + \\frac{\\gamma - 1}{|\\boldsymbol{\\beta}|^2} (\\boldsymbol{\\beta} \\cdot \\mathbf{x}) \\boldsymbol{\\beta}\n", "$$\n", "\n", "You can verify that this matches the \"boost matrix\" above. As a simple check, notice that this reduces to the textbook version if $\\boldsymbol{\\beta} = (\\beta_x, 0)$. For our purposes, we can either use the vector version, or we can write out the components to find:\n", "\n" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "cfcee45004df275446033a390590803b", "grade": false, "grade_id": "cell-bc701a22e47dd815", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "$$\n", "t' = \\gamma(t - \\beta_x x - \\beta_y y) \\\\\n", "x' = x - \\gamma c t \\beta_x + (\\gamma - 1) \\beta_x \\frac{\\beta_x x + \\beta_y y}{\\beta^2} \\\\\n", "y' = y - \\gamma c t \\beta_y + (\\gamma - 1) \\beta_y \\frac{\\beta_x x + \\beta_y y}{\\beta^2}\n", "$$\n", "\n", "Exactly the same equations hold for the energy-momentum vectors: just replace $(t, \\mathbf{x})$ with $(E, \\mathbf{p})$." ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "c4fc161fa4251d284c93163f30b4d8df", "grade": false, "grade_id": "cell-e81a415872c361e6", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "$$\n", "\\mathbf{p} = \\gamma m \\mathbf{v} = \\frac{m\\mathbf{v}}{\\sqrt{1 - v^2/c^2}},\n", "$$\n", "\n", "which we can invert first by noticing that\n", "\n", "$$\n", "p^2 = \\gamma^2 m^2 v^2 = \\gamma^2 m^2 c^2 (1 - \\frac{1}{\\gamma^2}) \\\\\n", "\\frac{p^2}{m^2 c^2} = \\gamma^2 - 1 \\\\\n", "\\gamma = \\sqrt{1 + (p/mc)^2}\n", "$$\n", "\n", "to get\n", "\n", "$$\n", "\\mathbf{v} = \\frac{\\mathbf{p}/m}{\\sqrt{1 + (p/mc)^2}}\n", "$$" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "f36384663aa48523882b3f09776fa06c", "grade": false, "grade_id": "cell-8b35cbd5976de3af", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "code", "source": [ "import numpy as np\n", "\n", "1.23 * np.exp(-15/7)\n", "NA = 6.022e23\n", "1.23 * NA / 28.9644e-3" ], "outputs": [], "execution_count": null, "metadata": {} }, { "cell_type": "markdown", "source": [ "## Appendix C: Finding the rate for proton downscattering\n", "\n", "In reality, the process we write as $p \\rightarrow p + \\pi + \\pi$ is really a two-to-four scattering process, $p + X \\rightarrow p + X + \\pi + \\pi$, where $X$ is a nucleus inside of an air molecule. The rate at which this will occur is described by a _cross section_, $\\sigma \\approx 114$ millibarn (very roughly, based on experiments on nitrogen nuclei and then corrected for the fact that we always scatter with enough energy to make pions.) A millibarn is equal to $10^{-31}$ square meters. The inverse scattering rate (average time to scatter) is related to the cross section as\n", "\n", "$$\n", "\\tau = 1/(\\sigma n v)\n", "$$\n", "\n", "where $v$ is the speed of the proton and $n$ is the number density of air. An important effect is that the density of air decreases rapidly further up from the Earth's surface. We will adopt [a simple isothermal model](https://www.spaceacademy.net.au/watch/debris/atmosmod.htm), in which the mass density of the atmosphere at height $h$ is equal to\n", "\n", "$$\n", "\\rho(h) \\approx 1.23\\ {\\rm kg/m}^3 \\exp \\left(-\\frac{h}{7\\ {\\rm km}} \\right)\n", "$$\n", "\n", "or translating to number density using $n = (N_A / M) \\rho$,\n", "$$\n", "n(h) \\approx 2.56 \\times 10^{25}\\ {\\rm m}^{-3} \\exp \\left(-\\frac{h}{7\\ {\\rm km}} \\right).\n", "$$\n", "\n", "Rewriting $v$ as $\\beta c$, we find that\n", "\n", "$$\n", "\\tau \\approx \\exp \\left( \\frac{h}{7\\ {\\rm km}} \\right) \\frac{1}{\\beta} 1.14 \\times 10^{-5}\\ \\rm{s}.\n", "$$\n", "\n", "Cosmic ray showers are most likely to be initiated in the stratosphere at an altitude of around 15 km.\n", "\n", "## Appendix D: Drawing random pion energies\n", "\n", "From quantum mechanics, an approximate model for the distribution of pions produced in proton downscattering gives the probability distribution\n", "\n", "$$\n", "p(E_\\pi) = b e^{-b(E_\\pi - m_\\pi c^2)}.\n", "$$\n", "\n", "To draw randomly from this distribution, we use the [inverse transform sampling](https://en.wikipedia.org/wiki/Inverse_transform_sampling) method. First, we integrate to find the cumulative distribution function:\n", "\n", "$$\n", "F(E_\\pi) = 1 - e^{-b(E_\\pi - m_\\pi c^2)}\n", "$$\n", "\n", "To draw randomly from this distribution, we set $y = F(E_\\pi)$ and invert to find $E_\\pi$ as a function of $y$:\n", "\n", "$$\n", "y = 1 - e^{-b(E_\\pi - m_\\pi c^2)} \\\\\n", "e^{-b(E_\\pi - m_\\pi c^2)} = 1 - y \\\\\n", "-b(E_\\pi - m_\\pi c^2) = \\log (1-y) \\\\\n", "E_\\pi = m_\\pi c^2 - \\frac{1}{b} \\log (1-y)\n", "$$\n", "\n", "Now drawing $y$ from the uniform distribution over $[0,1]$ gives the $E_\\pi$ distributed as we want. (Read the link above to see why this trick works, but the one-sentence explanation is that the CDF $F(E_\\pi)$ is also distributed over $[0,1]$ by definition.)" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "83471ee03bc30f530832003ce08b4896", "grade": false, "grade_id": "cell-697f10de55d07c71", "locked": true, "schema_version": 3, "solution": false } } } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.7.4", "mimetype": "text/x-python", "codemirror_mode": { "name": "ipython", "version": 3 }, "pygments_lexer": "ipython3", "nbconvert_exporter": "python", "file_extension": ".py" }, "toc": { "toc_position": {}, "skip_h1_title": false, "number_sections": true, "title_cell": "Table of Contents", "toc_window_display": true, "base_numbering": 1, "toc_section_display": true, "title_sidebar": "Contents", "toc_cell": false, "nav_menu": {}, "sideBar": true }, "nteract": { "version": "0.22.4" } }, "nbformat": 4, "nbformat_minor": 2 } __MACOSX/Final Proj/._cosmic_assign.ipynb Final Proj/.DS_Store __MACOSX/Final Proj/._.DS_Store Final Proj/cosmic_assign.pdf PHYS 2600 - Final Project 1: cosmic Monte Carlo Simulation of Cosmic Ray Showers Project Contact: LA Branden Esses (all LAs/Prof can answer questions about any project, but for fastest results, use the Project Contact) This is one possible final project; you must choose one (and only one) project to complete and hand in. The deadline for this project is 12:00 midnight (Boulder time) at the end of Thursday, 30 April. Since final grades must be submitted to the university shortly after the end of the final exam period, to allow for timely grading no late submissions will be accepted. You are permitted to discuss the final projects with other students in the class, including both the physics and computing aspects, but your submitted project must be entirely your own; no direct copying of Python code is permitted, and you must write up the presentation in your own words. When you submit your project, you must provide two (2) files: cosmic.py: a Python module (see lecture 22) which implements the functions listed in the Application Progamming Interface (API) below. cosmic_presentation.ipynb: a Jupyter notebook which uses the code from your Python module to answer the physics questions. Use of Markdown text, MathJax equations, and plots are all encouraged to make your presentation clear and compelling! The rubric for grading will be split 50/50 between the code (in your .py module and Jupyter notebook) and the response to the physics questions (provided in cosmic_presentation.ipynb): Basic functionality of code (passes provided tests, follows API specification): 30% Six (6) additional API tests you write (tests are provided and are correct and useful): 10% Documentation of code (docstrings and comments make the code clear): 10% Correct and detailed answers to physics questions: 40% Quality of presentation (clear and readable, appropriate use of Markdown, equations, and plots): 10% A bonus of up to 10% is available for the extra-credit "challenge" physics question at the end. (These challenge questions are meant to be really hard! But partial credit will be awarded if you make some progress on them.) Overview Every second, the Earth is bombarded by cosmic rays: high-energy particles flying at us from the cosmos. These naturally-accelerated particles have been a treasure trove for research, leading to discoveries of new elementary particles such as the positron and the muon in the early 20th century. The collisions of a cosmic proton with Earth's atmosphere lead to a spectacular multi-pronged trail of new particles called a cosmic ray shower. For this project, you will build a Monte Carlo simulation (in two dimensions to make things simpler) that will reproduce the cosmic-ray air shower effect. The set of interactions and particles you study will be limited to five: the proton, electron, pion, muon, and neutrino. You'll need to account for special relativity, which is a crucial ingredient since the cosmic ray particles are moving very fast! In addition to making nice simulated pictures of a spectacular physical phenomenon, you'll be able to directly investigate how the flux of cosmic rays in the upper atmosphere translates into detection of muons in ground-based observatories. Physics background and important equations For this project, we will adopt a simplified model of particle physics, which omits a lot of complicated details but includes enough physics to get the right qualitative behavior of cosmic ray showers. We will consider only the following set of particles, and we will ignore their electric charges completely: name symbol mass proton 938.3 MeV/c^2 electron 0.5110 MeV/c^2 pion 139.6 MeV/c^2 muon 105.7 MeV/c^2 neutrino 0 ("MeV" stands for "mega-electron-volt", a unit of energy; electron-volt based units are common in particle physics.) The proton and electron are familiar; the muon is a heavier (unstable) cousin of the electron, with the same charge. Neutrinos are ghostly particles with nearly zero mass and very weak interactions. Finally, the pion is a bound state of the strong nuclear force, like the proton. We're missing some details from the real Standard Model of particle physics, including the distinction between different types of neutrinos, the anti-particles for each of the above particles, and neutrons. But this is enough to get reasonable cosmic ray physics! We will consider only the following interactions in our model: 1. Proton downscattering: 2. Pion decay: 3. Muon decay: The first process is due to collisions of the cosmic ray proton with the nuclei inside of air molecules in the Earth's atmosphere; this inelastic collision produces a pair of pions (due to charge conservation: the final state is really .) The other two processes are decays of unstable particles, and would happen even in a vacuum. Strictly speaking, there should be two neutrinos in the decay of the muon - but we'll simplify by just considering one. The pion can also decay as , but this is sufficiently rare that we can just ignore it. We can break our understanding of all of these processes into three smaller questions" What is the probability of a process happening, per unit time? What does the process look like in the rest frame of the initial particle? Given the process in the rest frame, what does it look like in the lab frame where we observe it? Probability of decay/downscattering As we saw in detail when dealing with Monte Carlo simulations, for a decay process with lifetime , the probability that we observe a decay in time interval is In our (laboratory) frame of reference, the three lifetimes we will need are: where and are the usual relativistic factors, Note that is a vector since we're working in two dimensions. The variable is the height above the Earth's surface; the proton downscattering rate depends strongly on the density of the air, and thus on the height. The first process, , arises from scattering off of air molecules; as a result, the time between scatters is inversely proportional to the speed of the proton, . Our model assumes a simple model for the density for air consistent (see appendix C). The other two processes are decays; in the rest frame of a or the lifetime is constant, but when they move at high speeds their apparent lifetime is extended due to time dilation. One more complication: the proton scattering requires the proton to exceed a minimum "threshold energy" for the interaction to take place at all. (A lot of energy is used up to convert to the mass of the two pions that are created.) Assuming the nucleus being recoiled from is very heavy, the threshold for this interaction is For a proton with less energy than the threshold, the downscattering cannot occur; in our model we will just assume that such a proton doesn't interact anymore. Finding the decay/scattering products in the rest frame Let's start with decay (processes 2 and 3). In the rest frame of the mother particle, the system before and after decay looks like this: Conservation of momentum works for any value of the angle ; quantum mechanics predicts the probability of any particular . We will assume the decay is isotropic, which means that all s are equally likely. Once is chosen, conservation of momentum tells us the rest - a quick calculation (see appendix) gives where is the mass of the mother particle, and is the mass of the heavy daughter particle. (The other daughter particle is always a neutrino, so we take .) Now the other process. This time, there are three particles in the final state: In this case, it turns out that on top of the angles, the momentum of each pion is random and determined by quantum effects. (It has to be, since momentum conservation alone doesn't tell us what the pion momentum should be!) A simple approximate model for the pion momentum distribution is where the coefficient , and . We can draw from this distribution by drawing from the uniform distribution over , and then computing and then the pion momentum is Once the two pion momenta are known, then momentum conservation tells us that If you check the energy of the proton, you'll find that conservation of energy is violated! This is because we're ignoring the heavy nucleus that the proton is recoiling off of - so we don't expect energy conservation in the subsystem we're looking at. Boosting back to the lab frame Once we have momentum vectors in the rest frame of an interaction, we have to apply a Lorentz transformation (a "boost") to go back to the lab frame. In two dimensions, the lab components of the momentum vector are given by now defining two directional velocity factors , so that . The right-hand side requires , and relativistic energy in the rest frame of the decaying particle, and then are taken from the lab-frame velocity of the decaying particle. Finally, to get the speed from the momentum, starting from the relativistic formula for a massive particle we can show (see appendix A) that for a massive particle. (If we have a massless particle like a neutrino, then its speed is always just , pointing in the same direction as . Your code should treat this as a special case!) Computational Strategy and Algorithms There are a lot of parts to keep track of in a cosmic ray shower! The good news is that once we implement each small piece of physics above, we can use the magic of Monte Carlo simulation to chain them all together easily! Aside from the Monte Carlo, all we need are some basic kinematics (force-free motion) to keep track of particle positions. Here is the overall algorithm we want to use: 1. Keep track of the following info for all particles: identity (proton, muon, ...), current position, and current velocity. 2. At each timestep , we use Monte Carlo to test for random decay/interaction, and split into multiple particles if such an event happens. (For the proton scattering only, we check to see if is above the threshold first.) 3. Otherwise, we advance the positions of all particles as . All of the complications happen in that second step. Whenever one of the three interaction processes from above happens, we find the outcome with the following sub-algorithm: 1. Choose a random direction for one of the decay products (for the proton, choose two random directions and two random momenta for each of the two pions.) 2. The momentum vectors of all decay products are then determined using conservation of energy and momentum. 3. Using the known velocity of the decaying particle, apply a Lorentz boost to find the momentum vectors of all decay products in the lab frame. 4. Convert from momentum to velocity, and update the list of particles. Application Programming Interface (API) specification You must implement all of the below functions, according to the specifications given. Two commonly-used tuple structures in our API are the particle and the decay product: my_particle = (name, p_x, p_y, traj_x, traj_y) my_product = (name, p_x, p_y) Here name is a string, one of: 'proton', 'pion', 'muon', 'electron', or 'neutrino'. traj_x and traj_y are arrays of distances, corresponding to the history of a particle's position. Finally, p_x and p_y are the momentum components of a particle in some reference frame. (Why not just use "particle" tuples for everything? A design decision to avoid potential confusion. Product tuples will be created outside the lab frame, and they won't have any position assigned to them until we do some work.) You can use any set of units for this project, and things will work as long as you're self-consistent! I recommend "particle physics units", in which the speed of light is set to : this puts energy, momentum, and mass all in common units of energy. My code uses MeV for energy and s for time, which automatically requires that all distances are given in units of light-seconds (since is 1.) Special relativity: lorentz_beta_gamma(beta_x, beta_y): Computes and returns the dimensionless Lorentz factors (beta, gamma) from the components of velocity. boost_momentum(m, px, py, beta_x, beta_y): Computes and returns (px_lab, py_lab), the components of momentum in the new frame of reference after a boost by . momentum_to_velocity(m, px, py): Returns the dimensionless velocity (beta_x, beta_y) corresponding to the given momentum and mass. Particle interactions: check_for_interaction(particle, dt, h): Checks to see if a particle will interact within time dt at height h above the Earth's surface. Returns True if a decay or downscatter occurs, and False otherwise. (As a convention, we assume that corresponds to the Earth's surface.) particle_decay(name): Simulates the decay of a pion or muon in its rest frame. Returns a tuple (product_Y, product_Z) containing decay product tuples for each decay product. proton_downscatter(): Simulates the downscattering of a proton in its rest frame. Returns a tuple (product_pi_1, product_pi_2, product_p) containing decay product tuples for each downscattering product. Other functions: boost_product_to_lab(product, beta_x, beta_y): Given a decay product tuple (name, m, px, py) (see above) and the two components of velocity (beta_x, beta_y) of the parent particle, returns the tuple (name, beta_x_lab, beta_y_lab) giving the velocity in the lab frame, i.e. after boosting the momentum by (note the minus signs since we're boosting to the lab frame.) particle_mass(name): Given a particle name, return the corresponding mass in whatever units you have chosen. velocity_to_momentum(m, beta_x, beta_y): Returns the momentum vector components (p_x, p_y) corresponding to the given dimensionless velocity and mass. Not used in the main loop, but useful for initializing particles with a given speed. (Note: this function should give an error if is zero!) run_API_tests(): A custom function that should use assertions to test the other functions implemented in the API. If all tests are passed, the function should simply return the None object. You should implement at least six (6) tests inside this function on your own; additional tests will be provided after the checkpoint is due. Main loop The below code implements the "main loop", running the Monte Carlo given an initial state, timestep, and number of timesteps to evolve. Once you have implemented the API above, add this code to your cosmic.py module, then call it to run the Monte Carlo simulation! In [ ]: def run_cosmic_MC(particles, dt, Nsteps): """ Main loop for cosmic ray Monte Carlo project. Arguments: ===== * particles: a list of any length, containing "particle" tuples. A particle tuple has the form: (name, p_x, p_y, traj_x, traj_y) where p_x and p_y are momentum vector components, and traj_x and traj_y are NumPy arrays of distance. * dt: time between Monte Carlo steps, in seconds. * Nsteps: Number of steps to run the Monte Carlo before returning. Returns: ===== A list of particle tuples. Example usage: ===== >> init_p_x, init_p_y = velocity_to_momentum(particle_mass('muon'), 0.85, -0.24) # beta ~ 0.8, relativistic >> init_particles = [ ('muon', init_p_x, init_p_y, np.array(), np.array([1e-4]) ) ] # ~ 30 km height in light- sec >> particles = run_cosmic_MC(init_particles, 1e-5, 100) The 'particles' variable after running should contain three particle tuples: the initial muon, an electron, and a neutrino. (Even though the muon decays, we keep it in the final particle list for its trajectory.) """ stopped_particles = [] for step_i in range(Nsteps): updated_particles = [] for particle in particles: # Unpack particle tuple (name, p_x, p_y, traj_x, traj_y) = particle (beta_x, beta_y) = momentum_to_velocity(particle_mass(name), p_x, p_y) # Check for interaction does_interact = check_for_interaction(particle, dt, traj_y[-1]) if does_interact: if name == 'proton': decay_products = proton_downscatter() else: stopped_particles.append(particle) decay_products = particle_decay(name) # Transform products back to lab frame for product in decay_products: (product_name, product_p_x, product_p_y) = boost_product_to_lab(product, beta_x, beta_y) # If this was a proton scatter, then the "new" proton is # the same as the original, so keep track of its trajectory! if name == 'proton' and product_name == 'proton': product_traj_x = traj_x product_traj_y = traj_y else: product_traj_x = np.array([traj_x[-1]]) product_traj_y = np.array([traj_y[-1]]) # Make new particle tuple and append product_particle = (product_name, product_p_x, product_p_y, product_traj_x, product_traj_y) updated_particles.append( product_particle ) else: # Doesn't interact, so compute motion traj_x = np.append(traj_x, traj_x[-1] + beta_x * dt) traj_y = np.append(traj_y, traj_y[-1] + beta_y * dt) updated_particles.append( (name, p_x, p_y, traj_x, traj_y) ) # Run next timestep particles = updated_particles # Add stopped particles back to list and return particles.extend(stopped_particles) return particles The below appendices are included for your information, but you don't need to understand them in gory detail for the project - just use the resulting formulas as given above. (You may need some of them for the challenge question, though.) Appendix A: Kinematics for decay/scattering Here I show the full derivations for the equations governing decay and scattering in our simplified model. Let's start with the simpler case of decay for processes 2 and 3 above. With only two particles in the final state, momentum conservation requires them to be back-to-back. Once the angle is fixed (by randomly choosing it), the only quantity left to determine is the magnitude of the momentum vectors. For a relativistic system, the expression for energy is . The mass of the is zero, so applying conservation of energy to the initial and final states, we find which we can solve to find Now, proton downscattering. This is actually nice and simple: since we are taking both the pion angles and momenta to be determined randomly, then the proton momentum is immediately given by Conservation of energy will be violated unless we consider the full 4-object process - but that's fine. The one assumption we need to add is that (a nucleus from the Earth's atmosphere) is not traveling relativistically, so our incoming proton must have enough energy in the lab frame to produce a pair of pions. This minimum is called the threshold energy. For , the minimum possible energy occurs if all three products are at rest, which yields (If we work this out carefully including the air nucleus we're scattering off of, we find the same result in the limit that the nucleus is very heavy.) Appendix B: Lorentz boosts in two spatial dimensions Most intro textbooks that discuss special relativity just show the effects of a Lorentz boost only along one of the coordinate axes. For this project, we need the more complicated expression for a coordinate change in an arbitrary direction. The mixing between coordinates can be described by a boost matrix of the form Textbook version: for boost by speed along the x-axis, we have In general, for a boost by in a general direction, we can write a vector version of these equations: You can verify that this matches the "boost matrix" above. As a simple check, notice that this reduces to the textbook version if . For our purposes, we can either use the vector version, or we can write out the components to find: Exactly the same equations hold for the energy-momentum vectors: just replace with . which we can invert first by noticing that to get In [ ]: import numpy as np 1.23 * np.exp(-15/7) NA = 6.022e23 1.23 * NA / 28.9644e-3 Appendix C: Finding the rate for proton downscattering In reality, the process we write as is really a two-to-four scattering process, , where is a nucleus inside of an air molecule. The rate at which this will occur is described by a cross section, millibarn (very roughly, based on experiments on nitrogen nuclei and then corrected for the fact that we always scatter with enough energy to make pions.) A millibarn is equal to square meters. The inverse scattering rate (average time to scatter) is related to the cross section as where is the speed of the proton and is the number density of air. An important effect is that the density of air decreases rapidly further up from the Earth's surface. We will adopt a simple isothermal model, in which the mass density of the atmosphere at height is equal to or translating to number density using , Rewriting as , we find that Cosmic ray showers are most likely to be initiated in the stratosphere at an altitude of around 15 km. Appendix D: Drawing random pion energies From quantum mechanics, an approximate model for the distribution of pions produced in proton downscattering gives the probability distribution To draw randomly from this distribution, we use the inverse transform sampling method. First, we integrate to find the cumulative distribution function: To draw randomly from this distribution, we set and invert to find as a function of : Now drawing from the uniform distribution over gives the distributed as we want. (Read the link above to see why this trick works, but the one- sentence explanation is that the CDF is also distributed over by definition.) p e π μ ν ∼ p → p + π + π π → μ + ν μ → e + ν p + +π+ π− π → e + ν τ dt p(τ, dt) = 1 − .e−dt/τ = exp( ) × × (1.14 × s)τp h7 km 1β 10−5 = γ × (2.603 × s)τπ 10−8 = γ × (2.197 × s)τμ 10−6 β γ β = ,v c γ = .1 1 − /v2 c2‾ ‾‾‾‾‾‾‾‾√ β = ( , )βx βy h p → p + π + π β π μ ≥ ( + 2 ) ≈ 1218 MeV.Ep mp mπ c2 θ θ θ θ |p| = c( − )M 2X M 2Y 2MX MX MY = 0MZ p( ) = bEπ e−b( − )Eπ mπc 2 b = 1/(500 MeV) = ≥Eπ | +pπ |2 c2 m2π c4‾ ‾‾‾‾‾‾‾‾‾‾‾‾‾√ mπ c2 y [0, 1] = − log(1 − y),Eπ mπ c2 1 b = /cpπ −E2π m2π c4‾ ‾‾‾‾‾‾‾‾‾√ = − − .pp pπ,1 pπ,2 = − γ + (γ − 1)px,lab px E c βx βx +βxpx βypy β2 = − γ + (γ − 1)py,lab py E c βy βy +βxpx βypy β2 ( , ) = ( /c, /c)βx βy vx vy β = +β2x β2y‾ ‾‾‾‾‾‾√ ,px py E , , γβx βy p = γmv = ,mv 1 − /v2 c2‾ ‾‾‾‾‾‾‾‾√ v = p/m 1 + (p/mc)2‾ ‾‾‾‾‾‾‾‾‾‾√ c p dt Ep (Δx, Δy) = ( dt, dt)vx vy c = 1 c ( , )βx βy h = 0 (− , − )βx βy m X → Y + ν θ p E = +p2 c2 m2c4‾ ‾‾‾‾‾‾‾‾‾‾√ ν = +EX EY EZ = + pcmXc2 +p2 c2 m2Y c 4‾ ‾‾‾‾‾‾‾‾‾‾‾√ p = . c( − )m2X m2Y 2mX = ( cos , sin ),pπ,i pπ θi pπ θi = − − = (− [cos + cos ], − [sin + sin ]).pp pπ,1 pπ,2 pπ θ1 θ2 pπ θ1 θ2 p + X → p + X + π + π X p → p + π + π ≥ ( + 2 ) ≈ 1218 MeV.Ep mp mπ c2 = ⋅ ⎛ ⎝ ⎜⎜ ct ′ x′ y′ ⎞ ⎠ ⎟⎟ ⎛ ⎝ ⎜⎜⎜⎜ γ −γβx −γβy −γβx 1 + γ−1 |β|2 β2x γ−1 |β|2 βxβy −γβy γ−1 |β|2 βxβy 1 + γ−1 |β|2 β2y ⎞ ⎠ ⎟⎟⎟⎟ ⎛ ⎝ ⎜⎜ ct x y ⎞ ⎠ ⎟⎟ = /cβx vx = γ(t − x)t ′ βx = γ(x − ct)x′ βx = yy′ β = ( , )βx βy = γ(t − β ⋅ x)t ′ = x − γctβ + (β ⋅ x)βx′ γ − 1 |β|2 β = ( , 0)βx = γ(t − x − y)t ′ βx βy = x − γct + (γ − 1)x′ βx βx x + yβx βy β2 = y − γct + (γ − 1)y′ βy βy x + yβx βy β2 (t, x) (E, p) p = γmv = ,mv 1 − /v2 c2‾ ‾‾‾‾‾‾‾‾√ = = (1 − )p2 γ2 m2v2 γ2 m2c2 1 γ2 = − 1 p2 m2c2 γ2 γ = 1 + (p/mc)2‾ ‾‾‾‾‾‾‾‾‾‾√ v = p/m 1 + (p/mc)2‾ ‾‾‾‾‾‾‾‾‾‾√ p → p + π + π p + X → p + X + π + π X σ ≈ 114 10−31 τ = 1/(σnv) v n h ρ(h) ≈ 1.23 exp(− )kg/m3 h7 km n = ( /M)ρNA n(h) ≈ 2.56 × exp(− ).1025 m−3 h7 km v βc τ ≈ exp( ) 1.14 × s.h7 km 1β 10−5 p( ) = b .Eπ e−b( − )Eπ mπc 2 F( ) = 1 −Eπ e−b( − )Eπ mπc 2 y = F( )Eπ Eπ y y = 1 − e−b( − )Eπ mπc2 = 1 − ye−b( − )Eπ mπc2 − b( − ) = log(1 − y)Eπ mπ c2 = − log(1 − y)Eπ mπ c2 1 b y [0, 1] Eπ F( )Eπ [0, 1] __MACOSX/Final Proj/._cosmic_assign.pdf Final Proj/cosmic.py import numpy as np __MACOSX/Final Proj/._cosmic.py Final Proj/cosmic_presentation.ipynb { "cells": [ { "cell_type": "markdown", "source": [ "# Presentation Notebook - Final Project 1: cosmic\n", "\n", "This notebook contains physics questions to be answered as part of the cosmic final project - see companion assignment notebook for details." ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "1343f7093bcaee2f49538c1c93cbd91f", "grade": false, "grade_id": "cell-143e07fdb29fa8ab", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "_Type your answer here using Markdown._" ], "metadata": { "deletable": false, "nbgrader": { "cell_type": "markdown", "checksum": "273c3ad184d2c9a2c5b2119b01ce8ec3", "grade": true, "grade_id": "cell-a257c56d54421841", "locked": false, "points": 20, "schema_version": 3, "solution": true } } }, { "cell_type": "code", "source": [ "# Import cell\n", "\n", "from cosmic import *" ], "outputs": [], "execution_count": null, "metadata": {} }, { "cell_type": "markdown", "source": [ "# Physics Questions\n", "\n", "Add cells below as necessary to answer each of the physics questions.\n", "\n", "__Make sure your notebook runs with no errors by using the \"Restart & Run All\" command before you submit!__" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "81a4f7d861f6f32c20da5a197dede896", "grade": false, "grade_id": "cell-32a242f48cb1212c", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "### 1. Testing with simple showers\n", "\n", "As a simple test of your Monte Carlo code, find the \"shower\" induced by an initial electron, initial muon, and initial pion, each traveling at 0.999c towards the ground. (None of these will produce a complicated shower; the electron doesn't interact, and the other interactions can only happen once.) Make a plot or plots showing the particle trajectories, and explain them.\n", "\n", "Once you think your plots make sense, check conservation of energy and momentum for the muon decay (not true for downscatters since the air nuclei provide/absorb energy!)" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "e59450e78618ba2190a8f223de305f4b", "grade": false, "grade_id": "cell-654b697fc5411f33", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "### 2. Simulating a single proton shower\n", "\n", "Simulate a cosmic-ray shower with a single incoming proton with initial velocity $(\\beta_x, \\beta_y) = (0.95999, -0.28)$, and starting at $(x,y) = (0, 20)$ km. Make a two-dimensional figure showing the trajectories of all particles produced in the shower, and use the figure to explain cosmic rays at a level that a freshman physics major would understand. (How does it compare to a real cosmic ray shower?)" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "e129d1493da37ecaca34381678046339", "grade": false, "grade_id": "cell-403ceee6eee7f4da", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "### 3. Ground flux of cosmic muons\n", "\n", "For protons with an initial total energy of 30 GeV, the flux of cosmic rays (rate at which they bombard the Earth's atmosphere) is approximately equal to\n", "\n", "$$\n", "I(\\theta) = (1.0\\ \\rm{m}^{-2} \\rm{s}^{-1}) \\times \\cos^2(\\theta)\n", "$$\n", "\n", "where $\\theta$ is the angle from the vertical. Run your Monte Carlo simulation several times for protons incident on the upper atmosphere at an altitude of 20 km, incoming towards the point $(0,h)$ on the surface at height $h$ above sea level at an angle $\\theta$. Based on how many muons cross through ground level in your simulation results, what is the total flux of muons (number of muons observed per second) at ground level?\n", "\n", "If you had a detector which was 0.1 square meters in size, how long would you have to wait per muon on average at sea level ($h=0$)?" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "f63ecb99b761eacfecbeff07eac6b931", "grade": false, "grade_id": "cell-489d4fffab529961", "locked": true, "schema_version": 3, "solution": false } } }, { "cell_type": "markdown", "source": [ "### 4. [Extra-credit challenge!] Cosmic muons in underground laboratories\n", "\n", "Cosmic ray muons are a nuisance if you're trying to do a very sensitive search for rare particles, like hunting for dark matter. To shield against these muons, most dark matter detectors are placed deep underground. When traveling through the solid rock that makes up the Earth's crust, muons suffer energy losses [as described in section 30.4.1 of this document](http://pdg.lbl.gov/2017/reviews/rpp2017-rev-cosmic-rays.pdf) (don't worry about understanding the rest of the document!) Implement energy loss for muons passing below $y=0$ in your Monte Carlo simulation. How far underground would a detector have to be for the muon detection rate (in your model) to pass below 1 muon per year?" ], "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "d70a8a5ab0664639135287c13a4032db", "grade": false, "grade_id": "cell-0947c7f3cf39c562", "locked": true, "schema_version": 3, "solution": false } } } ], "metadata": { "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.8" }, "toc": { "toc_position": {}, "skip_h1_title": false, "number_sections": false, "title_cell": "Table of Contents", "toc_window_display": true, "base_numbering": 1, "toc_section_display": true, "title_sidebar": "Contents", "toc_cell": false, "nav_menu": {}, "sideBar": true }, "nteract": { "version": "0.22.4" } }, "nbformat": 4, "nbformat_minor": 2 } __MACOSX/Final Proj/._cosmic_presentation.ipynb Final Proj/cosmic_presentation.pdf Presentation Notebook - Final Project 1: cosmic This notebook contains physics questions to be answered as part of the cosmic final project - see companion assignment notebook for details. Honor Code Pledge Please type your name in the box below to agree to the Honor Code Pledge for this assignment: "On my honor, as a University of Colorado Boulder student, I have neither given nor received unauthorized assistance." Type your answer here using Markdown. In [ ]: # Import cell from cosmic import * Physics Questions Add cells below as necessary to answer each of the physics questions. Make sure your notebook runs with no errors by using the "Restart & Run All" command before you submit! 1. Testing with simple showers As a simple test of your Monte Carlo code, find the "shower" induced by an initial electron, initial muon, and initial pion, each traveling at 0.999c towards the ground. (None of these will produce a complicated shower; the electron doesn't interact, and the other interactions can only happen once.) Make a plot or plots showing the particle trajectories, and explain them. Once you think your plots make sense, check conservation of energy and momentum for the muon decay (not true for downscatters since the air nuclei provide/absorb energy!) 2. Simulating a single proton shower Simulate a cosmic-ray shower with a single incoming proton with initial velocity , and starting at km. Make a two- dimensional figure showing the trajectories of all particles produced in the shower, and use the figure to explain cosmic rays at a level that a freshman physics major would understand. (How does it compare to a real cosmic ray shower?) 3. Ground flux of cosmic muons For protons with an initial total energy of 30 GeV, the flux of cosmic rays (rate at which they bombard the Earth's atmosphere) is approximately equal to where is the angle from the vertical. Run your Monte Carlo simulation several times for protons incident on the upper atmosphere at an altitude of 20 km, incoming towards the point on the surface at height above sea level at an angle . Based on how many muons cross through ground level in your simulation results, what is the total flux of muons (number of muons observed per second) at ground level? If you had a detector which was 0.1 square meters in size, how long would you have to wait per muon on average at sea level ( )? 4. [Extra-credit challenge!] Cosmic muons in underground laboratories Cosmic ray muons are a nuisance if you're trying to do a very sensitive search for rare particles, like hunting for dark matter. To shield against these muons, most dark matter detectors are placed deep underground. When traveling through the solid rock that makes up the Earth's crust, muons suffer energy losses as described in section 30.4.1 of this document (don't worry about understanding the rest of the document!) Implement energy loss for muons passing below in your Monte Carlo simulation. How far underground would a detector have to be for the muon detection rate (in your model) to pass below 1 muon per year? ( , ) = (0.95999, −0.28)βx βy (x, y) = (0, 20) I(θ) = (1.0 ) × (θ)m−2s−1 cos2 θ (0, h) h θ h = 0 y = 0 __MACOSX/Final Proj/._cosmic_presentation.pdf Final Proj/cosmic_API_tests_full.pdf In [ ]: %matplotlib inline from cosmic import * import matplotlib.pyplot as plt 1. lorentz_beta_gamma Verify that lorentz_beta_gamma computes the correct values of and . Test input: (beta_x, beta_y) = (0.3, 0.4) Expected output: In [ ]: test_beta_xy = (0.30, 0.40) test_beta, test_gamma = lorentz_beta_gamma(*test_beta_xy) print(test_beta, test_gamma) assert np.abs(test_beta - 0.5) <= 1e-6 assert np.abs(test_gamma - 1.15470) <= 1e-6 2. boost_momentum Using the test beta/gamma from above, boost a known momentum to another frame of reference. Test input: (p_x, p_y) = (40, -40) MeV/c, for a pion (mass m = 139.6 MeV/c**2.) Same (beta_x, beta_y) = (0.3, 0.4) from above. Expected output: (calculated by hand) Energy: MeV In [ ]: test_p_xy = (40., -40.) test_p = np.sqrt(test_p_xy**2 + test_p_xy**2) test_m = 139.6 test_E = np.sqrt(test_p**2 + test_m**2 ) print(test_E) px_boost, py_boost = boost_momentum(test_m, test_p_xy, test_p_xy, test_beta_xy, test_beta_xy) print([px_boost, py_boost] ) assert np.abs(test_E - 150.63)/150.63 <= 1e-4 assert np.abs(px_boost + 12.92)/12.92 <= 1e-4 assert np.abs(py_boost + 110.56)/110.56 <= 1e-4 3. check_for_interaction Several test cases for interactions with very high/low probability, so the outcome should always be True/False. In [ ]: ## Unyt version h = 15 * 1000 / 3e8 # km --> light-seconds test_p1 = velocity_to_momentum(particle_mass('electron'), 0.3, 0.4) test_part_1 = ('electron', test_p1, test_p1, [], []) print("Electron interacts: (must be False)") print(check_for_interaction(test_part_1, 1e6, h)) test_p2 = velocity_to_momentum(particle_mass('muon'), 0.3, 0.4) test_part_2 = ('muon', test_p2, test_p2, [], []) print("Muon interacts in 1e-10 seconds: (False)") print(check_for_interaction(test_part_2, 1e-10, h)) print("Muon interacts in 100 seconds: (True)") print(check_for_interaction(test_part_2, 10, h)) test_p3 = velocity_to_momentum(particle_mass('proton'), 0.7, 0.4) test_part_3 = ('proton', test_p3, test_p3, [], []) print("Proton interacts in 1 second: (True)") print(check_for_interaction(test_part_3, 1 , h)) test_p4 = velocity_to_momentum(particle_mass('proton'), 0.3, 0.2) test_part_4 = ('proton', test_p4, test_p4, [], []) print("Slow proton doesn't interact (below threshold): (False)") print(check_for_interaction(test_part_4, 1, h )) 4. proton_downscatter Verify that proton downscattering returns decay-product tuples, and that total momentum is zero in the proton's rest frame. In [ ]: part1, part2, part3 = proton_downscatter() print(part1,'\n', part2, '\n', part3) # Should be a proton and two pions as decay product tuples # Check momentum conservation: should print 0 twice print(part1 + part2 + part3) print(part1 + part2 + part3) 5. Main loop test In [ ]: import unyt as u print(1*u.m) # Define light-second distance u.define_unit('ls', 1*u.s * u.physical_constants.c) (1*u.ls).to('km') In [ ]: # Proton shower px, py = velocity_to_momentum(particle_mass('proton'), 0.959999, -0.28) particles = [ ('proton', px, py, ( * u.m).to_value('ls'), ([20e3] * u.m).to_value('ls')), ] particles = run_cosmic_MC(particles, 1e-5, 50) pc = { 'proton': 'black', 'pion': 'red', 'muon': 'green', 'electron': 'blue', 'neutrino': 'gold', } #print(particles) for particle in particles: plt.plot(particle, particle, color=pc[particle]) In [ ]: # Muon shower px, py = velocity_to_momentum(particle_mass('muon'), 0.92, -0.24) init_particles = [ ('muon', px, py, np.array(), np.array([1e-4]) ) ] # ~ 30 km height in light-sec particles = run_cosmic_MC(init_particles, 2e-7, 100) for particle in particles: plt.plot(particle, particle, color=pc[particle]) β γ β = 0.5 γ = 1/ = 2/ ≈ 1.15471 − β2‾ ‾‾‾‾‾√ 3‾√ E = ≈ 150.63+m2 p2‾ ‾‾‾‾‾‾‾√ ( , ) = (−12.92, −110.56)p′x p′y __MACOSX/Final Proj/._cosmic_API_tests_full.pdf

## Related Questions

Similar orders to Simulating Cosmic Ray Showers using Monte Carlo in Python 3
12
Views
0
Lab 4 - Hooks law and springs
Please follow the instructions in the document below. All of it must be done exactly how it is written. In the part 1, there needs to be a data table with the data that will be gathered. No need to do part IV, V or VI. In the summary provide what have ...
75
Views
1
Physics I for Engineering and Science
3 questions will be assigned Dec 21st @4pm EST. I will send to you at the specified time and need worked sent back to me before 6pm EST. I do not know the question specifically. However, I do know what the 3 question will be about. The questions need to be...
38
Views
0
Physics( semiconductors ) with 4 homework questions
I need the solutions of these four questions which i have attached...
44
Views
0 