diff --git a/AdditionalPDFs.py b/AdditionalPDFs.py index 1c3885b5e698506a4e58163d84d673659f8943f5..2fed1c0aa498441681c930c8d7c05e32a2857771 100644 --- a/AdditionalPDFs.py +++ b/AdditionalPDFs.py @@ -1,30 +1,21 @@ import numpy as np from scipy.stats import rv_discrete, uniform import scipy.special as sc - - -from scipy.stats._distn_infrastructure import ( - rv_discrete, get_distribution_names) - +from scipy.stats._distn_infrastructure import (rv_discrete, get_distribution_names) class gpd_gen(rv_discrete): def _argcheck(self, mu, lbda): return mu >= 0.0 and lbda >= 0.0 and lbda <= 1.0 def _rvs(self, mu, lbda): - population = np.asarray( - self._random_state.poisson(mu, self._size) - ) + population = np.asarray(self._random_state.poisson(mu, self._size)) if population.shape == (): population = population.reshape(-1) offspring = population.copy() while np.any(offspring > 0): # probability dists are NOT ufuncs # print("offspring", offspring) - offspring[:] = [ - self._random_state.poisson(m) - for m in lbda*offspring - ] + offspring[:] = [self._random_state.poisson(m) for m in lbda*offspring] population += offspring return population @@ -41,9 +32,7 @@ class gpd_gen(rv_discrete): elif n == 2: return (mu/(1-lbda))**2+mu/(1-lbda)**3 - gpoisson = gpd_gen(name='gpoisson') - class borel_gen(rv_discrete): def _argcheck(self, mu): @@ -82,7 +71,6 @@ class borel_gen(rv_discrete): return _rnd - def _stats(self, mu): _mu = 1/(1-mu) _var = mu/(1-mu)**3 @@ -93,15 +81,9 @@ class borel_gen(rv_discrete): g2 = scipy._lib._util._lazywhere(mu_nonzero, (tmp,), lambda x: 3 + (1 + 8*x+6*x**2)/(x*(1-x)), np.inf) return _mu, _var, g1, g2 - borel= borel_gen(name='borel') - - - + class erlang_gen(rv_discrete): - - - def _pdf(self, x, a): # gamma.pdf(x, a) = x**(a-1) * exp(-x) / gamma(a) return np.exp(self._logpdf(x, a)) @@ -109,11 +91,6 @@ class erlang_gen(rv_discrete): def _logpdf(self, k, mu, nu): return sc.xlogy(a-1.0, x) - x - sc.gammaln(a) - - - - - # def _rvs(self, mu, nu, size=None, random_state=None): # u = scipy.stats.uniform.rvs(loc=0, scale = 1, size=size) # cum = np.cumsum([self._pmf(_k, mu, nu) for _k in range(0, 100)]) diff --git a/LightSimtastic.py b/LightSimtastic.py index e3c1651c82ff60d155ab513901dd422e18fcc069..9cdbb2d89f3195307d93d72063fd673a429848e0 100644 --- a/LightSimtastic.py +++ b/LightSimtastic.py @@ -1,3 +1,4 @@ +import warnings import time import numpy as np import pandas as pd @@ -7,293 +8,127 @@ import scipy.signal as signal from tqdm.auto import tqdm from itertools import chain, product, zip_longest from operator import add -from types import SimpleNamespace +from types import SimpleNamespace from AdditionalPDFs import borel import codecs np.set_printoptions(suppress=True) # evil but works: just suppress the numpy nested ragged array warning -import warnings if hasattr(np, 'VisibleDeprecationWarning'): warnings.filterwarnings("ignore", category=np.VisibleDeprecationWarning) else: warnings.filterwarnings("ignore", category=DeprecationWarning) -class SiPMSimulation: - - - +class SiPMSimulation: def __init__(self): - - '"Name_Sim"', '"mu"', '"ppXT"', '"pdXT"', '"taudXT"', '"rdXT"', '"pAp"', '"tauAp"', '"taur"', '"Sig0"', '"Sig1"', '"DCR"', '"Sigt"', '"GSig"', '"tslow"', '"tfast"', '"rfast"', '"start_tgate"', '"len_tgate"', '"t0"', '"tl0"', '"tl1"', '"tl2"', - - - - self.VariableDictionary ={ - "Name_Sim": { - "Description": "Name for the simulation.", - "Unit": None, - }, - - "mu": { - "Description": "Mean Number of Primary Geiger discharges in SiPM Cell", - "Unit": "Units of Charge, Normalized to Gain [Q/G]" - }, - "ppXT": { - "Description": "Characteristic Probability of Prompt Cross Talk in SiPM Cell", - "Unit": None - }, - "pdXT": { - "Description": "Characteristic Probability of Delayed Cross Talk in SiPM Cell", - "Unit": None - }, - "taudXT":{ - "Description": "Time Constant of Delayed Cross Talk in SiPM Cell", - "Unit": "Nanoseconds [ns]" - }, - "rdXT":{ - "Description": " Fraction of Delayed Cross Talk from Adjacent SiPM Cells", - "Unit": None - }, - "pAp":{ - "Description": " Probability of Afterpulse occurrence in SiPM Cell", - "Unit": None - }, - "tauAp":{ - "Description": " Time delay constant of Afterpulsing in SiPM Cell", - "Unit": None - }, - "taur":{ - "Description": " PDE Recovery constant in SiPM Cell", - "Unit": None - }, - "Sig0":{ - "Description": "Fluctuation in Gain of Geiger Discharge", - "Unit": "Units of Charge, Normalized to Gain [Q/G]" - }, - "Sig1":{ - "Description": " Pedestal width due to Noise from Electronic Components", - "Unit": "Units of Charge, Normalized to Gain [Q/G]" - }, - "DCR":{ - "Description": "Dark Count Rate", - "Unit": "Gigahertz [GHz]" - }, - "Sigt":{ - "Description": "Electronics Noise for Current Transient", - "Unit": "Nanoseconds [ns]" - }, - "GSig":{ - "Description": "Width of Electronics Transfer Function", - "Unit": "Nanoseconds [ns]" - }, - "tslow":{ - "Description": "Time constant of Slow Component of SiPM pulse shape.", - "Unit": "Nanoseconds [ns]" - - }, - "tfast":{ - "Description": "Time constant of Fast Component of SiPM pulse shape.", - "Unit": "Nanoseconds [ns]" - - }, - "rfast":{ - "Description": "Fraction of Fast Component contribution to SiPM pulse shape.", - "Unit": None - - }, - "start_tgate":{ - "Description": "Start time of SiPM Integration Window", - "Unit": "Nanoseconds [ns]" - - }, - "len_tgate":{ - "Description": "Duration of SiPM Integration Window", - "Unit": "Nanoseconds [ns]" - - }, - "t0":{ - "Description": "Duration before Dark Count Evaluation", - "Unit": "Nanoseconds [ns]" - - }, - "tl0":{ - "Description": "Free Parameter 0 for Light Pulse Distribution", - "Unit": "Nanoseconds [ns]" - - }, - - "tl1":{ - "Description": "Free Parameter 1 for Light Pulse Distribution", - "Unit": "Nanoseconds [ns]" - - }, - "tl2":{ - "Description": "Free Parameter 2 for Light Pulse Distribution", - "Unit": "Nanoseconds [ns]" - - }, - - "Gen_mu":{ - "Description": "Function for the random sampling of the number of Primary Gieger Discharges", - "Unit": None, - "Options": ["Poisson", "Gauss", "Fixed"] - - }, - - "Gen_tmu":{ - "Description":"Function for the random sampling of the time of primary Gieger discharges", - "Unit": None, - "Options": ["Gauss", "Exp"] - - }, - - "Gen_gain":{ - "Description":"Function for the random sampling of the SiPM Gain", - "Unit": None, - "Options": ["Gauss"] - - }, - - "Gen_npXT":{ - "Description":"Function for the random sampling of the prompt cross-talk", - "Unit": None, - "Options": ["Poisson", "Binom", "Borel"] - - }, - - "Gen_ndXT":{ - "Description":"Function for the random sampling of the number delayed cross-talk discharges", - "Unit": None, - "Options": ["Poisson", "Binom", "Borel"] - - }, - - "Gen_tdXT":{ - "Description":"Function for the random sampling of the time of delayed cross-talk discharges", - "Unit": None, - "Options": ["Exp"] - - } - + '"Name_Sim"', '"mu"', '"ppXT"', '"pdXT"', '"taudXT"', '"rdXT"', '"pAp"', '"tauAp"', '"taur"', '"Sig0"', '"Sig1"', '"DCR"', '"Sigt"', '"GSig"', '"tslow"', '"tfast"', '"rfast"', '"start_tgate"', '"len_tgate"', '"t0"', '"tl0"', '"tl1"', '"tl2"', + + self.VariableDictionary = { + "Name_Sim": {"Description": "Name for the simulation.", "Unit": None, }, + "mu": {"Description": "Mean Number of Primary Geiger discharges in SiPM Cell", "Unit": "Units of Charge, Normalized to Gain [Q/G]"}, + "ppXT": {"Description": "Characteristic Probability of Prompt Cross Talk in SiPM Cell", "Unit": None}, + "pdXT": {"Description": "Characteristic Probability of Delayed Cross Talk in SiPM Cell", "Unit": None}, + "taudXT": {"Description": "Time Constant of Delayed Cross Talk in SiPM Cell", "Unit": "Nanoseconds [ns]"}, + "rdXT": {"Description": " Fraction of Delayed Cross Talk from Adjacent SiPM Cells", "Unit": None}, + "pAp": {"Description": " Probability of Afterpulse occurrence in SiPM Cell", "Unit": None}, + "tauAp": {"Description": " Time delay constant of Afterpulsing in SiPM Cell", "Unit": None}, + "taur": {"Description": " PDE Recovery constant in SiPM Cell", "Unit": None}, + "Sig0": {"Description": "Fluctuation in Gain of Geiger Discharge", "Unit": "Units of Charge, Normalized to Gain [Q/G]"}, + "Sig1": {"Description": " Pedestal width due to Noise from Electronic Components", "Unit": "Units of Charge, Normalized to Gain [Q/G]"}, + "DCR": {"Description": "Dark Count Rate", "Unit": "Gigahertz [GHz]"}, + "Sigt": {"Description": "Electronics Noise for Current Transient", "Unit": "Nanoseconds [ns]"}, + "GSig": {"Description": "Width of Electronics Transfer Function", "Unit": "Nanoseconds [ns]"}, + "tslow": {"Description": "Time constant of Slow Component of SiPM pulse shape.", "Unit": "Nanoseconds [ns]"}, + "tfast": {"Description": "Time constant of Fast Component of SiPM pulse shape.", "Unit": "Nanoseconds [ns]"}, + "rfast": {"Description": "Fraction of Fast Component contribution to SiPM pulse shape.", "Unit": None}, + "start_tgate": {"Description": "Start time of SiPM Integration Window", "Unit": "Nanoseconds [ns]"}, + "len_tgate": {"Description": "Duration of SiPM Integration Window", "Unit": "Nanoseconds [ns]"}, + "t0": {"Description": "Duration before Dark Count Evaluation", "Unit": "Nanoseconds [ns]"}, + "tl0": {"Description": "Free Parameter 0 for Light Pulse Distribution", "Unit": "Nanoseconds [ns]"}, + "tl1": {"Description": "Free Parameter 1 for Light Pulse Distribution", "Unit": "Nanoseconds [ns]"}, + "tl2": {"Description": "Free Parameter 2 for Light Pulse Distribution", "Unit": "Nanoseconds [ns]"}, + "Gen_mu": {"Description": "Function for the random sampling of the number of Primary Gieger Discharges", "Unit": None, "Options": ["Poisson", "Gauss", "Fixed"]}, + "Gen_tmu": {"Description": "Function for the random sampling of the time of primary Gieger discharges", "Unit": None, "Options": ["Gauss", "Exp"]}, + "Gen_gain": {"Description": "Function for the random sampling of the SiPM Gain", "Unit": None, "Options": ["Gauss"]}, + "Gen_npXT": {"Description": "Function for the random sampling of the prompt cross-talk", "Unit": None, "Options": ["Poisson", "Binom", "Borel"]}, + "Gen_ndXT": {"Description": "Function for the random sampling of the number delayed cross-talk discharges", "Unit": None, "Options": ["Poisson", "Binom", "Borel"]}, + "Gen_tdXT": {"Description": "Function for the random sampling of the time of delayed cross-talk discharges", "Unit": None, "Options": ["Exp"]} } - - - - self.pb_ag_text=["Sampling Primary Geiger Distributions...", - "Sampling SiPM Noise Distributions...", - "Calculating Prompt Cross Talk...", - "Calculating Delayed Cross Talk...", - "Calculating Afterpulses...", - "Combining SiPM Noise...", - "All Geiger Discharges Calculated.", - ] - - self.pb_Q_text=["Sampling Charge Distributions...", - "Calculating Charge from Geiger Discharges...", - "Charge Spectrum Calculated." - - ] - - self.pb_I_text=["Sampling Current Distributions...", - "Calculating Noise from Geiger Discharges...", - "Convolving Electronic Noise Function..." - "Transient Calculated." - ] - - - self.out_columns = ["GeigerArray", - "ChargeSpectrum", - "Transients", - "TimeElapsed_GeigerArray", - "TimeElapsed_ChargeSpectrum", - "TimeElapsed_Transients", - "NSim" - ] - - self.pars = pd.DataFrame(columns = list(self.VariableDictionary.keys())+ self.out_columns) - - + self.pb_ag_text = ["Sampling Primary Geiger Distributions...", "Sampling SiPM Noise Distributions...", + "Calculating Prompt Cross Talk...", "Calculating Delayed Cross Talk...", + "Calculating Afterpulses...", "Combining SiPM Noise...", "All Geiger Discharges Calculated.",] + self.pb_Q_text = ["Sampling Charge Distributions...", + "Calculating Charge from Geiger Discharges...", "Charge Spectrum Calculated."] + self.pb_I_text = ["Sampling Current Distributions...", "Calculating Noise from Geiger Discharges...", + "Convolving Electronic Noise Function...", "Transient Calculated."] + self.out_columns = ["GeigerArray", "ChargeSpectrum", "Transients", "TimeElapsed_GeigerArray", + "TimeElapsed_ChargeSpectrum", "TimeElapsed_Transients", "NSim"] + self.pars = pd.DataFrame(columns=list( + self.VariableDictionary.keys()) + self.out_columns) self.pars["GeigerArray"] = self.pars["GeigerArray"].astype(object) self.pars["ChargeSpectrum"] = self.pars["ChargeSpectrum"].astype(object) self.pars["Transients"] = self.pars["Transients"].astype(object) - - - def ReadMathCad(self, inputfile): variables = {} with codecs.open(inputfile, encoding='utf-8') as file: def _SetFloating(string): try: - string=float(string) + string = float(string) return string except ValueError: return string - def _ToFloats(list_strings): return [_SetFloating(string) for string in list_strings] - - def _RemoveNonAscii(string): - return "".join(i for i in string if ord(i)<126 and ord(i)>31) - + def _RemoveNonAscii(string): + return "".join(i for i in string if ord(i) < 126 and ord(i) > 31) lines = file.read() - for line in lines.splitlines(True): line = _RemoveNonAscii(line) elems = line.replace('"', "").split(",") elems = _ToFloats(elems) - if elems[0] in self.VariableDictionary: variables[elems[0]] = elems[1::] - self.AddVariables(variables) - - - def ClearVariables(self): self.pars = self.pars.iloc[0:0] - + def AddVariables(self, variables): if not isinstance(variables, dict): - raise Exception("Expected a dictionary of parameters. Please enter a dictionary with all of the following elements: {0}.".format( - ", ".join(list(self.VariableDictionary.keys())))) - - - missing = np.setdiff1d(list(self.VariableDictionary.keys()),list(variables.keys())) - - if (len(missing)>0): - raise Exception("Variables are missing from the entered dictionary. Please add following elements: {0}.".format( - ", ".join(list(missing)))) - + raise Exception("Expected a dictionary of parameters. Please enter a dictionary with all of the following elements: {0}.".format(", ".join(list(self.VariableDictionary.keys())))) + + missing = np.setdiff1d( + list(self.VariableDictionary.keys()), list(variables.keys())) + + if (len(missing) > 0): + raise Exception("Variables are missing from the entered dictionary. Please add following elements: {0}.".format(", ".join(list(missing)))) + for key, value in variables.items(): if not isinstance(value, list): - raise Exception("Expected a list of simulation parameters, but a type {0} was found for variable {1}. Please enter your simulation parameters in a list.".format(type(value), key)) - + raise Exception("Expected a list of simulation parameters, but a type {0} was found for variable {1}. Please enter your simulation parameters in a list.".format(type(value), key)) + self.ClearVariables() - + for key, value in variables.items(): self.pars[key] = value - + def GetInformation(self, key): try: - temp = pd.DataFrame.from_dict(self.VariableDictionary[key],orient='index') + temp = pd.DataFrame.from_dict(self.VariableDictionary[key], orient='index') print(temp) except: - raise Exception("Parameter not found. Available parameters are: {0}".format( - ", ".join(list(self.VariableDictionary.keys())))) - + raise Exception("Parameter not found. Available parameters are: {0}".format(", ".join(list(self.VariableDictionary.keys())))) + def PrintVariables(self, latex=False): - if (latex): df_var = pd.DataFrame.from_dict(self.VariableDictionary) print(df_var.to_latex()) else: print(self.pars) - + def dXTtSplit(self, dXT_nt, dXT_n): output = [] partialSum = 0 @@ -305,485 +140,311 @@ class SiPMSimulation: output.append(temp) return output - - def AllGeiger(self, ns): - pbar = tqdm(range(len(self.pb_ag_text)-1)) pbar.set_description("{0}".format(self.pb_ag_text[0])) - nlight = np.asarray(ns.dist_flight.rvs(size=ns.n_sim)) nDC = np.asarray(ns.dist_fDC.rvs(size=ns.n_sim)) - - if(len(nlight)>len(nDC)): + if (len(nlight) > len(nDC)): nDC.resize(nlight.shape) - elif (len(nDC)>len(nlight)): + elif (len(nDC) > len(nlight)): nlight.resize(nDC.shape) - nPG = nDC + nlight - pbar.update(1) pbar.set_description("{0}".format(self.pb_ag_text[1])) - + nlight_sum = np.sum(nlight) nDC_sum = np.sum(nDC) nPG_sum = nlight_sum + nDC_sum - + npXT = np.array(ns.dist_fpXT.rvs(size=nPG_sum)) ndXT = np.array(ns.dist_fdXT.rvs(size=nPG_sum)) - + npXT_sum = np.sum(npXT) - - nPG_pXT_sum = nPG_sum + npXT_sum - + nPG_pXT_sum = nPG_sum + npXT_sum pg_l_d = np.ones(nlight_sum) pg_l_t = np.asarray(ns.dist_ftlight.rvs(size=nlight_sum)) pg_dc_d = np.ones(nDC_sum) pg_dc_t = np.asarray(ns.dist_ftDC.rvs(size=nDC_sum)) - + nlight_cumsum = np.cumsum(nlight) nDC_cumsum = np.cumsum(nDC) nPG_cumsum = np.cumsum(nPG) ndXT_cumsum = np.cumsum(ndXT) - + pg_l_d = np.split(pg_l_d, nlight_cumsum)[:-1] pg_l_t = np.split(pg_l_t, nlight_cumsum)[:-1] pg_dc_d = np.split(pg_dc_d, nDC_cumsum)[:-1] pg_dc_t = np.split(pg_dc_t, nDC_cumsum)[:-1] - + pXT_n = np.split(npXT, nPG_cumsum)[:-1] dXT_n = np.split(ndXT, nPG_cumsum)[:-1] - + pg_d = [np.concatenate(elem) for elem in zip(pg_l_d, pg_dc_d)] pg_t = [np.concatenate(elem) for elem in zip(pg_l_t, pg_dc_t)] - pbar.update(1) pbar.set_description("{0}".format(self.pb_ag_text[2])) - - pXT_d = [ - np.array( - list( - chain.from_iterable( - [ - [elem_pg_d]*elem_pXT_n - if(elem_pXT_n > 0 and elem_pg_d >0) + + pXT_d = [np.array(list(chain.from_iterable([[elem_pg_d]*elem_pXT_n + if (elem_pXT_n > 0 and elem_pg_d > 0) else [] - for elem_pg_d, elem_pXT_n in zip(list_pg_d, list_pXT_n) - ] - ) - ) - ) - for list_pg_d, list_pXT_n in zip(pg_d, pXT_n) - ] - - - pXT_t = [ - np.array( - list( - chain.from_iterable( - [ - [elem_pg_t]*elem_pXT_n - if(elem_pXT_n > 0 and elem_pg_d >0) + for elem_pg_d, elem_pXT_n in zip(list_pg_d, list_pXT_n)]))) + for list_pg_d, list_pXT_n in zip(pg_d, pXT_n)] + + pXT_t = [np.array(list(chain.from_iterable([[elem_pg_t]*elem_pXT_n + if (elem_pXT_n > 0 and elem_pg_d > 0) else [] - for elem_pg_d, elem_pg_t, elem_pXT_n in zip(list_pg_d, list_pg_t, list_pXT_n) - ] - ) - ) - ) - for list_pg_d, list_pg_t, list_pXT_n in zip(pg_d, pg_t, pXT_n) - ] - - + for elem_pg_d, elem_pg_t, elem_pXT_n in zip(list_pg_d, list_pg_t, list_pXT_n)]))) + for list_pg_d, list_pg_t, list_pXT_n in zip(pg_d, pg_t, pXT_n)] + pbar.update(1) pbar.set_description("{0}".format(self.pb_ag_text[3])) - - - dXT_d = [ - np.array( - list(chain.from_iterable( - [ - [elem_pg_d]*elem_dXT_n - if(elem_dXT_n > 0 and elem_pg_d >0) - else [] - for elem_pg_d, elem_dXT_n in zip(list_pg_d, list_dXT_n) - ] - ) - ) - ) - for list_pg_d, list_dXT_n in zip(pg_d, dXT_n) - ] - - dXT_t = [ - np.array( - list( - chain.from_iterable( - [ - list(map(add, [elem_pg_t]*elem_dXT_n, ns.dist_ftdXT.rvs(size=elem_dXT_n))) - if(elem_dXT_n > 0 and elem_pg_d >0) + + dXT_d = [np.array(list(chain.from_iterable([[elem_pg_d]*elem_dXT_n + if (elem_dXT_n > 0 and elem_pg_d > 0) + else [] + for elem_pg_d, elem_dXT_n in zip(list_pg_d, list_dXT_n)]))) + for list_pg_d, list_dXT_n in zip(pg_d, dXT_n)] + + dXT_t = [np.array(list(chain.from_iterable([list(map(add, [elem_pg_t]*elem_dXT_n,ns.dist_ftdXT.rvs(size=elem_dXT_n))) + if (elem_dXT_n > 0 and elem_pg_d > 0) else [] - for elem_pg_d, elem_pg_t, elem_dXT_n in zip(list_pg_d, list_pg_t, list_dXT_n) - ] - ) - ) - ) - for list_pg_d, list_pg_t, list_dXT_n in zip(pg_d, pg_t, dXT_n) - ] - - + for elem_pg_d, elem_pg_t, elem_dXT_n in zip(list_pg_d, list_pg_t, list_dXT_n)]))) + for list_pg_d, list_pg_t, list_dXT_n in zip(pg_d, pg_t, dXT_n)] + pbar.update(1) pbar.set_description("{0}".format(self.pb_ag_text[4])) - - apg_d = [np.concatenate(elem) for elem in zip(pg_d, pXT_n, dXT_n)] apg_t = [np.concatenate(elem) for elem in zip(pg_t, pXT_t, dXT_t)] - Ap_d = [] Ap_t = [] for list_apg_d, list_apg_t in zip(apg_d, apg_t): _Apd_list = [] _Apt_list = [] - for d, t in zip(list_apg_d, list_apg_t): - - if(d>0): - + if (d > 0): _tAp = ns.dist_ftAp.rvs(size=int(d)) - _nAp = stats.binom.rvs( - n=np.ones(int(d)).astype(int), - p=(ns.pAp*(1 - np.exp(-_tAp/ns.taur))), - size=int(d) - ) - - cond_sel = (_nAp>0) - - if(sum(cond_sel)>0): + n=np.ones(int(d)).astype(int), + p=(ns.pAp*(1 - np.exp(-_tAp/ns.taur))), + size=int(d)) + cond_sel = (_nAp > 0) + if (sum(cond_sel) > 0): _tAp = _tAp[cond_sel] _nAp = _nAp[cond_sel] - - _Apd = (1-np.exp(-_tAp/ns.tslow))*_nAp _Apt = t + _tAp else: _Apd = [] _Apt = [] - else: _Apd = [] _Apt = [] - _Apd_list.extend(list(_Apd)) _Apt_list.extend(list(_Apt)) - Ap_d.append(_Apd_list) Ap_t.append(_Apt_list) - Ap_d = np.array(Ap_d, dtype="object") Ap_t = np.array(Ap_t, dtype="object") - - pbar.update(1) pbar.set_description("{0}".format(self.pb_ag_text[5])) - - ag_d = [ - np.concatenate( (list_pg_d, list_pXT_d, list_dXT_d, list_Ap_d)) - for list_pg_d, list_pXT_d, list_dXT_d, list_Ap_d, in zip(pg_d, pXT_d, dXT_d, Ap_d) - ] - - ag_t = [ - np.concatenate( (list_pg_t, list_pXT_t, list_dXT_t, list_Ap_t)) - for list_pg_t, list_pXT_t, list_dXT_t, list_Ap_t, in zip(pg_t, pXT_t, dXT_t, Ap_t) - ] - + + ag_d = [np.concatenate((list_pg_d, list_pXT_d, list_dXT_d, list_Ap_d)) + for list_pg_d, list_pXT_d, list_dXT_d, list_Ap_d, in zip(pg_d, pXT_d, dXT_d, Ap_d)] + + ag_t = [np.concatenate((list_pg_t, list_pXT_t, list_dXT_t, list_Ap_t)) + for list_pg_t, list_pXT_t, list_dXT_t, list_Ap_t, in zip(pg_t, pXT_t, dXT_t, Ap_t)] + ns.ag_d = np.asarray(ag_d, dtype="object") ns.ag_t = np.asarray(ag_t, dtype="object") - pbar.update(1) pbar.set_description("{0}".format(self.pb_ag_text[6])) - - + def QsG(self, ns, t): - tslow = ns.tslow tfast = ns.tfast rfast = ns.rfast tstart = ns.start_tgate tend = ns.start_tgate + ns.len_tgate - - - if(t<tstart): - QsG_slow = (1 - rfast)*(np.exp(-(tstart - t)/tslow) - np.exp(-(tend - t)/tslow)) - QsG_fast = (rfast)*(np.exp(-(tstart - t)/tfast) - np.exp(-(tend - t)/tfast)) + if (t < tstart): + QsG_slow = (1 - rfast)*(np.exp(-(tstart - t)/tslow) - np.exp(-(tend - t)/tslow)) + QsG_fast = (rfast)*(np.exp(-(tstart - t)/tfast) - np.exp(-(tend - t)/tfast)) QsG = QsG_slow + QsG_fast - - elif(t>=tstart and t<tend): + elif (t >= tstart and t < tend): QsG_slow = (1 - rfast)*np.exp(-(tend - t)/tslow) QsG_fast = (rfast)*np.exp(-(tend - t)/tfast) QsG = 1 - QsG_slow - QsG_fast - else: QsG = 0 - - return QsG - - + return QsG + def IsG(self, ns, t): - Islow= ((1 - ns.rfast)/ns.tslow)*(np.exp(-t/ns.tslow)) + Islow = ((1 - ns.rfast)/ns.tslow)*(np.exp(-t/ns.tslow)) Ifast = (ns.rfast/ns.tfast)*(np.exp(-t/ns.tfast)) IsG = Islow + Ifast return IsG - - + def QSiPM(self, ns): - pbar = tqdm(range(len(self.pb_Q_text)-1)) - pbar.set_description("{0}".format(self.pb_Q_text[0])) - Q = np.array(ns.dist_Q.rvs(ns.n_sim)) - n_ag_d = np.asarray([len(list_ag_d) for list_ag_d in ns.ag_d]) Sig1 = ns.dist_Sig1.rvs(np.sum(n_ag_d)) n_ag_d_cumsum = np.cumsum(n_ag_d)[:-1] ns.Sig1 = np.split(Sig1, n_ag_d_cumsum) - pbar.update(1) pbar.set_description("{0}".format(self.pb_Q_text[1])) - - Q_SiPM = np.array([ - np.sum( - np.asarray( - [ - elem_Sig1*elem_ag_d*self.QsG(ns, elem_ag_t) - if(elem_ag_d>0) - else 0 - for elem_Sig1, elem_ag_d, elem_ag_t in zip(list_Sig1, list_ag_d, list_ag_t) - ] - ) - ) - - for list_Sig1, list_ag_d, list_ag_t in zip(ns.Sig1, ns.ag_d, ns.ag_t) - - ]) - + Q_SiPM = np.array([np.sum(np.asarray([elem_Sig1*elem_ag_d*self.QsG(ns, elem_ag_t) + if (elem_ag_d > 0) + else 0 + for elem_Sig1, elem_ag_d, elem_ag_t in zip(list_Sig1, list_ag_d, list_ag_t)])) + for list_Sig1, list_ag_d, list_ag_t in zip(ns.Sig1, ns.ag_d, ns.ag_t)]) Q += Q_SiPM - - pbar.update(1) pbar.set_description("{0}".format(self.pb_Q_text[2])) pbar.close() - return Q - - def ISiPM(self, ns, ntim=1000, n_transients=200, convolve=True): - _range = [-ns.t0, ns.start_tgate+ns.len_tgate] - dtim = int(np.ceil((_range[1] - _range[0])/ntim-1)) tim = np.linspace(_range[0], _range[1], ntim) I = ns.dist_I.rvs(n_transients*len(tim)) I = np.array(np.split(I, len(tim))).T Gt = ns.dist_Gt.pdf(tim) - - n_ag_d = np.asarray([len(list_ag_d) for list_ag_d in ns.ag_d]) Sig1 = ns.dist_Sig1.rvs(np.sum(n_ag_d)) n_ag_d_cumsum = np.cumsum(n_ag_d)[:-1] Sig1 = np.split(Sig1, n_ag_d_cumsum) ns.Sig1 = Sig1 - - ag_d_temp = ns.ag_d[0:int(n_transients)] ag_t_temp = ns.ag_t[0:int(n_transients)] Sig1_temp = ns.Sig1[0:int(n_transients)] - - - I_SiPM = np.array([ - np.sum( - np.asarray( - [ - elem_Sig1*elem_ag_d*self.IsG(ns, elem_tim - elem_ag_t) - if(elem_ag_d>0 and elem_ag_t<elem_tim) - else 0 - for (elem_Sig1, elem_ag_d, elem_ag_t), elem_tim in product( - zip(list_Sig1, - list_ag_d, - list_ag_t), tim) - ] - ).reshape(len(list_ag_d), ntim), axis=0 - ) - - for list_Sig1, list_ag_d, list_ag_t in zip(tqdm(Sig1_temp), - ag_d_temp, - ag_t_temp) - - ]) - - I+=I_SiPM - + I_SiPM = np.array([np.sum(np.asarray([elem_Sig1*elem_ag_d * self.IsG(ns, elem_tim - elem_ag_t) + if (elem_ag_d > 0 and elem_ag_t < elem_tim) + else 0 + for (elem_Sig1, elem_ag_d, elem_ag_t), elem_tim in product(zip(list_Sig1,list_ag_d,list_ag_t), tim)]).reshape(len(list_ag_d), ntim), axis=0) + for list_Sig1, list_ag_d, list_ag_t in zip(tqdm(Sig1_temp), ag_d_temp, ag_t_temp)]) + I += I_SiPM I = np.apply_along_axis(signal.fftconvolve, 1, I, Gt, 'same') - + return I, tim - - - - def Simulate(self, - n_sim, - transients=False, - n_transients=None, - verbose=False - ): - - - + + def Simulate(self,n_sim,transients=False,n_transients=None,verbose=False): n_sim = int(n_sim) - for index, row in self.pars.iterrows(): - - dictionary = row.to_dict() - - if(verbose): + if (verbose): print("\n-------------------\n") print("\nSIMULATION {0}\n - {1}".format(index, row["Name_Sim"])) print("\n-------------------\n") for key, value in dictionary.items(): - if(key!="Name_Sim" and key not in self.out_columns): + if (key != "Name_Sim" and key not in self.out_columns): print("{0} \t\t {1} {2}".format(key, value, self.VariableDictionary[key]["Unit"])) print("\n-------------------\n") - - ns = SimpleNamespace(**dictionary) - + ns = SimpleNamespace(**dictionary) - if(ns.Gen_mu == "Poisson"): + if (ns.Gen_mu == "Poisson"): ns.dist_flight = scipy.stats.poisson(mu=ns.mu) - elif(ns.Gen_mu == "Fixed"): - ns.dist_flight = scipy.stats.uniform(loc = ns.mu, scale=0) - - if(ns.Gen_tmu == "Gauss"): + elif (ns.Gen_mu == "Fixed"): + ns.dist_flight = scipy.stats.uniform(loc=ns.mu, scale=0) + if (ns.Gen_tmu == "Gauss"): ns.dist_ftlight = scipy.stats.norm(loc=ns.tl0, scale=ns.tl1) - elif(ns.Gen_tmu == "Exp"): + elif (ns.Gen_tmu == "Exp"): ns.dist_ftlight = scipy.stats.expon(loc=ns.tl0, scale=ns.tl1) else: - raise Exception("The Time for Primary Geiger Discharge distribution mode that has been selected is invalid. Please choose from the following options: {0}".format( - ", ".join(self.VariableDictionary["Gen_tmu"]["options"]))) - + raise Exception("The Time for Primary Geiger Discharge distribution mode that has been selected is invalid. Please choose from the following options: {0}".format(", ".join(self.VariableDictionary["Gen_tmu"]["options"]))) + ns.dist_fDC = scipy.stats.poisson(mu=(ns.t0 + ns.start_tgate + ns.len_tgate)*ns.DCR) ns.dist_ftDC = scipy.stats.uniform(loc=-ns.t0, scale=ns.t0+ns.start_tgate+ns.len_tgate) - - if(ns.Gen_npXT == "Poisson"): - ns.dist_fpXT = scipy.stats.poisson(mu = ns.ppXT) - elif(ns.Gen_npXT == "Binomial"): - ns.dist_fpXT = scipy.stats.binom(n=1, p = ns.ppXT) - elif(ns.Gen_npXT == "Borel"): + + if (ns.Gen_npXT == "Poisson"): + ns.dist_fpXT = scipy.stats.poisson(mu=ns.ppXT) + elif (ns.Gen_npXT == "Binomial"): + ns.dist_fpXT = scipy.stats.binom(n=1, p=ns.ppXT) + elif (ns.Gen_npXT == "Borel"): ns.dist_fpXT = borel(mu=ns.ppXT) else: - raise Exception("The Prompt Cross-talk Discharge distribution mode that has been selected is invalid. Please choose from the following options: {0}".format( - ", ".join(self.VariableDictionary["Gen_npXT"]["options"]))) - - if(ns.Gen_ndXT == "Poisson"): - ns.dist_fdXT = scipy.stats.poisson(mu = ns.pdXT) - elif(ns.Gen_ndXT== "Binomial"): - ns.dist_fdXT = scipy.stats.binom(n=1, p = ns.pdXT) - elif(ns.Gen_ndXT == "Borel"): + raise Exception("The Prompt Cross-talk Discharge distribution mode that has been selected is invalid. Please choose from the following options: {0}".format(", ".join(self.VariableDictionary["Gen_npXT"]["options"]))) + + if (ns.Gen_ndXT == "Poisson"): + ns.dist_fdXT = scipy.stats.poisson(mu=ns.pdXT) + elif (ns.Gen_ndXT == "Binomial"): + ns.dist_fdXT = scipy.stats.binom(n=1, p=ns.pdXT) + elif (ns.Gen_ndXT == "Borel"): ns.dist_fdXT = borel(mu=ns.pdXT) else: - raise Exception("The Delayed Cross-talk Discharge mode that has been selected is invalid. Please choose from the following options: {0}".format( - ", ".join(self.VariableDictionary["Gen_ndXT"]["options"]))) + raise Exception("The Delayed Cross-talk Discharge mode that has been selected is invalid. Please choose from the following options: {0}".format(", ".join(self.VariableDictionary["Gen_ndXT"]["options"]))) - if(ns.Gen_tdXT == "Exp"): + if (ns.Gen_tdXT == "Exp"): ns.dist_ftdXT = scipy.stats.expon(scale=(ns.taudXT)) else: - raise Exception("The Time of Delayed Cross-talk Discharge mode that has been selected is invalid. Please choose from the following options: {0}".format( - ", ".join(self.VariableDictionary["Gen_tdXT"]["options"]))) - - + raise Exception("The Time of Delayed Cross-talk Discharge mode that has been selected is invalid. Please choose from the following options: {0}".format(", ".join(self.VariableDictionary["Gen_tdXT"]["options"]))) + ns.dist_ftAp = scipy.stats.expon(scale=(ns.tauAp)) ns.dist_Q = scipy.stats.norm(loc=0, scale=ns.Sig0) ns.dist_Sig1 = scipy.stats.norm(loc=1, scale=ns.Sig1) ns.dist_I = scipy.stats.norm(loc=0, scale=ns.Sigt) ns.dist_Gt = scipy.stats.norm(loc=0, scale=ns.GSig) ns.n_sim = n_sim - - - self.pars.at[index, "NSim"] = n_sim - + self.pars.at[index, "NSim"] = n_sim + start_time = time.time() self.AllGeiger(ns) end_time = time.time() - + start_time = time.time() Qs = self.QSiPM(ns) end_time = time.time() - - - self.pars.at[index, "GeigerArray"] = np.vstack([ns.ag_d, ns.ag_t]) -# self.pars.at[index, "TimeElapsed_GeigerArray"] = end_time-start_time + self.pars.at[index, "GeigerArray"] = np.vstack([ns.ag_d, ns.ag_t]) +# self.pars.at[index, "TimeElapsed_GeigerArray"] = end_time-start_time - self.pars.at[index, "ChargeSpectrum"] = Qs + self.pars.at[index, "ChargeSpectrum"] = Qs self.pars.at[index, "TimeElapsed_ChargeSpectrum"] = end_time-start_time - if(transients): + if (transients): start_time = time.time() I, t = self.ISiPM(ns, n_transients=n_transients) - #I, t = self.ISiPM(ns, ntim=6001, n_transients=n_transients) + # I, t = self.ISiPM(ns, ntim=6001, n_transients=n_transients) end_time = time.time() - - - Trans = { - "I":I, - "t":t - } - - self.pars.at[index, "Transients"] = Trans + Trans = {"I": I, "t": t} + self.pars.at[index, "Transients"] = Trans self.pars.at[index, "TimeElapsed_Transients"] = end_time-start_time - - def WriteDataFrame(self, name, directory, ext=None): - - if (self.pars.empty): raise Exception("The simulation data is empty. Please simulate an SiPM before writing a file.") - - if(ext is None): + if (ext is None): outdir = directory + name + ".pkl" self.pars.to_pickle(outdir) else: outdir = directory + name + ext - if(ext == ".pkl"): + if (ext == ".pkl"): self.pars.to_pickle(outdir) - elif(ext == ".hdf"): + elif (ext == ".hdf"): self.pars.to_hdf(outdir) - elif(ext == ".csv"): + elif (ext == ".csv"): self.pars.to_csv(outdir) - elif(ext == ".json"): + elif (ext == ".json"): self.pars.to_json(outdir) else: - raise Exception("The file format was not recognized. Please choose an output file extension .pkl, .hdf, .csv or .json " ) + raise Exception("The file format was not recognized. Please choose an output file extension .pkl, .hdf, .csv or .json ") def WriteASCII(self, name, directory): - if (self.pars.empty): raise Exception("The simulation data is empty. Please simulate an SiPM before writing a file.") - variables = self.pars.drop(['ChargeSpectrum'], axis=1) Qs = self.pars["ChargeSpectrum"] - with open(directory + name + ".head", 'a') as f: - f.write(variables.to_string(header = True, index = False)) - + f.write(variables.to_string(header=True, index=False)) with open(directory + name + ".dat", 'a') as f: for _elem in zip_longest(*Qs): formatted = (["{:.3f}".format(__elem) if __elem is not None else " " for __elem in _elem]) formatted = ','.join(formatted) f.write(formatted + '\n') -