Basics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic concepts\n",
    "\n",
    "By now, CRPropa is installed together with the Juypter interface and we are ready to begin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Candidates and modules\n",
    "\n",
    "To use CRPropa commands, tools and classes within Pyhton, the first step is to import it into the main namespace:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from crpropa import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or, alternatively, if one needs better control over different commands and modules `import crpropa as crp` which will require adding the `crp.` prefix before every command or class initiated from CRPropa.\n",
    "\n",
    "The principal object of interest, a cosmic ray, within the framework is called a **candidate**. The candidate is a particle described by an associated physical state. The particle state can be defined with many different properties, but two of the most fundamental ones are the particle's **position** and the particle's movement **direction**. The actual state of a candidate is referred as \"current\" one.\n",
    "Here is an example where Vector3d is a CRPropa object representing a 3D vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = Candidate()\n",
    "c.current.setPosition(Vector3d(10, 0, 0))\n",
    "c.current.setDirection(Vector3d(-1, 0, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The candidate itself is a passive object and only external mechanisms, called **modules**, can change its state.\n",
    "The first example of such a module in the context of the framework application is the rectilinear propagational module - ``SimplePropagation``. Its purpose, as the name suggests, is to propagate the candidate in the simplest way, i.e., to change the candidate's position following a straight line in the defined direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "module = SimplePropagation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To describe how modules act on candidates, one should recall those numerical simulations necessarily discretise the continuity of the real world with computational **steps**. The candidate state is changed by modules in steps. In every step, the module process candidate with the **process method**.\n",
    "\n",
    "Before continuing with processing, it is worth mentioning that modules can or even must take arguments during the initialization stage to set-up or change their default behaviour and it depends on the specific module which settings can be defined or modified. The modules' settings are described in [CRPropa documentation](http:https://crpropa.github.io/CRPropa3/annotated.html).\n",
    "\n",
    "In the already introduced module, ``SimplePropagation``, one could change the minimum and the maximum distance for which the candidate can jump within one step. These changes can also be done later, after the initialization, with so-called **setters** (set methods) or just checked with **getters** (get methods)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0\n"
     ]
    }
   ],
   "source": [
    "module = SimplePropagation(1, 3) # (minStep, maxStep)\n",
    "module.setMaximumStep(2)\n",
    "print(module.getMaximumStep())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The process method can be used in a conjunction with any kind of loop, e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In step 1. position is Vector(9, 0, 0)\n",
      "In step 2. position is Vector(7, 0, 0)\n",
      "In step 3. position is Vector(5, 0, 0)\n",
      "In step 4. position is Vector(3, 0, 0)\n",
      "In step 5. position is Vector(1, 0, 0)\n"
     ]
    }
   ],
   "source": [
    "for i in range(5):\n",
    "    module.process(c)\n",
    "    print(\"In step {}. position is {}\".format(i+1, c.current.getPosition()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After five iterations, the position is changed five times in the negative direction of the x-axis which is defined by ``setDirection()``. What needs an additional explanation is the propagated distance in these steps. To avoid the clash between different active modules in determining the step size, modules modify the candidate's current step with the minimum prescribed step while setting the maximum one for the next round. In that way, the next step would be bounded by the module which required the smallest maximum step. In case of the current example, the first step was 1 (minStep argument in ``SimplePropagation``), while all the later steps are 2 (given by ``setMaximumStep`` above).\n",
    "\n",
    "At this point, before introducing a new physical content, a slight digression in a form physical units is needed. All physical quantities and parameters used in the code are in expressed the SI units. For example, those numbers, 1, 2 and 3, in the above example are 1, 2 and 3 meter(s). From now on, the explicit units will be used to avoid possible confusion, and it is advisable to do the same in every CRPropa code. If some units are not defined within the framework, one could use external python library such as [pint](http:https://pint.readthedocs.io) to carry out the conversions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 Mpc is 3.085677581491367e+22 meters in the international system of units, while 1 GeV is 1.602176487e-10 joules\n"
     ]
    }
   ],
   "source": [
    "print(\"1 Mpc is {} meters in the international \\\n",
    "system of units, while 1 GeV is {} joules\".format(1*Mpc, 1*GeV))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the next example how to use pint:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flux: 1 eV^-1 m^-2 s^-1 sr^-1 = 3.1556925974699995e+22 GeV^-1 km^-2 year^-1 sr^-1\n"
     ]
    }
   ],
   "source": [
    "from pint import UnitRegistry\n",
    "ureg = UnitRegistry()\n",
    "\n",
    "def convertunits(unitstring, convertto):\n",
    "        unitfrom = ureg.parse_expression(unitstring)\n",
    "        unitto = unitfrom.to(ureg.parse_expression(convertto))\n",
    "        return unitto.magnitude\n",
    "\n",
    "fromunits = \"eV^-1 m^-2 s^-1 sr^-1\"\n",
    "tounits = \"GeV^-1 km^-2 year^-1 sr^-1\"\n",
    "print(\"Flux: 1 {1} = {0} {2}\".format(\n",
    "    1*convertunits(fromunits, tounits), fromunits, tounits))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The candidate can virtually be any classified particle or nucleus with the assumption that its energy is dominating its mass, i.e., $pc=\\sqrt{E^2-(mc^2)^2}\\approx E$, so that it travels at the speed of light in the observer reference frame. The reason behind it is that the framework by design does not explicitly take into account time. The time parameter is merely recovered from the propagated length. Thus, a propagation module is always required since any other physical module stays inactive when time does not pass.\n",
    "\n",
    "The particles are distinguished by the particle ID (PID) which inherits a numbering scheme defined by the Particle Data Group. CRPropa uses the [HepPID](http:https://lcgapp.cern.ch/project/simu/HepPDT/HepPDT.2.06.01/HepPID.pdf) implementation of the scheme. For the nuclei identification, CRPropa has its own implementation which is also specified by [PDG](http:https://pdg.lbl.gov/2008/reviews/rpp2008-rev-monte-carlo-numbering.pdf).\n",
    "\n",
    "After demonstrating ``SimplePropagation`` module, ``NuclearDecay`` is yet another module which is relatively self-contained, hence, can be presented without too much overhead. ``NuclearDecay`` is a module which represents the decay of a nucleus. It will act only on unstable nuclei, of course. Let's examine the following lines which employ the mentioned concepts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean length until the decay is 1.64 pc\n"
     ]
    }
   ],
   "source": [
    "def decay_distance(pid, E):\n",
    "    c = Candidate(pid, E)\n",
    "    c.current.setPosition(Vector3d(0, 0, 0)*Mpc)\n",
    "    c.current.setDirection(Vector3d(1, 0, 0))\n",
    "    init_id = c.current.getId()\n",
    "\n",
    "    module1 = SimplePropagation(0.01*pc, 10*pc)\n",
    "    module2 = NuclearDecay()\n",
    "\n",
    "    for step in range(200):\n",
    "        module1.process(c)\n",
    "        module2.process(c)\n",
    "        if c.current.getId() != init_id:\n",
    "            return c.current.getPosition().getR()\n",
    "        \n",
    "    return c.current.getPosition().getR()\n",
    "            \n",
    "# Li-8, half-life = 839.9 ms, E = 1 EeV\n",
    "pid = nucleusId(8, 3)\n",
    "E = 1*EeV\n",
    "\n",
    "# make sample\n",
    "dists = [decay_distance(pid, E)/pc for i in range(3000)]\n",
    "mean_dist = sum(dists)/len(dists)\n",
    "print(\"The mean length until the decay is {:.3g} pc\".format(mean_dist))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Depending on the particle energy, unstable Li-8 will travel a certain distance until it decays. The result can also easily be checked analytically because Li-8 has only one decay channel and then compared visually with the ``matplotlib`` plotting library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.57537201026\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XmcVNWd/vHP090sjawKGmxMAMEF\nXCK2iqLRuGKikBjXiVnUiZmMSzQmM1l+mkxWk6gTExfiOEbHGJdAEomjQaNMoqiRBhVBgyCiQEAR\nkE2Q7fv741Y3RdPd1Q19+1Z3Pe/Xq1636tatul8R+ulzzr3nKCIwMzMDKMu6ADMzKx4OBTMzq+NQ\nMDOzOg4FMzOr41AwM7M6DgUzM6vjUDAzszoOBTMzq+NQMDOzOhVZF9BSffv2jYEDB2ZdhplZuzJt\n2rR3IqJfoePaXSgMHDiQmpqarMswM2tXJL3RnOPcfWRmZnUcCmZmVsehYGZmdVILBUl3SHpb0sxG\n3pekn0uaK2mGpBFp1WJmZs2TZkvhTmB0E++fCgzNPS4Gbk2xFjMza4bUQiEi/gosb+KQscD/ROJZ\noLek/mnVY2ZmhWU5plAFLMh7vTC3z8zMMtIuBpolXSypRlLN0qVLd+g7Zs6ECRNg2bJWLs7MrAPJ\nMhQWAXvlvR6Q27ediLgtIqojorpfv4I35DXo8svhzDNh2rQd+riZWUnIMhQmAp/NXYU0ElgZEYvT\nOlntzBjz56d1BjOz9i+1aS4k3QscB/SVtBD4NtAJICLGAQ8DHwPmAu8BF6RVCzgUzMyaI7VQiIjz\nCrwfwCVpnb++D30o2b7RrNk/zMxKU7sYaG4NbimYmRXmUDAzszolEwpVVVBeDosXw/vvZ12NmVlx\nKplQqKiAAQMgAhYsKHy8mVkpKplQAHchmZkVUlKhUHsFkkPBzKxhJRUKtS0FX5ZqZtawkgwFtxTM\nzBpWUqHg7iMzs6aVVCi4+8jMrGklFQoDBkBZGSxaBBs2ZF2NmVnxKalQ6NwZ9twTtmyBhQuzrsbM\nrPiUVCiAB5vNzJpSsqHgcQUzs+2VXCj4CiQzs8aVXCi4+8jMrHElGwruPjIz217JhoJbCmZm2yu5\nUNhrr2S7cCFs2pRtLWZmxabkQqFLl+Rehc2bk5vYzMxsq5ILBXAXkplZY0oyFHxZqplZw0oyFNxS\nMDNrWEmHgi9LNTPbVkmGgruPzMwaVpKh4O4jM7OGlWQofPCDyXbBguTSVDMzS5RkKFRWwgc+kNy8\n5nUVzMy2KslQABgyJNnOmZNtHWZmxaRkQ2Ho0GTrUDAz26pkQ2GffZKtQ8HMbKuSDQW3FMzMtlfy\nofDqq9nWYWZWTFINBUmjJc2WNFfS1xt4/4OSJkt6XtIMSR9Ls558tQPN8+Z5Cm0zs1qphYKkcuBm\n4FRgGHCepGH1Dvt/wAMRcQhwLnBLWvXU160bVFUlgeDpLszMEmm2FA4H5kbEvIjYANwHjK13TAA9\nc897Af9IsZ7teLDZzGxbaYZCFbAg7/XC3L583wHOl7QQeBi4LMV6tuPBZjOzbWU90HwecGdEDAA+\nBtwtabuaJF0sqUZSzdKlS1vt5B5sNjPbVpqhsAjYK+/1gNy+fBcBDwBExDNAV6Bv/S+KiNsiojoi\nqvv169dqBbqlYGa2rTRDYSowVNIgSZ1JBpIn1jvmTeAEAEn7k4RC6zUFCvCYgpnZtlILhYjYBFwK\nTAJeIbnKaJak70oakzvsKuALkl4E7gU+HxGRVk31DR4MZWXJFNobNrTVWc3MildFml8eEQ+TDCDn\n77sm7/nLwKg0a2hKly7JNNrz5yf3K+y3X1aVmJkVh6wHmjPncQUzs61KPhQ8rmBmtlXJh4JbCmZm\nWzkUfK+CmVkdh4JbCmZmdUo+FAYOhIoKWLAA1q3Luhozs2yVfCh06gSDBiXPX3st21rMzLJW8qEA\n7kIyM6vlUMCDzWZmtRwKuKVgZlbLoYBvYDMzq+VQYGtLYfbsbOswM8uaQ4FkUrxu3eCtt2DZsqyr\nMTPLjkOBZPrs4cOT57NmZVuLmVmWHAo5BxyQbGfOzLYOM7MsORRy3FIwM3Mo1HFLwczMoVAnPxTa\nbkFQM7Pi4lDI2XNP6N0bli+HJUuyrsbMLBsOhRzJ4wpmZg6FPB5XMLNS51DI41Aws1LnUMjj7iMz\nK3UOhTy+AsnMSp1DIU+/frD77rBmDbz5ZtbVmJm1PYdCPbVdSB5XMLNS5FCop7YLyeMKZlaKHAr1\n+AokMytlDoV6HApmVsqaFQqSfifp45I6fIgMG5ZsX3kFNm/OthYzs7bW3B/ytwD/BMyRdK2kfVOs\nKVO9e8OAAbB+Pcybl3U1ZmZtq1mhEBF/johPAyOA+cCfJT0t6QJJndIsMAvuQjKzUtXs7iBJuwGf\nB/4ZeB64kSQkHkulsgz5slQzK1XNHVP4PfAk0A04PSLGRMT9EXEZ0L2Jz42WNFvSXElfb+SYsyW9\nLGmWpN/syH9Ea3NLwcxKVUUzj/uviHg4f4ekLhHxfkRUN/QBSeXAzcBJwEJgqqSJEfFy3jFDgW8A\noyJihaTdd+i/opUdfHCynT492zrMzNpac7uPvt/AvmcKfOZwYG5EzIuIDcB9wNh6x3wBuDkiVgBE\nxNvNrCdVw4dD584wdy6sXJl1NWZmbafJUJD0AUmHApWSDpE0Ivc4jqQrqSlVwIK81wtz+/LtA+wj\naYqkZyWNbqSOiyXVSKpZunRpgdPuvM6d3Vows9JUqPvoFJLB5QHADXn7VwPfbKXzDwWOy53jr5IO\njIh38w+KiNuA2wCqq6vbZP7SQw+FqVNh2jT46Efb4oxmZtlrMhQi4i7gLkmfiogJLfzuRcBeea8H\n5PblWwj8LSI2Aq9LepUkJKa28Fyt7tBDk21NTbZ1mJm1pSZDQdL5EfFrYKCkr9R/PyJuaOBjtaYC\nQyUNIgmDc0lugMv3B+A84FeS+pJ0J6V3y9if6/3Kf+LkRg+tzg2fT5uWWjVmZkWnUPfRLrlto5ed\nNiYiNkm6FJgElAN3RMQsSd8FaiJiYu69kyW9DGwGvhYRy1p6rjQMHw5duiSDze++m9zpbGbW0Sna\n2RJj1dXVUbOjfTotaCkAHHEEPPccPP44HH/8jp3SzKwYSJrW2C0E+Qp1H/28qfcj4vKWFtaeHHpo\nEgrTpjkUzKw0FOo+Kuke9drBZo8rmFmpaM7VRyWrdrDZVyCZWako1H30s4i4QtIfge0GHyJiTGqV\nFYFhw5LB5tde82CzmZWGQt1Hd+e216VdSDHq1Cm5s/m555I7mz2uYGYdXZPTXETEtNz2LyRzHa0A\nlgPP5PZ1eO5CMrNS0qxZUiV9HBgHvAYIGCTpixHxSJrFFQMPNptZKWnu1NnXAx+NiLkAkvYG/hdw\nKJiZdSDNnTp7dW0g5MwjmRSvwxs2DLp2TQabV6zIuhozs3QVmjr7DElnADWSHpb0eUmfA/5IEUxa\n1xZqB5vB02ibWcdXqKVweu7RFXgLOJZkmuulQGWqlRURz5hqZqWi0M1rF7RVIcVs5Ei45RZ4+ums\nKzEzS1dzrz7qClwEDCdpNQAQERemVFdROfroZDtlCmzZAmXNHYkxM2tnmvvj7W7gAyQrsf2FZMGc\nkhhoBhg4EPbcE5Ytg9mzs67GzCw9zQ2FIRFxNbA2Nx/Sx4Ej0iuruEhwzDHJ8yefzLYWM7M0NTcU\nNua270o6AOgF7J5OScWptgvpqaeyrcPMLE3NvXntNkl9gKuBiSQrsV2dWlVFyKFgZqWgWaEQEbfn\nnv4FGJxeOW2sBSuxHXgg9OwJr78OixZBVVXKtZmZZaBZ3UeSdpP0C0nTJU2T9DNJu6VdXDEpL4ej\njkqeu7VgZh1Vc8cU7gPeBj4FnAm8A9yfVlHFyl1IZtbRNXdMoX9EfC/v9fclnZNGQcXMoWBmHV1z\nWwqPSjpXUlnucTYwKc3CitHhhydzIb34IqxcmXU1Zmatr9BynKtJluEUcAXw69xbZcAa4KupVtfW\nCgw8V1Ymi+4880zyGD26DWszM2sDhVZe6xERPXPbsoioyD3KIqJnWxVZTNyFZGYdWbNn8ZE0RtJ1\nucdpaRZVzGpDwXc2m1lH1NxLUq8Fvgy8nHt8WdKP0iysWI0alWyfew7efz/bWszMWltzWwofA06K\niDsi4g5gNMn8RyVnt92S1djWr/cSnWbW8bRkEujeec97tXYh7cmxxybbxx/Ptg4zs9bW3FD4EfC8\npDsl3QVMA36QXlnF7eSTk+2kllyU++ePbvswMytCBW9ekyTgKWAkcFhu979HxJI0Cytmxx8PFRXw\n7LPJ/Qq9SrrdZGYdScGWQkQE8HBELI6IiblHyQYCJBPjHXkkbN7sLiQz61ia2300XdJhhQ8rHaec\nkmwffTTbOszMWlNzQ+EI4FlJr0maIeklSTPSLKzY5Y8rRGRbi5lZa2luKJxCso7C8cDpwGm5bZMk\njZY0W9JcSV9v4rhPSQpJ1c2sJ3MjRiSXp86fD3PmZF2NmVnraDIUJHWVdAXwNZJ7ExZFxBu1jwKf\nLQduBk4FhgHnSRrWwHE9SG6M+9sO/jdkorwcTjoped6iq5DMzIpYoZbCXUA18BLJD/frW/DdhwNz\nI2JeRGwgWZNhbAPHfQ/4MbC+Bd9dFDyuYGYdTaFQGBYR50fEL0kW1zmmBd9dBSzIe70wt6+OpBHA\nXhHxv019kaSLJdVIqlm6dGkLSkhXbUth8mTYsCHbWszMWkOhUNhY+yQiNrXmiSWVATcAVxU6NiJu\ni4jqiKju169fa5axU6qq4IADYO1amDIl62rMzHZeoVA4WNKq3GM1cFDtc0mrCnx2EbBX3usBuX21\negAHAP8naT7JzXET29NgM7gLycw6lkLrKZTn1lOoXVOhIu95ofUUpgJDJQ2S1Bk4F5iY990rI6Jv\nRAyMiIHAs8CYiKjZyf+mNlUbCh5sNrOOoCUT4rVIrrvpUpJlO18BHoiIWZK+K2lMWudta8ccA127\nwvPPw5KSvs/bzDqC1EIBICIejoh9ImLviPhBbt81ETGxgWOPa2+tBEgCoXbA+Q9/yLYWM7OdlWoo\nlIpPfSrZTpiQbR1mZjvLodAKxoxJZk2dPBmWLcu6GjOzHedQaEoz10Do0yeZTnvzZnjwwTasz8ys\nlTkUWqKJgHAXkpl1BA6FVvKJT0BZGTz2WLLwjplZe+RQaCW7755cnrpxIzz0UNbVmJntGIdCK3IX\nkpm1dw6FVnTGGcn2T39K5kMyM2tvKrIuoN2qP9h84mSqqmDkSHj2WXjkETjzzGxKMzPbUW4ptDJ3\nIZlZe+ZQaGW1ofDHP8KaNdnWYmbWUg6FVjZoEIwalYwpjB+fdTVmZi3jUEjBhRcm2zvuyLYOM7OW\n8kBza8kbeD6rbyWXdXuYJ5+EOXNg6NAM6zIzawG3FFLQo9s6zj47eX7nnZmWYmbWIg6FlFxwQbK9\n665kojwzs/bAoZCSY46BIUNg0aJkPiQzs/bAoZASaWtrwQPOZtZeOBRS9NnPJjOnPvggLFvVM+ty\nzMwKciikaMAAOPlk2LABfvPECVmXY2ZWkEMhZbX3LIx7aAwR2dZiZlaIQyFlY8dCVRW8/MZAJtUc\nlnU5ZmZNciikrHNnuOyy5PkN48/KthgzswIcCm3g4othl67reGz6YcyYNzjrcszMGuVQaAN9+sCF\nox8B3Fows+LmUGgjV5wxHmkLv5l8AouX7Zp1OWZmDfKEeGmptzLb4P7wyVFP8bunPsJND36SH5yT\nUV1mZk1wS6ENXXXmAwCMe+h01v7x1CQ46i/raWaWIYdCGzpq+CxG7j+L5at78atJo7Mux8xsOw6F\nNva1s+4H4If3fpr31nfJuBozs205FNrYJ49+kkOHzmbx8r7c9OAnsy7HzGwbDoU2JsEPL7wdgGvv\nP4931+yScUVmZlulGgqSRkuaLWmupK838P5XJL0saYakxyV9KM16isVJh9Zw7EEvsGJ1T64ff3bW\n5ZiZ1UktFCSVAzcDpwLDgPMkDat32PNAdUQcBIwHfpJWPcUkv7XwnxPO4q23Mi7IzCwnzZbC4cDc\niJgXERuA+4Cx+QdExOSIeC/38llgQIr1FJWjhs/itJFPs3Z9JT/6UdbVmJkl0gyFKmBB3uuFuX2N\nuQh4JMV6is4PLvhvAG69FebPz7YWMzMokoFmSecD1cBPG3n/Ykk1kmqWLl3atsWl6KDB8/in4//M\nhg1w5ZVZV2Nmlm4oLAL2yns9ILdvG5JOBL4FjImI9xv6ooi4LSKqI6K6X79+qRSblZ984Zd07w5/\n+ANMnJh1NWZW6tIMhanAUEmDJHUGzgW2+bEn6RDglySB8HaKtRStqr7v8P3vJ88vuwzWrs22HjMr\nbamFQkRsAi4FJgGvAA9ExCxJ35U0JnfYT4HuwG8lvSCpJH9XvmS/EzhkyKu8+SZ898J7sy7HzEqY\nop0tHFxdXR01NTU79uEinnzuub/vx8jLb6a8bAvTb72YAwe9vvXNEydnV5iZdQiSpkVEdaHjimKg\n2eDw/f7Ol06byKbNFXzpxivZskVZl2RmJcihUER+cOHt7NFnOVNmHch1v/WCC2bW9hwKRaR397X8\n91XJVbnf+tVFTJ29b8YVmVmpcSgUmY8f8SyXf2ICmzZXcN4Pr2b1e5VZl2RmJcTLcRahH3/hl/xl\nxsG8OG8Il/ziCv6nW70Bcg88m1lK3FIoQl07b+Teb36Pyi7rufvPJ3PP4ydmXZKZlQiHQpHa/0Nv\n8vN//QUAX/zZV3hh7t6t88W160J7fWgza4BDoYhddOrDfObER1m7vpLTrv4R/3hnt6xLMrMOzqFQ\nxCT4ryuv4+gDZrDonX6cfs0PWbuua9ZlmVkH5oHmItel80Z+/52rGXn5LUyfsw+fvvZbTOh8AuXl\nW7Ye5IFnM2slbim0A317reKh732D3t1X8+DTR3PluEtoZ7OTmFk74VBoJ/b74AImXPNtOlVs5Bd/\nOIOv/vJLDgYza3UOhXbk+EOe57dX/wedKjZyw4SzHQxm1uocCu3M2KOmbBMMX7vtXxwMZtZqHArt\nUH4wXD/+HP71X2HTpqyrMrOOwKHQTo09agrjr/4OXTptYNw4OO3I51j14Md9c5qZ7RSHQjs25qin\nmXzdlfTrvYJJNYcz6opf8MZbe2Rdlpm1Yw6Fdu7IYS/z7M8vYf8Pzmfm/MEcfuktTH7hw1mXZWbt\nlEOhAxjcfzFP33gpJ46o4e13d+WEf7ueq++8gE2by9ydZGYt4lDoIHp3X8sjP/x3rjn/LqTg+/d8\nlo985UbmL3F3kpk1n0OhA6ko38J/fO5OnvjJVQzo9zbPvHwAB3/xdm5+8BNs3uz/1WZWmH9SdEDH\nHvwiL477Z844+q+seq87l970ZY788k1MnzM069LMrMgp2tmdT9XV1VFTU7NjHy6xPvUIePDpUVx2\n8+UsXLo7ZWWb+dJpE7n6/LvZo8+Khj+UP7le/T8vT7xn1m5JmhYR1QWPcyh0fGvWdeXbd13Ajb//\nFJu3lLNL13VcccZ4vnrW/fTuvrb5X+RQMGu3mhsK7j4qAd0r13P9v9zK8+O+wOkjp7B2fSU/+M1n\nGPzZ3/D9e85n2aqeWZdoZkXCLYUS9PSs4Xzzjn/mLzOS+xkqu6zn8ydP4sozfsvQAYsa/2ChloK7\nm8yKVnNbCl5kpwQdNXwWk6+7kieeH8F148/mT1OP4NY/jmXcQ6dz0ohpXDj6EcYe9RRdO29MrwgH\niFlRciiUKAlOGDGdE0ZMZ9b8gfznhDP59eMn8ei0w3h02mH06bGK8457gnOOm8yo4TO3Xemtllte\nZh2Ou4+szvJVPbh38gncMelUps/Zp27/Hn2W88lRT3LG0U9yzIEzdqwFUb8l4JaCWZvy1UcNcSg0\n2wtz9+aeJ05kwpMf4fUle9btr+yynuMOeoFTqms4ccQ09v/gG5SVNePvkEPBLFMOhYY4FFosAl54\nbQjj/3osD/3tSGbM23ub93ftsZJRw2dyzIEvMXL/lzlkyBy6V65v+YkcGmapcig0xKGw0xYv25VH\npx3GpJrD+MuMg/jHsn7bvC9tYb+93uTQoa9y8N6vMfxD8zlg4OsM6LcUqYkvLhQKTR1rZgX56iNL\nRf/dlvO5kyfxuZMnEQHzl3yAp2YeyFOzDmTq7H156fXBvPLmQF55cyC/fnzr53p2W8PQqkUMrVrI\nPgMWMmTPRQzqv5iBeyyh/67LKW/NIlvSynCLxGwbqYaCpNHAjUA5cHtEXFvv/S7A/wCHAsuAcyJi\nfpo1WeuRYFD/JQzqv4TPnPQYAOs3dOKl1wczfc4+zJw/kJnzB/HS64NZtqoX0+bsy7Q5+273PRXl\nmxjQbzFVu71DVd/k0X/Xc9ijzwr26L2CPfosp1/vlezWc2XrXCbrFqNZo1ILBUnlwM3AScBCYKqk\niRHxct5hFwErImKIpHOBHwPnpFWTpa9r540ctu9sDtt3dt2+CFj6bm/mLBrAnEVVzFk0gLn/qOKN\nt/bgjbf3YMny3Zi/pD/zl/Qv+P27dF3Hbj2X0Kf7Gvr0WE2f7qvp3X0tvXZZQ89u79Gz21p6dDuN\nHpXv0b1yHd0r17HL379Ety7r6db1fbp1WU9ll12o7Pw+nSo2Nd2lBYUDZGdaFm6lWBFKs6VwODA3\nIuYBSLoPGAvkh8JY4Du55+OBmyQp2ttAhzVJgt37vMvufd5l1AEzt3t/3fudWfROXxa904+F7/Rj\n0Tt9WbJiV95a0afusWxVT95Z1Yu16ytZu76SN9/e+brKyjZT2XkDXTtvoEunDXTptJGunefTuWIT\nnTttpEunjXSuuJ5OFZvoXLGRThWb6VSxiU7lm+hUsZmK8s10+j1UVCSP8vKtj/zXZWXbvle379XT\nkYIybaGsLChbmOwvK0v+zBrb5j/y99X+WTf2qH2/0Lax9/L/fzb3daFj62vq/YIBXsDOfr617Ewd\nffrA4MGtV0tD0gyFKmBB3uuFwBGNHRMRmyStBHYD3kmxLisylV02MKTqHwyp+keTx0XA6ve6sWxV\nT1as6cGK1T1YsaY7767pzup13Vi5dhdWrt2FNesqk8f6Sla/14333u/C2vWVuW1X1r3fhfUbOrNp\nc0VdyGTjKxmd19qrc86B++5L9xztYqBZ0sXAxbmXayTNbur4JvSl/QROe6oV2le97alWaF/1utb0\n9L3/ft65//4d/vyHmnNQmqGwCNgr7/WA3L6GjlkoqQLoRTLgvI2IuA24bWcLklTTnEuyikF7qhXa\nV73tqVZoX/W61vS0Vb1pTp09FRgqaZCkzsC5wMR6x0wEPpd7fibwhMcTzMyyk1pLITdGcCkwieSS\n1DsiYpak7wI1ETER+G/gbklzgeUkwWFmZhlJdUwhIh4GHq6375q85+uBs9KsoZ6d7oJqQ+2pVmhf\n9banWqF91eta09Mm9ba7aS7MzCw9Xo7TzMzqlEwoSBotabakuZK+nnU9jZG0l6TJkl6WNEvSl7Ou\nqRBJ5ZKel/RQ1rUUIqm3pPGS/i7pFUlHZl1TYyRdmfs7MFPSvZK6Zl1TPkl3SHpb0sy8fbtKekzS\nnNy2T5Y11mqk1p/m/h7MkPR7Sb2zrDFfQ/XmvXeVpJDUN41zl0Qo5E25cSowDDhP0rBsq2rUJuCq\niBgGjAQuKeJaa30ZeCXrIprpRuBPEbEfcDBFWrekKuByoDoiDiC5WKPYLsS4Exhdb9/XgccjYijw\neO51MbiT7Wt9DDggIg4CXgW+0dZFNeFOtq8XSXsBJwNvpnXikggF8qbciIgNQO2UG0UnIhZHxPTc\n89UkP7Sqsq2qcZIGAB8Hbs+6lkIk9QI+QnLVGxGxISLezbaqJlUAlbl7eLoBTd/y3cYi4q8kVw3m\nGwvclXt+F/CJNi2qEQ3VGhGPRsSm3MtnSe6lKgqN/NkC/Cfwb0Bqg8GlEgoNTblRtD9oa0kaCBwC\n/C3bSpr0M5K/pA0s4lx0BgFLgV/lurtul7RL1kU1JCIWAdeR/Ea4GFgZEY9mW1Wz7BERi3PPlwB7\nZFlMC1wIPJJ1EU2RNBZYFBEvpnmeUgmFdkdSd2ACcEVErMq6noZIOg14OyKmZV1LM1UAI4BbI+IQ\nYC3F072xjVxf/FiSINsT2EXS+dlW1TK5G1GL/vJGSd8i6ba9J+taGiOpG/BN4JpCx+6sUgmF5ky5\nUTQkdSIJhHsi4ndZ19OEUcAYSfNJuuSOl/TrbEtq0kJgYUTUtrzGk4REMToReD0ilkbERuB3wFEZ\n19Qcb0nqD5DbtsJ8tumR9HngNODTRT6bwt4kvyC8mPv3NgCYLukDrX2iUgmF5ky5URQkiaTP+5WI\nuCHrepoSEd+IiAERMZDkz/SJiCja32YjYgmwQFLtSj8nsO1U7sXkTWCkpG65vxMnUKSD4vXkT13z\nOeDBDGtpUm4RsH8DxkTEe1nX05SIeCkido+Igbl/bwuBEbm/062qJEIhN5hUO+XGK8ADETEr26oa\nNQr4DMlv3S/kHh/LuqgO5DLgHkkzgA8DP8y4ngblWjPjgenASyT/VovqDlxJ9wLPAPtKWijpIuBa\n4CRJc0haO9c29R1tpZFabwJ6AI/l/p2Ny7TIPI3U2zbnLu4Wk5mZtaWSaCmYmVnzOBTMzKyOQ8HM\nzOo4FMzMrI5DwczM6jgUrEORtDl3eeEsSS/mZpQsy71XLennTXx2oKR/artqG6xhnKRRWdZgpc2X\npFqHImlNRHTPPd8d+A0wJSK+3YzPHgd8NSJOS7fKJmt4ATg0IjZnVYOVNrcUrMOKiLeBi4FLlTiu\nds0HScfm3Rz4vKQeJDdaHZPbd2Wu5fCkpOm5x1G5zx4n6f/y1mW4J3fXMZIOk/R0rpXynKQeStab\n+Kmkqbm5+7/YUL2S9gderR8Iku7MtSBqJL2am3Oqdh2L65SstzBD0mWp/WFayUh1jWazrEXEvNx6\nGrvXe+urwCURMSU3+eB6ksnx6loKuUnIToqI9ZKGAvcC1bnPHwIMJ5nOegowStJzwP3AORExVVJP\nYB1wEcksp4dJ6gJMkfRoRLxer6ZTgT818p8ykGQK+L2ByZKGABfk9n84IjZJ2rXFf0Bm9TgUrFRN\nAW6QdA/wu4hYmPtlP18n4CaVn8mAAAABs0lEQVRJHwY2A/vkvfdcRCyEui6fgcBKYHFETAWond1W\n0snAQZLOzH22FzAUqB8Kp5D8oG/IAxGxBZgjaR6wH8k0EuNq1wSIiIbm3zdrEYeCdWiSBpP8QH8b\n2L92f0RcK+l/gY+R/OZ+SgMfvxJ4i2SFtjKS1kSt9/Oeb6bpf0sCLouISU3U2Q3oHRGNLaRTf/DP\ng4GWCo8pWIclqR8wDrip/rTIkvbOzTz5Y5JZdPcDVpNMkFarF8lv/ltIJiksL3DK2UB/SYflztFD\nyappk4Av5aZER9I+2n5xn48Ck5v47rMklUnaGxicO9djwBdz58DdR9Ya3FKwjqYy153TiWThlLuB\nhqYgv0LSR0lWjJtFsurWFmCzpBdJ1si9BZgg6bMkff1rmzpxRGyQdA7wC0mVJOMJJ5IsVTqQZP57\nkaz+Vn+ZylNJZkVtzJvAc0BP4F9y4xy3k3RpzZC0Efgvkpk/zXaYL0k1KwKSpgNH5BbUqf/encBD\nEdFUaJi1CrcUzIpARBTrCnBWYtxSMDOzOh5oNjOzOg4FMzOr41AwM7M6DgUzM6vjUDAzszoOBTMz\nq/P/AfXZgieAXSY+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f539704fef0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def p_decay(dist):\n",
    "    \"\"\"distance should be in pc\"\"\"\n",
    "    halflife = 0.8399 # in s\n",
    "    mass = 8.02248736*amu\n",
    "    E = 1*EeV\n",
    "    \n",
    "    lifetime = halflife #* np.log(2)\n",
    "    gamma = E / (mass * c_squared)\n",
    "    mean_dist = lifetime*gamma*c_light / pc\n",
    "    \n",
    "    return math.exp(-dist/mean_dist)\n",
    "\n",
    "print(0.8399 / np.log(2) *1*EeV / (8.02248736*amu * c_squared)*c_light / pc)\n",
    "\n",
    "n, bins, patches = plt.hist(dists, 80, normed=1, facecolor='orange', alpha=0.75)\n",
    "y = [p_decay(b) for b in bins]\n",
    "_ = plt.plot(bins, y, 'b', linewidth=2)\n",
    "_ = plt.xlabel('Distance / pc')\n",
    "_ = plt.ylabel('Probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some modules have an additional functionality besides `process`. For example, there is a possibility to plot the dependence on energy of the abovementioned decay using the same module but without candidate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XecVOXZ//HPRS/CIl16B5eiwAKi\nxl5QwYZJLNFYIponRpP8HgUVo8YYscXE2IKKhEQxhiYgil0siQIqCyx9aQvIUpcOW67fHzPk2ZDd\nZQZm5szMft+v177cuc/ZmeuW8uU+5Trm7oiIiESqStAFiIhIalFwiIhIVBQcIiISFQWHiIhERcEh\nIiJRUXCIiEhUFBwiIhIVBYeIiERFwSEiIlFRcIiISFSqBV1APDRu3NjbtWsXdBkiIilj7ty5m929\nSST7pmVwtGvXjjlz5gRdhohIyjCz1ZHuq0NVIiISFQWHiIhERcEhIiJRSavgMLMhZja6oKAg6FJE\nRNJWWgWHu09z92EZGRlBlyIikrbSKjhERCT+FBwiImlgx75C3srekJDPSsv7OEREKpP3czZy75T5\nbN19gD5tG3BcRu24fp6CQ0QkRW3ZtZ8HpuUwbd56ujWvx4vXZcU9NEDBISKSctydqfPW88DUheza\nX8Svzu3Crad3pEa1xJx9UHCIiKSQ9dv3MnLKAj5cnE/vNg14bGgvOjerl9AaUiI4zKwu8AnwgLtP\nD7oeEZFEKylxXvtqDaPeXkxxifPrwZn8+OR2VK1iCa8lkOAwszHAYCDf3XuUGh8E/BGoCrzk7qPC\nm4YDbyS8UBGRJLBy825GTMzmy5VbObVTYx65vCetG9YJrJ6gVhxjgWeAcQcHzKwq8CxwLpAHzDaz\nqUBLIAeolfgyRUSCU1RcwsufreT37y2lRrUqPDa0F9/PaoVZ4lcZpQUSHO4+y8zaHTLcH1ju7rkA\nZvY6cAlwDFAXyAT2mtkMdy9JYLkiIgm3aMMO7pqQzfx1BZyX2YyHLu1Bs/rJ8e/nZDrH0RJYW+p1\nHjDA3W8DMLPrgc3lhYaZDQOGAbRp0ya+lYqIxMn+omKe+XA5z3+8ggZ1qvPs1X24sGfzwFcZpSVT\ncFTI3cceZvtoYDRAVlaWJ6ImEZFYmrt6G8MnZrM8fxeX927JfYMzObZujaDL+i/JFBzrgNalXrcK\nj0XMzIYAQzp16hTLukRE4mr3/iKeeHcJY79YxXH1a/HKDf04s2vToMsqVzIFx2ygs5m1JxQYVwJX\nR/MG7j4NmJaVlXVzHOoTEYm5z5ZtZsSkbPK27eW6gW25a1A3jqmZTH81/7egLscdD5wBNDazPOB+\nd3/ZzG4DZhK6HHeMuy+M8n214hCRlFCwt5CH38rhjTl5tG9clzduGUj/9g2DLisi5p5+pwOysrJ8\nzpw5QZchIlKmdxd+x8gpC9iy+wDDTuvAHWd3plb1qoHWZGZz3T0rkn2Tez0UJa04RCSZbdq5nwem\nLeSt7A1kHlefMdf3o0fL1HvwnFYcIiJx5u5M/mYdv5mew579xdx+diduOb0j1asmzyORKu2KQ0Qk\n2azbvpd7J8/n4yWb6Nv2WB4d2otOTY8JuqyjklbBoUNVIpIsSkqcV79czai3F+PAA0MyuXZgME0J\nYy2tgkOX44pIMlixaRcjJmYze9U2vte5Mb+7LNimhLGWVsEhIhKkwuISXvw0lz+8v4xa1arw+BW9\nuKJv8E0JYy2tgkOHqkQkKAvWFTB8YjYL1+9gUPfm/ObS7jStlxxNCWMtrYJDh6pEJNH2FRbzpw+X\n8cInuRxbpwbPXdOHC3seF3RZcZVWwSEikkhzVm3lronZ5G7azdA+rbhv8PE0qJN8TQljTcEhIhKl\n3fuLeHzmEv7yz1W0yKjNuBv7c1qXJkGXlTBpFRw6xyEi8fbJ0k3cM2k+6wv2ct1JoaaEdZO8KWGs\npdVsdY5DROJl+54DPDR9ERO/zqNDk7r845aBZLVLjaaEsZZWwSEiEg9vz9/AfW8uZNueA/zszI78\n/KzgmxIGScEhIlKO/B37+PWbC3ln4Xd0b1Gfv9zYj+4tUq8pYawpOEREDuHuTJibx0PTc9hXVMLw\nQd24+XvtqZZETQmDlFbBoZPjInK01m7dwz2T5/Ppss30a3cso4b2omOT1G5KGGtqqy4iQqgp4bh/\nruKxmUswYMQF3bhmQFuqpEFTwkiorbqISBSW5+9i+MRs5q7exuldmvC7y3vSskHtoMtKWgoOEam0\nCotLGD0rlz++v4w6Navy+x+cwGW9W6ZdU8JYU3CISKW0YF0Bd07IZtGGHVzU6zgeGNKdJvVqBl1W\nSlBwiEilsq+wmD+8v4wXP82lYd0a/PnavpzfvXnQZaUUBYeIVBpfrdzKiInZ5G7ezQ+zWnPPhceT\nUad60GWlnLQKDl2OKyJl2bW/iMfeWcy4f66m1bG1+dtNAzi1c+Ogy0pZaRUc6lUlIof6aEk+906a\nz4Yd+7jxlPb87/ldqFMjrf7qSzj93xORtLRt9wEemp7DpG/W0anpMUy49WT6tj026LLSgoJDRNKK\nu/PW/A3c/+ZCCvYWcvtZnfjZWZ2oWa3yNiWMNQWHiKSN/B37GDllAe/mbKRnywz+etMAMlvUD7qs\ntKPgEJGU5+78Y04eD72Vw4GiEu6+oBs3naqmhPGi4BCRlLZmyx7unpzN58u30L99Qx4d2ov2jesG\nXVZaU3CISEoqLnHGfrGKJ2YuoWoV47eX9uDq/m0qTVPCICV9cJjZ8cAdQGPgA3d/PuCSRCRgyzbu\n5K6J2XyzZjtndm3Cw5f1pIWaEiZMIAcAzWyMmeWb2YJDxgeZ2RIzW25mIwDcfZG73wr8ADgliHpF\nJDkcKCrh6Q+WceHTn7Jq827+eOWJjLm+n0IjwYJacYwFngHGHRwws6rAs8C5QB4w28ymunuOmV0M\n/BT4awC1ikgSmLd2O8MnZrP4u50MOaEFDwzJpNExakoYhECCw91nmVm7Q4b7A8vdPRfAzF4HLgFy\n3H0qMNXM3gJeS2StIhKsvQeKeer9pbz0aS5N6tXkxeuyODezWdBlVWrJdI6jJbC21Os8YICZnQFc\nDtQEZpT3w2Y2DBgG0KZNm/hVKSIJ888VW7h7Ujartuzhqv6tufvC46lfS00Jg5ZMwVEmd/8Y+DiC\n/UYDoyH06Nj4ViUi8bRjXyGj3l7Ma1+uoW2jOrx28wBO7qimhMkimYJjHdC61OtW4bGIqTuuSOr7\ncPFG7pm0gPyd+/jJqe35f+d1pXYNtQtJJuUGh5k1jODnS9x9e4xqmQ10NrP2hALjSuDqaN5A3XFF\nUtfW3Qd4cNpC3vx2PV2b1eOFa/tyYusGQZclZahoxbE+/FXR3TRVgahPKJjZeOAMoLGZ5QH3u/vL\nZnYbMDP8vmPcfWGU76sVh0iKcXemzlvPg9Ny2LmvkDvO7szPzuxEjWpqF5KszL3s0wFm9o27967w\nhyPYJwhZWVk+Z86coMsQkcPYULCX+6Ys4P1F+ZzQKoNHr+hFt+ZqShgEM5vr7lmR7FvRimNgBD8f\nyT4JoxWHSGpwd8Z/tZZHZiyisKSEkRcdzw2ntKeq2oWkhHLXgu6+D8DMTjKzegfHzay+mQ0ovU+y\ncPdp7j4sIyMj6FJEpByrt+zm6he/5J7J8+nRMoOZvziNn3yvg0IjhURyVdXzQJ9Sr3eVMSYiUqHi\nEmfMZyt58r0lVK9ShVGX9+SH/VpjpsBINZEEh3mpEyHuXmJmyXQZ77/pUJVIclry3U7umjCPeXkF\nnHN8U357aU+aZ9QKuiw5QpFctpBrZrebWfXw1x1AbrwLOxI6VCWSXPYXFfPUe0sZ/KdPydu2lz9d\n1ZsXr8tSaKS4SFYOtwJPAyPDr98n3NpDRKQ836zZxvCJ2SzduIvLerfkvsGZNKxbI+iyJAYOGxzu\nnk/oZrykp0NVIsHbc6CIJ99dypjPV9K8fi3GXJ/FWd3UlDCdHPZQlZl1MLNpZrYp/AyNN82sQyKK\ni5YOVYkE64vlmxn0h095+bOVXDOgDe/+8jSFRhqK5FDVa4Sek3FZ+PWVwHhgQLyKEpHUUrC3kEdm\nLOL12Wtp16gOfx92EgM6NAq6LImTSIKjjruXfoDS38zszngVJCKp5b2cjYycMp9NO/dzy+kd+OU5\nXahVXU0J01kkwfF2+DGurwMO/BCYcbAJortvjWN9UdE5DpHE2bxrPw9MXcj07A10a16PF6/Lolcr\nNSWsDMrtVfXvHcxWVrDZ3T3pzneoV5VI/Lg7b367ngenLWT3/mJ+flYnbjm9o5oSprhY9aoCwN3b\nH31JIpIO1m/fy8gpC/hwcT692zTgsaG96Nys3uF/UNJKRc/j6OPuX1f0w5HsIyKpr6TEee2rNYx6\nezHFJc6vB2fy45Pbqb9UJVXRiuOV8PO+K/qd8TKQNG3VdY5DJPZWbt7N8InZfLVyK6d0asSoy3vR\numGdoMuSAFUUHBnAXCoOjk2xLefo6AmAIrFTVFzCS5+t5Kn3llKjWhUeG9qL72e1UlNCKT843L1d\nAusQkSSSs34HwydmM39dAedmNuO3l/agWX31l5KQpOxyKyLB2F9UzLMfLue5j1fQoE51nr26Dxf2\nbK5VhvwHBYeIADB3dagp4fL8XVwebkp4rJoSShkUHCKV3J4DRTw+cwljv1hFi4zajL2hH2d0bRp0\nWZLEIgoOM2sJtC29v7vPildRIpIYny3bzIhJ2eRt28u1J7Vl+AXdOKam/j0pFTvs7xAze5RQm5Ec\noDg87EDSBYcuxxWJTMGeQh6ekcMbc/Lo0Lgub9wykP7tGwZdlqSISFqOLAF6ufv+xJR09NRyRKR8\n7yz4jvveXMDW3Qe45bQO3H52ZzUllNi2HCH0mNjqQMoEh4j8t007Q00J35q/gczj6vPK9f3o0VLP\nrpHoVdRy5E+EDkntAb41sw8oFR7ufnv8yxORo+XuTPp6Hb+ZnsPeA8XceX5Xhp3WgepV1ZRQjkxF\nK46Dx3rmAlMP2Vbx8S0RSQrrtu/lnknz+WTpJvq2PZZHh/aiU9Njgi5LUlxFd47/BcDM7nD3P5be\nZmZ3xLswETlyJSXO375czaNvL8aBB4Zkct3AdlRRU0KJgUjOcfwY+OMhY9eXMSYiSWDFpl2MmJjN\n7FXb+F7nxvzusp5qSigxVdE5jquAq4H2Zlb6UFU9IGme+iciIYXFJbz4aS5/eH8ZtatX5Ynvn8DQ\nPi3VLkRirqIVxxfABqAx8GSp8Z1AdjyLKs3MLgUuAuoDL7v7u4n6bJFUsWBdAcMnZrNw/Q4u7Nmc\nBy7uTtN6akoo8VHROY7VwGpgYKw/1MzGAIOBfHfvUWp8EKFDYFWBl9x9lLtPAaaY2bHAE4CCQyRs\nX2ExT3+wjD/PyuXYOjV44Ud9GNTjuKDLkjQXyZ3jJwF/Ao4HahD6S323u9c/is8dCzwDjCv1OVWB\nZ4FzgTxgtplNdfec8C4jw9tFBJizait3Tcwmd9Nuvt+3FSMvyiSjTvWgy5JKIJKT488AVwL/ALKA\n64AuR/Oh7j7LzNodMtwfWO7uuQBm9jpwiZktAkYBb+sxtSKwa38Rj7+zmHH/Wk2LjNqMu7E/p3Vp\nEnRZUolE1M3M3ZebWVV3Lyb0SNlvgLtjXEtLYG2p13nAAODnwDlAhpl1cvcXyvphMxsGDANo06ZN\njEsTSQ6zlm7i7knzWV+wlx8PbMed53elrpoSSoJF8jtuj5nVIHT3+GOETpgn7JZTd38aeDqC/UYD\noyHUqyredYkk0vY9B3ho+iImfp1HxyZ1mXDrQPq2VVNCCUYkwXEtoaC4Dfgl0BoYGoda1oXf+6BW\n4bGIqTuupKO352/gvjcXsn3PAW47sxO3ndVJTQklUIcNDndfHV5xtAMmAUvc/UAcapkNdDaz9oQC\n40pC95FEzN2nAdOysrJujkN9IgmVv2Mfv35zIe8s/I4eLevzlxv70b2FmhJK8CK5quoi4AVgBWCE\nbgi8xd3fPtIPNbPxwBlAYzPLA+5395fN7DZgJqErt8a4+8Io31crDkl57s6EuXk8ND2HfUUljLig\nGz85tT3V1JRQkkQkz+NYDAx29+Xh1x2Bt9y9WwLqOyJ6HoekqrVb93DP5Pl8umwz/dqFmhJ2aKKm\nhBJ/sX4ex86DoRGWS+ju8aSjFYekquISZ9w/V/H4zCUY8NAl3blmQFs1JZSkFMmK43lCzxt/g1A7\n9e8Da4D3Adx9UpxrjJpWHJJKlufv5K4J2Xy9Zjund2nC7y7vScsGtYMuSyqZWK84agEbgdPDrzcB\ntYEhhIIk6YJDJBUUFpfw509W8PQHy6lTsyq//8EJXNZbTQkl+UVyVdUNiSgkFnSoSlLF/LwC7pww\nj8Xf7eSinsfxwMXdaVKvZtBliUQkrW451eW4kuz2FRbzh/eX8eKnuTSqW4M/X9uX87s3D7oskaik\nVXCIJLMvc7cwYtJ8Vm7ezQ+zWnPPRceTUVtNCSX1pFVw6FCVJKOd+wp59J3F/O1fa2jdsDav/mQA\np3RqHHRZIkfssHcUmVkzM3vZzN4Ov840s5viX1r03H2auw/LyNDdtZIcPlqcz/lPzeLVL9dw4ynt\nmfmL0xQakvIiWXGMBV4B7g2/Xgr8HXg5TjWJpLytuw/w0PQcJn+zjs5Nj2HiT0+mT5tjgy5LJCYi\nCY7G7v6Gmd0N4O5FZlYc57pEUpK789b8Ddz/5kIK9hZy+9md+dmZHalZTU0JJX1EEhy7zawRoXs2\nDj4RsCCuVR0hneOQIG3csY+RUxbwXs5GerXK4NWbB9Ct+dE8KFMkOUUSHL8CpgIdzexzoAlwRVyr\nOkK6HFeC4O68MWctv31rEQeKSrj7gm7cpKaEksYiuQHwazM7HehKqDvuEncvjHtlIilgzZY9jJiU\nzRcrtjCgfUNGDe1F+8Z1gy5LJK4iaateh9Cqo62732xmnc2sq7tPj395IsmpuMQZ+8Uqnpi5hKpV\njIcv68FV/dqoKaFUCpEcqnoFmAsMDL9eB/wDUHBIpbR0Y6gp4bdrt3NWt6Y8fFkPjstQU0KpPCIJ\njo7u/kMzuwrA3fdYknZh08lxiacDRSU8//EKnvloGfVqVeePV57IxSe0UFNCqXQiCY4DZlab/7uq\nqiOwP65VHSGdHJd4mbd2O8MnZrP4u51cfEIL7h+SSaNj1JRQKqdIguN+4B2gtZm9CpwCXB/PokSS\nxd4Dxfz+vSW8/NlKmtarxUvXZXFOZrOgyxIJVIXBET4ktRi4HDiJ0FVVd7j75gTUJhKof67YwohJ\n2azesoer+rfh7gu7Ub+WmhKKVBgc7u5mNsPdewJvJagmkUDt2FfIIzMWM/6rNbRtVIfXbh7AyR3V\nX0rkoEgOVX1tZv3cfXbcqxEJ2AeLNnLv5AXk79zHsNM68MtzulC7htqFiJQWSXAMAK4xs9XAbkKH\nq9zde8W1MpEE2rJrPw9Oy2HqvPV0bVaPF67ty4mtGwRdlkhSKjc4zKy9u68Ezk9gPUdFl+NKtNyd\nqfPW8+C0HHbuK+SX53Thp2d0pEY1tQsRKY+5e9kbzOa6e18z+8Ddz05wXUclKyvL58yZE3QZkuQ2\nFOxl5OQFfLA4nxNaN+Cxob3o2rxe0GWJBCL8d35WJPtWdKiqipndA3Qxs18dutHdf3+kBYoEqaTE\neX32Wh6ZsYjCkhJGXnQ8N5zSnqpqFyISkYqC40rg0vA++meYpIVVm3czYlI2/8rdysAOjRg1tCdt\nG6kpoUg0yg0Od18CPGpm2e7+dgJrEom54hJnzGcrefK9JVSvUoVRl/fkh/1aq12IyBGIpK26QkNS\n2uLvdjB8Qjbz8go45/hm/PbSHjTPqBV0WSIpK5LLcUVS0v6iYp79aAXPfbScjNrVefqq3gzpdZxW\nGSJHScEhaembNdu4a0I2y/J3cemJLfj1kO40rFsj6LJE0kJEwWFmJwPtSu/v7uPiVNOhn90BuBfI\ncPekfGStJI89B4p48t2ljPl8Jc3r1+KV6/txZremQZclklYieQLgX4GOwLdAcXjYgSMODjMbAwwG\n8t29R6nxQcAfgarAS+4+yt1zgZvMbMKRfp5UDl8s38yISfNZs3UPPzqpDcMHdaOemhKKxFwkK44s\nINPLu1PwyIwFnqFU+JhZVeBZ4FwgD5htZlPdPSeGnytpqGBvIY/MWMTrs9fSvnFd/j7sJAZ0aBR0\nWSJpK5LgWAA0BzbE6kPdfZaZtTtkuD+wPLzCwMxeBy4BFBxSrncXfsfIKQvYvGs/t5weakpYq7qa\nEorEUyTB0RjIMbOvKPXkP3e/OMa1tATWlnqdBwwws0bAw0BvM7vb3R8p64fNbBgwDKBNmzYxLk2S\nzeZd+3lg6kKmZ2+gW/N6vPTjLHq1UlNCkUSIJDgeiHcRFXH3LcCtEew3GhgNoV5V8a5LguHuTPl2\nHQ9Oy2HP/mL+97wu3HJ6R6pXVVNCkUSJ5AbATxJRCLAOaF3qdavwWMTUHTe9rd++l3snz+ejJZvo\n3SbUlLBzM3XDEUm0w/4zzcxOMrPZZrbLzA6YWbGZ7YhDLbOBzmbW3sxqEOqVNTWaN3D3ae4+LCMj\nIw7lSVBKSpy//ms15z01i3/lbuXXgzOZcOvJCg2RgERyqOoZQn+J/4PQFVbXAV2O5kPNbDxwBtDY\nzPKA+939ZTO7DZhJ6HLcMe6+MMr31YojzeRu2sWISfP5auVWTu3UmEcu70nrhnWCLkukUiv3eRz/\n3sFsjrtnhZsd9gqPfePuvRNS4RHQ8zhSX1FxCS99tpKn3ltKzWpVGDk4k+/3baV2ISJxEqvncRy0\nJ3zo6Fsze4zQZblJeSZSK470kLN+B8MnZjN/XQHnZYaaEjatr6aEIskikhVHW2AjUAP4JZABPOfu\ny+Nf3pHRiiM17S8q5pkPl/P8xytoUKc6v7mkBxf0aK5VhkgCxHTF4e6rzaw2cJy7P3jU1YmUYe7q\nbQyfmM3y/F0M7dOK+wYfT4M6akookowi6VU1BHiC0IqjvZmdCPwmDjcAHjUdqko9u/cX8cS7Sxj7\nxSpaZNRm7A39OKOrmhKKJLNIbwDsD3wM4O7fmln7ONZ0xNx9GjAtKyvr5qBrkcP7bNlmRkzKJm/b\nXn48sC13DurGMTXV6V8k2UXyp7TQ3QsOOc6sO7PliBXsKeThGTm8MSePDk3q8o9bB9KvXcOgyxKR\nCEUSHAvN7Gqgqpl1Bm4HvohvWUdGh6qS38xwU8Ktuw/wP2d05PazO6spoUiKieSy2p8D3Qk1OBwP\n7AB+Ec+ijpTuHE9e+Tv38T+vzuWWv86lyTE1efNnp3DXoG4KDZEUFMlVVXsIPYHv3viXI+nG3Zn4\n9Toemp7D3sJi7jy/K8NO66CmhCIprNzgMLMK+0Ql41VVklzytu3hnskLmLV0E33bHsujQ3vRqekx\nQZclIkepohXHQELPxxgPfAkk/V1YOseRHA42JXz0ncUAPHhxd649qS1VqiT9byERiUC5d46HH+V6\nLnAV0At4CxgfbePBIOjO8eCs2LSL4ROymbN6G9/r3JjfXaamhCKpICZ3jrt7MfAO8I6Z1SQUIB+b\n2YPu/kxsSpV0UVhcwuhZufzxg2XUrl6VJ75/AkP7tFS7EJE0VOHJ8XBgXEQoNNoBTwOT41+WpJIF\n6wq4a0I2ORt2cEGP5jx4SXea1lNTQpF0VdHJ8XFAD2AG8KC7L0hYVZIS9hUW8/QHy/jzrFwa1q3B\nCz/qw6AexwVdlojEWUUrjh8Bu4E7gNtLHXIwwN29fpxri5pOjifO7FVbGT4hm9zNu7mibyvuuyiT\njDrVgy5LRBKgonMcKXehvXpVxd+u/UU8/s5ixv1rNS0b1Gbcjf05rUuToMsSkQRSRzmJ2MdL8rl3\n8gLWF+zlxwPbcef5XamrpoQilY7+1Mthbdt9gIfeymHS1+vo2KQuE24dSN+2akooUlkpOKRc7s7b\nC77j128uYPueQm47sxM/P7sTNaupv5RIZabgkDLl79jHfW8uYObCjfRsmcG4GweQ2SLprocQkQCk\nVXDoqqqj5+78Y24ev52ew/6iEkZc0I2fnNqeampKKCJhaRUcuqrq6Kzduod7Js/n02Wb6d+uIaOG\n9qRDEzUlFJH/lFbBIUemuMQZ989VPPbOEqoYPHRpD67p30ZNCUWkTAqOSm55/k7umpDN12u2c0bX\nJjx8WU9aNqgddFkiksQUHJVUYXEJf/5kBU9/sJw6Navy1A9P4NIT1ZRQRA5PwVEJzc8r4M4J81j8\n3U4G9zqOBy7uTuNjagZdloikCAVHJbKvsJin3l/Ki7NyaXxMTUZf25fzujcPuiwRSTEKjkriy9wt\njJg0n5Wbd3Nlv9bcfeHxZNRWU0IRiV7SB4eZ1QWeAw4AH7v7qwGXlFJ27ivk0XcW87d/raF1w9q8\n+pMBnNKpcdBliUgKC+SuLjMbY2b5ZrbgkPFBZrbEzJab2Yjw8OXABHe/Gbg44cWmsI8W53P+U7N4\n7cs1/OTU9sz8xWkKDRE5akGtOMYCzwDjDg6En3H+LKHnnOcBs81sKtAKmB/erTixZaamrbsP8Jtp\nC5ny7Xo6Nz2GiT89md5tjg26LBFJE4EEh7vPMrN2hwz3B5a7ey6Amb0OXEIoRFoB3xLQCilVuDvT\nszfwwNSFFOwt5I6zO/M/Z3ZUU0IRialkOsfRElhb6nUeMIDQc86fMbOLgGnl/bCZDQOGAbRp0yaO\nZSanjTv2ce/kBby/aCMntMrg1ZsH0K25mhKKSOwlU3CUyd13AzdEsN9oYDRAVlaWx7uuZOHu/H32\nWh6esYjC4hLuvfB4bjy1PVXVLkRE4iSZgmMd0LrU61bhsYhVtu64a7bsYcSkbL5YsYWTOjRk1OW9\naNe4btBliUiaS6bgmA10NrP2hALjSuDqaN6gsnTHLS5xXvl8JU+8u4TqVarw8GU9uKqfmhKKSGIE\nEhxmNh44A2hsZnnA/e7+spndBswEqgJj3H1hlO+b9iuOpRt3cueEbOat3c7Z3Zry28t6cFyGmhKK\nSOKYe/qdDsjKyvI5c+YEXUZMHSgq4fmPV/DMR8uoV6s69w/J5OITWqgpoYjEhJnNdfesSPZNpkNV\nUo55a7dz14RslmzcycUntOCXl7kxAAAKmElEQVT+IZk0UlNCEQlIWgVHuh2q2nugmN+/t4SXP1tJ\n03q1eOm6LM7JbBZ0WSJSyaVVcKTTyfEvVmzm7knzWb1lD1f1b8PdF3ajfi01JRSR4KVVcKTDimPH\nvkIembGY8V+toW2jOrx28wBO7qj+UiKSPNIqOFJ9xfHBoo3cO3kB+Tv3Mey0DvzynC7UrqF2ISKS\nXNIqOFLVll37+c30HN78dj3dmtfjz9f25YTWDYIuS0SkTGkVHKl2qMrdmTpvPQ9Oy2HnvkJ+eU4X\nfnpGR2pUUy9HEUleaRUcqXSoakPBXkZOXsAHi/M5sXUDHruiF12a1Qu6LBGRw0qr4EgFJSXO+Nlr\neGTGYopKShh50fHccIqaEopI6lBwJNCqzbsZMSmbf+Vu5ZROjXjksl60aVQn6LJERKKSVsGRrOc4\niopLGPP5Sp58dyk1qlXh0aE9+UFWa7ULEZGUlFbBkYznOBZ/t4PhE7KZl1fAOcc347eX9qB5Rq2g\nyxIROWJpFRzJZH9RMc9+tILnPlpORu3q/Omq3gzudZxWGSKS8hQccfDNmm0Mn5jN0o27uKx3S+4b\nnEnDujWCLktEJCYUHDG050ART767lDGfr6R5/Vq8cn0/zuzWNOiyRERiKq2CI8iT458v38yISdms\n3bqXa09qy12DulJPTQlFJA2lVXAEcXK8YG8hj8xYxOuz19K+cV3+PuwkBnRolKiPFxFJuLQKjkR7\nL2cjI6fMZ/OuA9x6ekd+cU5nalVXU0IRSW8KjiOwedd+Hpi6kOnZG+jWvB4vXdePnq0ygi5LRCQh\nFBxRcHemfLuOB6flsGd/Mf97XhduOb0j1auqKaGIVB4Kjgit276XeyfP5+Mlm+jdpgGPDe1FZzUl\nFJFKSMFxGCUlzqtfrmbU24spcbh/SCbXDWynpoQiUmmlVXDE+nLc3E27GDFxPl+t2sqpnRrzyOU9\nad1QTQlFpHJLq+CI1eW4RcUlvPjpSp56fym1qlXhsSt68f2+rdQuRESENAuOWFi4voDhE7NZsG4H\n52WGmhI2ra+mhCIiByk4wgqLS/jD+0t54ZNcjq1Tneeu6cMFPZprlSEicggFR1i1Ksa3a7dz6Ykt\nuW/w8TSoo6aEIiJlUXCEmRljru9HzWq681tEpCK6c60UhYaIyOEpOEREJCpJHxxm1sHMXjazCUHX\nIiIicQ4OMxtjZvlmtuCQ8UFmtsTMlpvZiIrew91z3f2meNYpIiKRi/fJ8bHAM8C4gwNmVhV4FjgX\nyANmm9lUoCrwyCE/f6O758e5RhERiUJcg8PdZ5lZu0OG+wPL3T0XwMxeBy5x90eAwfGsR0REjl4Q\n5zhaAmtLvc4Lj5XJzBqZ2QtAbzO7u4L9hpnZHDObs2nTpthVKyIi/yHp7+Nw9y3ArRHsNxoYDZCV\nleXxrktEpLIKIjjWAa1LvW4VHjtqB7vjAjvMLB8oKLU5o4LXpb9vDGyOQTmHft7R7FvW9kjGUnXO\n5W3TnMse05zLnnOs5lteTUeyX6zmHI9f47YR7+nucf0C2gELSr2uBuQC7YEawDygexw+d3Skrw/5\nfk48Pv9o9i1reyRjqTrn8rZpzppzNHOO1XyjmfOR/Fk+kjnH+9f4cF/xvhx3PPBPoKuZ5ZnZTe5e\nBNwGzAQWAW+4+8I4fPy0KF4fui0en380+5a1PZKxVJ1zeds057LHNOfkmfOR/FkubzzSOcZjvhWy\ncEpJmJnNcfesoOtIJM25cqhsc65s84XEzTnp7xwPwOigCwiA5lw5VLY5V7b5QoLmrBWHiIhERSsO\nERGJioJDRESiouAQEZGoKDjKUFYrdzP7npm9YGYvmdkXQdYXD+XMuY2ZTQl3Oa6wi3EqKmfOmWb2\nhpk9b2ZXBFlfPJjZpWb2opn93czOC4/VNbO/hMevCbrGWCtnzmn9uIZy5vxfY0csETeLJMMXMAbI\np9TNiOHxQcASYDkw4pBtE8p4n0uBW4KeTyLmDFwE/Cj8/d+Dnk+C5vz/gO+Fv58a9HziOOdjgZfD\n318LDKkEv87/nnNZv/bJ/hXDOf/XWNS1BP0/I4H/008D+vCfd7FXBVYAHfi/u9gzK/pNBbwB1At6\nPomYM9AI+Aj4ELgh6PkkaM5NCbX9fxz4POj5xHHOTwJ9wt/fDZwY/v61oOeTiDmX9Wuf7F8xnPN/\njUX7VWkOVbn7LGDrIcP/bvHu7geA14FLynsPM2sDFLj7zvhVGjsxmPMNwP3ufhah1UfSO9o5u3u+\nu/8MGEHs+hzFVTRztpBHgbfd/evwvnmEesZBihy+jsGcU87RzjmW/x9S4jdJHJXZ4r2CVu43Aa8k\nssA4iGbO7wC3h8dXJbbMmIp4zmbWzsxGE3r42OOJLzVmynt8wc+Bc4ArzOxg1+lJwFAze54A2lfE\nUMRzjvRxDSkgml/nssaOSNK3VQ+Cl9PK3d3vD6CchChrzu6+AEi7E8QHlTPnVcCwQApKAHd/Gnj6\nkLHdhFaXaamcOUf0uIZUVc6c/2vsSFX2FUfcWrwnMc1Zc05XmnOC5lzZg2M20NnM2ptZDeBKYGrA\nNcWb5qw5pyvNOUFzrjTBEXCL90BozpozmrPmHI9awpdniYiIRKTSrDhERCQ2FBwiIhIVBYeIiERF\nwSEiIlFRcIiISFQUHCIiEhUFh1RqZlZsZt+W+kqK546Y2Sozm29mWRXs8+Pwtf2lxxqb2SYzq2lm\nr5rZ1nR8rogES72qpLLb6+4nxvINzaxa+Maso3Wmu1fUoXcy8KSZ1XH3PeGxK4Bp7r4fuMbMxsag\nDpH/oBWHSBnC/+J/0My+Dv/Lv1t4vK6Fnoj4lZl9Y2aXhMevN7OpZvYh8IGZVTGz58xssZm9Z2Yz\nzOwKMzvLzKaU+pxzzWxyBPX0NbNPzGyumc00s+PcfQfwCTCk1K5XAuPLfheR2FBwSGVX+5BDVT8s\ntW2zu/cBngf+Nzx2L/Chu/cHzgQeN7O64W19gCvc/XTgcqAdkEnoCXsDw/t8BHQzsybh1zcQerJb\nucysOvCn8Hv3De//cHjzeEJhgZm1ALoQevCWSNzoUJVUdhUdqpoU/u9cQkEAcB5wsZkdDJJaQJvw\n9++5+8EH7ZwK/MPdS4DvzOwjAHd3M/sr8CMze4VQoFx3mBq7Aj2A98wMQk992xDe9hbwnJnVB34A\nTHT34sNNWuRoKDhEyrc//N9i/u/PigFD3X1J6R3NbACwO8L3fYXQA5P2EQqXw50PMWChuw88dIO7\n7zWzd4DLCK08fhVhDSJHTIeqRKIzE/i5hf/pb2a9y9nvc0JP1atiZs2AMw5ucPf1wHpgJJE9UXIJ\n0MTMBoY/s7qZdS+1fTyhwGhGqHuqSFwpOKSyO/Qcx6jD7P8QUB3INrOF4ddlmUjoMZ45wN+Ar4GC\nUttfBda6+6LDFRh+lvQVwKNmNg/4Fji51C7vAS2Av7vaXUsCqK26SJyY2THuvsvMGgFfAae4+3fh\nbc8A37j7y+X87Cog6zCX40ZSw1hgurtPOJr3ESlNKw6R+JluZt8CnwIPlQqNuUAvQiuR8mwidFlv\nuTcAHo6ZvQqcTuhcikjMaMUhIiJR0YpDRESiouAQEZGoKDhERCQqCg4REYmKgkNERKKi4BARkaj8\nfxZn7TLOjWHRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f53971810f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pid = nucleusId(8, 3)\n",
    "mass = nuclearMass(pid)\n",
    "energies = 10**np.linspace(17, 22, 100)*eV\n",
    "meanpath = np.zeros(100)\n",
    "module = NuclearDecay()\n",
    "\n",
    "for i, E in enumerate(energies):\n",
    "    gamma = E / (mass * c_squared)\n",
    "    meanpath[i] = module.meanFreePath(pid, gamma)\n",
    "\n",
    "_ = plt.plot(energies / eV, meanpath / pc)\n",
    "_ = plt.xlabel('Energy [eV]')\n",
    "_ = plt.ylabel('Mean free path [pc]')\n",
    "_ = plt.loglog()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ModuleList and sources\n",
    "\n",
    "For the sake of transparency and simplicity the first two simulation examples use the explicit for-loop to iterate modules, but it is better to delegate this main loop to CRPropa itself due to performance reasons and some additional features. The **module iterator**, named ``ModuleList``, operates basically like above for-loops: desired modules are added to ``ModuleList`` and then a simulation is being run in its run method.\n",
    "\n",
    "Here is an example with several other available modules: ``PhotoPionProduction``, ``ElectronPairProduction`` and ``MaximumTrajectoryLength``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Charge 2.0 e, Energy: 0.49889470221750665 EeV\n"
     ]
    }
   ],
   "source": [
    "pid = nucleusId(8, 3)\n",
    "E = 1*EeV\n",
    "c = Candidate(pid, E)\n",
    "c.current.setPosition(Vector3d(0, 0, 0)*Mpc)\n",
    "c.current.setDirection(Vector3d(1, 0, 0))\n",
    "\n",
    "sim = ModuleList()\n",
    "sim.add(SimplePropagation())\n",
    "sim.add(NuclearDecay())\n",
    "sim.add(PhotoPionProduction(CMB))\n",
    "sim.add(ElectronPairProduction(CMB))\n",
    "sim.add(MaximumTrajectoryLength(10*Mpc))\n",
    "\n",
    "sim.run(c)\n",
    "\n",
    "print(\"Charge {} e, Energy: {} EeV\"\n",
    "      .format(c.current.getCharge() / eplus, c.current.getEnergy() / EeV))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One should take care that there is at least one stopping condition, like ``MaximumTrajectoryLength`` module here, otherwise the loop would never end.\n",
    "\n",
    "Manually preparing and inserting individual candidates into the simulation box, like in above examples, is sometimes useful for specific purposes, but generally an automatic insertion of candidates according to specified properties is a more common scenario. This mechanism is provided by **sources of particles** which are described by ``Source`` class and its derivated classes. ``Source`` accepts different properties, refered as ``SourceFeature``, like a source position, composition and energy spectrum of injected particles, redshift, directions of the injection and so on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "source = Source()\n",
    "source.add(SourceParticleType(nucleusId(1, 1))) # protons only\n",
    "source.add(SourcePowerLawSpectrum(1*EeV, 1000*EeV, -1))\n",
    "source.add(SourceIsotropicEmission())\n",
    "source.add(SourcePosition(Vector3d(5, 5, 5)*Mpc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Source object is supplied to the ModuleList as an argument of the run method together with the total number of candidates which would be injected. The syntax is the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sim.run(source, 10000, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The third argument means that the ModuleList should propagate children candidates created in physical interactions - recursively. Processing of children candidates starts after finishing propagation of the original particle.\n",
    "\n",
    "CRPropa optionally utilises [OpenMP](http:https://openmp.org), a shared-memory parallel processing framework, which enables it to process multiple candidates in parallel if the underlying system and hardware unit supports multithreading. Usually, OpenMP speeds up the simulation almost linearly (up to 8 cores), but in cost of an indeterministic ordering of results in output files. Hence if the ordering is important, like with trajectory output, or in a case of debugging the code to simplify the execution, the parallel processing can be disabled by specifying the environment variable ``OMP_NUM_THREADS=1`` before running the script."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Observer and output methods\n",
    "\n",
    "With the source object CRPropa offers systematical injection of candidates, but after reaching a stopping condition they are just discared. Thus, there should also exist antoher side of the pipline where candidates are systematically detected. An object which fulfill the role of detector is called **observer** and it is embodied in the Observer class. The observer is a module serving as an artificial geometrical boundary where different actions can be triggered when candidate cross the boundary, the most common action is to save the state of the detected candidate into a file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}