Simulation_new_william.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install scipy\n",
    "#!pip install matplotlib\n",
    "\n",
    "import ROOT\n",
    "import numpy as np\n",
    "#from root_numpy import random_sample\n",
    "from scipy.ndimage import uniform_filter\n",
    "import time\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import multiprocessing as mp\n",
    "\n",
    "# This import registers the 3D projection, but is otherwise unused.\n",
    "from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 unused import\n",
    "from matplotlib import cm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Simulation Model</h1>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "[4.87993747 5.99305597 4.60918068 5.38785233 3.23975128 6.15557717\n",
      " 3.68560935 2.05695411 1.37657942 1.31518549 2.50045579 2.0289163\n",
      " 2.00388827 1.8140459  1.28681887 2.098115   2.78698301 1.13799802\n",
      " 1.22317833 1.48158871 1.49344488]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5aUlEQVR4nO3dd3hc5ZX48e8ZjXofaVRsy5YtF8k2bsim2/TQyaZCQghppJAASTbJhuSXzW7KbgJLAqQtkAoBspBmwHRMDQHce5FsWbIlWZJVrF5m3t8fc8cIWWUkzZ2m83meeTyauXPv8cUcvfPec88rxhiUUkrFHke4A1BKKWUPTfBKKRWjNMErpVSM0gSvlFIxShO8UkrFKE3wSikVozTBq6ghIjNFpENE4ib4+dtE5P5gx6VUpNIEr4JKRKpEpFtE2kWkVUT+ISKfE5FJ/1szxlQbY9KMMZ4A4jhXRA4P+fwPjTGfnmwcoSYiM0TkzyLSJCJtIrJDRG6w3isWESMi64Z85kER+a71/FwR8Vq/HNtFZK+IfCL0fxMVaprglR2uNMakA7OA/wa+Afw6vCFFtQeAGnznMwf4GHB0yDaniciZo+yj1hiTBmQAXwbuE5EFdgSrIocmeGUbY0ybMWYt8GHg4yKyGEBEEkXkDhGpFpGjIvIrEUm23tstIlf49yEiThFpFJEVg0arTuu9T1jbt4vIARH5rPV6KvAUMM0atXaIyDQR+a6IPDho31eJyE7rm8ZLIlI26L0qEflXEdlmjZr/JCJJ1nu5IvKE9blmEXl1uG8oIvJLEbljyGt/F5GvWM+/ISJHBo2qLxjhVK4EfmeM6TTGDBhjNhtjnhqyzY+BHwTw38QYY9YBzcCSsbZX0U0TvLKdMeYt4DBwjvXSfwPzgWXAXGA68B3rvYeBawd9/D1AkzFm0zC7bgCuwDcq/QTwExFZYYzpBC7FGrVaj9rBHxSR+daxbgXcwDrgcRFJGLTZh4BLgNn4kuEN1utftf4+biAfuA0YrufHw8CHRUSsY2YDFwOPWKPnLwIrrW877wGqhtkHwD+Bn4vINSIyc4RtfgHMF5ELR3jf//d2iMhVQC5QMdq2KvppglehUgu4rGR3I/BlY0yzMaYd+CFwjbXdQ8BVIpJi/fwRfInyJMaYJ40xldao9GXgWd75JTKWDwNPGmOeM8b0A3cAycDgaY67jTG1xphm4HF8v5AA+oFCYJYxpt8Y86oZvqnTq/gSvz+mDwBvWL9sPEAisFBE4o0xVcaYyhFi/aC1r/8HHBSRLSKycsg23fhG8N8fYR/TRKTV2u6vwFeMMZtH2FbFCE3wKlSm45sWcAMpwEZriqMVeNp6HWNMBbAbuNJK8lfhS/onEZFLReSf1jRJK3AZvpFpIKYBh/w/GGO8+Oa5pw/apn7Q8y4gzXp+O77R77PW1NC/DXcAK+k/wjvfSD4C/HHQ3/NW4LtAg4g8IiLTRthPizHm34wxi/B9Y9gC/M3/zWCQ+4F8EblymN3UGmOy8H3buRs4f7hjqdiiCV7ZzhptTgdeA5rwjSIXGWOyrEemdQHQzz9NczWwy0qGQ/eZCPwZ38g730pe6wB/0hurTWotvouW/v0JUAQcGevvY4xpN8Z81RgzB98voK+MMn/+MPABEZkFnGbF7N/PQ8aYs604DPCjAI7dhO/vPA1wDXmvD/gP4Hu8cx6Gfr4X30XvU0TkvWMdT0U3TfDKNiKSYV0wfQR40Biz3Rop34dvvjzP2m66iLxn0EcfwTdX/XlGGL0DCfimOBqBARG51PqM31EgR0QyR/j8/wGXi8gFIhKPb169F/hHAH+vK0RkrvVLoQ3fdIt3uG2taZAmfKPrZ4wxrdY+FojI+dYvqh58v/SG3YeI/EhEFlsXnNPxnZcKY8yxYTZ/AEjCd+1gWNYvgv/hneseKkZpgld2eFxE2vFNeXwLuBPfRVC/b+Cb4viniBwHngdOlOwZY+qAN/DNh/9puANYc/c340vULfimP9YOen8PvtHzAWsqaNqQz+8FrgPuwZeAr8RX3tkXwN9vnhVzhxXnL4wx60fZ/iHgQt79yyoR38XmJnxTQXnAN0f4fAq+efNW4AC+Ef9Vw21o3SPwHYaM7ofxG2DmCNM5KkaILvihlFKxSUfwSikVozTBK6VUjNIEr5RSMUoTvFJKxShnuAMYLDc31xQXF4c7DKWUihobN25sMsa4h3svohJ8cXExGzZsCHcYSikVNUTk0Ejv6RSNUkrFKE3wSikVozTBK6VUjNIEr5RSMUoTvFJKxShN8EopFaM0wSulVIzSBD8FvV7RREVDe7jDUErZTBP8FHTLI1u4/Zm94Q5DKWUzTfBTTGfvAE0dvVQ2doY7FKWUzTTBTzGHW7oBOHSsk37PsCvEKaVihK0JXkSyROQxEdkjIrtF5Aw7j6fGVt3cBUC/x1BjPVdKxSa7R/B3AU8bY0qBpcBum4+nxlA9KKlXNHSEMRKllN1sS/DWavargV+DbyV3/4rykWTA48XjnTrr0tY0d5Ho9P1n13l4pWKbnSP42UAj8FsR2Swi94tI6tCNRORGEdkgIhsaGxttDGd4X/jjJj77wNRpUVzd3EWJO4289EQqG3UEr1QsszPBO4EVwC+NMcuBTuDfhm5kjLnXGFNujCl3u4ftWW+bhvYent99lJf2NtLROxDSY4dLdXMXM10plLjTNMErFePsTPCHgcPGmDetnx/Dl/AjxrptdXgNDHgN/6hoCnc4tvN6fRdWi1zJlOSlUtnQgTFTZ3pKqanGtgRvjKkHakRkgfXSBcAuu443EWu31jIvL43UhDhe2R/66aFQa+zopXfAe2IEf7xngKaOvnCHpZSyid1L9n0J+KOIJAAHgE/YfLyA1TR3sam6la9fsoBNh1p4eV8jxhhEJNyh2cZfQVPkSsFh/T0rGztwpyeGMyyllE1sTfDGmC1AuZ3HmKjHt9UCcOWSaaQnOnl+dwNVx7qYnXvSdeCYUX3Ml+BnulJIjI8DfAn+9Dk54QxLKWWTKXsn69ottayYmUWRK4XV830Xd1/ZF9vTNNXNXYjA9OxkCjOSSI6Po7JBSyWVilVTMsHvP9rOnvp2rlo6DYBZOakU56Twcown+JqWLgozkkh0xuFwCHPcqVpJo1QMm5IJ/vGttTgELl8y7cRrq+e7eaPyGL0DnjBGZi9fBU3KiZ+1VFKp2DblErwxhrVbazmzJPddFxfXzHfT3e9hQ1VLGKOzl78G3q/EncaR1m66+2L3l5pSU9mUS/Dbj7RRdazrxPSM3+lzcoiPk5idh+/p93D0eO+7E3xeKsbAwSadh1cqFk25BL92Sy3xccJ7FhW86/XURCfls1wxOw9/uMWqoMl59wge0GkapWLUlErwXq/hiW11rJmfR2ZK/Envr1ngZk99O0eP94QhOnsNroH3m52biogmeKVi1ZRK8G9VNVN/vIerlk0b9v3V83zlkrE4ih9cA++XFB/HjOxk7SqpVIyaUgn+8a21JMfHcWFZ3rDvlxWm405PjMl5+OrmbpLj48hJTXjX6yXuNCq1L7xSMWnKJPh+j5d12+u4aGE+KQnD38ArIqye5+a1iqaY6xHvr6AZ2oqhxJ3GgaYOvDH291VKTaEE/1pFEy1d/SdVzwy1ZoGb1q5+th1uDU1gITK0Bt5vbl4aPf1eatu6wxCVUspOUybBP76llowkJ+fMzx11u3Pm5iICr+yLnfbBxpiTauD93qmk0Xl4pWLNlEjwPf0entlZz6WLC0l0xo26bXZqAktmZPHyvoYQRWe/Y519dPd7mOlKPum9ErevuZrOwysVe6ZEgn9xTwOdfZ4Rq2eGWjMvly01rbR19dscWWj4SyQH18D7uVITyEqJ11JJpWLQlEjwa7fU4k5PDLgt7poFbrzGN28fC2qaTy6R9BMR7UmjVIyK+QTf3tPPi3sbuPyUQuIcgS3msXRGFulJzpgpl/TXwM/IPjnBg2+aRufglYo9MZ/gn915lL4Bb8DTMwDOOAfnzMs9scpTtKtu7iI/I5Gk+OGvP5S402hs76WtOzampJRSPjGf4NdurWVGdjLLi7LG9bnV89zUH+9hv00XH40xVITowuZIFTR+2pNGqdgU0wn+WEcvr1U0ceXSaeNea9W/ytPLe+2Zpvm/DTVceOfL7Ko9bsv+BxupBt6vJM9K8FpJo1RMiekEv25HPR6vGfPmpuFMy0pmXl4ar+wPfoLvG/By9wsVALbfUNU74KHueA9FI8y/AxRlJxMfJzoPr1SMiekE//iWWublpVFakD6hz6+e7+bNg81BXxDjL5sOc6TVd+fonvr2oO57qCMt3RgzfAWNnzPOQXGOLt+nVKyJ2QRf29rNW1XNXDWB6Rm/NfPd9A14+efBY0GLq9/j5WfrK1g6I5NlRVnsrrN3ima0GvjBtFRSqdgTswn+yW11AFw5gekZv1WzXSQ6HUEtl/zrpiMcbunmlgvnUVaYwZ76dlsrdUargR+sJC+V6mNd9Hu8tsWilAqtmE3wa7fWsnRGJsW5qRPeR1J8HKfPyQlaf3j/6H3JjEzOW5BHWWE6bd391Nu4wEh1cxeJTgfutMRRtytxpzHgNRyyauaVUtEvJhP8waZOth9pm9To3W/1fDcHGjtPjIQn42+bj1Dd3MXN589DRCgtyABgT5198/A1zd0UuVJwjHGTl5ZKKhV7YjLBr91SiwhcsWTyCX6NVS452WqaAWv0vnh6BhdYC44ssC7+7q63bx5+rBp4vzn+pmOa4JWKGTGX4I0xrN16hFXFLgoykya9vxJ3KtOzkic9D//3LbUcOvbO6B0gMzme6VnJto3gjTHUBJjg05Piyc9IpLJBSyWVihUxl+B31R2nsrFzXK0JRiMirJ7v5vWKYxO+AOkfvS8szOCihfnveq+0IJ09No3gW7v6ae8dGPUmp8G0kkap2BJzCf7xrXU4HcJliwuDts8183Pp6B1gc3XrxGLaVsvBpk5uvmDeSSWbpYXpVDZ20tMf3Fp7GFQiOc4EHwv9d5RSMZbgjTE8vrWWc+blkj1kcenJOHNuLnEOmdAiIB6v4Z4XKygtSOfiIaN3gLLCDDxee/rSjD/Bp9LeM0BjR2/QY1FKhV5MJfhN1S0cae0O2vSMX0ZSPCtmZk1oGb8nttVyoLGTWy6YN2wly4lKGhvuaPUn+BnZJ6/kNJx3etLoPLxSscDWBC8iVSKyXUS2iMgGO48FvuqZRKeDixYWBH3fq+e52X6kjaZxjG49XsPdL+xnQX4671k0fEzFOSkkOh3sseGO1prmLnLTEkhNdAa0vZZKKhVbQjGCP88Ys8wYU27nQQY8Xp7cXscFZXmkBZjQxmPNAl+55Gv7Ax/FP7m9jspG39z7SHXozjgH8/PTbRvBB3qBFaAgI4mUhDhN8ErFiJiZonnjwDGaOvom1DkyEIunZeJKTQi4XNLrNdzzwn7m5aVx6eLRv1HYVUkTaA28n8MhzNHVnZSKGXYneAM8KyIbReTG4TYQkRtFZIOIbGhsnHit+dottaQnOjl3Qd6E9zEah0M4Z14ur+xvxOsdu8pk3Y469jd08KVRRu9+pYUZNHX00dgevIub/R4vta3d40rwAHPdadoXXqkYYXeCP9sYswK4FLhJRFYP3cAYc68xptwYU+52uyd0kN4BD0/vrOfiRQUjLksXDKvnuWnq6GPXGPPlXmvuvcSdyuWnjF2uWWbd0RrMUXxdaw9ew7imaMA3D3+ktTvoLZKVUqFna4I3xhyx/mwA/gqssuM4L+9tpL1nIOjVM0OdMz/Xd7wxpmme3lnPvqMd3HzBvIAW+va3LAjmHa3jLZH081fSHGjSUbxS0c62BC8iqSKS7n8OXAzssONYa7fW4kpN4KySHDt2f0JeehILCzNGnYf3j97nuFMD7oWTk5ZIXnpiUHvSTDjBn6ik0Xl4paKdnSP4fOA1EdkKvAU8aYx5OtgH6eob4PndR7nslAKccfZfM16zwM3GQy209/QP+/6zu+rZU9/Ol86fG9Do3a+0MCPoI/iEOAf5GePrxzMrJwWH6PqsSsUC2zKiMeaAMWap9VhkjPmBHcdJjo/jsc+dyafPnmPH7k+yep6bAa/hjcqTV3nyeg13vVDB7NxUrhxnJ8uygnQqGjqCtuBGTXMXM7KTx/VLBnw98ItcKVoqqVQMiPoySRFh8fTJLewxHqfOyiY1IW7Y9sHP7T7K7rrjfPG8ueP+NlFamE6fx8vBpuBMjYy3Bn4wX08anaJRKtpFfYIPtQSngzNKcnl5X+O7mnIZ45t7L85J4eoJXOz1tywI1hqt462BH6zEncqBxg48AZSDKqUilyb4CVizwE1NczdVg5a3e353Aztrj3PTBEbv4Bs1Ox0SlDta27r6aevup8gVWA+a4WLpHfDV0Sulopcm+AlYM89Xr//yXl93SWMMd72wj5muFP5l+fQJ7TPB6WBuXlpQetLUtEysgsbPXypZofPwSkU1TfATMDMnheKcFF6x+tK8uKeBHUcmNvc+mK9lweRH8P4SycnMwYNW0igV7TTBT9Ca+W7eqDxGT7+Hu17YT5ErmX9ZMbHRu19pYQZ1bT20dvVNaj+TTfCu1ASyU+L1QqtSUU4T/AStnu+mu9/DHc/sZdvhNm46dy7xk6zDLz3RsmByo/ia5i6yU+LJSIqf8D50+T6lop8m+Ak6fU4OCXEO7n/tINOzknnfihmT3mdZobX4xyTn4SdTQeNX4k7jgCZ4paKaJvgJSk10Ul6cDcBN580lwTn5U5mXnkh2SnxQRvATnZ7xK8lLpamjb9LTRUqp8NEEPwnXrJrJmSU5fODUyY/ewXfTVmlBBrsnkeA9XsPhlvG3CR5Ke9IoFf00wU/CVUun8dBnTg/K6N2vrDCDffXtE77JqK6tmwGvCWKC12kapaKVJvgIU1qYTne/50QlzHhNtIvkUDOyk0mIc2iCVyqKaYKPMGUFk7vQWjPJEkk/Z5yD4twUKht0ikapaKUJPsLMy0/DIUx4Hr66uYs4h1CYOb42wcPRShqlopsm+AiTFB/H7NzUCY/gq5u7mZ6VHJTe+CXuNA41d9E3EJwWxkqp0NIEH4FKCzMmXCoZjBp4v5K8VDxeQ3WzTtMoFY00wUegsoJ0qpu76OgdGPdng1ED7+evpKnQeXilopIm+Ajk7w2/d5yj+I7eAZo7+4I2gp+jpZJKRTVN8BGotNDXk2a8i3/UBKlE0i8t0UlhZpImeKWilCb4CDQ9K5n0RCd76seX4INVAz+YLt+nVPTSBB+BRITSwnT21I1viibYI3iwlu9r6HjX8oRKqeigCT5ClRb4KmnGk1irm7vISHKSmTLxNsFDleSl0d47QGN7b9D2qZQKDU3wEaq0MJ2O3gEOtwS+Lmp1cxczc4I3eodBlTQ6D69U1AkowYvIWSKSaj2/TkTuFJFZ9oY2tfkracZTDx/MGng/7SqpVPQKdAT/S6BLRJYCXwUqgT/YFpVigX91pwArabxew+Hm7qDVwPvlZySSmhCn67MqFYUCTfADxjcZfDXwM2PMz4F0+8JSaYlOZrpSAh7BH23voc/jpSg7uAleRCjJ0+X7lIpGgSb4dhH5JnAd8KSIOIDgXclTwyotSGd3gKWS1ceCX0HjV+JO0xG8UlEo0AT/YaAX+JQxph6YAdxuW1QK8PWkqWrqpLvPM+a2dtTA+5W4U6lt66FzAq0TlFLhE1CCN8bUG2PuNMa8av1cbYzROXiblRWk4zWwv2HsaZqa5i4cAtOykoMeh/9C68EmvdCqVDQJtIrmfSKyX0TaROS4iLSLyMT62aqAlRb6F/8IIMG3dFOYmRzU5QP9SvK0J41S0SjQbPBj4CpjTKYxJsMYk26MyQjkgyISJyKbReSJiYc5Nc10pZAcHxfQPLwdJZJ+s3JScAg6D69UlAk0wR81xuye4DFuASb62SktziHMLwisZYGdCT7RGcdMV4rWwisVZQJN8BtE5E8icq01XfM+EXnfWB8SkRnA5cD9k4pyCltYmM6e+uOjtizo7vPQ2N4b9LtYB/M1HdMRvFLRJNAEnwF0ARcDV1qPKwL43E+BrwMjrvkmIjeKyAYR2dDY2BhgOFNHaUEGLV39NIzSC6amJTgLbY+mJC+NA02deLzadEypaOEMZCNjzCfGu2MRuQJoMMZsFJFzR9n3vcC9AOXl5Zo9higteKc3fH7G8Atp21kD71fiTqVvwMuRlm5bvykopYIn0CqaGSLyVxFpsB5/tqZfRnMWcJWIVAGPAOeLyIOTjHfKCaQnjZ018H4lurqTUlEn0Cma3wJrgWnW43HrtREZY75pjJlhjCkGrgFeNMZcN4lYp6TMlHimZSaN2pOmurmL1IQ4soPYJngoTfBKRZ9AE7zbGPNbY8yA9fgd4LYxLjVIaWHGqCN4/0LbImJbDNmpCbhSEzTBKxVFAk3wx6w2wXHW4zrgWKAHMca8ZIwJ5KKsGkZpQToVDR30DQx/rdrOEsnBStypVDZoqaRS0SLQBP9J4ENAPVAHfAAY94VXNTGlhRkMeM2wo2djDDUtoUrwWiqpVDQJtBfNIWPMVcYYtzEmzxjzXmNMtd3BKZ8yf2/4Ye5obezopaffG5LKlhJ3Gsc6+2jp7LP9WEqpyRu1TFJEvm6M+bGI3AOcVMJojLnZtsjUCbNzU0mIc/juaF3+7vf8C23bWQPvV5KXCsCBpg5OTXXZfjyl1OSMVQfvbzGwwe5A1MiccQ7m5aexe5gLraEokfSb6/Z9k6hs6OTUWZrglYp0oyZ4Y8zj1tMuY8yjg98TkQ/aFpU6SWlBBq/uP/lO3+pj3YjAdBvaBA81PdvXrVLn4ZWKDoFeZP1mgK8pm5QVptPQ3suxjne3LKhu7qIgI4mk+DjbY4hzCHNyUzXBKxUlxpqDvxS4DJguIncPeisD0OV9Qsh/R+ve+nbOnJt44nV/DXyolOSlsbWmNWTHU0pN3Fgj+Fp88+89wMZBj7XAe+wNTQ1WWmj1pBkyDx+qGni/5UVZHG7p5ujxnpAdUyk1MWPNwW8FtorIQ8aY/hDFpIaRm5ZIblriu1oW9PR7qD/eE9IEv7LYd3H17apmrlgyLWTHVUqNX6Bz8MUi8piI7BKRA/6HrZGpk5QVpr9rdafDLd0AFLnsv8Dqt3BaBsnxcbx9sDlkx1RKTcx4mo39Et+8+3nAHwDtDBlipQXp7DvawYDH17KgJoQlkn7xcQ5WzMri7aqWkB1TKTUxgSb4ZGPMC4BYd7V+F99KTSqESgsy6BvwUnXM1w+mOoQ3OQ22stjF7vrjHO/RWTulIlmgCb5XRBzAfhH5ooj8C5BmY1xqGCcutFprtFY3d5EU78Cdljjax4JuVbELY2DjIR3FKxXJAk3wtwApwM3AqcDHgI/bFZQa3ty8NJwOOdGTpsaqoLGzTfBwls3MwukQnYdXKsIFumTf29bTDrSLZNgkOuMocaf5etIQ+hJJv5QEJ4umZ7JB5+GVimhj3ej0U2PMrSLyOMM3G7vKtsjUsEoL09lQ1eJrE9zcxRklOWGJY1VxNr9/4xC9Ax4SnfbfRauUGr+xRvAPWH/eYXcgKjClBRn8fUstVce66OzzhGUED1Be7OK+Vw+y7XDbidp4pVRkGetGp43W0xzgSWNM72jbK/v5L7Q+t6seCG2J5GCDb3jSBK9UZAr0IuuVwD4ReUBErhCRgObuVfCVWT1pnt15FAhfgnelJjA3L00vtCoVwQJd0ekTwFzgUeBaoFJE7rczMDW8/IxEslLi2Vjtu8A5Izs8CR58o/gNh1rweE+6PKOUigCBjuCxetE8BTyCr+HYe22KSY1CRCgtSMcYcKcnkpwQvgucK4uzae8ZYO8wC5EopcIvoAQvIpeKyO+A/cD7gfuBAhvjUqPwtw4O1/SMn3/ufcMhnaZRKhIFOoK/HvgbsMAYc4MxZp0xRvvBh0mZdaE13Al+RnYyhZlJvKXz8EpFpEDn4K8FNgPnAIhIsoik2xmYGpl/BB/qHjRDiQjlxS7ermrGGJ2HVyrSBDpF8xngMeB/rZdm4BvRqzAoLUzn3AVuzlvgDncorCrO5ujxXmqau8MdilJqiECnaG4CzgKOAxhj9gN5dgWlRpfojON3n1jF8pnZ4Q6FlbPfqYdXSkWWgLtJGmP6/D9YdfD6nVwxPy+djCSnJnilIlCgCf5lEbkNSBaRi/DVwz9uX1gqWjgcvnn4tzTBKxVxAk3w/wY0AtuBzwLrgG/bFZSKLiuLXRxo7ORYh3ayUCqSBNou2CsifwP+ZoxptDckFW1WzfZdC3i7qoVLFuvtEUpFilFH8OLzXRFpAvYCe0WkUUS+M9aORSRJRN4Ska0islNE/iNYQavIsnh6JglOh87DKxVhxpqi+TK+6pmVxhiXMcYFnAacJSJfHuOzvcD5xpilwDLgEhE5fbIBq8iT6IxjWVGWJnilIsxYCf5jwLXGmIP+F4wxB4Dr8N3dOiLj02H9GG89tPImRq0qdrGz9jidvXqDs1KRYqwEH2+MaRr6ojUPHz/WzkUkTkS2AA3Ac8aYNycUpYp45cXZeLyGzdWt4Q5FKWUZK8H3TfA9AIwxHmPMMnx3vq4SkcVDtxGRG0Vkg4hsaGzU67fR6tRZ2TgELZdUKoKMleCXisjxYR7twCmBHsQY0wqsBy4Z5r17jTHlxphytzv8t96riUlPiqesMIMNmuCVihijJnhjTJwxJmOYR7oxZtQpGhFxi0iW9TwZuAjYE7TIVcRZWexic3Ur/R5vuENRSjGOBT8moBBYLyLbgLfxzcE/YePxVJitLHbR3e9hx5G2cIeilCLAG50mwhizDVhu1/5V5Fl54oan5ohohKbUVGfnCF5NMXnpSRTnpPB2VUu4Q1FKoQleBVl5sYsNVc14dSFupcJOE7wKqlXFLlq6+qls7Bh7Y6WG+EdFEzc/vBmPDhCCQhO8Cqp3FgDRaRo1Pr0DHv7tL9tZu7WWvfXt4Q4nJmiCV0FVnJNCblqi9qVR4/bAG4eobu4CYMMh/fcTDJrgVVCJCCuLs3nroP4PqgLX3NnHXS/sZ818N/kZifoNMEg0waugW1ns4khrN7WtuhC3CszdL+yns3eAb11eRnmxi436DTAoNMGroFulC3Grcahs7ODBfx7i2lUzmZ+fzspZ2dS29XBEBwiTpgleBV1pQTqpCXGa4FVA/mvdHpLi4/jyRfMBX6ktoH2NgkATvAo6Z5yDFbOyefugzqOq0f2jsonndx/lC+eVkJuWCPgGCGmJTh0gBIEmeGWLVcUu9h5tp7VrzK7Saoryeg0/eHI307OS+eRZs0+87oxzsHxmFhv0QuukaYJXtvB/zd54SP8nVcP7y+Yj7Kw9ztcvWUBSfNy73ltpDRDauvrDFF1s0ASvbLF8ZhbxcaILgKhhdfUNcPsze1halMVVS6ed9H55cTbGwKZqHSBMhiZ4ZYuk+DhOmZ7J2zFQD7+puoUtNa3hDiOm3PfKQY4e7+X/XV6GiJz0/rKiLJwO0Xn4SdIEr2yzcraL7Ufa6On3hDuUCfN6DZ9/cCPX3PuG9rkPkqPHe/jVy5VcdkrBiam8oVISnCyanqnz8JOkCV7ZZuUsF/0eE9Wj3801LRw93ovXC5/5wwYa23vDHVLUu+OZvXi8hm9cUjrqditnZbPlcCu9A9E7QAg3TfDKNuXF1gIgUTxNs257PQlxDv7wqVW0dPXxuQc3asKZhJ21bTy26TA3nFXMrJzUUbctL3bRN+DVb06ToAle2SYrJYEF+elRe6HVGMPTO+o5Z14up8/J4X8+uIyNh1r49l93YIy2sx0vY3xlkVnJ8dx03twxtz8xQNBpmgnTBK9sVV6czaZDLQxE4ULc2w63caS1m0sWFwBw+ZJCvnT+XB7deJjfvl4V3uCi0It7GvhH5TFuvXA+mcnxY26fm5bInNxUvaN1EjTBK1utmu2is8/Dnijs7/3UjnqcDuGihfknXvvyhfO5aGE+339yF6/ubwxjdNGl3+PlB+t2M8edykdOmxnw58qLs9lwqEVXCJsgTfDKViutKoloax/sm56p44ySHLJSEk687nAIP/nwMublpfPFhzZzsKkzjFFGj4ffquZAYye3XVpGfFzgaae82EWrrhA2YZrgla2mZSUzPSs56hZw2F3XTtWxLi47pfCk99ISndz/8XIc4qusae/Ruy1H09bdz0+e28cZc3K4oCxvXJ/1DxB0Hn5iNMEr2/kWAGmJqguTT++owyFw8aDpmcGKXCn8/KMrONjUyS2PbNE1REfx8/UVtHb3860RbmoajW+FsASdh58gTfDKditnu2jq6KXqWFe4QwnYuh31nDY7hxyrw+FwzizJ5btXLuTFPQ3c8ezeEEYXPaqPdfG716t4/4oZLJ6eOe7Piwjls1y8HWXfACOFJnhlu1X+r9lRMg+//2g7FQ0dXHpKwZjbXnf6LD5y2kx++VIlf99yJATRRZcfPb2HOIfwrxcvmPA+youzqWnupr6tJ4iRTQ2a4JXtStxpZKXER01fkad21APwnkVjJ3gR4btXLmLVbBdff2wb2w632hxd9Nh4qJknt9dx4+o5FGQmTXg//nn4aLuOEwk0wSvbORzW1+woSvDls7LJzwgsKSU4HfzyoyvITUvkxj9spOG4jjSNMXzvid3kpSfy2TVzJrWvhdMySI6P0740E6AJXoXEqtnZVB3roqE9spNfVVMnu+uOn7i5KVA5aYncd305bd393PjAxqhusBYMj2+rY0tNK197zwJSEpyT2le8tQBItAwQIokmeBUSJ75mR/gozD89c+kw5ZFjWTgtgzs/tJQtNa18awq3M+jp9/Cjp/awsDCD96+YEZR9lhe72F13XEtSx0kTvAqJRdMySYp3RPwNT0/vqGPpjEymZyVP6POXnlLIrRfO48+bDvPr1w4GObro8NvXqzjS2s23Ly/D4RhfWeRIVhZn4zWwubo1KPubKjTBq5BIcDpYXpQd0V+zD7d0sfVw24RG74PdfP48Ll1cwA/X7eblfVOrnUFX3wD3vlLJuQvcnDk3N2j7XT4zG4eg9fDjZFuCF5EiEVkvIrtEZKeI3GLXsVR0WDnb9zW7rq073KEM62n/9Mw459+HcjiEOz64lPn56XzxoU0cmEK32T/8Vg0tXf186fx5Qd1vWqKThdMy9I7WcZrc1Y/RDQBfNcZsEpF0YKOIPGeM2WXjMVUE++CpM/jVS5Xc/sxe7vzQsnCHc5KndtRTVpgxZp/yQKQmOrnv+nKu/vnrXHHPa+RnJJGR5CQjOZ6M5HgyrUdGkvVnsvOk19KTnDjH0bcl3HoHPNz3ygFOn+Pi1FnZQd9/+SwXf3q7hn6Pd1z9bKYy2xK8MaYOqLOet4vIbmA6oAl+iipypfDJs2fzq5cr+fgZxSwtygp3SCccPd7DxkMtfPWi+UHbZ5ErhQc+tYqH36qmrXuAtu5+jnf3c6Slm7buftq6+xkYo8VBbloCHz1tFp88azaZKWO32A2nv246Qv3xHm7/4BJb9r+y2MXv/lHFrtrjEfVvJ5LZOYI/QUSKgeXAm8O8dyNwI8DMmYG3EVXR6abzSnhsYw3fe2IXj37ujHH3JrHLMzv91TOTm54ZatG0TL7/3lOGfc8YQ3e/h+NW8vf/Amjr7ud4j+/PHUfauOuF/fz6tYPccGYxnzp7NtmpCcPuL5wGPF5++XIlS2ZkcnYQ594He2cBkGZN8AGyPcGLSBrwZ+BWY8zxoe8bY+4F7gUoLy+fmnVlU0h6UjxfvXgB3/zLdtZtr+fyJZO7oBks67bXMS8vjbl56SE7poiQkuAkJcE56p2eu+uO87MXK/j5SxX85vWDfOyMWXzmnDnkjtInJ9TW7ajn0LEufnXdqbb90s7PSGKmK4UNVS18+hxbDhFzbJ3IEpF4fMn9j8aYv9h5LBU9PlReRGlBOv/11O6IuCGoqaOXtw42T/riql3KCjP4+UdX8Mytq7mwLJ97XznA2T96ke8/sSsi7po1xvCL9RXMzUsbsftmsPgWAGmesvcYjJedVTQC/BrYbYy5067jqOgT5xC+c8VCDrd085vXw18r/uzOo3jNxG5uCqX5+encfe1ynvvyGi5bXMhvXj/IOT9ez3fX7gxrI64X9zSwp76dz68pCVrd+0hWFrto6uiLqs6k4WTnCP4s4GPA+SKyxXpcZuPxVBQ5c24uF5bl84v1lWFvX/DUjjqKc1IoLQjd9MxkzM1L484PL+PFr57L1cum8eA/D7H6x+v5f3/bwZHW0JagGmP4+foKpmclc9WyabYfb+WgeXg1NtsSvDHmNWOMGGOWGGOWWY91dh1PRZ/bLiulp9/Dnc/uC1sMrV19vFF5jEsWF0bMBd9AFeem8uMPLGX9v57L+0+dwSNvV3Pu7ev55l+2UdMcmhHumweb2VTdyufWzAlJ6WKJO43slHi94SlAWkyqwmaOO42Pn1nMnzbUsKv2pOvvIfHcrqMMeA2XBbl6JpSKXCn81/tO4aWvnce1q2by541HOPeOl/jao1tp7eqz9dg/X19BbloiHywvsvU4fiLCqbNcEd/TKFJogldhdfP588hKjud7T+wKy4Wzp3fUMz0rmVMmsNpQpJmelcx/Xr2YV75+HtefMYu/bTnCLY9swWvTcoLbDrfy6v4mPn3ObJLi42w5xnBWFmdzoKmTpo7ekB0zWmmCV2GVmRLPrRfO540Dx3hu19GQHru9p59X9zdx6eKCqJueGU1BZhL/fuUivnPlIl7e12hb07NfrK8kI8nJR08L7f0r5VHSmTQSaIJXYfeR02YyNy+NH67bTd+AN2THfXFPA30eb9BvbooU1502k/csyudHT+9ha01rUPe9/2g7T++s54Yzi0lPCu0dtounZ5DodOg8fAA0wauwi49z8K3Ly6g61sUf3qgK2XHXba8jPyOR5UXB75sSCUSEH79/KfkZSXzp4c1B7aX+y5crSY6P44azZgdtn4FKdMaxtCiLtw/pCH4smuBVRDhvQR6r57u564X9NHfae2EQoLN3gJf2NnLJogLba7fDKTMlnruuWcaR1u6gLUJS09zF37fUcu2qmbjC1DZhZXE2O4+00dU3EJbjRwtN8CpifPvyMrr6PPz0efvLJl/a20jvgJdLFkf2zU3BUF7s4ssXzmPt1loe3Xh40vu779UDOAQ+szr0o3e/8mIXA17DliBPPcUaTfAqYszPT+cjq2byxzer2X+03dZjPbWjjpzUBFbNdtl6nEjx+XPncsacHP797zupaJj4uW1o7+GRt2t4/4oZFGZObNWrYFgxMxsRvdA6Fk3wKqJ8+aL5pCTE8f0nd9t2jJ5+Dy/uaeDiRQXExfD0zGBxDuGn1ywjOSGOLz60ecI9gH7zWhUDHi+fXVMS5AjHJzM5ngX56XpH6xg0wauI4kpN4JYL5vHyvkZe2ttgyzFe2ddIV58nqm9umoj8jCT+54NL2VPfzg/Xjf8XaFtXPw/+8xCXL5nG7NzJL4oyWSuLXWw61MKAJ3SVV9FGE7yKONefUUxxTgrff3K3Lf/zPr2jnszkeE6fkxP0fUe680rz+PTZs/nDG4dOLFEYqD+8UUVH7wBfODe8o3e/8uJsOvs87Km3dzovmmmCVxEnwengtsvKqGjo4KG3qoO6774BL8/tPsrFC/On7LJvX7+klFOmZ/KNP28LuDlZV98Av3n9IBeU5lFWmGFzhIFZeeKGJ52mGcnU/BeuIt5FC/M5Y04OP3luH21dwavffr2yifaegZi9uSkQCU4H91y7nAGPl1se3hzQtyT/YtpfOG9uCCIMzLSsZKZnJWs9/ChCsmSfUuMlInz7ijKuuOc17nlxP9++YmFQ9vvU9jrSE52cZdOyctGiODeVH77vFG55ZAt3v7Cfr1y8YMRt7V5MezLKi7P554FjGGPC2m5i++E2dta20TvgpXfAQ2+/lz6P1/dzv8d63fden/95v7XtgJf0JCePfu7MoMelCV5FrEXTMvlweRG/f6OKj54+a9IX9vo9Xp7ddZQLyvJIdIauOVakunrZdF7d38Q96ys4vSSHM0uG/6X3t832LqY9GeXFLv6+pZbDLd0UuVJCfvzuPg8/fmYPv/tHFUPvIROBRKeDRGec78/4Qc+dDhKcDrJSEkh0OshJs+eGMU3wKqJ95eL5PL61lh+u281915dPal9vHmimtat/StzcFKj/uGoRm6pbuPWRLTx1yznkDFnn1eM1/PKlSk6Zbt9i2pMxeAGQUCf4t6ua+dqjW6k61sX11jq5KQlxJMb7krjTIWFvYqdz8Cqi5aUn8YXz5vLcrqP8o6JpUvt6akcdKQlxnLvAHaTool9qopN7rl1Oa3c///ro1pNaGazbXkfVsS5uOq8k7MlqOPPz0klPcvJ2CG946u7z8L0ndvGh/32DAa/hoc+cxn9evZgiVwo5aYmkJTqJj3NExPnSBK8i3qfOnu3rdf7ELhraeybUT8XjNTyz8yjnLcgLae/yaLBoWibfuqyM9Xvf3VrYvxxfiTuVixdG5kVph0M4dVY2Gw+FppJm46FmLrv7VX792kGuO20Wz9y6esSprUigUzQq4iXFx3HbZWXc9NAmVv3gBTKSnJTkpTHXncbcPN+jxJ1GkStlxDtTN1Q109TRO6WrZ0Zz/RmzeK2iiR89vYfTZudwyoxM1u/1Lab9Px9cGtEN2VYWu7j9mb20dvWRlWLPXHZPv4c7ntnLr18/yPSsZB769GmcGYFTVkNpgldR4fIlhRRknsG2w21UNnZQ0dDB+r2N72qeleB0MCc3lRJ3mu8XQF4aJW7fz0/tqCfR6eC8BXlh/FtELhHh9g8s4dK7XuVLD2/iiZvP4Wcvhm4x7ckotyp7Nh5q4YKy/KDvf+OhFr726FYONHXy0dNm8s3LykhLjI7UGR1RKgWcOsvFqbPe3RysraufisYOKhs6TiT+HbVtPLWjDv9KdSLgEOH80jxSo+R/zHDISkngrmuWc829b/DR+99ka00r37t6UcTfELa0KIv4OOHtquAm+J5+D3c+t4/7Xz1AYWYyf/z0aVFXXqv/2lVUy0yJ59RZ2SfVZ/f0e6g61klFQweVDZ0cau7k+jOKwxNkFFk128UtF8znJ8/vC+li2pORFB/HKdMzg3pH66Zq36i9srGTj5w2k9uiaNQ+WPRFrFQAkuLjKC3IoLQgMm6rjyZfPH8uTR29nFGSEzUXpFcWu/jt61X09HsmFXNPv4efPL+P+145QEFGEg98ahXnzIveqitN8Eqpd4lzCN977+JwhzEu5cUu/veVA2w/0naiR814eL2GDYdauO2v26lo6ODaVUXcdllZyNebDTZN8EqpqOefonu7qjmgBN/W1c/mmhY2V7eyuaaVLdUtHO8ZoDAzid9/chVr5kfvqH0wTfBKqajnSk1gbl7asCs8ebyG/Q3tbDrUyubqFjZVt1DZ2AmAQ3wriV2+ZBorZmZxyeKCqB+1D6YJXikVE1YWZ/PktjqaOnrZWtPK5upWNlW3sLWmlc4+3wpWrtQElhdl8b4VM1helMWSoqyovHgaqNj9mymlppTyWS4efquG8u8/D/iuJZQVpvO+FTNYMSuL5UXZzMpJiYgWAqGiCV4pFRMuWpTPR6pnUpSdwoqZWSyZkUVyQnRUAdlFE7xSKiZkJMXzw385JdxhRBTbblETkd+ISIOI7LDrGEoppUZm5z3IvwMusXH/SimlRmFbgjfGvALoarhKKRUmYe8iJCI3isgGEdnQ2NgY7nCUUipmhD3BG2PuNcaUG2PK3e7YuHtMKaUiQdgTvFJKKXtogldKqRhlZ5nkw8AbwAIROSwin7LrWEoppU4mE1nA2C4i0g7sDXccw8gFmsIdxDA0rvHRuMZH4xqfcMU1yxgz7AXMSLuTda8xpjzcQQwlIhs0rsBpXOOjcY2PxhU4nYNXSqkYpQleKaViVKQl+HvDHcAINK7x0bjGR+MaH40rQBF1kVUppVTwRNoIXimlVJBogldKqRgV8gQ/Vp948blbRCpEZJuIrIiQuM4VkTYR2WI9vhOCmIpEZL2I7BKRnSJyyzDbhPx8BRhXyM+XddwkEXlLRLZasf3HMNskisifrHP2pogUR0hcN4hI46Bz9mm747KOGycim0XkiWHeC/m5CjCusJwr69hVIrLdOu6GYd4PSw4bljEmpA9gNbAC2DHC+5cBTwECnA68GSFxnQs8EeJzVQissJ6nA/uAheE+XwHGFfLzZR1XgDTreTzwJnD6kG2+APzKen4N8KcIiesG4GdhOGdfAR4a7r9XOM5VgHGF5VxZx64Cckd5Pyw5bLhHyEfwZuw+8VcDfzA+/wSyRKQwAuIKOWNMnTFmk/W8HdgNTB+yWcjPV4BxhYV1HjqsH+Otx9BKgquB31vPHwMuEJtXYg4wrpATkRnA5cD9I2wS8nMVYFyRLCw5bDiROAc/HagZ9PNhIiR5AGdYX7GfEpFFoTyw9dV4Ob6R32BhPV+jxAVhOl/WV/stQAPwnDFmxHNmjBkA2oCcCIgL4P3W1/rHRKTI7piAnwJfB7wjvB+WcxVAXBD6c+VngGdFZKOI3DjM+xGTwyIxwUeqTfh6PiwF7gH+FqoDi0ga8GfgVmPM8VAddyxjxBW282WM8RhjlgEzgFUisjhUxx5NAHE9DhQbY5YAz/HOyNkWInIF0GCM2WjnccYrwLhCeq6GONsYswK4FLhJRFaH8NjjEokJ/ggw+LfxDOu1sDLGHPd/xTbGrAPiRSTX7uOKSDy+JPpHY8xfhtkkLOdrrLjCdb6GxNAKrOfktYFPnDMRcQKZwLFwx2WMOWaM6bV+vB841eZQzgKuEpEq4BHgfBF5cMg24ThXY8YVhnM1+NhHrD8bgL8Cq4ZsEjE5LBIT/FrgeutK9OlAmzGmLtxBiUiBf+5RRFbhO3e2/kO3jvdrYLcx5s4RNgv5+QokrnCcL+tYbhHJsp4nAxcBe4Zsthb4uPX8A8CLxro6Fs64hszTXoXv2oZtjDHfNMbMMMYU47uA+qIx5rohm4X8XAUSV6jP1aDjpopIuv85cDEwtPIuYnJYyLtJiq9P/LlArogcBv4d3wUnjDG/AtbhuwpdAXQBn4iQuD4AfF5EBoBu4Bq7/6HjG8l8DNhuzd0C3AbMHBRXOM5XIHGF43yBr8Ln9yISh++Xyv8ZY54Qkf8ENhhj1uL75fSAiFTgu7B+TYTEdbOIXAUMWHHdEIK4ThIB5yqQuMJ1rvKBv1pjFyfwkDHmaRH5HIQ3hw1HWxUopVSMisQpGqWUUkGgCV4ppWKUJnillIpRmuCVUipGaYJXSqkYpQleTXki8i3xdXjcZnUIPE1EXhrcKVBEykXkJev54E6Ze0TkjrAFr9QoNMGrKU1EzgCuwNcdcwlwIe/0EckTkUtH+OirVtuB5cAVInKW7cEqNU6a4NVUVwg0+W97N8Y0GWNqrfduB7412oeNMd3AFiKnIZ5SJ2iCV1Pds0CRiOwTkV+IyJpB770B9InIeSN9WESygXnAKzbHqdS4aYJXU5rVEO1U4EagEfiTiNwwaJPvA98e5qPniMhWfE2knjHG1Nsdq1LjpQleTXlWG9+XjDH/DnwReP+g914EkvGtzDPYq1Yr5EXAp0RkWajiVSpQmuDVlCYiC0Rk3qCXlgGHhmz2fXyLT5zEGHMQ+G/gG7YEqNQkaIJXU10avi6Pu0RkG7AQ+O7gDax+9o2j7ONXwGoJ4YLUSgVCu0kqpVSM0hG8UkrFKE3wSikVozTBK6VUjNIEr5RSMUoTvFJKxShN8EopFaM0wSulVIz6/xphEQMh0t01AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulation took 186.2 seconds\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-19-58174f5c5a92>:1057: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n",
      "  ax = fig.add_subplot(111)\n"
     ]
    }
   ],
   "source": [
    "# the constants of the constants. Pretty much choose once and leave it. Do not change.\n",
    "C_total_channels = int(3.6 * 1e6)\n",
    "C_channels_per_module = int(15*30)\n",
    "C_total_modules = int(C_total_channels/C_channels_per_module)\n",
    "\n",
    "# skip\n",
    "#C_max_radius = 1\n",
    "\n",
    "### Skipping radius dependency for now...\n",
    "\n",
    "\n",
    "#==============================================================================\n",
    "\n",
    "### Code for the simulation:\n",
    "\n",
    "class parameter_object:\n",
    "    # This might look cleaner, but it is less understood exactly what type of parameters and how many each function needs.\n",
    "    # I don't think I will implement the usage of this, (yet at least)\n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "class noise_generator:\n",
    "    def __init__(self, threshold, channel_count): # channel_count = C_total_channels\n",
    "        global global_counter\n",
    "        low = 0\n",
    "        high = 1000\n",
    "        noise = ROOT.TF1(\"noise_gauss\", \"TMath::Gaus(x, 0, 5)\", low, high)\n",
    "        self._mean = channel_count * noise.Integral(threshold, high)/noise.Integral(low, high)/2\n",
    "        #self._prng = ROOT.TRandom3(time.time_ns())\n",
    "        global_counter += 1\n",
    "        self._prng = ROOT.TRandom3(int(time.time()) + global_counter)\n",
    "\n",
    "    def generate(self):\n",
    "        \"\"\"The first return value is the exact number of noise activated channels in the small window and\n",
    "            the second return value is the estimated number of noise activated channels in the small window. (estimated from large window)\n",
    "            (both results are in small window, both results are noise above threshold)\n",
    "            noise_exact, noise_est = noise.generate()\"\"\"\n",
    "        # Large time window 25ns\n",
    "        # Small time window 3.125ns\n",
    "        prng = self._prng\n",
    "        total = prng.Poisson(self._mean) # number of noise above threshold in large window in this run.\n",
    "        narrow_exact = 0\n",
    "        if total <= 5:\n",
    "            for i in range(total):\n",
    "                if prng.Uniform(25) < 3.125:\n",
    "                    narrow_exact += 1\n",
    "        else:\n",
    "            # faster to use numpy\n",
    "            #narrow_exact = (np.random.uniform(low=0, high=25, size=total) < 3.125).sum()\n",
    "            narrow_exact = (np.random.rand(total) < 0.125).sum()  # 3.125/25 = 0.125\n",
    "        \n",
    "        noise_per_time = (total - narrow_exact)/(25 - 3.125)\n",
    "        estimate = 3.125 * noise_per_time\n",
    "        return narrow_exact, estimate\n",
    "\n",
    "    \n",
    "class signal_generator:\n",
    "    def __init__(self, threshold, real_expected_voltage, guess_expected_voltage): # expected_voltage = expected peak voltage  (mu)\n",
    "        global global_counter\n",
    "        low = -100\n",
    "        high = 1000\n",
    "        self._threshold = threshold\n",
    "        self._expected_voltage = real_expected_voltage\n",
    "        signal_func = ROOT.TF1(\"sth\", \"TMath::Landau(x, %f, %f)\" % (guess_expected_voltage, guess_expected_voltage/4), low, high)\n",
    "        self.fraction_above = signal_func.Integral(threshold, high)/signal_func.Integral(low, high)\n",
    "        #self._prng = ROOT.TRandom3(time.time_ns())\n",
    "        global_counter += 1\n",
    "        self._prng = ROOT.TRandom3(int(time.time()) + global_counter)\n",
    "    \n",
    "    def generate(self, expected_count):\n",
    "        \"\"\"The paramater is the expected number of particles that reach the detector.\n",
    "            The first return value is the number of particles that reached the detector and\n",
    "            the second return value is the number of activated channels.\n",
    "            signal_real_total, signal_above = signal.generate(30)\"\"\"\n",
    "        prng = self._prng\n",
    "        expected_voltage = self._expected_voltage\n",
    "        threshold = self._threshold\n",
    "        \n",
    "        signal_count = prng.Poisson(expected_count)\n",
    "        above_threshold_count = 0\n",
    "        sigma = expected_voltage/4\n",
    "        for i in range(signal_count):\n",
    "            if prng.Landau(expected_voltage, sigma) > threshold: # 1/4 is just the refernce 25/100\n",
    "                above_threshold_count += 1\n",
    "        return signal_count, above_threshold_count\n",
    "    \n",
    "\n",
    "\n",
    "def simulation_get_counts(bunch_crossing_count , ratio_of_channels_to_sample, signal_to_channel_count_ratio, threshold, real_signal_mu, guess_signal_mu):\n",
    "    \"\"\"Runs a simulation over specified number of BC, and returns two lists:\n",
    "        total_actual_signal_count and total_estimated_signal_count for each BC\"\"\"\n",
    "\n",
    "    channel_count = get_number_of_channel_readouts(ratio_of_channels_to_sample)\n",
    "    avg_particles_reaching_detector = channel_count * signal_to_channel_count_ratio\n",
    "    \n",
    "    noise_gen = noise_generator(threshold, channel_count)\n",
    "    signal_gen = signal_generator(threshold, real_signal_mu, guess_signal_mu)\n",
    "    \n",
    "    result_actual_signal_count = np.zeros(bunch_crossing_count)\n",
    "    result_estimated_signal_count = np.zeros(bunch_crossing_count)\n",
    "    \n",
    "    signal_fraction_above = signal_gen.fraction_above\n",
    "    \n",
    "    for BC in range(bunch_crossing_count):\n",
    "        # The first return value is the exact number of noise activated channels in the small window and\n",
    "        # the second return value is the estimated number of noise activated channels in the small window (estimated from large window).\n",
    "        # (both results are in small window, both results are noise above threshold)\n",
    "        noise_exact, noise_est = noise_gen.generate()\n",
    "        \n",
    "        # The paramater is the expected number of particles that reach the detector.\n",
    "        # The first return value is the number of particles that reached the detector and \n",
    "        # the second return value is the number of activated channels.\n",
    "        total_actual_signal, signal_above = signal_gen.generate(avg_particles_reaching_detector)\n",
    "        \n",
    "        # from the total, estimate how much total signal there was\n",
    "        total_activations = signal_above + noise_exact  # total activations in small window (i.e. above threshold)\n",
    "        total_signal_est = (total_activations - noise_est)/signal_fraction_above\n",
    "\n",
    "        # asign the results\n",
    "        result_actual_signal_count[BC] = total_actual_signal\n",
    "        result_estimated_signal_count[BC] = total_signal_est\n",
    "\n",
    "    return result_actual_signal_count, result_estimated_signal_count\n",
    "# we don't need to know the noise. Therefore not returing the noise values\n",
    "\n",
    "def simulation_parallel_get_counts(number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, threshold, real_signal_mu, guess_signal_mu):\n",
    "    \"\"\"Runs the simulations in parallel inorder to faster simulate. It then combines each result.\n",
    "        and returns two lists: total_actual_signal_count and total_estimated_signal_count for each BC\"\"\"\n",
    "    global global_counter\n",
    "   # parameters for the simulation\n",
    "    bins_of_BC = eq_div_strict(number_of_BC, nprocs)\n",
    "    param_list = [ (bins_of_BC[i],  # number of BC for each parallel simulation\n",
    "                    ratio_of_channels_to_sample,\n",
    "                    signal_to_channel_count_ratio,\n",
    "                    threshold,\n",
    "                    real_signal_mu, guess_signal_mu)\n",
    "                 for i in range(nprocs)]\n",
    "   # parallel simulation\n",
    "    results = np.asarray(pool.starmap(simulation_get_counts, param_list))\n",
    "    results = np.hstack(results).squeeze() # reshape\n",
    "    \n",
    "    result_actual_signal_count = results[0]\n",
    "    result_estimated_signal_count = results[1]\n",
    "\n",
    "   # list of count,  count for each bunch crossing event\n",
    "    return result_actual_signal_count, result_estimated_signal_count\n",
    "\n",
    "\n",
    "\n",
    "def get_listed_deviation_from_expected(result_actual_signal_count, result_estimated_signal_count):\n",
    "    \"\"\"Get the deviation from the simulated signals vs the expected signals as a list for each BC.\n",
    "        Returns the deviation for each BC as a list of positive numbers\"\"\"\n",
    "    deviation_list = np.abs( np.divide(result_estimated_signal_count, result_actual_signal_count) - 1 )\n",
    "    return deviation_list\n",
    "    \n",
    "def get_total_deviation_from_expected(result_actual_signal_count, result_estimated_signal_count):\n",
    "    \"\"\"Get the total deviation from the total simulated signals vs the expected signals.\n",
    "        Returns both the total deviation as a positive number\"\"\"\n",
    "    ra = result_actual_signal_count.sum()\n",
    "    if ra == 0:\n",
    "        return 0\n",
    "    total_deviation = abs( result_estimated_signal_count.sum() / ra - 1 )\n",
    "    return total_deviation\n",
    "\n",
    "def get_absolute_total_deviation_from_expected(result_actual_signal_count, result_estimated_signal_count):\n",
    "    \"\"\"Get the absolute deviation as a difference between the [total actual signal] and [total estimated signal]\"\"\"\n",
    "    total_absolute_deviation = abs(result_estimated_signal_count.sum() - result_actual_signal_count.sum())\n",
    "    return total_absolute_deviation\n",
    "\n",
    "def get_total_deviation_from_expected_with_simulated_results(number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, threshold, real_signal_mu, guess_signal_mu):\n",
    "    \"\"\"Small helper method to pipe the simulation results into the deviation calculation directly.\"\"\"\n",
    "    result_actual_signal_count, result_estimated_signal_count = simulation_parallel_get_counts(number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, threshold, real_signal_mu, guess_signal_mu)\n",
    "    #result_actual_signal_count, result_estimated_signal_count = simulation_get_counts(number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, threshold, signal_mu)\n",
    "    return get_total_deviation_from_expected(result_actual_signal_count, result_estimated_signal_count)\n",
    "    \n",
    "\n",
    "def get_the_best_threshold_and_deviation_in_terms_of_least_deviation(list_of_thresholds_to_try, number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, real_signal_mu, guess_signal_mu):\n",
    "    \"\"\"Simulates for each threshold provided in list, tests which threshold has the least deviation in computed vs expected.\n",
    "        and returns the best threshold it has found, also returns the deviation.\n",
    "        ASLO SEE THE FUNCTION \"get_the_best_threshold_and_deviation_in_terms_of_least_average_deviation\"\n",
    "        \"\"\"\n",
    "    # should not contain too large steps.\n",
    "    list_of_thresholds_to_try = list(list_of_thresholds_to_try) # enforce list\n",
    "    # calculate deviation for each simulation for each threshold\n",
    "    deviation_iter = ([\n",
    "        get_total_deviation_from_expected_with_simulated_results(\n",
    "            number_of_BC,\n",
    "            ratio_of_channels_to_sample,\n",
    "            signal_to_channel_count_ratio,\n",
    "            threshold,\n",
    "            real_signal_mu, guess_signal_mu)\n",
    "        for threshold in list_of_thresholds_to_try])\n",
    "    deviations = np.fromiter(deviation_iter, dtype=np.float, count=len(list_of_thresholds_to_try))\n",
    "    \n",
    "    # Blur the deviations,  assumes deviations follow a smooth curve.\n",
    "    # The actuall minimum should not be discarded because of a random deviation, therefore blur.\n",
    "    size = 5  # default is 3,  might want to change to 3\n",
    "    deviations_blur = uniform_filter(deviations, size=size)#[1:-1] # nah no need to slice, the boundary error doesn't bother us\n",
    "    \n",
    "    # get the minimum\n",
    "    index = np.argmin(deviations_blur)\n",
    "    threshold = list_of_thresholds_to_try[index]\n",
    "    deviation = deviations[index]\n",
    "    \n",
    "    return threshold, deviation\n",
    "\n",
    "def get_the_best_threshold_in_terms_of_least_deviation(list_of_thresholds_to_try, number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, real_signal_mu, guess_signal_mu):\n",
    "    \"\"\"Same as above, except it only returns the threshold\"\"\"\n",
    "    threshold, deviation = get_the_best_threshold_and_deviation_in_terms_of_least_deviation(list_of_thresholds_to_try, number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, real_signal_mu, guess_signal_mu)\n",
    "    return threshold\n",
    "\n",
    "\n",
    "def caluclate_SB(expected_signal_voltage, noise_sigma):\n",
    "    \"\"\"Formula for calculating S/B\"\"\"\n",
    "    return expected_signal_voltage / noise_sigma\n",
    "\n",
    "def calculate_expected_signal_voltage(signal_mu):\n",
    "    \"\"\"Some method to calculate the expected signal voltage for provided signal_mu and sigma\"\"\"\n",
    "    \n",
    "    # INTEGRATION METHOD:  (not used)\n",
    "        # Use the definition:  ∫(x*f(x))dx\n",
    "    \"\"\"\n",
    "        signal_sigma = signal_mu/4  # get_signal_sigma_for_mu(signal_mu)\n",
    "        low = 0  # There are no particles hitting the detector with negative energies.\n",
    "        high = 500\n",
    "        fSignalEnergy = ROOT.TF1(\"Signal\", \"TMath::Landau(x, %f, %f)*x\" % (signal_mu, signal_sigma), low, high)\n",
    "        scale = 1  # 0.001  # mili, no, sigma is also in mV so it doesn't matter\n",
    "        return fSignalEnergy.Integral(low, high) * scale\n",
    "    \"\"\"\n",
    "    # ASSUME MU IS THE EXPECTED VALUE:  (currently used)\n",
    "    return signal_mu\n",
    "\n",
    "def calculate_list_of_SB_from_list_of_signal_mu(noise_sigma, iterator_mu, steps):\n",
    "    \"\"\"Returns a list of calculated S/B for a list of signal mu.\"\"\"\n",
    "    return np.fromiter(\n",
    "                map(lambda signal_mu : caluclate_SB(\n",
    "                        calculate_expected_signal_voltage(signal_mu),\n",
    "                        noise_sigma) ,\n",
    "                    iterator_mu)\n",
    "             , dtype=np.float, count=steps)\n",
    "\n",
    "def compute_list_of_least_deviation_from_list_of_signal_mu(number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, linspace_data_thresholds, linspace_data_signal_mu, number_of_mu_steps):\n",
    "    \"\"\"It will try a lot of thresholds provided for each signal_mu provided,\n",
    "        and get the best result from the threasholds and put that least deviation result in a list for each signal mu.\n",
    "        ALSO: SEE THE FUNCTION NAMED  \"compute_list_of_least_average_deviation_from_list_of_signal_mu\"\n",
    "        \"\"\"\n",
    "    return np.fromiter(\n",
    "                 map(lambda signal_mu :\n",
    "                         get_the_best_threshold_and_deviation_in_terms_of_least_deviation(\n",
    "                             linspace(*linspace_data_thresholds),  # list_of_thresholds_to_try\n",
    "                             number_of_BC,\n",
    "                             ratio_of_channels_to_sample,\n",
    "                             signal_to_channel_count_ratio,\n",
    "                             signal_mu, signal_mu)[1] ,  # deviation is at index 1\n",
    "                         linspace(*linspace_data_signal_mu))\n",
    "                 , dtype=np.float, count=number_of_mu_steps)\n",
    "\n",
    "def compute_list_of_best_thresholds_from_list_of_signal_mu(number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, linspace_data_thresholds, linspace_data_signal_mu, number_of_mu_steps):\n",
    "    \"\"\"It will try a lot of thresholds provided for each signal_mu provided,\n",
    "        and get the best result threasholds from the least deviation,\n",
    "        and put that threasholds in a list for each signal mu.\"\"\"\n",
    "    return np.fromiter(\n",
    "                 map(lambda signal_mu :\n",
    "                         get_the_best_threshold_in_terms_of_least_deviation(\n",
    "                             linspace(*linspace_data_thresholds),  # list_of_thresholds_to_try\n",
    "                             number_of_BC,\n",
    "                             ratio_of_channels_to_sample,\n",
    "                             signal_to_channel_count_ratio,\n",
    "                             signal_mu, signal_mu) ,  \n",
    "                         linspace(*linspace_data_signal_mu))\n",
    "                 , dtype=np.float, count=number_of_mu_steps)\n",
    "\n",
    "\n",
    "def compute_average_deviation_from_simulations(simulation_samples, number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, real_signal_mu, guess_signal_mu, threshold):\n",
    "    \"\"\"This function will compute a more statistical accurate deviation from the parameters given.\n",
    "       It does this by simulating multiple times and takes the average deviation from those simulations.\n",
    "        simulation_samples = how many simulations to run for the statistics\n",
    "        number_of_BC = how many BC to run for each single simulation\n",
    "        ratio_of_channels_to_sample = the scale of the simulation,  1.0 = full scale.\n",
    "        average_collisions_mu = average number of collisions per bunch crossing.\n",
    "        signal_mu = average energy of the signal\n",
    "        threshold = threshold of the detector in for the simulations.\n",
    "    \"\"\"\n",
    "    \n",
    "    # makes diviation list\n",
    "    deviations = np.zeros(simulation_samples)\n",
    "    \n",
    "    # makes multiple simulations and stores the deviation for each simulation\n",
    "    for i in range(simulation_samples):\n",
    "        deviations[i] = get_total_deviation_from_expected_with_simulated_results(\n",
    "                number_of_BC,\n",
    "                ratio_of_channels_to_sample,\n",
    "                signal_to_channel_count_ratio,\n",
    "                threshold,\n",
    "                real_signal_mu, guess_signal_mu\n",
    "            )\n",
    "    \n",
    "    # takes the mean deviation found,  other intesting statistics can be made: e.g.  variance:  deviations.var()\n",
    "    the_mean_deviation = deviations.mean()\n",
    "    \n",
    "    # print statistics\n",
    "    #print_count_statistics(deviations)  # commented because not needed.\n",
    "    \n",
    "    return the_mean_deviation\n",
    "\n",
    "def compute_list_of_least_average_deviation_from_list_of_signal_mu(simulation_samples, number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, linspace_data_thresholds, linspace_data_signal_mu, number_of_mu_steps):\n",
    "    \"\"\"It will try a lot of thresholds provided for each signal_mu provided,\n",
    "        and get the best result from the threasholds and put that least deviation result in a list for each signal mu.\"\"\"\n",
    "    return np.fromiter(\n",
    "                 map(lambda signal_mu :\n",
    "                         get_the_best_threshold_and_deviation_in_terms_of_least_average_deviation(\n",
    "                             simulation_samples,\n",
    "                             linspace(*linspace_data_thresholds),  # list_of_thresholds_to_try\n",
    "                             number_of_BC,\n",
    "                             ratio_of_channels_to_sample,\n",
    "                             signal_to_channel_count_ratio,\n",
    "                             signal_mu)[1] ,  # deviation is at index 1\n",
    "                         linspace(*linspace_data_signal_mu))\n",
    "                     , dtype=np.float, count=number_of_mu_steps)\n",
    "\n",
    "def get_the_best_threshold_and_deviation_in_terms_of_least_average_deviation(simulation_samples, list_of_thresholds_to_try, number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, signal_mu):\n",
    "    \"\"\"Simulates for each threshold provided in list, tests which threshold has the least deviation in computed vs expected.\n",
    "        and returns the best threshold it has found, also returns the deviation.\"\"\"\n",
    "    # should not contain too large steps.\n",
    "    list_of_thresholds_to_try = list(list_of_thresholds_to_try) # enforce list\n",
    "    # calculate deviation for each simulation for each threshold\n",
    "    deviation_iter = ([\n",
    "        compute_average_deviation_from_simulations(\n",
    "            simulation_samples,\n",
    "            number_of_BC,\n",
    "            ratio_of_channels_to_sample,\n",
    "            signal_to_channel_count_ratio,\n",
    "            signal_mu, signal_mu,\n",
    "            threshold)\n",
    "        for threshold in list_of_thresholds_to_try])\n",
    "    deviations = np.fromiter(deviation_iter, dtype=np.float, count=len(list_of_thresholds_to_try))\n",
    "    \n",
    "    # Blur the deviations,  assumes deviations follow a smooth curve.\n",
    "    # The actuall minimum should not be discarded because of a random deviation, therefore blur.\n",
    "    size = 5  # default is 3,  might want to change to 3\n",
    "    deviations_blur = uniform_filter(deviations, size=size)#[1:-1] # nah no need to slice, the boundary error doesn't bother us\n",
    "    \n",
    "    # get the minimum\n",
    "    index = np.argmin(deviations_blur)\n",
    "    threshold = list_of_thresholds_to_try[index]\n",
    "    deviation = deviations[index]\n",
    "    \n",
    "    return threshold, deviation\n",
    "        \n",
    "def compute_best_statistical_threshold_with_least_deviation_from_collision_and_signal_mu(simulation_samples, number_of_BC, ratio_of_channels_to_sample, average_collisions_mu, real_signal_mu, guess_signal_mu, max_threshold, min_threshold, steps_threshold):\n",
    "    \"\"\"This function will compute the statistical best threshold from [simulation_samples] number of simulations.\n",
    "       That way you get higher accuracy for the threshold and deviation.\n",
    "        simulation_samples = how many simulations to run for the statistics\n",
    "        number_of_BC = how many BC to run for each single simulation\n",
    "        ratio_of_channels_to_sample = the scale of the simulation,  1.0 = full scale.\n",
    "        average_collisions_mu = average number of collisions per bunch crossing.\n",
    "        signal_mu = average energy of the signal\n",
    "        max_threshold = highest threshold to test.\n",
    "        min_threshold = lowest threshold to test.\n",
    "        steps_threshold = how many steps inbetween [min_threshold] and [max_threshold].\n",
    "        RETURNS: the best threshold found from the [simulation_samples], and it's average deviation.\n",
    "    \"\"\"\n",
    "    \n",
    "    # set average_collisions_mu\n",
    "    signal_to_channel_count_ratio = get_signal_to_channel_count_ratio(\n",
    "                                        collisions_per_BC = average_collisions_mu, # here we set it.\n",
    "                                        particles_per_collision = 500,\n",
    "                                        ratio_of_particles_reaching_HGTD=0.1) # some example value, idk if it matters\n",
    "    \n",
    "    # threshold to try, should not contain too large steps.\n",
    "    list_of_thresholds_to_try = list(linspace(min_threshold, max_threshold, steps_threshold))\n",
    "    \n",
    "    # makes diviation list\n",
    "    deviations = np.zeros(len(list_of_thresholds_to_try))\n",
    "    \n",
    "    # simulates multiple times and sums the deviations\n",
    "    for i in range(simulation_samples):\n",
    "        # calculate deviation for each threshold by simulating with the given parameters\n",
    "        deviation_iter = ([\n",
    "            get_total_deviation_from_expected_with_simulated_results(\n",
    "                number_of_BC,\n",
    "                ratio_of_channels_to_sample,\n",
    "                signal_to_channel_count_ratio,\n",
    "                threshold,\n",
    "                real_signal_mu, guess_signal_mu)\n",
    "            for threshold in list_of_thresholds_to_try])\n",
    "        # adds all deviations together.\n",
    "        deviations += np.fromiter(deviation_iter, dtype=np.float, count=len(list_of_thresholds_to_try))\n",
    "    \n",
    "    # now we have summed up deviations from multiple simulations [simulation_samples]\n",
    "    # and it is time to find the element with least deviaiton\n",
    "    \n",
    "    # Blur the deviations,  assumes deviations follow a smooth curve.\n",
    "    # The actuall minimum should not be discarded because of a random deviation, therefore blur.\n",
    "    size = 5  # default is 3,  might want to change to 3\n",
    "    deviations_blur = uniform_filter(deviations, size=size)#[1:-1] # nah no need to slice, the boundary error doesn't bother us\n",
    "    #deviations_blur = deviations  # turns off the blur. No we want blur.\n",
    "    \n",
    "    # get the minimum / least deviaiton\n",
    "    index = np.argmin(deviations_blur)\n",
    "    threshold = list_of_thresholds_to_try[index]\n",
    "    deviation = deviations[index] / simulation_samples  # we have to divide it to get the average, since we summed it before.\n",
    "    \n",
    "    # returns the threshold found with least deviation, and it's deviation\n",
    "    return threshold, deviation\n",
    "\n",
    "\n",
    "\n",
    "# -End of simulation code-\n",
    "\n",
    "#==============================================================================\n",
    "\n",
    "### Code for the plots:\n",
    "\n",
    "def plot_deviation_vs_SB():\n",
    "    \"\"\"Plots deviation vs S/B. See settings below to calibrate the simulation\"\"\"\n",
    "    ## Simulator Settings:\n",
    "    ratio_of_channels_to_sample = 1.0  # scale of simulation,   0.001 = small fast,   1.0 = full scale slow.\n",
    "    number_of_BC = 80 # * 120\n",
    "    \n",
    "    simulation_samples = 3 # 10\n",
    "    \n",
    "    #signal_to_channel_count_ratio = get_signal_to_channel_count_ratio()\n",
    "    signal_to_channel_count_ratio = get_signal_to_channel_count_ratio(collisions_per_BC = 0.01)\n",
    "    noise_sigma = 5  #mv\n",
    "    \n",
    "    min_signal_mu = 6  #mv  # cannot be below 0\n",
    "    max_signal_mu = 27 #mv\n",
    "    steps_signal_mu = (max_signal_mu - min_signal_mu) * 1   # resolution\n",
    "    #steps_signal_mu = 20\n",
    "    \n",
    "    min_threshold = 6  #mv  # cannot be below 0\n",
    "    max_threshold = 26 #mv\n",
    "    steps_threshold = (max_threshold - min_threshold) * 1   # resolution\n",
    "\n",
    "    ## Simulate, get data:\n",
    "    #signal_mu_list = list(linspace(min_signal_mu, max_signal_mu, steps_signal_mu)) # can be used to plot\n",
    "    SB = calculate_list_of_SB_from_list_of_signal_mu(\n",
    "            noise_sigma,\n",
    "            linspace(min_signal_mu, max_signal_mu, steps_signal_mu),\n",
    "            steps_signal_mu\n",
    "        )\n",
    "    \n",
    "    linspace_data_signal_mu = (min_signal_mu, max_signal_mu, steps_signal_mu)\n",
    "    linspace_data_thresholds = (min_threshold, max_threshold, steps_threshold)\n",
    "    \n",
    "    deviations = compute_list_of_least_average_deviation_from_list_of_signal_mu(\n",
    "            simulation_samples,\n",
    "            number_of_BC,\n",
    "            ratio_of_channels_to_sample,\n",
    "            signal_to_channel_count_ratio,\n",
    "            linspace_data_thresholds,\n",
    "            linspace_data_signal_mu,\n",
    "            steps_signal_mu\n",
    "        )\n",
    "    runs = 5\n",
    "    n = 1 #For seeing how quick the program is\n",
    "    for _ in range(runs-1):\n",
    "        deviations += compute_list_of_least_average_deviation_from_list_of_signal_mu(\n",
    "            simulation_samples,\n",
    "            number_of_BC,\n",
    "            ratio_of_channels_to_sample,\n",
    "            signal_to_channel_count_ratio,\n",
    "            linspace_data_thresholds,\n",
    "            linspace_data_signal_mu,\n",
    "            steps_signal_mu\n",
    "        )\n",
    "        n +=1\n",
    "        print(n)\n",
    "    deviations = np.divide(deviations, runs)\n",
    "    \n",
    "    #deviations_log = np.log(deviations)\n",
    "    \n",
    "    # print data\n",
    "    print(deviations)\n",
    "    ## Plot:\n",
    "    x_list = SB\n",
    "    y_list = deviations\n",
    "    title = \"Deviations vs SNR\"\n",
    "    prefix_filename = \"Deviations_vs_SB_with_newNoise_BC\" + str(number_of_BC)\n",
    "    xlabel = \"SNR\"\n",
    "    ylabel = \"Deviations\"\n",
    "    make_integer_axis = False\n",
    "    make_margin = False\n",
    "    plot_the_data(title, prefix_filename, xlabel, ylabel, make_integer_axis, make_margin, x_list, y_list)\n",
    "    \n",
    "\n",
    "def plot3D_threshold_and_SB_vs_deviation():\n",
    "    \"\"\"Creates a 3D plot of threshold and S/B vs deviaiton\"\"\"\n",
    "    ## Simulator Settings:\n",
    "    ratio_of_channels_to_sample = 0.01  # scale of simulation,   0.001 = small fast,   1.0 = full scale slow.\n",
    "    number_of_BC = mp.cpu_count() * 1\n",
    "    \n",
    "    signal_to_channel_count_ratio = get_signal_to_channel_count_ratio()\n",
    "    noise_sigma = 5  #mv\n",
    "    \n",
    "    min_signal_mu = 15  #mv  # cannot be below 0\n",
    "    max_signal_mu = 40 #mv\n",
    "    steps_signal_mu = (max_signal_mu - min_signal_mu) * 1   # resolution\n",
    "    \n",
    "    min_threshold = 10  #mv  # cannot be below 0\n",
    "    max_threshold = 80 #mv\n",
    "    steps_threshold = (max_threshold - min_threshold) * 1   # resolution\n",
    "\n",
    "    canvas = ROOT.TCanvas(\"c\",\"Graph2D example\",0,0,600,400)\n",
    "    plot_threshold_SB_deviation = ROOT.TGraph2D()\n",
    "    point_index = 0\n",
    "    ## Simulate, get data:\n",
    "    for loop_signal_mu in linspace(min_signal_mu, max_signal_mu, steps_signal_mu):\n",
    "        for loop_treshold in linspace(min_threshold, max_threshold, steps_threshold):\n",
    "            \n",
    "            deviation = get_total_deviation_from_expected_with_simulated_results(\n",
    "                            number_of_BC,\n",
    "                            ratio_of_channels_to_sample,\n",
    "                            signal_to_channel_count_ratio,\n",
    "                            loop_treshold,\n",
    "                            loop_signal_mu\n",
    "                        )\n",
    "            SB = caluclate_SB(calculate_expected_signal_voltage(loop_signal_mu), noise_sigma)\n",
    "            \n",
    "            plot_threshold_SB_deviation.SetPoint(point_index, SB, loop_treshold, deviation) # i, x, y, z\n",
    "            point_index += 1\n",
    "    \n",
    "    \n",
    "    plot_threshold_SB_deviation.SetTitle(\"Graph title; X axis title; Y axis title; Z axis title\");\n",
    "    plot_threshold_SB_deviation.Draw(\"surf1\")\n",
    "    canvas.Draw()\n",
    "    print(\"done\")\n",
    "\n",
    "    \n",
    "def plot_particle_intensity_mu_vs_best_threshold():    \n",
    "    ## Simulator Settings:\n",
    "    ratio_of_channels_to_sample = 1.0  # scale of simulation,   0.001 = small fast,   1.0 = full scale slow.\n",
    "    number_of_BC = mp.cpu_count() * 8\n",
    "    \n",
    "    signal_mu = 50 # mv  # might want to change this\n",
    "    \n",
    "    min_threshold = 2  #mv  # cannot be below 0\n",
    "    max_threshold = 80 #mv\n",
    "    steps_threshold = (max_threshold - min_threshold) * 1   # resolution\n",
    "    \n",
    "    # important setting for this simulation\n",
    "    min_intensity_mu = 0.008\n",
    "    max_intensity_mu = 0.015\n",
    "    steps_intensity_mu = (15-8)*3\n",
    "    \n",
    "    intensity_mu_list = list(linspace(min_intensity_mu, max_intensity_mu, steps_intensity_mu))\n",
    "    threshold_list = np.zeros(len(intensity_mu_list))\n",
    "    \n",
    "    ## Simulate, get data:\n",
    "    for index, collisions_per_BC in enumerate(intensity_mu_list):\n",
    "        signal_to_channel_count_ratio = get_signal_to_channel_count_ratio(collisions_per_BC=collisions_per_BC)\n",
    "        \n",
    "        threshold, deviation = get_the_best_threshold_and_deviation_in_terms_of_least_deviation(\n",
    "                        linspace(min_threshold, max_threshold, steps_threshold), #list_of_thresholds_to_try,\n",
    "                        number_of_BC,\n",
    "                        ratio_of_channels_to_sample,\n",
    "                        signal_to_channel_count_ratio,\n",
    "                        signal_mu\n",
    "                    )\n",
    "        threshold_list[index] = threshold\n",
    "    \n",
    "    # print result in list form\n",
    "    print(threshold_list)\n",
    "    \n",
    "    ## Plot:\n",
    "    x_list = intensity_mu_list\n",
    "    y_list = threshold_list\n",
    "    title = \"threshold vs intensity_mu\"\n",
    "    prefix_filename = \"threshold_vs_intensity_mu_with_BC\" + str(number_of_BC)\n",
    "    xlabel = \"intensity_mu\"\n",
    "    ylabel = \"threshold\"\n",
    "    make_integer_axis = True\n",
    "    make_margin = False\n",
    "    plot_the_data(title, prefix_filename, xlabel, ylabel, make_integer_axis, make_margin, x_list, y_list)\n",
    "    \n",
    "\n",
    "def plot_signal_mu_vs_best_threshold():\n",
    "    ## Simulator Settings:\n",
    "    ratio_of_channels_to_sample = 1.0  # scale of simulation,   0.001 = small fast,   1.0 = full scale slow.\n",
    "    number_of_BC = mp.cpu_count() * 4\n",
    "    \n",
    "    signal_to_channel_count_ratio = get_signal_to_channel_count_ratio()\n",
    "    \n",
    "    noise_sigma = 5  #mv\n",
    "    \n",
    "    min_threshold = 2  #mv  # cannot be below 0\n",
    "    max_threshold = 80 #mv\n",
    "    steps_threshold = (max_threshold - min_threshold) * 1   # resolution\n",
    "    \n",
    "    min_signal_mu = 12  #mv  # cannot be below 0\n",
    "    max_signal_mu = 40 #mv\n",
    "    steps_signal_mu = (max_signal_mu - min_signal_mu) * 2   # resolution\n",
    "    \n",
    "    linspace_data_thresholds = (min_threshold, max_threshold, steps_threshold)\n",
    "    linspace_data_signal_mu = (min_signal_mu, max_signal_mu, steps_signal_mu)\n",
    "    \n",
    "    ## Simulate, get data:\n",
    "    signal_mu_list = list(linspace(*linspace_data_signal_mu))\n",
    "    \n",
    "    threshold_list = compute_list_of_best_thresholds_from_list_of_signal_mu(\n",
    "        number_of_BC,\n",
    "        ratio_of_channels_to_sample,\n",
    "        signal_to_channel_count_ratio,\n",
    "        noise_sigma,\n",
    "        linspace_data_thresholds,\n",
    "        linspace_data_signal_mu,\n",
    "        steps_signal_mu\n",
    "    )\n",
    "    \n",
    "    ## Plot:\n",
    "    x_list = signal_mu_list\n",
    "    y_list = threshold_list\n",
    "    title = \"threshold vs signal_mu\"\n",
    "    prefix_filename = \"threshold_vs_signal_mu_with_BC\" + str(number_of_BC)\n",
    "    xlabel = \"signal_mu\"\n",
    "    ylabel = \"threshold\"\n",
    "    make_integer_axis = False\n",
    "    make_margin = False\n",
    "    plot_the_data(title, prefix_filename, xlabel, ylabel, make_integer_axis, make_margin, x_list, y_list)\n",
    "\n",
    "    \n",
    "def plot_deviation_vs_threshold(signal_mu, min_threshold, max_threshold, number_of_thresholds_to_try):\n",
    "    \"\"\"Plots deviation vs threshold for given mu.\"\"\"\n",
    "    ## Simulator Settings:\n",
    "    ratio_of_channels_to_sample = 0.001  # scale of simulation,   0.001 = small fast,   1.0 = full scale slow.\n",
    "    number_of_BC = mp.cpu_count() * 4\n",
    "    \n",
    "    signal_to_channel_count_ratio = get_signal_to_channel_count_ratio()\n",
    "\n",
    "    steps_threshold = number_of_thresholds_to_try  # resolution, number of steps\n",
    "    \n",
    "    linspace_data_thresholds = np.linspace(min_threshold, max_threshold, steps_threshold)\n",
    "    deviation_list = []    # Är append ett ineffektivt sätt att koda på?\n",
    "    \n",
    "    for threshold in linspace_data_thresholds:\n",
    "        deviation_list.append(get_total_deviation_from_expected_with_simulated_results(\n",
    "            number_of_BC, \n",
    "            ratio_of_channels_to_sample, \n",
    "            signal_to_channel_count_ratio, \n",
    "            threshold, \n",
    "            signal_mu\n",
    "        ))\n",
    "    \n",
    "#     print(linspace_data_thresholds, deviation_list)\n",
    "    \n",
    "    ## Plot:\n",
    "    x_list = linspace_data_thresholds\n",
    "    y_list = deviation_list\n",
    "    title = \"deviation vs threshold for mu=\"+str(signal_mu)\n",
    "    prefix_filename = \"deviation_vs_threshold_with_BC\" + str(number_of_BC)+\"_mu\"+str(signal_mu)\n",
    "    xlabel = \"threshold\"\n",
    "    ylabel = \"deviation\"\n",
    "    make_integer_axis = True\n",
    "    make_margin = False\n",
    "    plot_the_data(title, prefix_filename, xlabel, ylabel, make_integer_axis, make_margin, x_list, y_list)\n",
    "\n",
    "    \n",
    "def plot_deviation_vs_threshold_2():\n",
    "    # Rough estimate of calculation time:\n",
    "    # It takes roughly:  40 seconds\n",
    "        # with following settings:\n",
    "        #  ratio_of_channels_to_sample = 1.0,   simulation_samples = 10,   number_of_BC = nprocs * 4\n",
    "        #  steps_threshold = 30\n",
    "    \n",
    "    ratio_of_channels_to_sample = 1.0  # scale of simulation,   0.001 = small fast,   1.0 = full scale slow.\n",
    "    \n",
    "    # how many simulations to sample for the deviation. - larger = more accurate\n",
    "    simulation_samples = 10\n",
    "    \n",
    "    # number of BC for each single simulation.\n",
    "    number_of_BC = nprocs * 4   # (works best with a multiple of [nprocs], nprocs is about equal to 4 or 8)\n",
    "    \n",
    "    # intensity of collisions, average collisions per bunch crossing\n",
    "    average_collisions_mu = 200\n",
    "    # set average_collisions_mu\n",
    "    signal_to_channel_count_ratio = get_signal_to_channel_count_ratio(\n",
    "                                        collisions_per_BC = average_collisions_mu, # here we set it.\n",
    "                                        particles_per_collision = 500,\n",
    "                                        ratio_of_particles_reaching_HGTD=0.1) # some example value, idk if it matters\n",
    "    \n",
    "    # average energy of the signal\n",
    "    signal_mu = 100\n",
    "    \n",
    "    # thresholds to use\n",
    "    min_threshold = 10\n",
    "    max_threshold = 100\n",
    "    steps_threshold = 30   # resolution,  it takes longer time if you increase resolution. but the graph is better\n",
    "    \n",
    "    # make deviation list\n",
    "    list_deviation = []\n",
    "    list_threshold = list(linspace(min_threshold, max_threshold, steps_threshold))\n",
    "    \n",
    "    # here we loop over signal_mu\n",
    "    # could also edit the code to have a fixed value singnal_mu and instead loop average_collisions_mu\n",
    "    for threshold in list_threshold:\n",
    "        # calculate the threshold and deviaiton\n",
    "        deviation = compute_average_deviation_from_simulations(\n",
    "                simulation_samples,\n",
    "                number_of_BC,\n",
    "                ratio_of_channels_to_sample,\n",
    "                signal_to_channel_count_ratio,\n",
    "                signal_mu,\n",
    "                threshold\n",
    "        )\n",
    "        list_deviation.append(deviation)\n",
    "        \n",
    "    # print data\n",
    "    #print(list_deviation)\n",
    "    \n",
    "    ## Plot:\n",
    "    x_list = list_threshold\n",
    "    y_list = list_deviation\n",
    "    title = \"deviation vs threshold,  for collision_mu=\" + str(average_collisions_mu) + \", signal_mu=\" + str(signal_mu)\n",
    "    prefix_filename = \"Deviation vs threshold_nr2_with_BC\" + str(number_of_BC)\n",
    "    xlabel = \"threshold\"\n",
    "    ylabel = \"deviation\"\n",
    "    make_integer_axis = False\n",
    "    make_margin = False\n",
    "    plot_the_data(title, prefix_filename, xlabel, ylabel, make_integer_axis, make_margin, x_list, y_list)\n",
    "\n",
    "    \n",
    "def plot_best_threshold_deviation_vs_signal_mu():\n",
    "    # Rough estimate of calculation time:\n",
    "    # It takes roughly:  900 seconds = 15 minuites\n",
    "        # with following settings:\n",
    "        #  ratio_of_channels_to_sample = 1.0,   simulation_samples = 10,   number_of_BC = nprocs * 4\n",
    "        #  steps_threshold = 30,   steps_signal_mu = 25\n",
    "    \n",
    "    ratio_of_channels_to_sample = 1.0  # scale of simulation,   0.001 = small fast,   1.0 = full scale slow.\n",
    "    \n",
    "    # how many simulations to sample for the deviation. - larger = more accurate\n",
    "    simulation_samples = 10 \n",
    "    \n",
    "    # number of BC for each single simulation.\n",
    "    number_of_BC = nprocs * 4   # (works best with a multiple of [nprocs], nprocs is about equal to 4 or 8)\n",
    "    \n",
    "    # intensity of collisions, average collisions per bunch crossing\n",
    "    average_collisions_mu = 200\n",
    "    \n",
    "    # threshold searches in the simulation (what thresholds to test deviation)\n",
    "    min_threshold = 10\n",
    "    max_threshold = 40\n",
    "    steps_threshold = 30   # resolution,  it takes longer time if you increase resolution. but the graph is better\n",
    "    \n",
    "    # signal mu to loop,  (could also edit for fixed signal mu and loop average_collisions_mu instead)\n",
    "    min_signal_mu = 5\n",
    "    max_signal_mu = 100\n",
    "    steps_signal_mu = 25\n",
    "    list_signal_mu = list(linspace(min_signal_mu, max_signal_mu, steps_signal_mu))\n",
    "    \n",
    "    # create lists\n",
    "    list_threshold = []\n",
    "    \n",
    "    # here we loop over signal_mu\n",
    "    # could also edit the code to have a fixed value singnal_mu and instead loop average_collisions_mu\n",
    "    for signal_mu in linspace(min_signal_mu, max_signal_mu, steps_signal_mu):\n",
    "        # calculate the threshold and deviaiton\n",
    "        threshold, deviation = compute_best_statistical_threshold_with_least_deviation_from_collision_and_signal_mu(\n",
    "                simulation_samples,\n",
    "                number_of_BC,\n",
    "                ratio_of_channels_to_sample,\n",
    "                average_collisions_mu,\n",
    "                signal_mu,\n",
    "                max_threshold,\n",
    "                min_threshold,\n",
    "                steps_threshold\n",
    "            )\n",
    "        list_threshold.append(threshold)\n",
    "        \n",
    "    # print data\n",
    "    #print(list_threshold)\n",
    "    \n",
    "    ## Plot:\n",
    "    x_list = list_signal_mu\n",
    "    y_list = list_threshold\n",
    "    title = \"best threshold vs signal mu\"\n",
    "    prefix_filename = \"Threshold_vs_signal_mu_with_BC\" + str(number_of_BC)\n",
    "    xlabel = \"signal mu\"\n",
    "    ylabel = \"threshold\"\n",
    "    make_integer_axis = False\n",
    "    make_margin = False\n",
    "    plot_the_data(title, prefix_filename, xlabel, ylabel, make_integer_axis, make_margin, x_list, y_list)\n",
    "    \n",
    "\n",
    "def plot_best_threshold_deviation_vs_BC_SLOW():\n",
    "    \n",
    "    ratio_of_channels_to_sample = 1.0  # scale of simulation,   0.001 = small fast,   1.0 = full scale slow.\n",
    "    \n",
    "    # how many simulations to sample for the deviation. - larger = more accurate\n",
    "    simulation_samples = 1\n",
    "    \n",
    "    # initial number and steps\n",
    "    number_of_BC = nprocs * 1   # (works best with a multiple of [nprocs], nprocs is about equal to 4 or 8)\n",
    "    BC_steps = 3\n",
    "    \n",
    "    # intensity of collisions, average collisions per bunch crossing\n",
    "    average_collisions_mu = 200\n",
    "    \n",
    "    # threshold searches in the simulation (what thresholds to test deviation)\n",
    "    min_threshold = 20\n",
    "    max_threshold = 30\n",
    "    steps_threshold = 10   # resolution,  it takes longer time if you increase resolution. but the graph is better\n",
    "    \n",
    "    # signal mu to loop,  (could also edit for fixed signal mu and loop average_collisions_mu instead)\n",
    "    signal_mu = 100\n",
    "    \n",
    "    # create lists\n",
    "    list_threshold = []\n",
    "    list_BC = []\n",
    "    \n",
    "    # loop for the number of BC to use (numBC)  and get the threshold for each BC\n",
    "    # NOTE: THIS IS THE SLOW METHOD!  SEE THE OTHER FUNCTION BELOW \n",
    "    for i_numBC in range(BC_steps):\n",
    "        numBC = number_of_BC + 48*i_numBC\n",
    "        list_BC.append(numBC)\n",
    "        \n",
    "        # calculate the threshold and deviaiton\n",
    "        threshold, deviation = compute_best_statistical_threshold_with_least_deviation_from_collision_and_signal_mu(\n",
    "                simulation_samples,\n",
    "                numBC,\n",
    "                ratio_of_channels_to_sample,\n",
    "                average_collisions_mu,\n",
    "                signal_mu,\n",
    "                max_threshold,\n",
    "                min_threshold,\n",
    "                steps_threshold\n",
    "            )\n",
    "        list_threshold.append(threshold)\n",
    "        \n",
    "    # print data\n",
    "    print(list_threshold)\n",
    "    \n",
    "    ## Plot:\n",
    "    x_list = list_BC\n",
    "    y_list = list_threshold\n",
    "    title = \"best threshold vs number of bunch crossings\"\n",
    "    prefix_filename = \"Threshold_vs_numbunch_slow\"\n",
    "    xlabel = \"number of bunch crossings\"\n",
    "    ylabel = \"threshold\"\n",
    "    make_integer_axis = False\n",
    "    make_margin = False\n",
    "    plot_the_data(title, prefix_filename, xlabel, ylabel, make_integer_axis, make_margin, x_list, y_list)\n",
    "    \n",
    "    \n",
    "    #----------------------------------------------------------------------------------------------\n",
    "def plot_best_threshold_deviation_vs_BC_custom_FAST():\n",
    "    \n",
    "    ratio_of_channels_to_sample = 1.0  # scale of simulation,   0.001 = small fast,   1.0 = full scale slow.\n",
    "    \n",
    "    # how many simulations to sample for the deviation. - larger = more accurate\n",
    "    simulation_samples = 5\n",
    "    \n",
    "    # intensity of collisions, average collisions per bunch crossing\n",
    "    average_collisions_mu = 200\n",
    "    signal_to_channel_count_ratio = get_signal_to_channel_count_ratio(collisions_per_BC = average_collisions_mu)\n",
    "    \n",
    "    # threshold searches in the simulation (what thresholds to test deviation)\n",
    "    min_threshold = 21\n",
    "    max_threshold = 29\n",
    "    steps_threshold = 8*4   # resolution,  it takes longer time if you increase resolution. but the graph is better\n",
    "    \n",
    "    # signal mu to loop,  (could also edit for fixed signal mu and loop average_collisions_mu instead)\n",
    "    real_signal_mu = 100\n",
    "    guess_signal_mu = 100\n",
    "    \n",
    "    # number of BC (should be multiple of 8)\n",
    "    #number_of_BC = 16*2\n",
    "    #BC_steps = 15+1\n",
    "    number_of_BC = 8000\n",
    "    extra_BC = 800\n",
    "    BC_steps = 10\n",
    "    \n",
    "    # the thresholds to use\n",
    "    list_of_thresholds_to_try = list(linspace(min_threshold, max_threshold, steps_threshold))\n",
    "    \n",
    "    # long cache list of simulations\n",
    "    list_of_results_i_threshold_i_samples = np.zeros(len(list_of_thresholds_to_try), dtype=np.ndarray)\n",
    "    \n",
    "    # simulate for each parameter and store in the chache\n",
    "    for index, threshold in enumerate(list_of_thresholds_to_try):\n",
    "        \n",
    "        # makes diviation list\n",
    "        simulations = np.zeros(simulation_samples, dtype=np.ndarray)\n",
    "\n",
    "        # makes multiple simulations and stores the deviation for each simulation\n",
    "        for i in range(simulation_samples):\n",
    "            \n",
    "            result_actual_signal_count, result_estimated_signal_count = simulation_parallel_get_counts(number_of_BC + extra_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, threshold, real_signal_mu, guess_signal_mu)\n",
    "            simulations[i] = np.array([result_actual_signal_count, result_estimated_signal_count])\n",
    "            \n",
    "        list_of_results_i_threshold_i_samples[index] = simulations\n",
    "        \n",
    "    \n",
    "    #print(list_of_results_i_threshold_i_samples)\n",
    "    \n",
    "    list_BC = []\n",
    "    list_threshold = []\n",
    "    \n",
    "    # extract data from the cache simulations (list_of_results_i_threshold_i_samples)\n",
    "    for i_BC in range(1, BC_steps+1):\n",
    "        BC = int(number_of_BC/BC_steps * i_BC)\n",
    "        list_BC.append(BC)\n",
    "        \n",
    "        deviations_per_t = np.zeros(len(list_of_thresholds_to_try))\n",
    "        \n",
    "        # make random table of index to extract random BC from\n",
    "        N = number_of_BC + extra_BC\n",
    "        K = BC # K ones, N-K zeros\n",
    "        arr = np.array([1] * K + [0] * (N-K))\n",
    "        \n",
    "        # for each threshold extract\n",
    "        for index, threshold in enumerate(list_of_thresholds_to_try):\n",
    "        \n",
    "            deviations = np.zeros(simulation_samples)\n",
    "            \n",
    "            for i in range(simulation_samples):\n",
    "                # extract BCs\n",
    "                np.random.shuffle(arr)\n",
    "                result_actual_signal_count = list_of_results_i_threshold_i_samples[index][i][0][arr]\n",
    "                result_estimated_signal_count = list_of_results_i_threshold_i_samples[index][i][1][arr]\n",
    "                \n",
    "                deviations[i] = get_total_deviation_from_expected(result_actual_signal_count, result_estimated_signal_count)\n",
    "\n",
    "            the_mean_deviation = deviations.mean()\n",
    "            \n",
    "            deviations_per_t[index] = the_mean_deviation\n",
    "            \n",
    "        size = 5  # default is 3,  might want to change to 3\n",
    "        deviations_t_blur = uniform_filter(deviations_per_t, size=size)#[1:-1] # nah no need to slice, the boundary error doesn't bother us\n",
    "\n",
    "        # get the minimum\n",
    "        index = np.argmin(deviations_t_blur)\n",
    "        threshold = list_of_thresholds_to_try[index]\n",
    "        deviation = deviations_per_t[index]\n",
    "        \n",
    "        # store the result\n",
    "        list_threshold.append(threshold)\n",
    "        \n",
    "    \n",
    "    \n",
    "        \n",
    "    # print data\n",
    "    print(list_threshold)\n",
    "    \n",
    "    ## Plot:\n",
    "    x_list = list_BC\n",
    "    y_list = list_threshold\n",
    "    title = \"best threshold vs number of bunch crossings\"\n",
    "    prefix_filename = \"Threshold_vs_numbunch_fast\"\n",
    "    xlabel = \"number of bunch crossings\"\n",
    "    ylabel = \"threshold\"\n",
    "    make_integer_axis = False\n",
    "    make_margin = False\n",
    "    plot_the_data(title, prefix_filename, xlabel, ylabel, make_integer_axis, make_margin, x_list, y_list)\n",
    "\n",
    "    \n",
    "def plot_model_error_VS_threshold_VS_deviation():\n",
    "    ratio_of_channels_to_sample = 1.0  # scale of simulation,   0.001 = small fast,   1.0 = full scale slow.\n",
    "    \n",
    "    # how many simulations to sample for the deviation. - larger = more accurate\n",
    "    simulation_samples = 1 * 20\n",
    "    \n",
    "    # initial number and steps\n",
    "    number_of_BC = nprocs * 2 * 20   # (works best with a multiple of [nprocs], nprocs is about equal to 4 or 8)\n",
    "    \n",
    "    # intensity of collisions, average collisions per bunch crossing\n",
    "    average_collisions_mu = 0.01\n",
    "    signal_to_channel_count_ratio = get_signal_to_channel_count_ratio(collisions_per_BC = average_collisions_mu)\n",
    "    \n",
    "    # threshold searches in the simulation (what thresholds to test deviation)\n",
    "    min_threshold = 15\n",
    "    max_threshold = 35\n",
    "    steps_threshold = 10*2 +1   # resolution,  it takes longer time if you increase resolution. but the graph is better\n",
    "    \n",
    "    # signal mu to loop,  (could also edit for fixed signal mu and loop average_collisions_mu instead)\n",
    "    real_signal_mu = 100\n",
    "    guess_pm_signal_mu = 5*2\n",
    "    \n",
    "    \n",
    "    fig = plt.figure(figsize=(10, 8))\n",
    "    ax = fig.gca(projection='3d')\n",
    "    \n",
    "    div_list_list = []\n",
    "    \n",
    "    for guess_signal_mu in range(real_signal_mu-guess_pm_signal_mu, real_signal_mu+guess_pm_signal_mu+1):\n",
    "        div_list = []\n",
    "        for threshold in linspace(min_threshold, max_threshold, steps_threshold):\n",
    "            deviation = compute_average_deviation_from_simulations(\n",
    "                    simulation_samples,\n",
    "                    number_of_BC,\n",
    "                    ratio_of_channels_to_sample,\n",
    "                    signal_to_channel_count_ratio,\n",
    "                    real_signal_mu,\n",
    "                    guess_signal_mu,\n",
    "                    threshold\n",
    "                )\n",
    "            #print(deviation)\n",
    "            deviation += 0.000001\n",
    "            div_list.append(math.log(deviation))\n",
    "        div_list_list.append(div_list)\n",
    "            \n",
    "            #some_plot.SetPoint(point_index, guess_signal_mu, threshold, deviation) # i, x, y, z\n",
    "            #point_index += 1\n",
    "    \n",
    "    \n",
    "    #print(point_index)\n",
    "    #canvas = ROOT.TCanvas(\"I dont even know where this will show up\", \"not this either\", 600, 400)\n",
    "    #some_plot.Draw(\"TRI1\")\n",
    "    #some_plot.SetTitle(\"Real sigma mu = \"+ str(real_signal_mu) + \"; model sigma mu; treshold; logarithmic deviation\");\n",
    "    #canvas.Draw()\n",
    "    \n",
    "    #X = np.arange(-5, 5, 1)\n",
    "    #Y = np.arange(-5, 5, 1)\n",
    "    #X, Y = np.meshgrid(X, Y)\n",
    "    #R = np.sqrt(X**2 + Y**2)\n",
    "    #Z = np.sin(R)\n",
    "    Z = np.array(div_list_list)\n",
    "    X = np.array(list(range(real_signal_mu-guess_pm_signal_mu, real_signal_mu+guess_pm_signal_mu+1)))\n",
    "    Y = np.array(list(linspace(min_threshold, max_threshold, steps_threshold)))\n",
    "    X, Y = np.meshgrid(X, Y)\n",
    "    \n",
    "    print(\"X\", X.shape, X)\n",
    "    print(\"Y\", Y.shape, Y)\n",
    "    print(\"Z\", Z.shape, Z)\n",
    "    \n",
    "    surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm, linewidth=0, antialiased=False)\n",
    "    \n",
    "    # Customize the z axis.\n",
    "    #ax.set_zlim(0, 0.2)\n",
    "    #ax.zaxis.set_major_locator(LinearLocator(10))\n",
    "    #ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))\n",
    "    \n",
    "    fig.colorbar(surf, shrink=0.5, aspect=5)\n",
    "    \n",
    "    plt.show()\n",
    "    \n",
    "    print(\"Done\")\n",
    "    \n",
    "    \n",
    "    \n",
    "def plot_the_data(title, prefix_filename, xlabel, ylabel, make_integer_axis, make_margin, x_list, y_list):\n",
    "    timestr = time.strftime(\"%Y%m%d-%H%M%S\")\n",
    "    FILENAME = prefix_filename + \"_(\" + timestr + \").pdf\"\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    plt.title(title)\n",
    "    plt.xlabel(xlabel)\n",
    "    plt.ylabel(ylabel)\n",
    "    if make_integer_axis:\n",
    "        plt.locator_params(axis=\"both\", integer=True) # integer axis\n",
    "    if make_margin:\n",
    "        plt.margins(y=1)  # add margin to y\n",
    "    \n",
    "    plt.plot(x_list, y_list)\n",
    "    plt.show()\n",
    "    ax = fig.add_subplot(111)\n",
    "    fig.savefig(FILENAME)\n",
    "    #fig.savefig(FILENAME, dpi=200)  # use for png\n",
    "    \n",
    "\n",
    "    \n",
    "# -End of plot code-\n",
    "    \n",
    "#==============================================================================\n",
    "\n",
    "### Helper functions:\n",
    "\n",
    "def eq_div(N, i):\n",
    "    return int([N // i + 1] * (N % i) + [N // i] * (i - N % i))\n",
    "def eq_div_strict(N, i):\n",
    "    return [N // i] * i\n",
    "def np_stats(np_list):\n",
    "    \"\"\"Use a list from e.g. simulaiton_combined_get_counts and get statistics about that list\"\"\"\n",
    "    return np.mean(np_list), np.var(np_list), np.std(np_list)\n",
    "def print_count_statistics(np_list):\n",
    "    \"\"\"Use a list from e.g. simulaiton_combined_get_counts and get statistics about that list\"\"\"\n",
    "    mean, variance, standard_deviation = np_stats(np_list)\n",
    "    print(\"mean:\", mean)\n",
    "    print(\"variance:\", variance)\n",
    "    print(\"standard_deviation:\", standard_deviation)\n",
    "def get_signal_sigma_for_mu(signal_mu):\n",
    "    \"\"\"Using references from the short paper, to calculate sigma\"\"\"\n",
    "    reference_singal_mu = 100  # from s.15  (at threshold 40)\n",
    "    reference_signal_sigma = 25  #mv  25 from reference for signal_mu = 100 (s.24)\n",
    "    minimum = 0.01\n",
    "    signal_sigma = minimum + reference_signal_sigma * (signal_mu/reference_singal_mu) \n",
    "    return signal_sigma\n",
    "def get_signal_to_channel_count_ratio(collisions_per_BC = 100, particles_per_collision = 75, ratio_of_particles_reaching_HGTD=0.1):\n",
    "    \"\"\"Determines how many signals there will be,\n",
    "        based from  collisions_per_BC,  particles_per_collision,  ratio_of_particles_reaching_HGTD\"\"\"\n",
    "    #collisions_per_BC = 100\n",
    "    #particles_per_collision = 75\n",
    "    #ratio_of_particles_reaching_HGTD = 0.1\n",
    "    signal_to_channel_count_ratio = collisions_per_BC * particles_per_collision * ratio_of_particles_reaching_HGTD / C_total_channels\n",
    "    return signal_to_channel_count_ratio\n",
    "def get_number_of_channel_readouts(ratio_of_channels_to_sample):\n",
    "    \"\"\"Returns the number of channels to read from - used to downscale simulation\"\"\"\n",
    "    return int(C_total_channels * ratio_of_channels_to_sample)\n",
    "\n",
    "def linspace(start, end, steps):\n",
    "    scale = (end-start)/(steps)\n",
    "    return map(lambda x : start + x*scale, range(steps)) # does not include the last value, since steps = length\n",
    "\n",
    "\n",
    "# -End of helper functions-\n",
    "\n",
    "#==============================================================================\n",
    "\n",
    "###  Main code - what to run, what to plot, what do calculate\n",
    "\n",
    "def main():\n",
    "    # code... what to plot, what do calculate, ...\n",
    "    \n",
    "    #plot_model_error_VS_threshold_VS_deviation()\n",
    "    \n",
    "    #plot_best_threshold_deviation_vs_BC_custom_FAST()\n",
    "    \n",
    "    plot_deviation_vs_SB()\n",
    "    #plot3D_threshold_and_SB_vs_deviation()\n",
    "    #plot_particle_intensity_mu_vs_best_threshold()\n",
    "    #plot_signal_mu_vs_best_threshold()\n",
    "    #plot_deviation_vs_threshold(0.01,0,0.03,100)\n",
    "    #plot_best_threshold_deviation_vs_signal_mu()\n",
    "    #plot_deviation_vs_threshold_2()\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    # -code end-\n",
    "    return\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### Not important to understand ... it's just a startup script, see main() above instead.\n",
    "if __name__ == '__main__':\n",
    "    mp.freeze_support()\n",
    "    \n",
    "    # profiles to determine compute resources\n",
    "    profiles = {\"normal\": 2, \"fast\": 4, \"fastest\": 99}\n",
    "    max_processes = profiles[\"fastest\"]  # edit this\n",
    "    \n",
    "    # setup pool\n",
    "    nprocs = min(mp.cpu_count(), max_processes)\n",
    "    pool = None\n",
    "    with mp.Pool(processes=nprocs) as poool:\n",
    "        pool = poool\n",
    "    \n",
    "        global_counter = 0\n",
    "\n",
    "        # run code\n",
    "        t0 = time.time()\n",
    "        main()\n",
    "        t1 = time.time()\n",
    "        print(\"Simulation took\", round(t1-t0, 1), \"seconds\") # real time not the actuall process_time()\n",
    "\n",
    "        # clean up\n",
    "        pool.close()\n",
    "        pool.join()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "====================================================\n",
    "\n",
    "# Notes\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vi vill mäta antalet proton-proton-kollisioner - men vi kan inte se själva kollisionerna!\n",
    "I varje kollision skapas allt ifrån 5 till 500 partiklar - variationen är stor, men antalet\n",
    "följer en stabil distribution. En del av dessa partiklar flyger i väg mot HGTD och kan mätas där.\n",
    "\n",
    "Om vi medelvärdesbildar över en tid är antalet partiklar som når HGTD direkt\n",
    "proportionellt mot luminositeten!\n",
    "Vi behöver alltså bara räkna antalet träffar i vår detektor\n",
    "för att se hur luminositeten varierar under en körning.\n",
    "\n",
    "\n",
    "Att mäta luminositet\n",
    "\n",
    "Luminositeten behöver mätas både vid max (ca 200) och min (ca 0.01) antal\n",
    "kollisioner per “bunch crossing” (dvs då klungor av 10^11 protoner korsar\n",
    "varandra). Fördelningen av antalet kollisioner i en bunch crossing (BC) definieras\n",
    "av en Poisson-fördelning:\n",
    "N ~Po(u)\n",
    "Här motsvarar μ det genomsnittliga antalet kollisioner per BC. Var 25:e\n",
    "nanosekund sker en BC, så vid μ=200 sker kollisioner med en frekvens av 40MHz.\n",
    "\n",
    "\n",
    "Alltså\n",
    "1. Sätt collisions_mu till 200 (alt 0.01)\n",
    "2. Välj ett lämpligt expected voltage signal_mu = 100\n",
    "3. Räkna ut bästa threshold.  (se till att det används ett mycket bra val för många simulationer / deviaitons)\n",
    "4. Medelvärdes bilda med number_of_BC = [stort tal]\n",
    "    simulation(number_of_BC, collisions_mu, signal_mu, threshold)\n",
    "5. Räkna  total_signal_estimate / ratio_reaching_HGTD / time  =  estimated_luminosity\n",
    "(luminositeten borde ha enheten signals per second)\n",
    "6. Jämför estimated_luminosity  med  collisions_mu/(time_BC)\n",
    "    eller  estimated_luminosity*time_BC  med  collisions_mu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "====================================================\n",
    "\n",
    "\n",
    "<h1>Testing</h1>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}