Skip to content
Snippets Groups Projects
Commit d52e2e5e authored by Pavel Parygin's avatar Pavel Parygin
Browse files

1-cell developments incl. fixed AP, recursive AP, custom electronics noise, script parameters etc

parent 2a069b32
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Example.py 100644 → 100755
+ 76
− 29
View file @ d52e2e5e
#!/usr/bin/python3
from LightSimtastic import SiPMSimulation
import matplotlib.pyplot as plt
from matplotlib.cm import plasma
import numpy as np
import sys
#################################
# THE SIMULATION TOOL TAKES A DICTIONARY AS AN INPUT, WITH EACH OF THE VARIABLES AS KEYS.
# FOR EACH ELEMENT YOU ADD TO THE LISTS, ANOTHER SIMULATION WILL BE PERFORMED.
# THIS ALLOWS SCANS THROUGH PARAMETERS.
#################################
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('--tauap', help='tau AP')
parser.add_argument('--tauap2', default=None, help='second tau AP')
parser.add_argument('--pap', help='PAP')
parser.add_argument('--nevents', help='N events')
parser.add_argument('--dcr', help='DCR rate (GHz)')
parser.add_argument('--noise', default=False, help='Generate electronics noise (True/False)')
parser.add_argument('--folder', help='Folder name to store results')
args = parser.parse_args()
tauAp = float(args.tauap)
tauAp2 = float(args.tauap2) if args.tauap2 else None
pap = float(args.pap)
nevn = int(args.nevents)
dcr = float(args.dcr)
noise = bool(args.noise)
folder = args.folder
variables={
"Name_Sim":[0], #THIS CAN BE ANYTHING YOU LIKE
"mu":[7], # MEAN NUMBER OF GEIGER DISCHARGES
"mu":[1], # MEAN NUMBER OF GEIGER DISCHARGES
"ppXT":[0.0], # PROBABILITY OF PROMPT CROSS-TALK
"pdXT":[0.2], #PROBABILITY OF DELAYED CROSS-TALK
"taudXT":[25], #TIME CONSTANT FOR DELAYED CROSS-TALK (NS)
"rdXT":[0.5], #FRACTION OF DELAYED CROSS-TALK FROM ADJACENT CELLS
"pAp":[0.15], #PROBABILITY OF AFTERPULSE
"tauAp":[7.5], #TIME CONSTANT OF AFTERPULSE
"taur":[20], # SIPM RECOVERY TIME CONSTANT
"Sig0":[0.075], #WIDTH OF PEDESTAL [NORMALISED ADC]
"Sig1":[0.02], # INCREASE IN PEAK WIDTH PER ADDITIONAL DISCHARGE [NORMALISED ADC]
"DCR":[0.0], # DARK COUNT RATE [GHZ]
"Sigt":[0.02], # ELECTRONICS NOISE FOR CURRENT TRANSIENT (FOR TRANSIENT)
"GSig":[0.75], # WIDTH OF ELECTRONICS NOISE TRANSFER FUNCTION (FOR TRANSIENT)
"tslow":[20], # TIME CONSTANT FOR SLOW PART OF PULSE
"tfast":[1.5], # TIME CONSTANT FOR FAST PART OF PULSE
"rfast":[0.2], # PROPORTION OF SLOW/FAST
"start_tgate":[-5], # GATE START TIME [NS]
"len_tgate":[100], # LENGTH OF GATE
"t0":[100,], # STARTING POINT OF GATE
"tl0":[0], #GAUSSIAN MEAN OF PRIMARY GEIGER DICHARGE TIME DISTRIBUTION
"tl1":[0.1], #GAUSSIAN WIDTH OF PRIMARY GEIGER DICHARGE TIME DISTRIBUTION
"pdXT":[0.0], #PROBABILITY OF DELAYED CROSS-TALK
"taudXT":[0], #TIME CONSTANT FOR DELAYED CROSS-TALK (NS)
"rdXT":[0.0], #FRACTION OF DELAYED CROSS-TALK FROM ADJACENT CELLS
"pAp":[pap], #PROBABILITY OF AFTERPULSE
"tauAp":[tauAp], #TIME CONSTANT OF AFTERPULSE
"tauAp2" : [tauAp2],
"taur":[15.], # SIPM RECOVERY TIME CONSTANT
"Sig0":[0.9], #WIDTH OF PEDESTAL [NORMALISED ADC]
"Sig1":[0.], # INCREASE IN PEAK WIDTH PER ADDITIONAL DISCHARGE [NORMALISED ADC]
"DCR":[dcr], # DARK COUNT RATE [GHZ]
"Sigt":[1e-3], # ELECTRONICS NOISE FOR CURRENT TRANSIENT (FOR TRANSIENT)
"GSig":[1.], # WIDTH OF ELECTRONICS NOISE TRANSFER FUNCTION (FOR TRANSIENT)
"tslow":[15.], # TIME CONSTANT FOR SLOW PART OF PULSE
"tfast":[0.01], # TIME CONSTANT FOR FAST PART OF PULSE
"rfast":[.0001], # PROPORTION OF SLOW/FAST
"start_tgate":[-250], # GATE START TIME [NS]
"len_tgate":[500], # LENGTH OF GATE
"t0":[250], # STARTING POINT OF GATE
"tl0":[1], #GAUSSIAN MEAN OF PRIMARY GEIGER DICHARGE TIME DISTRIBUTION
"tl1":[0.], #GAUSSIAN WIDTH OF PRIMARY GEIGER DICHARGE TIME DISTRIBUTION
"tl2":[0], #FREE PARAMETER
"Gen_mu":["Poisson"], # NUMBER PRIMARY GEIGER DISCHARGE PDF
"Gen_mu":["Fixed"], # NUMBER PRIMARY GEIGER DISCHARGE PDF
"Gen_tmu":["Gauss"], # TIME OF PRIMARY GEIGER DISCHARGE PDF
"Gen_gain":["Gauss"], # GAIN PDF (FOR TRANSIENT)
"Gen_npXT":["Binomial"], #NUMBER PROMPT X-TALK DISCHARGE DISTRIBUTION
"Gen_ndXT":["Binomial"], #NUMBER PROMPT X-TALK DISCHARGE DISTRIBUTION
"Gen_tdXT":["Exp"], #TIME DELAYED X-TALK DISTRIBUTION
"Gen_tAP":["Exp"], #AFTERPULSE TIME DISTRIBUTION} #fixed
"gennoise":[noise],
"folder": [args.folder]
}
'''
"Gen_nAP":["Binom"], #NUMBER AFTER DISTRIBUTION
"Gen_tAP":["Exp"], #AFTERPULSE TIME DISTRIBUTION
"Gen_noise":["Gauss"] #ELECTRONIC NOISE DISTRIBUTION (FOR TRANSIENT)
}
'''
#################################
# CREATE A SIPM SIMULATION CLASS OBJECT
#################################
......@@ -55,14 +84,14 @@ s = SiPMSimulation()
s.AddVariables(variables)
n_events = int(1e5)
n_events = nevn
#################################
# SIMULATE
#################################
s.Simulate(n_events, transients=False)
n_transients = n_events
s.Simulate(n_events, transients=True, n_transients=n_transients, dumpTransients=True)
#################################
# EXTRACT SIMULATION DATAFRAME
......@@ -80,14 +109,32 @@ Qs = df.iloc[0]["ChargeSpectrum"]
# PLOT
#################################
'''
plt.figure(figsize=(15.5,10.5))
H, edges = np.histogram(Qs, bins=1000)
edges = edges[:-1]+(edges[1]-edges[0])
plt.plot(edges, H)
plt.title("Simulated Charge Spectrum")
plt.yscale("log")
#plt.yscale("log")
plt.xticks(fontsize=25)
plt.yticks(fontsize=25)
plt.xlabel("# GD", fontsize=25),
plt.savefig("./Example.png")
#plt.show()
plt.savefig("/eos/user/p/pparygin/www/jsroot7/charge.png")
'''
cmap = plasma
if n_transients >100: n_transients = 100
plt.figure(figsize=(15.5,10.5))
for n in range(n_transients):
color = cmap(10+ n % n_transients)
#plt.plot(s.pars.at[0,"Transients"]['t'][4000:6000],s.pars.at[0,"Transients"]['I'][n][4000:6000], color = color, alpha=0.2)
plt.plot(s.pars.at[0,"Transients"]['t'],s.pars.at[0,"Transients"]['I'][n], color = color, alpha=0.2)
plt.title(f"{n_transients} simulated waveform(s)")
plt.xticks(fontsize=25)
plt.yticks(fontsize=25)
plt.xlabel("Time,ns", fontsize=25),
#plt.savefig("/eos/user/p/pparygin/www/jsroot7/transients_tauAp_{}.png".format(variables["tauAp"][0]))
plt.show()
GeNoise_upd.py 0 → 100644
+ 55
− 0
View file @ d52e2e5e
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.ar_model import AutoReg
from joblib import Parallel, delayed
def generate_ar_noise(model_params, initial_noise, length, noise_std):
noise = np.zeros(length)
noise[:len(initial_noise)] = initial_noise
# Vectorized computation for each step after initial noise
for i in range(len(initial_noise), length):
noise[i] = model_params[0] + np.dot(model_params[1:], noise[i - len(model_params[1:]):i][::-1])
noise[i] += np.random.normal(scale=noise_std) # Adding random noise
return noise
def get_random_initial_noise(original_noise, lags):
start_idx = np.random.randint(0, len(original_noise) - lags)
return original_noise[start_idx:start_idx + lags]
def noise(noise_sample_length=1000, num_transients=10, plot=False):
# Load waveform and prepare noise segments
waveform = np.genfromtxt('waveform_data.txt', delimiter='\n') + 6.5
noise_segments_indices = [(0, 3000), (4000, 6000)]
noise_segments = [waveform[start:end] for start, end in noise_segments_indices]
all_noise = np.concatenate(noise_segments)
noise_std = np.std(all_noise) / 2
lags = np.random.randint(1, 30, 1)[0]
# Fit AR model
ar_model = AutoReg(all_noise, lags=lags).fit()
model_params = ar_model.params
# Generate noise in parallel for each transient
noises = Parallel(n_jobs=-1)(
delayed(generate_ar_noise)(
model_params,
get_random_initial_noise(all_noise, lags),
noise_sample_length,
noise_std
) for _ in range(num_transients)
)
# Convert list to array
noises = np.array(noises)
if plot:
plt.figure(figsize=(10, 4))
plt.plot(noises[0]) # Plotting only the first transient
plt.title('Generated Noise Using AR Model')
plt.savefig("./noise.png")
return noises
LightSimtastic.py
+ 229
− 89
View file @ d52e2e5e
import time
import math
import numpy as np
import pandas as pd
import scipy
......@@ -10,6 +11,10 @@ from operator import add
from types import SimpleNamespace
from AdditionalPDFs import borel
import codecs
import GeNoise_upd
import h5py
import matplotlib.pyplot as plt
np.set_printoptions(suppress=True)
......@@ -22,12 +27,9 @@ else:
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"',
'"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"' '"gennoise"',
......@@ -65,6 +67,10 @@ class SiPMSimulation:
"Description": " Time delay constant of Afterpulsing in SiPM Cell",
"Unit": None
},
"tauAp2":{
"Description": " Second time delay constant of Afterpulsing in SiPM Cell",
"Unit": None
},
"taur":{
"Description": " PDE Recovery constant in SiPM Cell",
"Unit": None
......@@ -176,6 +182,16 @@ class SiPMSimulation:
"Unit": None,
"Options": ["Exp"]
},
"gennoise":{
"Description" : "Generate electronics noise from real transient",
"Unit" : None,
"Options": [True, False]
},
"folder":{
"Description":"Directory name for output",
"Unit" : None
}
}
......@@ -219,7 +235,103 @@ class SiPMSimulation:
self.pars["GeigerArray"] = self.pars["GeigerArray"].astype(object)
self.pars["ChargeSpectrum"] = self.pars["ChargeSpectrum"].astype(object)
self.pars["Transients"] = self.pars["Transients"].astype(object)
#self.napcounts = None
self.napcounts = np.array([0])
self.ndcrcounts = None
def generate_afterpulses(self, primary_times, primary_amplitudes, ns, max_time, recursive=False):
all_Ap_d = []
all_Ap_t = []
all_count_ap = [] # To store count_ap for each waveform
# Process each primary pulse (each waveform)
for list_apg_d, list_apg_t in zip(primary_amplitudes, primary_times):
Ap_d = []
Ap_t = []
count_ap = np.array([0]) # Initialize count_ap for this waveform
def recursive_afterpulse(t, d):
nonlocal count_ap # This makes count_ap accessible inside this function
if t > max_time or d <= 0:
return [], []
_Apd_list = []
_Apt_list = []
# Generate afterpulse times
#_tAp = ns.dist_ftAp.rvs(size=int(d))
n_afterpulses = math.ceil(d)
p1 = 0.7
p2 = 1 - p1
#_tAp = ns.dist_ftAp.rvs(size=math.ceil(d))
tAp1 = ns.dist_ftAp.rvs(size=math.ceil(d))
tAp2 = ns.dist_ftAp2.rvs(size=math.ceil(d)) if ns.dist_ftAp2 else 0
use_tau1 = np.random.rand(n_afterpulses) < p1
_tAp = np.where(use_tau1, tAp1, tAp2) if ns.dist_ftAp2 else ns.dist_ftAp.rvs(size=int(d))
#for i, is_tau1 in enumerate(use_tau1):
# tau_used = "tau1" if is_tau1 else "tau2"
# print(f"Afterpulse {i}: Using {tau_used} with _tAp = {_tAp[i]}")
# Calculate the number of afterpulses based on the probability
#if ns.Gen_tAP == "Exp":
# print('probability: ', ns.pAp * (1 - np.exp(-_tAp / ns.taur)))
_nAp = stats.binom.rvs(
n=np.ones(math.ceil(d)).astype(int),
p=(ns.pAp * (1 - np.exp(-_tAp / ns.taur))) if ns.Gen_tAP == "Exp" else ns.pAp,
size=math.ceil(d)
)
# Select only those that generated afterpulses
cond_sel = (_nAp > 0)
if sum(cond_sel) > 0:
_tAp = _tAp[cond_sel]
_nAp = _nAp[cond_sel]
count_ap = np.append(count_ap, _nAp) # Appending afterpulse counts
# Calculate afterpulse amplitudes
_Apd = (1 - np.exp(-_tAp / ns.tslow)) * _nAp
_Apt = t + _tAp
# Add to list
_Apd_list.extend(list(_Apd))
_Apt_list.extend(list(_Apt))
# Recursively generate afterpulses for each afterpulse
if recursive:
for new_t, new_d in zip(_Apt, _Apd):
if new_t <= max_time: # Ensure it stays within the time window
next_Apd, next_Apt = recursive_afterpulse(new_t, new_d)
_Apd_list.extend(next_Apd)
_Apt_list.extend(next_Apt)
#print(f'recurse \t {new_t}, {new_d}')
return _Apd_list, _Apt_list
for d, t in zip(list_apg_d, list_apg_t):
if d > 0:
Apd, Apt = recursive_afterpulse(t, d)
Ap_d.extend(Apd)
Ap_t.extend(Apt)
all_Ap_d.append(Ap_d)
all_Ap_t.append(Ap_t)
all_count_ap.append(count_ap) # Store count_ap for this waveform
# Convert to numpy arrays
all_Ap_d = np.array(all_Ap_d, dtype="object")
all_Ap_t = np.array(all_Ap_t, dtype="object")
return all_Ap_d, all_Ap_t, all_count_ap
......@@ -313,7 +425,12 @@ class SiPMSimulation:
pbar.set_description("{0}".format(self.pb_ag_text[0]))
nlight = np.asarray(ns.dist_flight.rvs(size=ns.n_sim))
if (ns.Gen_mu == "Fixed"): nlight = np.asarray(ns.dist_flight.rvs(size=ns.n_sim), dtype=int)
else: nlight = np.asarray(ns.dist_flight.rvs(size=ns.n_sim))
nDC = np.asarray(ns.dist_fDC.rvs(size=ns.n_sim))
#nDC[nDC>1] = 1
self.ndcrcounts = np.sum(nDC)
#print(nDC)
if(len(nlight)>len(nDC)):
nDC.resize(nlight.shape)
......@@ -321,25 +438,26 @@ class SiPMSimulation:
nlight.resize(nDC.shape)
nPG = nDC + nlight
#nPG[nPG>1] = 1
pbar.update(1)
pbar.set_description("{0}".format(self.pb_ag_text[1]))
nlight_sum = np.sum(nlight)
nlight_sum = int(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 = np.array(ns.dist_fpXT.rvs(size=int(nPG_sum)))
ndXT = np.array(ns.dist_fdXT.rvs(size=int(nPG_sum)))
npXT_sum = np.sum(npXT)
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_l_t = np.array([1,20,1,20,1,20,1,20,1,20])
pg_dc_d = np.ones(nDC_sum)
pg_dc_t = np.asarray(ns.dist_ftDC.rvs(size=nDC_sum))
......@@ -348,6 +466,7 @@ class SiPMSimulation:
nPG_cumsum = np.cumsum(nPG)
ndXT_cumsum = np.cumsum(ndXT)
#print("\nDEBUG ", pg_l_d, nlight, nlight_cumsum, "\n")
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]
......@@ -359,6 +478,7 @@ 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)]
#print(pg_d)
pbar.update(1)
pbar.set_description("{0}".format(self.pb_ag_text[2]))
......@@ -440,52 +560,11 @@ class SiPMSimulation:
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)
)
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")
Ap_d, Ap_t, _napcounts = self.generate_afterpulses(apg_t, apg_d, ns, 250.)
self.napcounts = np.concatenate(_napcounts, axis=0).sum()
#print(self.napcounts)
#print(Ap_d)
pbar.update(1)
pbar.set_description("{0}".format(self.pb_ag_text[5]))
......@@ -495,6 +574,8 @@ class SiPMSimulation:
for list_pg_d, list_pXT_d, list_dXT_d, list_Ap_d, in zip(pg_d, pXT_d, dXT_d, Ap_d)
]
#print('Light/DC and AP:\n', pg_t, '\n', 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)
......@@ -502,10 +583,24 @@ class SiPMSimulation:
ns.ag_d = np.asarray(ag_d, dtype="object")
ns.ag_t = np.asarray(ag_t, dtype="object")
#print(ns.ag_d)
pbar.update(1)
pbar.set_description("{0}".format(self.pb_ag_text[6]))
flattened_data = []
ii=0
for t_row, a_row in zip(pg_t, pg_d):
ii+=1
for t, a in zip(t_row, a_row):
flattened_data.append([t, a, ii, 'light' if t == 1.0 else 'DC'])
ii=0
for t_row, a_row in zip(Ap_t, Ap_d):
ii+=1
for t, a in zip(t_row, a_row):
flattened_data.append([t, a, ii, 'AP'])
_f = f'./rawSim_r{ns.DCR}_dcr{self.ndcrcounts}_Nap_{self.napcounts}_APdT{ns.tauAp}.csv'
np.savetxt(_f,np.array(flattened_data), delimiter='\t', fmt='%s', comments='')
print('Raw simulations info stored at:', _f)
def QsG(self, ns, t):
......@@ -533,10 +628,13 @@ class SiPMSimulation:
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
t0 = 1.
A = 7
Irise = (1-np.exp(-t/0.1))/30. if t<=t0 else 0
Islow = ((1 - ns.rfast)/ns.tslow)*(np.exp(-(t)/ns.tslow)) if t>t0 else 0
Ifast = (ns.rfast/ns.tfast)*(np.exp(-(t)/ns.tfast)) if t>t0 else 0
IsG = Irise + Islow + Ifast
return IsG*A
def QSiPM(self, ns):
......@@ -581,13 +679,12 @@ class SiPMSimulation:
return Q
def ISiPM(self, ns, ntim=1000, n_transients=200, convolve=True):
_range = [-ns.t0, ns.start_tgate+ns.len_tgate]
ntim = (ns.start_tgate+ns.len_tgate + ns.t0) *10
dtim = int(np.ceil((_range[1] - _range[0])/ntim-1))
tim = np.linspace(_range[0], _range[1], ntim)
tim = np.linspace(_range[0], _range[1], ntim, endpoint=False)
I = ns.dist_I.rvs(n_transients*len(tim))
I = np.array(np.split(I, len(tim))).T
Gt = ns.dist_Gt.pdf(tim)
......@@ -595,16 +692,30 @@ class SiPMSimulation:
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))
#Sig1 = np.ones(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
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)]
#print('\n',ag_d_temp,'\n',ag_t_temp)
for k in range(len(ag_d_temp)):
sorted_indices = np.argsort(ag_t_temp[k])
ag_t_temp[k] = np.array(ag_t_temp[k])[sorted_indices]
ag_d_temp[k] = np.array(ag_d_temp[k])[sorted_indices]
for k in range(len(ag_d_temp)):
prev_amplitude = 0
for i in range(1, len(ag_d_temp[k])):
deltaT = ag_t_temp[k][i] - ag_t_temp[k][i - 1]
if ag_d_temp[k][i] > (1 - np.exp(-deltaT / ns.tslow)):
ag_d_temp[k][i] = (1 - np.exp(-deltaT / ns.tslow))
#print(ag_d_temp, '\n')
I_SiPM = np.array([
np.sum(
np.asarray(
......@@ -626,19 +737,35 @@ class SiPMSimulation:
])
'''
flattened_data = []
ii=0
for t_row, a_row in zip(ag_t_temp, ag_d_temp):
ii+=1
for t, a in zip(t_row, a_row):
flattened_data.append([t, a, ii])
_f = f'/eos/cms/store/user/pparygin/LightSim/noNoise/rawSim_r{ns.DCR}_dcr{self.ndcrcounts}_Nap_{self.napcounts}_APdT{ns.tauAp}.csv'
np.savetxt(_f,np.array(flattened_data), delimiter='\t', fmt='%s', comments='')
'''
print("Transients...")
noise = GeNoise_upd.noise(len(I_SiPM[0]), len(I_SiPM), plot=True) if ns.gennoise else 0
I+=I_SiPM
I = np.apply_along_axis(signal.fftconvolve, 1, I, Gt, 'same')
I = I+noise
return I, tim
print("Transients are done")
return I, tim
def Simulate(self,
n_sim,
transients=False,
n_transients=None,
verbose=False
verbose=False,
dumpTransients=False
):
......@@ -704,8 +831,15 @@ class SiPMSimulation:
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"])))
if (ns.Gen_tAP == 'Exp'):
ns.dist_ftAp = scipy.stats.expon(scale=(ns.tauAp))
ns.dist_ftAp2 = scipy.stats.expon(scale=(ns.tauAp2)) if ns.tauAp2 else None
elif(ns.Gen_tAP == 'Fixed'):
ns.dist_ftAp = scipy.stats.uniform(ns.tauAp, scale=0)
else:
raise Exception("The Time of Afterpulsing Discharge mode that has been selected is invalid. Please choose from the following options: {0}".format(
", ".join(self.VariableDictionary["Gen_tAP"]["options"])))
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)
......@@ -731,6 +865,7 @@ class SiPMSimulation:
self.pars.at[index, "ChargeSpectrum"] = Qs
self.pars.at[index, "TimeElapsed_ChargeSpectrum"] = end_time-start_time
if(transients):
start_time = time.time()
I, t = self.ISiPM(ns, n_transients=n_transients)
......@@ -745,6 +880,11 @@ class SiPMSimulation:
self.pars.at[index, "Transients"] = Trans
self.pars.at[index, "TimeElapsed_Transients"] = end_time-start_time
if dumpTransients == True:
h5f = h5py.File(f'./Sim_r{ns.DCR}_dcr{self.ndcrcounts}_Nap_{self.napcounts}_APdT{ns.tauAp}_AP2dT{ns.tauAp2}.h5', 'w')
h5f.create_dataset('simulation', data=np.vstack((t.reshape(1,len(t)),I)))
print("File with transients: {}".format(h5f.filename))
h5f.close()
def WriteDataFrame(self, name, directory, ext=None):
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment