Simulation_new_hanna.ipynb

{
 "cells": [
  {
   "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 numpy as np\n",
    "#from root_numpy import random_sample\n",
    "from scipy.ndimage import uniform_filter\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import multiprocessing as mp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Simulation Model</h1>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'mp' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-48f794396bb3>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m    986\u001b[0m \u001b[0;31m### Not important to understand ... it's just a startup script, see main() above instead.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    987\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0m__name__\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'__main__'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 988\u001b[0;31m     \u001b[0mmp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfreeze_support\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    989\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    990\u001b[0m     \u001b[0;31m# profiles to determine compute resources\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'mp' is not defined"
     ]
    }
   ],
   "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",
    "        #for i in range(total):\n",
    "        #    if prng.Uniform(25) < 3.125:\n",
    "        #        narrow_exact += 1\n",
    "        # faster to use numpy\n",
    "        narrow_exact = (np.random.uniform(low=0, high=25, size=total) < 3.125).sum()\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, 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 = expected_voltage\n",
    "        signal_func = ROOT.TF1(\"sth\", \"TMath::Landau(x, %f, %f)\" % (expected_voltage, 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, 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",
    "    peak_signal_voltage = signal_mu\n",
    "    \n",
    "    noise_gen = noise_generator(threshold, channel_count)\n",
    "    signal_gen = signal_generator(threshold, peak_signal_voltage)\n",
    "    \n",
    "    result_actual_signal_count = np.zeros(bunch_crossing_count)\n",
    "    result_estimated_signal_count = np.zeros(bunch_crossing_count)\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_gen.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, 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",
    "                    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",
    "\n",
    "### CALCULATIONS\n",
    "\n",
    "# ------------------------------------------------------------------------------------------------------------\n",
    "\n",
    "### Deviation\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",
    "    \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 the total deviation as a positive number\"\"\"\n",
    "    total_deviation = abs( result_estimated_signal_count.sum() / result_actual_signal_count.sum() - 1 )\n",
    "    return total_deviation\n",
    "    \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, 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, 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_listed_deviation_from_expected_with_simulated_results():\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",
    "    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, signal_mu)\n",
    "    return get_listed_deviation_from_expected(result_actual_signal_count, result_estimated_signal_count)\n",
    "\n",
    "\n",
    "def get_list_of_best_threshold_in_terms_of_least_deviation(list_of_thresholds_to_try, max_number_of_BC, ratio_of_channels_to_sample, signal_to_channel_count_ratio, signal_mu):\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_listed_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",
    "        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",
    "    \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, 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",
    "        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",
    "        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, 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, signal_mu)\n",
    "    return threshold\n",
    "\n",
    "\n",
    "# -------------------------------------------------------------------------------------------------------------\n",
    "\n",
    "### SB\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_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",
    "# -------------------------------------------------------------------------------------------------------------\n",
    "\n",
    "### Signal voltage\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",
    "# =============================================================================================================\n",
    "\n",
    "### Computations\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",
    "    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)[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) ,  \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, average_collisions_mu, 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 absolute deviation from those \n",
    "        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",
    "    # 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",
    "    # 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",
    "                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",
    "    \n",
    "def compute_maximum_and_mean_deviation_from_simulations(simulation_samples, number_of_BC, ratio_of_channels_to_sample, average_collisions_mu, signal_mu, threshold):\n",
    "    \"\"\"Returns the maximum deviation, mean deviation and fraction of deviations smaller than one percent\"\"\"\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",
    "    # 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",
    "                signal_mu\n",
    "            )\n",
    "        print(i)\n",
    "    \n",
    "    # takes the mean deviation found,  other intesting statistics can be made: e.g.  variance:  deviations.var()\n",
    "    the_maximum_deviation = deviations.max()   # maximum deviation\n",
    "    the_mean_deviation = deviations.mean()     # mean deviation\n",
    "    fraction_of_deviations_smaller_than_one_percent = (deviations < 0.01).sum()/simulation_samples\n",
    "    \n",
    "    # print statistics\n",
    "    #print_count_statistics(deviations)  # commented because not needed.\n",
    "    \n",
    "    return the_maximum_deviation, the_mean_deviation, fraction_of_deviations_smaller_than_one_percent\n",
    "    \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, 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",
    "                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",
    "        print(i)  # Not needed, but makes it easier to see how far the simulation has come\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",
    "\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 = mp.cpu_count() * 4# * 120\n",
    "    \n",
    "    signal_to_channel_count_ratio = get_signal_to_channel_count_ratio()\n",
    "    noise_sigma = 5  #mv\n",
    "    \n",
    "    min_signal_mu = 8  #mv  # cannot be below 0\n",
    "    max_signal_mu = 35 #mv\n",
    "    steps_signal_mu = (max_signal_mu - min_signal_mu) * 2   # resolution\n",
    "    #steps_signal_mu = 20\n",
    "    \n",
    "    min_threshold = 3  #mv  # cannot be below 0\n",
    "    max_threshold = 80 #mv\n",
    "    steps_threshold = (max_threshold - min_threshold) * 2   # 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_deviation_from_list_of_signal_mu(\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 = 10*10\n",
    "    for _ in range(runs-1):\n",
    "        deviations += compute_list_of_least_deviation_from_list_of_signal_mu(\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",
    "    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_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_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",
    "    \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",
    "                average_collisions_mu,\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_vs_number_of_BC():\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 = 20\n",
    "    \n",
    "    # number of BC for each single simulation.\n",
    "    min_number_of_BC = nprocs * 100   # (works best with a multiple of [nprocs], nprocs is about equal to 4 or 8)\n",
    "    max_number_of_BC = nprocs * 1001\n",
    "    step_size_number_of_BC = nprocs * 400 # (works best with a multiple of [nprocs], nprocs is about equal to 4 or 8)\n",
    "    list_number_of_BC = list(range(min_number_of_BC, max_number_of_BC, step_size_number_of_BC))\n",
    "    \n",
    "    # Mean energy\n",
    "    signal_mu = 100\n",
    "    \n",
    "    # intensity of collisions, average collisions per bunch crossing\n",
    "    average_collisions_mu = 0.01\n",
    "    \n",
    "    # threshold searches in the simulation (what thresholds to test deviation)\n",
    "    min_threshold = 20\n",
    "    max_threshold = 40\n",
    "    steps_threshold = 80   # 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",
    "    list_deviation = []\n",
    "    \n",
    "    # here we loop over number_of_BC\n",
    "    # could also edit the code to have a fixed value singnal_mu and instead loop average_collisions_mu\n",
    "    for number_of_BC in list_number_of_BC:\n",
    "        # calculate the threshold and deviation\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",
    "        print(number_of_BC)\n",
    "        list_threshold.append(threshold)\n",
    "        list_deviation.append(deviation)\n",
    "        \n",
    "    # print data\n",
    "    print(list_threshold)\n",
    "    print(list_deviation)\n",
    "    \n",
    "    ## Plot threshold:\n",
    "    x_list = list_number_of_BC\n",
    "    y_list = list_threshold\n",
    "    title = \"best threshold vs number of BC, mu=\" + str(average_collisions_mu)\n",
    "    prefix_filename = \"Threshold_vs_number_of_BC_with_BC\" + str(number_of_BC) + \"_sim_sampl\" + str(simulation_samples)\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 print_best_threshold():\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 = 100\n",
    "    \n",
    "    # number of BC for each single simulation.\n",
    "    number_of_BC = nprocs * 400   # (works best with a multiple of [nprocs], nprocs is about equal to 4 or 8), at least 400 for average_collisions_mu = 0.01, enough with 40 for average_collisions_mu = 200\n",
    "    \n",
    "    # Mean energy\n",
    "    signal_mu = 100\n",
    "    \n",
    "    # intensity of collisions, average collisions per bunch crossing\n",
    "    average_collisions_mu = 0.01\n",
    "    \n",
    "    # threshold searches in the simulation (what thresholds to test deviation)\n",
    "    min_threshold = 25\n",
    "    max_threshold = 40\n",
    "    steps_threshold = 80   # resolution,  it takes longer time if you increase resolution. but the graph is better\n",
    "\n",
    "    # calculate the threshold and deviation\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",
    "        \n",
    "    print('Optimal threshold for average_collisions_mu = '+str(average_collisions_mu)+' is '+str(threshold))\n",
    "    print('with a mean absolute deviation of '+str(deviation))\n",
    "    \n",
    "\n",
    "def print_maximum_mean_and_fraction_deviation():\n",
    "    \"\"\"Print maximum deviation, mean deviation and fraction of deviations smaller than 1%\"\"\"\n",
    "    simulation_samples = 10000\n",
    "    number_of_BC = nprocs * 400\n",
    "    ratio_of_channels_to_sample = 1.0\n",
    "    average_collisions_mu = 0.01\n",
    "    signal_mu = 100\n",
    "    threshold = 32\n",
    "    \n",
    "    maximum_deviation, mean_deviation, fraction = compute_maximum_and_mean_deviation_from_simulations(simulation_samples, number_of_BC, ratio_of_channels_to_sample, average_collisions_mu, signal_mu, threshold)\n",
    "    \n",
    "    print('For average_collisions_mu='+str(average_collisions_mu)+', threshold='+str(threshold)+' mV, '+str(simulation_samples)+' simulations, signal_mu='+str(signal_mu)+', BC='+str(number_of_BC))\n",
    "    print('maximum deviation = '+str(maximum_deviation*100)+'%')\n",
    "    print('mean deviation = '+str(mean_deviation*100)+'%')\n",
    "    print('fraction of deviations smaller than 1%: '+str(fraction))\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 = 50, particles_per_collision = 500, 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 = 50\n",
    "    #particles_per_collision = 500\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_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_best_threshold_deviation_vs_signal_mu()\n",
    "#     plot_deviation_vs_threshold_2()\n",
    "#     plot_best_threshold_vs_number_of_BC()\n",
    "    print_maximum_mean_and_fraction_deviation()\n",
    "#     print_best_threshold()\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 = mp.Pool(processes=nprocs)\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": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "====================================================\n",
    "\n",
    "\n",
    "<h1>Testing</h1>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(signal_count):\n",
    "            if prng.Landau(expected_voltage, expected_voltage/4) > threshold: # 1/4 is just the refernce 25/100\n",
    "                above_threshold_count += 1"
   ]
  },
  {
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}