PHYS 2600代做、Python编程语言调试、代做Python、代写data

日期：2020-04-30 11:22

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.)