From d52e2e5ec4c41678def5d38d37fbe1f0421dc767 Mon Sep 17 00:00:00 2001
From: pparygin <pavel.parygin@cern.ch>
Date: Tue, 19 Nov 2024 15:13:49 +0100
Subject: [PATCH] 1-cell developments incl. fixed AP, recursive AP, custom
 electronics noise, script parameters etc

---
 Example.py        | 105 ++++++++++-----
 GeNoise_upd.py    |  55 ++++++++
 LightSimtastic.py | 318 +++++++++++++++++++++++++++++++++-------------
 3 files changed, 360 insertions(+), 118 deletions(-)
 mode change 100644 => 100755 Example.py
 create mode 100644 GeNoise_upd.py

diff --git a/Example.py b/Example.py
old mode 100644
new mode 100755
index 69f0d7c..ce87627
--- a/Example.py
+++ b/Example.py
@@ -1,48 +1,77 @@
+#!/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()
+
diff --git a/GeNoise_upd.py b/GeNoise_upd.py
new file mode 100644
index 0000000..04a1ad2
--- /dev/null
+++ b/GeNoise_upd.py
@@ -0,0 +1,55 @@
+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
diff --git a/LightSimtastic.py b/LightSimtastic.py
index a4c72ca..04f38ad 100644
--- a/LightSimtastic.py
+++ b/LightSimtastic.py
@@ -1,4 +1,5 @@
 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,10 +235,106 @@ 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
+
+
+    
     def ReadMathCad(self, inputfile):
         variables = {}
         with codecs.open(inputfile, encoding='utf-8') as file:
@@ -313,33 +425,39 @@ 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)
         elif (len(nDC)>len(nlight)):
             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))
         
@@ -347,19 +465,21 @@ class SiPMSimulation:
         nDC_cumsum = np.cumsum(nDC)
         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]
         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]
+        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)]
     
-
+        #print(pg_d)
+        
         pbar.update(1)
         pbar.set_description("{0}".format(self.pb_ag_text[2]))
         
@@ -439,53 +559,12 @@ 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]))
@@ -494,19 +573,35 @@ class SiPMSimulation:
             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)  
         ]
+
+        #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)  
         ]
-        
+
         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):
                
         tslow = ns.tslow
@@ -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]
-       
-        dtim = int(np.ceil((_range[1] - _range[0])/ntim-1))
-        tim = np.linspace(_range[0], _range[1], ntim)
+        ntim = (ns.start_tgate+ns.len_tgate + ns.t0) *10 
+        
+        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
-
         
         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(
@@ -619,26 +730,42 @@ class SiPMSimulation:
                             ]
                         ).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)
-                    
-                ])
             
+                ])
+
+        '''
+        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
+        
+        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"])))
             
-                
-            ns.dist_ftAp = scipy.stats.expon(scale=(ns.tauAp))
+            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)
@@ -722,7 +856,7 @@ class SiPMSimulation:
             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
@@ -730,7 +864,8 @@ 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)
@@ -744,7 +879,12 @@ 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):
-- 
GitLab