Skip to content
Snippets Groups Projects
Commit c13ce820 authored by Antonello, Dr. Massimiliano's avatar Antonello, Dr. Massimiliano
Browse files

Formatting

parent 3064373b
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
AdditionalPDFs.py
+ 4
− 27
View file @ c13ce820
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,10 +32,8 @@ 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):
return ((mu > 0) & (mu<1))
......@@ -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)])
......
LightSimtastic.py
+ 179
− 518
View file @ c13ce820
import warnings
import time
import numpy as np
import pandas as pd
......@@ -14,215 +15,62 @@ 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": {"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_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)
"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:
......@@ -232,42 +80,31 @@ class SiPMSimulation:
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)
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()))))
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()))
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("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):
......@@ -283,11 +120,9 @@ class SiPMSimulation:
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())
......@@ -305,23 +140,16 @@ 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)):
nDC.resize(nlight.shape)
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]))
......@@ -336,8 +164,6 @@ class SiPMSimulation:
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)
......@@ -359,297 +185,160 @@ class SiPMSimulation:
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
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
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
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)))
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):
_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)
)
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_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_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))
QsG = QsG_slow + QsG_fast
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
def IsG(self, ns, t):
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)
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)
])
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)
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)
])
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):
print("\n-------------------\n")
print("\nSIMULATION {0}\n - {1}".format(index, row["Name_Sim"]))
......@@ -658,22 +347,18 @@ class SiPMSimulation:
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)
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"):
ns.dist_ftlight = scipy.stats.norm(loc=ns.tl0, scale=ns.tl1)
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)
......@@ -685,8 +370,7 @@ class SiPMSimulation:
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"])))
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)
......@@ -695,15 +379,12 @@ class SiPMSimulation:
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"):
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)
......@@ -711,8 +392,6 @@ class SiPMSimulation:
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
start_time = time.time()
......@@ -723,11 +402,9 @@ class SiPMSimulation:
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, "ChargeSpectrum"] = Qs
self.pars.at[index, "TimeElapsed_ChargeSpectrum"] = end_time-start_time
......@@ -736,24 +413,13 @@ class SiPMSimulation:
I, t = self.ISiPM(ns, n_transients=n_transients)
# I, t = self.ISiPM(ns, ntim=6001, n_transients=n_transients)
end_time = time.time()
Trans = {
"I":I,
"t":t
}
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):
outdir = directory + name + ".pkl"
self.pars.to_pickle(outdir)
......@@ -771,19 +437,14 @@ class SiPMSimulation:
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))
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')
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment