联系方式

  • QQ:99515681
  • 邮箱:99515681@qq.com
  • 工作时间:8:00-23:00
  • 微信:codehelp

您当前位置:首页 >> R语言程序R语言程序

日期:2020-05-01 09:14

2020/4/28 cosmic_assign
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 24) 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 p 938.3 MeV/c
2
electron e 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.
(For the kinematics of the proton's interaction, you will also need the mass of the nitrogen-14 nucleus, which is approximately MN
≈ 13040 MeV/ c
2
in
the same units.)
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!
We will consider only the following interactions in our model:
2020/4/28 cosmic_assign
https://hepjyp01.colorado.edu/user/zizo9379/notebooks/final_cosmic_fix/cosmic_assign.ipynb 2/7
1. Proton downscattering: p → p + π + π
2. Pion decay: π → μ + ν
3. Muon decay: μ → e + ν
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 + π
+ + π
?
.) 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
π → e + ν, 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 dt is
p(τ, dt) = 1 ? e? dt /τ.
In our (laboratory) frame of reference, the three lifetimes we will need are:
τp = exph7 km×1β× (1.14 × 10 ? 5s)τπ = γ × (2.603 × 10 ? 8s)τμ = γ × (2.197 × 10 ? 6s)
where β and γ are the usual relativistic factors,
β =vc, γ =1√1 ? v2/ c2.
Note that β = (βx, β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.
The first process, p → p + π + π, 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 (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
Ep ≥ (mp+ 2mπ)c2 ≈ 1218 MeV.
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
| p | =c(M2X? M2Y)2MX
where MX
is the mass of the mother particle, and MY
is the mass of the heavy daughter particle. (The other daughter particle is always a neutrino, so
we take MZ = 0.)
Now the other process. To avoid running into problems with energy conservation, we need to include the heavy nitrogen-14 nucleus in the scattering.
We take it to be at rest in the lab frame, which means it will be incoming at high speed when we switch to the proton rest frame:
( )
2020/4/28 cosmic_assign
https://hepjyp01.colorado.edu/user/zizo9379/notebooks/final_cosmic_fix/cosmic_assign.ipynb 3/7
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
p(Eπ) = be? b (Eπ? mπc2)
where the coefficient b = 1/(500 MeV), and Eπ = | pπ|2c2 + m2πc4 ≥ mπc2. We can draw from this distribution by drawing y from the uniform
distribution over [0, 1], and then computing,
and then the pion momentum is
The creation of two pions from nothing requires a large amount of energy, at least mπc2
for each pion. Assuming the extra energy of the recoiling
proton is small by comparison, we have the relation
are the energy of the nitrogen-14 nucleus before and after the collision, respectively. We can use this to compute the magnitude of
the nitrogen-14 momentum after collision as
Making one more assumption, that the direction of →p′N
is the same as the direction of the incident momentum pN
(which is true in the limit that the
nitrogen-14 mass is very heavy), we can compute the components of p′N
by rescaling them using the ratio of the magnitudes | p′N| / | pN| . Then we
can determine the proton's momentum by adding up initial and final states:
pp = pN? p′N? pπ , 1 ? pπ , 2.
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
px , lab = px ? γEc
βx + (γ ? 1)βxβxpx + βypyβ2
py , lab = py ? γ
now defining two directional velocity factors (βx, βy) = (vx/ c, vy/ c), so that β = β2x + β2y. The right-hand side requires px, py, and relativistic energy E in
the rest frame of the decaying particle, and then βx, βy, γ 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 c, pointing in the same direction as p. 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 (forcefree
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.
2020/4/28 cosmic_assign
https://hepjyp01.colorado.edu/user/zizo9379/notebooks/final_cosmic_fix/cosmic_assign.ipynb 4/7
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 Ep
is above the threshold first.)
3. Otherwise, we advance the positions of all particles as (Δx, Δy) = (v
xdt, vydt).
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, rest_mass, 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. rest_mass is the rest mass of the particle. 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 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.)
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 (βx
, βy
).
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 h = 0 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(beta_x , beta_y) : Simulates the downscattering of a proton in its rest frame. The inputs are the current components of
relativistic velocity of the proton, which are used to determine the momentum of the recoiling nitrogen-14 nucleus in the proton 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, p_x_lab, p_y_lab) giving the product's momentum in the lab frame,
i.e. after boosting the momentum by ( ? βx
, ? βy
) (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 m 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!
2020/4/28 cosmic_assign
https://hepjyp01.colorado.edu/user/zizo9379/notebooks/final_cosmic_fix/cosmic_assign.ipynb 5/7
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([0]), 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(beta_x, beta_y)
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
2020/4/28 cosmic_assign
https://hepjyp01.colorado.edu/user/zizo9379/notebooks/final_cosmic_fix/cosmic_assign.ipynb 6/7
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 X → Y + ν
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 p of the momentum vectors.
For a relativistic system, the expression for energy is E = √p. The mass of the ν is zero, so applying conservation of energy to the initial and
final states, we find
EX = EY + EZmX
which we can solve to find.
Now, proton downscattering. The details for these were laid out above, but an important assumption we made is that the nucleus scattered off of in
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 p → p + π + π, the minimum possible energy occurs if all three products are at rest, which yields
Ep ≥ (mp
+ 2mπ
)c
2 ≈ 1218 MeV.
(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;
nitrogen-14 is heavy enough that we don't have to modify this estimate for our purposes.)

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
You can verify that this matches the "boost matrix" above. As a simple check, notice that this reduces to the textbook version if β = (βx
, 0). 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 (t, x) with (E, p).
which we can invert first by noticing that
2020/4/28 cosmic_assign
https://hepjyp01.colorado.edu/user/zizo9379/notebooks/final_cosmic_fix/cosmic_assign.ipynb 7/7
Appendix C: Finding the rate for proton downscattering
In reality, the process we write as p → p + π + π is really a two-to-four scattering process, p + X → p + X + π + π, where X is a nucleus inside of an air
molecule. The rate at which this will occur is described by a cross section, σ ≈ 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
τ = 1/(σnv)
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
ρ(h) ≈ 1.23 kg/m3exp ?h7 km
or translating to number density using n = (NA
/M)ρ,n(h) ≈ 2.56 × 1025 m? 3exp ?h7 km.
Rewriting v as βc, 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 (https://en.wikipedia.org/wiki/Inverse_transform_sampling) method.
First, we integrate to find the cumulative distribution function:
To draw randomly from this distribution, we set y = F(Eπ) and invert to find Eπ
as a function of y:
Now drawing y from the uniform distribution over [0, 1] gives the Eπ
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π) is also distributed over [0, 1] by definition.)

版权所有:留学生编程辅导网 2021,All Rights Reserved 联系方式:QQ:99515681 电子信箱:99515681@qq.com
免责声明:本站部分内容从网络整理而来,只供参考!如有版权问题可联系本站删除。