How to simulate noise using both time windows.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Large time window 25ns\n",
    "\n",
    "small time window 3.125ns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to JupyROOT 6.22/02\n"
     ]
    }
   ],
   "source": [
    "import ROOT\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Class definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class noise_generator:\n",
    "    def __init__(self, threshold, channel_count):\n",
    "        noise = ROOT.TF1(\"noise_gauss\", \"TMath::Gaus(x, 0, 5)\", 0, 1000)\n",
    "        self._mean = channel_count * noise.Integral(threshold, 1000)/noise.Integral(0, 1000)/2\n",
    "        self._prng = ROOT.TRandom3(time.time_ns())\n",
    "\n",
    "    def generate(self):\n",
    "        total = self._prng.Poisson(self._mean)\n",
    "        narrow = 0\n",
    "        for i in range(total):\n",
    "            if self._prng.Uniform(25) < 3.125:\n",
    "                narrow += 1\n",
    "        estimate = 3.125 * (total - narrow)/(25 - 3.125)\n",
    "        return narrow, estimate\n",
    "\n",
    "    \n",
    "class signal_generator:\n",
    "    def __init__(self, threshold, expected_voltage):\n",
    "        self._threshold = threshold\n",
    "        self._expected_voltage = expected_voltage\n",
    "        signal_func = ROOT.TF1(\"sth\", \"TMath::Landau(x, %f, %f)\" % (expected_voltage, expected_voltage/4), -100, 1000)\n",
    "        self.fraction_above = signal_func.Integral(threshold, 1000)/signal_func.Integral(-100, 1000)\n",
    "        self._prng = ROOT.TRandom3(time.time_ns())\n",
    "    \n",
    "    def generate(self, expected_count):\n",
    "        signal_count = self._prng.Poisson(expected_count)\n",
    "        above_threshold_count = 0\n",
    "        for i in range(signal_count):\n",
    "            if self._prng.Landau(self._expected_voltage, self._expected_voltage/4) > self._threshold:\n",
    "                above_threshold_count += 1\n",
    "        return signal_count, above_threshold_count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Paramaters and return values explained"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "17 14.0\n",
      "36 36\n"
     ]
    }
   ],
   "source": [
    "# Create a noise generator by passing threshold and channel count as paramaters\n",
    "my_noise_generator = noise_generator(20, 3.6 * 1e6)\n",
    "\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.\n",
    "noise_exact, noise_est = my_noise_generator.generate()\n",
    "\n",
    "print(noise_exact, noise_est)\n",
    "\n",
    "# Create a signal generator, the first paramater is the threshold, the second paramater is the\n",
    "# expected peak voltage\n",
    "my_signal_generator = signal_generator(20, 100)\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 the second return\n",
    "# value is the number of activated channels.\n",
    "signal_total, signal_above = my_signal_generator.generate(30)\n",
    "\n",
    "print(signal_total, signal_above)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Minimum working example\n",
    "This is how I imagine this will be used for a simulation of bunch crossings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9994578844940053\n"
     ]
    }
   ],
   "source": [
    "channel_count = 3.6 * 1e6\n",
    "bunch_crossing_count = 10000\n",
    "threshold = 20\n",
    "peak_signal_voltage = 100\n",
    "avg_particles_reaching_decector = 30\n",
    "\n",
    "noise = noise_generator(threshold, channel_count)\n",
    "signal = signal_generator(threshold, peak_signal_voltage)\n",
    "\n",
    "actual_signal_sum = 0\n",
    "estimated_signal_sum = 0\n",
    "\n",
    "for i in range(bunch_crossing_count):\n",
    "    noise_exact, noise_est = noise.generate()\n",
    "    signal_total, signal_above = signal.generate(avg_particles_reaching_decector)\n",
    "    \n",
    "    total_activations = signal_above + noise_exact\n",
    "    signal_est = (total_activations - noise_est)/signal.fraction_above\n",
    "    \n",
    "    actual_signal_sum += signal_total\n",
    "    estimated_signal_sum += signal_est\n",
    "    \n",
    "print(estimated_signal_sum/actual_signal_sum)"
   ]
  }
 ],
 "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": 4
}