Instructions: Create a new notebook called MonteCarloExercises in your MonteCarlo directory and solve the following problems inside it. Be sure to include the problem statements in a markdown cell above your solution. You don't need to put the "helper" code in the markdown cell, just implement the helper code in your code cell with your solution.
Preliminaries: At the top of your notebook, include a "Heading 1" cell with the title Monte Carlo Sampling and Simulation Exercises. Then include the inline functions and libraries by adding a code cell that invokes the %pylab inline magic and imports the needed packages.
The Weibull distribution, with parameters $\alpha > 0$ and $\beta > 0$ is described by the density function
$$f(x) = \begin{cases} \alpha\beta^{-\alpha}x^{\alpha -1}e^{-(x/\beta)^{\alpha}} & x > 0 \\ 0 & \text{otherwise} \end{cases}$$(a) Solve for the CDF of the Weibull distribution analytically, using
$$ F(x) = \int_{-\infty}^x f(t)dt$$Write the solution as a $\LaTeX$ formatted equation in a markdown cell.
In [ ]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a MARKDOWN cell
(b) Find the inverse of the CDF $F^{-1}(x)$ analytically. Write the solution as a $\LaTeX$ formatted equation in a markdown cell.
In [ ]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a MARKDOWN cell
(c) Write a program to generate 10,000 samples from the Weibull distribution and plot the true distribution, the CDF and the sampled distribution as a histogram, all on the same graph. Let $\alpha$ = 1.5, $\beta$ = 6.
Note: Do not use the built-in Weibull function from numpy.random.
In [ ]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(d) Run the four statistical tests of randomness (mean, variance, chi-square, K-S) on your results from question 1(c). Does it pass all of them?
In [ ]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
The Beta(4,3) distribution has density
$$ f(x) = \begin{cases} 60 x^3 (1 - x)^2 & 0 \le x \le 1 \\ 0 & \text{otherwise} \end{cases} $$The maximum of this distribution is $f(0.6) = 2.0736$ (exactly).
(a) Write a program that uses the Rejection Sampling Method to generate a set of 10,000 samples from this distribution and plot the density distribution, the CDF, and a histogram of the sampled distribution on the same graph.
Note: Do not use the built-in Beta function from numpy.random.
In [ ]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(b) Run the four statistical tests of randomness (mean, variance, chi-square, K-S) on your results from question 2(a). Does it pass all of them?
In [ ]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
Nuclei are composed of an ensemble of tightly bound protons and neutrons, which are in turn, tightly bound ensembles of quarks and gluons. If two nuclei are made to collide at very high energies, such as at the CERN Large Hadron Collider (LHC), they can vaporize into quarks and gluons and form a Quark Gluon Plasma.
The left-hand figure below shows two Au (gold) nuclei just after a grazing collision in which only a fraction of the 197 protons and neutrons (collectively called nucleons) actually interact. The right-hand side of the figure shows the substructure of the nucleons, which are composed of three valence quarks and the gluons that hold them together.
Only some of the nucleons interact when the nuclei collide unless they hit perfectly head-on. We define the offset of the two nuclei by an impact parameter, $b$, defined as the distance between the centers of the two nuclei as seen in a beam's-eye view. This figure illustrates the two-dimensional geometric overlap region shaded green.
The distribution of nucleons within a nucleus is not uniform. The radial distribution for spherical nuclei is generally described by the Woods-Saxon density profile, given by
$$ \rho(r) = \frac{\rho_0 (1 + wr^2/R^2)}{1 + \exp((r-R)/a)} $$where $R$ is the average radius of the nucleus, and $a$, $w$ are density parameters. These parameters come from empirical observations of electron scattering off various nuclei. The three parameter Woods-Saxon distribution describes a central nucleon density suppression to minimize the Coulombic potential, a maximum nucleon density radius, and then a fall off to zero density at infinite radius. Typical values for an assortment of nuclei is given in the Table.
| Nucleus | A | R (fm) | a (fm$^{-1}$) | w |
|---|---|---|---|---|
| C | 12 | 2.47 | 0 | 0 |
| O | 16 | 2.608 | 0.513 | -0.051 |
| Al | 27 | 3.07 | 0.519 | 0 |
| S | 32 | 3.458 | 0.61 | 0 |
| Ca | 40 | 3.76 | 0.586 | -0.161 |
| Ni | 58 | 4.309 | 0.516 | -0.1308 |
| Cu | 63 | 4.2 | 0.596 | 0 |
| W | 186 | 6.51 | 0.535 | 0 |
| Au | 197 | 6.38 | 0.535 | 0 |
| Pb | 208 | 6.68 | 0.546 | 0 |
| U | 238 | 6.68 | 0.6 | 0 |
(a) Plot the Woods-Saxon distribution $\rho(r)$ and $r^2\rho(r)$ for gold from $r$ = 0 fm to $r$ = 18 fm.
In [ ]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(b) Let’s create a realistic distribution for two gold ions, A, and B.
First, distribute 197 nucleons for each nucleus using the Woods-Saxon distribution multiplied by the spherical coordinate weighting factor, $r^2$. Use the distribute1d function from the tour to sample $r$ values numerically.
Then, use uniform sampling of a number $u$ from [0,1) scaled by $2\pi$ to obtain the azimuthal angle $\phi$ = 2$\pi u$
Followed by uniform sampling of another number $v$ = [0,1), transformed to $\theta = \cos^{-1}(2v - 1)$ to obtain the polar angle $\theta$.
(Note that this is the physics definition of azimuthal and polar angles, not the math definition.)
The reason for these transformations is that the volume element $dV = r^2\sin\theta dr d\theta d\phi$ has both radial and polar angle dependence. If you do not sample in this way, your results will be bunched up at the poles and not uniformly distributed inside of the sphere.
In [ ]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(c) Once you have the spherical coordinates for each of the nucleons in each nucleus, convert to cartesian coordinates. In nuclear collisions, the beam axis along which the particles travel/collide is the $z$ axis. The polar angle $\theta$ is measured up from the positive-$z$ axis and the azimuthal angle $\phi$ is measured up from the positive $x$-axis.
Plot the nuclei in 2d for two different planar projections: the beam's eye view ($x$-$y$) and the side view ($x$-$z$) side-by-side.
Let the nucleons have radius of 1 fm, make the nucleons from each nucleus a different color, and displace them by a randomly chosen impact parameter $b$ between 0 and 18 fm. For example, shift ion A by $b/2$ to the left along the $x$-axis and B by $b/2$ to the right along the x-axis.
In [ ]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(d) Since we cannot measure the impact parameter directly for each collision, we need a method to estimate the geometric overlap so that we can categorize the collisions based on the number of nucleons participating in each collision and how many nucleon-nucleon collisions occurred. Nucleons that pass close enough to each other will experience an inelastic collision that will rip the quarks right out of the nucleon. The distance that dictates when an inelastic collision occurs is governed by the interaction cross-section, which has been measured very carefully at proton colliders. The inelastic cross-section at collision energy of 200 GeV is $\sigma_{inel}$ = 42 mb (1 mb = 1 millibarn = 10$^{-24}$ cm$^2$ = 10$^{-31}$ m$^2$), defining an overlap area between two nucleons.
Make an interact object for your two plots in part (c) that let's you vary the impact parameter from 0 to 18 fm. On the plots, change the color of the nucleons that "participate" in interactions to something else. (Let the participating nucleons in nucleus A be one color and those in nucleus B be another color.)
To do this, you’ll need to find all the nucleons from nucleus A that are within a distance, $D$, of a nucleon from nucleus B using $D = \sqrt{\sigma_{inel}/\pi}$, where $\sigma_{inel}$ is the cross sectional area measured in millibarns. Don’t forget to convert the millibarns to fm$^{2}$.
In [ ]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell