From e72d4cb3b6f16bcdc4f9ed49b70a1b7200864744 Mon Sep 17 00:00:00 2001
From: Frank Sauerburger <frank@sauerburger.com>
Date: Tue, 8 Aug 2017 18:36:31 +0200
Subject: [PATCH] Full review

Change typos and minor details after a full review. The review also lead to a
new section about further reading.
---
 README.md | 294 +++++++++++++++++++++++++++++++++++++-----------------
 1 file changed, 200 insertions(+), 94 deletions(-)

diff --git a/README.md b/README.md
index 99fb0c0..9df20fa 100644
--- a/README.md
+++ b/README.md
@@ -2,8 +2,7 @@ This repository consists of a collection of python examples intended as an
 introduction on the usage of python in data analysis, especially for the
 advanced laboratories in physics at the University of Freiburg. In previous
 years code examples for [ROOT](https://root.cern.ch/) have been provided.
-Material on the usage of python was missing. The code examples shown in this
-repository follow the examples shown in the ROOT introductions.
+Material showing how to use python for the same task was missing. 
 
 # Table of Contents
 1. [Installation](#installation)
@@ -13,13 +12,14 @@ repository follow the examples shown in the ROOT introductions.
 5. [Plotting Functions](#numpy-arrays)
 6. [Plotting Data Points](#plotting-data-points)
 7. [Reading, Plotting and Fitting Experimental Data](#reading-plotting-and-fitting-experimental-data)
+8. [Further Reading](#furhter-reading)
 
 # Installation
 To get started with python for data analysis in the advanced laboratories you
 need the python interpreter. In this document we will use `python3`. The
 additional packages `numpy`, `scipy` and `matplotlib` are useful for data
 analysis and data presentation. To install all the packages on Ubuntu, you
-can run the following command line.
+can run the following commands.
 
 <!-- console -->
 ```bash
@@ -28,35 +28,43 @@ $ sudo apt-get install -y python3 python3-pip
 $ pip3 install --user numpy scipy matplotlib
 ```
 
+The `--user` argument for `pip` installs the python package in your home
+directory, which hides potentially older packages installed with `apt-get`.
+
 # Prerequisites and Structure
 
-This tutorial assumes that you have a little experience in python, which
+This tutorial assumes that you have some experience with python, which
 includes variable assignment, function calling and function definition, if statements and
 for loops. The tutorial uses only the very basics, such as variable assignments
 and function calls, but it is certainly advisable to know about control
 structures.
 
 To catch up on these aspects, you can refer to the [python
-documentation](https://docs.python.org/3/tutorial/).
+documentation](https://docs.python.org/3/tutorial/). If your are already
+familiar with another programming language, it should be quite intuitive to
+switch to python.
 
-The tutorial is structured into several examples, which might depend on each
+This tutorial is structured into several examples, which might depend on each
 other. The examples show code snippets which you are supposed to copy to a text
-editor. The scripts can then be executed in a terminal. You are invited to use
-the interactive mode of python or ipython instead of copying the code to a text
-editor. Recently jupyter notebooks have been
+editor. The scripts can then be executed in a terminal. Besides this modus
+operandi, you are invited to use
+the interactive mode of python or ipython instead and copying the code directly
+to the python interpreter. Recently jupyter notebooks have been
 become very popular. I recommend you to try out these different platforms and
 choose the one most suited for you.
 
-This repository does not contain ready-made python example scripts or plots.The
-only resource with code is this README file. The idea behind this is, that there
+This repository does not contain ready-made python example scripts or plots. The
+only code resource is this README file. The idea behind this is, that there
 shouldn't be duplications of code snippets, which will be out-of-sync eventually.
 The repository is set up, such that
-each commit triggers continues integration tasks, which parses the examples from
+each commit triggers continues integration tasks on the server, which parses the examples from
 the README and executes them with the
 [doxec](https://srv.sauerburger.com/esel/doxec) package. This means, you can
-download [read-made scripts and
+download [ready-made scripts and
 plots](https://srv.sauerburger.com/esel/FP-python-examples/-/jobs/artifacts/master/download?job=doxec_test)
- produced by the continues integration task.
+produced by the continues integration task.
+
+Let's get going!
 
 # 'Hello World' Example
 The first example is basically a 'Hello World' script, to check whether python
@@ -65,7 +73,7 @@ content.
 
 <!-- write hello_world.py -->
 ```python
-# load math library with sqrt function
+# Load math library with sqrt function.
 import math
 
 print("Example 1:")
@@ -94,23 +102,26 @@ arrays are defined in the numpy package, and are implemented in a very
 efficient way. 
 
 To get stared with numpy arrays create a file `np_arrays.py` and add all lines
-listed in this chapter. The first line should be an import statement.
+listed in this section. The first line should be an import statement.
 <!-- write np_arrays.py -->
 ```python
+# Import the numpy library.
 import numpy as np
 ```
 In this example we create a numpy array `numbers` containing my favorite numbers from
-a python list.
+the python list `[4, 9, 16, 36, 49]`.
 <!-- append np_arrays.py -->
 ```python
+# Create a numpy array by passing a list to np.array.
 numbers = np.array([4, 9, 16, 36, 49])
 ```
 Having all these numbers in a numpy array makes bulk computations very efficient.
-Assume, we want to calculate the square root of all these numbers, wen can
+Assume, we want to calculate the square root of all these numbers. We can
 simply use numpy's `sqrt` method do perform the same operation on all elements
 of the array at the same time.
 <!-- append np_arrays.py -->
 ```pyton
+# Calculte the square root for each item in the array numbers.
 roots = np.sqrt(numbers)
 ```
 
@@ -118,6 +129,7 @@ Since the resulting variable `roots` is also a numpy array, we can perform
 similar operations on this variable.
 <!-- append np_arrays.py -->
 ```pyton
+# Perform other calucaltions for each element separately.
 something_else = 1.5 * roots - 4
 ```
 Numpy arrays overload the typical arithmetic operations, such that the above
@@ -127,13 +139,15 @@ use such vectorized statements, and try to avoid manually looping over all the
 values. Using a python loop to run over $`10^3`$ values is probably fine, but you
 don't want to wait for a python loop iterating over 10^6 or 10^9 values.
 
-Finally add a print statement to check that all the caluclations are as expected. 
+Finally add a print statement to check that all the calculations have been
+carried out as expected. 
 <!-- append np_arrays.py -->
 ```python
 print("The result is", something_else)
 ```
 
-When executed you should get the following printout.
+When executed, you should get the following printout. You can convince yourself,
+that the calculations are correct.
 <!-- console_output -->
 ```bash
 $ python3 np_arrays.py
@@ -145,19 +159,19 @@ introduction. It is definetely worth glancing at the
 [documentation](https://docs.scipy.org/doc/numpy/index.html).
 
 # Plotting Functions
-One major aspect of data analysis is also data presentation. This includes the
-geneation of diagrams and plots. You can use the powerful library matplotlib to
-create publication-quality plots from python. The goal of this example is to plot
-the cropped parabola f(x), which is limited y=4 for x>=2. 
+One major aspect of data analysis is the visualization of data. This includes the
+generation of diagrams and plots. You can use the powerful library matplotlib to
+create high-quality plots in a python environment. The goal of this example is to plot
+the cropped parabola
 
 ```math
     f(x) = \left\{\begin{array}{lr}
               x^2, & \text{for } x < 2\\
               4, & \text{for } 2 \leq x
-            \end{array}\right.
+            \end{array}\right.,
 ```
 
-The final plot should look like this.
+which looks like this:
 ![Plot of cropped parabola](https://srv.sauerburger.com/esel/FP-python-examples/-/jobs/artifacts/master/raw/cropped_parabola.png?job=doxec_test)
 
 Create the file `func_plot.py` and add the following lines.
@@ -171,39 +185,59 @@ matplotlib.use('Agg')
 -->
 <!-- append func_plot.py -->
 ```python
+# Import the numpy library.
 import numpy as np
+
+# Import the powerfull plotting library.
 import matplotlib.pyplot as plt
 ```
 
 Plotting a function with matplotlib means plotting many points connected by a
-line. First we create an array with 200 equidistant values in the interval
-[-2.5, 3]. This array functions as a grid of x-values, for which we calculate
-the y values. 
+line. First we create an array with 200 equidistant $`x`$-values in the interval
+$`[-2.5, 3]`$. This array functions as a grid, for which we calculate
+the $`y`$ values. 
 <!-- append func_plot.py -->
 ```python
 x = np.linspace(-2.5, 3, 200)
 ```
 
 We can easily calculate the square of all these values with `x**2`. Cropping the
-right part is a bit more complex. First we create an index array of 1's and 0's, which
-indicate whether x >= 2. This index array can be used to select a subset of
-y-values. Finally we can assign the value 4 to this subset, and therefore
-effectively cropping the parabola. The full example reads:
+right part is a bit more complex. First we create an index array of `1`'s and
+`0`'s, which
+indicate whether $`x \geq 2`$. This index array has the same length as our
+$`x`$-grid. The first elements of the index array are `0`'s, since the
+corresponding $`x`$-value is below two. At some point in the array, the value
+changes to `1`, since then the corresponding $`x`$ values satisfy $`x\geq2`$.
+The index array can be used to select a subset of
+$`y`$-values, namely all $`y`$-values, for which $`x\geq 2`$. Finally we can
+assign the value $`4`$ to this subset, and therefore
+effectively cropping the parabola. The implementation in python of the algorithm outlined
+above  is rather short.
 <!-- append func_plot.py -->
 ```python
+# Calculate the regualar parabola.
 y = x**2
+
+# Create index array.
 idx = (x >= 2)
+
+# Set all y-values to 4, for which x >= 2.
 y[idx] = 4
 ```
 
-The final step of this example is to plot the points and connect the with a
-line by using matplolib's plot method. We can also add axis labels and save the
+The final step of this example is to call matplotlib, which plots the points and
+connects the with a line. Additionally, We can add axis labels and save the
 resulting figure.
 <!-- append func_plot.py -->
 ```python
+# Plot a line specified by x- and y-arrays.
 plt.plot(x, y)
-plt.xlabel("$x$")  # latex synatx can be used
+
+# Set axis label. Latex expression can be used.
+plt.xlabel("$x$") 
 plt.ylabel("cropped parabola")
+
+# Save the figure. Various different output formats are available.
 plt.savefig("cropped_parabola.eps")
 ```
 <!-- append func_plot.py
@@ -213,7 +247,7 @@ plt.savefig("cropped_parabola.png")
 ```
 -->
 
-Run your script and check the file `cropped_parabola.eps` is created.
+Run your script and check that the file `cropped_parabola.eps` is created.
 <!-- console -->
 ```bash
 $ python3 func_plot.py
@@ -231,14 +265,14 @@ cropped_parabola.png
 # Plotting Data Points
 A typical task in the advanced laboratory might be to compare measured data
 points to an expected function. Lets assume the expected function is the cropped
-parabola $`f(x)`$ from the previous examples. We will use random data points in
+parabola $`f(x)`$ from the previous example. We will use random data points in
 this example, since we haven't actually measured real data, which is expected to
 follow $`f(x)`$.
 
 This example is based on the code from the previous example. Copy the file from
-the previous examples to `data_plot.py`, such that we an append to it and keep
-the function plotting code. Please add the code snippets in this section to the
-file `data_plot.py`.
+the previous examples to `data_plot.py`, such that we an append the following
+code snippets to `data_plot.py` and keep
+the plotting code from the previous example as it is.
 <!-- console
 ```bash
 $ cp func_plot.py data_plot.py
@@ -246,39 +280,51 @@ $ cp func_plot.py data_plot.py
 -->
 
 We generate the pseudo data points by adding random deviations to the expected
-y-values. Lets assume we measured data points for all quarter-integer x-values.
+$`y`$-values. Lets pretend we have measured data points for all half-integer
+$`x`$-values in the interval $`[-2.5, 3]`$.
 <!-- append data_plot.py -->
 ```python
+# Create x-value grid for the measured data.
 x_data = np.array([-2.5, -2, -1.5, -1, -.5, 0, .5, 1, 1.5, 2, 2.5, 3])
 ```
 
-We evaluate the function $`f(x)`$ again for the `x_data` values. The random
-deviations are drawn from a centered normal distribution with a standard
-deviation of 0.3. The third argument of `numpy.random.normal` specifies, how
+We evaluate the function $`f(x)`$ again for the `x_data` values. The resulting
+array `y_data` matches the curve from the previous example perfectly. We draw random
+deviations from a centered normal distribution with a standard
+deviation of 0.3 and add them to `y_data`. The third argument of `numpy.random.normal` specifies, how
 many random samples we want to draw. We use `len(x_data)`, since we want to draw
-a random deviation for each x-value.
+a random deviation for each $`x`$-value.
 <!-- append data_plot.py -->
 ```python
+# Calculate square of x_data points.
 y_data = x_data**2
+
+# Crop y_data points if x_data >= 2.
 y_data[x_data >= 2] = 4
+
+# Add random fluctuations to y_data.
 y_data += np.random.normal(0, 0.3, len(y_data))
 ```
 
 Finally we can add this to our plot. Since our data points are subject to
 statistical fluctuations, we would like to use matplotlib's `errorbar` method,
-which draw our data points with errorbars.
+which draws our data points with as dots with error bars.
 
 <!-- append data_plot.py -->
 ```python
+# Draw with error bars, similar to plot().
 plt.errorbar(x_data, y_data, 0.3, fmt="ko", capsize=0)
+
+# Save figure.
 plt.savefig("measurement.eps")
 ```
 
-The character `k` in the format parameter `fmt` sets the color to blac*k*, the
+The character `k` in the format parameter `fmt` sets the color to *black*, the
 `o` in `fmt` changes the style to *large dots*. The optional parameter `capsize`
-modifies the style of the error bars. You can play with these options and see
+modifies the style of the error bars. You can play with these options to see
 what happens or have a look at the
-[documentation](https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.errorbar.html).
+[documentation](https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.errorbar.html)
+for more information about the options.
 
 <!-- append data_plot.py
 ```python
@@ -300,18 +346,18 @@ measurement.png
 After running `data_plot.py` you should have a plot similar to this.
 ![Plot of cropped parabola with data points](https://srv.sauerburger.com/esel/FP-python-examples/-/jobs/artifacts/master/raw/measurement.png?job=doxec_test)
 
-# Reading, Plotting and Fitting Experimental data
+# Reading, Plotting and Fitting Experimental Data
 We are given with experimental data from a radioactive decay in this example.
 The experimental setup consisted of a radioactive probe, a detector and a
 multi-channel-analyzer. The recorded data in `decay.txt` consist of two
 tab-separated columns. The first column is called channel. Each channel
 corresponds to a certain energy range. The multi-channel-analyzer maintains a
-counter for each channel. For each recorded decay the internal counter which
-corresponds to the measured energy is incremented. The seconds column stored
+counter for each channel. Decay causes the multi-channel-analyzer to increment the internal counter which
+corresponds to the energy of the measured decay. The seconds column stored
 these counts. Open the file with out favorite text editor and have a look at the
 data.
 
-As usual, create the file `decay.py` and the import statements, which we need
+As usual, create the file `decay.py` and add the import statements, which we need
 for this example.
 
 <!-- Add additional files for non-X11 environment in CI -->
@@ -323,27 +369,41 @@ matplotlib.use('Agg')
 -->
 <!-- append decay.py -->
 ```python
+# Import the numpy library.
 import numpy as np
-import scipy.optimize  # provides fit routines
-import scipy.stats  # provides statistical tests and distributions
+
+# Import the scipy library with fit routines.
+import scipy.optimize
+
+# Import the scipy library with probability distributions and statistical tests.
+import scipy.stats
+
+# Import the plotting library.
 import matplotlib.pyplot as plt
 ```
 
 To inspect the provided data, we can plot the raw data points first. Numpy
-provides the function `loadtxt`, which read a whitespace-separated file into a
+provides the function `loadtxt`, which reads a whitespace-separated file into a
 numpy array. The function returns a two dimensional array. The outer array has
-one entry for each line, and the inner two entries in our case, one for the
+one entry for each line in the text file. The inner array has two entries in our
+case, one for the
 channel and the other one for the event count. We can use `transpose()` to flip
-the matrix, such that the outer array has two entries, one with an array of
+the matrix, such that the outer array has tswo entries, one with an array of
 channel values and the other one with event counts. Since our measured event
 counts stem from radioactive decay, we know that the event counts follow a
-Poisson distribution. Therefore, the uncertainty of the events counts is simply
-the square root of the number of events.
+Poisson distribution. Therefore, the uncertainties of the event counts are simply
+the square roots of the number of events.
 
 <!-- append decay.py -->
 ```python
+# Read both columns from the text file.
 channel, count = np.loadtxt("decay.txt").transpose()
+
+# Calculate the uncertaintiy on the number of events per channel.
 s_count = np.sqrt(count)
+
+# Create and save a raw version of the plot with data points.
+# The label will be used later to identify the curves in a legend.
 plt.errorbar(channel, count, s_count, fmt='.k', capsize=0, label="Data")
 plt.savefig("decay_raw.eps")
 ```
@@ -359,23 +419,23 @@ $ python3 decay.py
 -->
 
 The plot of the raw data is shown below. The plot shows the channel on the
-x-axis and the number of events per bin on the y-axis. From the experimental setup we expect a
+$`x`$-axis and the number of events per bin on the $`y`$-axis. From the experimental setup we expect a
 linearly rising background plus a Gaussian peak.
 ![Data only for the decay](https://srv.sauerburger.com/esel/FP-python-examples/-/jobs/artifacts/master/raw/decay_raw.png?job=doxec_test)
 
-Judging from the plot, it looks like the assumed model describes the data. We
-would like this model to our data to determine the optimal value of the model
-parameter and their uncertainties (and the covariance matrix). Lets give a more
+Judging from the plot, it looks like the assumed model could describe the data. We
+would like to fit this model to our data to determine the optimal values of the model
+parameters and their uncertainties (and the covariance matrix). Lets give a more
 formal version of the expected model 
 ```math
-  n(c) = A exp\left(-\frac{(c-m)^2}{2 s^2}\right) + y_0 + bc,
+  n(c) = A \exp\left(-\frac{(c-m)^2}{2 s^2}\right) + y_0 + bc,
 ```
-where $`n(c)`$ is the expected number of events in channel $`c`$; $`A, m`$
-and  $`s`$ are the height, center and width of the Gaussian respectively. The
+where $`n(c)`$ is the expected number of events in channel $`c`$ ; $`A, m`$
+and  $`s`$ are the height, center and width of the Gaussian, respectively. The
 parameters $`y_0`$ and $`b`$ are the usual parameter of a linear curve which is
 assumed to describe our background. The model can be implemented in python as a
-function. The first parameter should be the x value, all following arguments are
-free parameters of the model. The return value corresponds to the y value, in
+function. The first parameter should be the $`x`$-value, all following arguments are
+free parameters of the model. The return value corresponds to the $`y`$-value, in
 our case the expected number of events.
 
 <!-- append decay.py -->
@@ -384,48 +444,66 @@ def model(channel, m, s, A, y0, b):
   return  A * np.exp(-0.5 * (channel - m)**2 / s**2) + y0 + b * channel
 ```
 
-Please note that we are going to make an approximation with this definiton.
-Strictly speaking, comparing thhe return values of our model to the measured
+Please note that we are making an approximation with this definition.
+Strictly speaking, comparing the return values of our model to the measured
 count is not correct. The variable channel corresponds to the radiation energy
 measured with the setup. Lets assume channel $`c_i`$ corresponds to energy
 $`E_i`$. If we measure $`n_i`$ events in channel $`c_i`$, this means that we
 have measured $`n_i`$ in the energy interval $`[\frac{1}{2}(E_{i-1} + E_i),
-\frac{1}{2}(E_i} + E_{i+}]`$. The proper way is to integrate our continuous
-function $`n(c)`$ in each bin $`[c_{i} - \frac{1}{2}, c_{i} } \frac{1}{2}]`$. The
-prodcedure show here is a good approximation, if the function can be considered
-is linear within each bin. However, the parameter $`A`$ and $`b`$ are not normalized to
+\frac{1}{2}(E_i + E_{i+1})]`$. The proper way is to integrate our continuous
+function $`n(c)`$ in each bin $`[c_{i} - \frac{1}{2}, c_{i} + \frac{1}{2}]`$ and
+compare these bin-wise integrals to the measured data. The
+procedure shown here is a good approximation, if the function can be considered
+to be linear within each bin. However, the parameter $`A`$ and $`b`$ are not normalized to
 the bin width in this case.
 
 
 To fit this model to our experimental data, we can use the function `curve_fit`
 provided by the scipy package. The function `curve_fit` performs a least square
 fit and returns the optimal parameters and the covariance matrix. The fit might
-not converge on its one. We can guide optimization procedure by providing
+not converge on its one. We can guide the optimization procedure by providing
 suitable start values of the free parameters. From the plot I read off a height
 $`A=50`$, a center $`m=60`$ and a width $`s = 10`$ for the Gaussian part and
-$`y_0 = 20`$, $`b = 1`$ for the linear part. These values don't have to be
+$`y_0 = 20`$ and $`b = 1`$ for the linear part. These values don't have to be
 accurate. They should be a rough estimation, this is usually enough to get a
-stable fit result.
+stable fit result. More information on the fitting method can be found in the
+[documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html).
 
 <!-- append decay.py -->
 ```
+# Define the intial values of the free parameters.
+# Remember, that we defined our model as n(c; m, s, A, y0, b)
 p0 = (60, 10, 50, 20, 1)
-popt, pcov = scipy.optimize.curve_fit(model, channel, count, p0, np.sqrt(count))
+
+# Perform the actual fit. The parameters are
+#   (1) Model to fit
+#   (2) Array of x-values
+#   (3) Array of y-values to which the model shold be fitted
+#   (4) Array with inital values for the free parameters
+#   (5) Array with uncertainties on the y-values.
+popt, pcov = scipy.optimize.curve_fit(model, channel, count, p0, s_count)
 ```
 
 To visualize the fitted model, we need to evaluate our model with the optimized
-parameters `popt`. We are now ready to add the fitted curve to the plot and save
-it.
+parameters `popt`. 
 
 <!-- append decay.py -->
 ```
+# Evaluate the model with the optimized parameter.
 fit_count = model(channel, *popt)
-plt.plot(channel, model(channel, *popt), label="Linear + Gauss")
 
+# Plot a curve representing the fitted model.
+plt.plot(channel, fit_count, label="Linear + Gauss")
+
+# Add axis labels.
 plt.xlabel("Channel")
 plt.ylabel("Counts")
-plt.ylim(0, 1.1 * max(count))
+
+# Add a legend to identify data and our fit. This method uses values passed to
+# the the optional arguemnt 'label' of plot() and errorbar().
 plt.legend()
+
+# Save the figure.
 plt.savefig("decay.eps")
 ```
 <!-- append decay.py
@@ -441,10 +519,11 @@ $ python3 decay.py
 The result should look like this.
 ![Data and fit for the decay](https://srv.sauerburger.com/esel/FP-python-examples/-/jobs/artifacts/master/raw/decay.png?job=doxec_test)
 
-Usually we want to measure some quantity with the experimental setup. For this
+Usually we want to measure some quantity with an experimental setup. For this
 we need the optimized parameters and the covariance matrix returned by the
 fit. Lets assume we are interested in the best fit value of the parameters and
-their uncertainties. We can add the following print outs, to display this kind
+their uncertainties. The uncertainties are the square roots of the diagonal of
+the covariance matrix. We can add the following print statements, to display this kind
 of information.
 
 <!-- append decay.py -->
@@ -455,32 +534,42 @@ print("  s = %g +- %g"  % (popt[1], np.sqrt(pcov[1][1])))
 print("  A = %g +- %g"  % (popt[2], np.sqrt(pcov[2][2])))
 print("  y0 = %g +- %g" % (popt[3], np.sqrt(pcov[3][3])))
 print("  b = %g +- %g"  % (popt[4], np.sqrt(pcov[4][4])))
-print()
+print()  # print blank line
 ```
 
 A $`\chi^2`$-test can also be performed, to assess the goodness of this fit. In
 a counting experiment like this one, we can rely on scipy's `chisquare`, which
 returns the $`\chi^2`$ and the $`p`$-value. The `chisquare` method assumes, that
-the uncertainties are the square root of the expectation. If this is not the
+the uncertainties are the square root of the expected number of events. If this is not the
 case, we have to compute the $`\chi^2`$ manually. The following example shows
-both examples. The print statements produce the same output. Please note, that
+both, the usage of `chisquare` and the manual computation. The print statements
+for each method produce the same output. Please note, that.
 we have five degrees of freedom, since we have five free parameters in our
 model.
 
 <!-- append decay.py -->
 ```python
+# Degrees of freedom in our model.
+dof = 5
+
 print("chi^2 from scipy:")
-chi2, p = scipy.stats.chisquare(count, fit_count, ddof=5)
-print("  chi2 / ndf = %g / %d" % (chi2, len(count) - 6))
+
+# Call chisquare and print.
+chi2, p = scipy.stats.chisquare(count, fit_count, ddof=dof)
+print("  chi2 / ndf = %g / %d" % (chi2, len(count) - dof))
 print("  p-value = %g" % p)
 
-print()
+print()  # print blank line
 
 print("Manual chi^2 test:")
+
+# Calculate chi^2
 uncertainty = np.sqrt(fit_count)
 chi2 = ((count - fit_count)**2 / uncertainty**2).sum()
-p = scipy.stats.distributions.chi2.sf(chi2, len(count) - 6)
-print("  chi2 / ndf = %g / %d" % (chi2, len(count) - 6))
+
+# Calculate p-value and print
+p = scipy.stats.distributions.chi2.sf(chi2, len(count) - 1 - dof)
+print("  chi2 / ndf = %g / %d" % (chi2, len(count) - dof))
 print("  p-value = %g" % p)
 
 ```
@@ -504,3 +593,20 @@ Manual chi^2 test:
   chi2 / ndf = 120.577 / 122
   p-value = 0.519418
 ```
+
+Congratulations! You have mastered the first steps to analysis experimental data
+with python.
+
+# Further Reading
+This tutorial can not cover all topics which can be relevant for the advanced
+laboratories. Here is a list with online resources, which might be useful.
+
+ - [Python](https://docs.python.org/3/) documentation
+ - [Numpy and Scipy](https://docs.scipy.org/doc/) documentation
+ - [Matplotlib](http://matplotlib.org/contents.html) documentaiton
+ - [Orthogonal distance regression](https://docs.scipy.org/doc/scipy/reference/odr.html) - least square fit, which considers
+   uncertainties in $`x`$ and $`y`$ directions.
+ - Histograms in
+   [numpy](https://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram.html)
+   and
+   [matplotlib](https://matplotlib.org/api/pyplot_api.html?highlight=matplotlib%20pyplot%20hist#matplotlib.pyplot.hist).
-- 
GitLab