Update LightSimtastic.py - Corrected Error in Afterpulse Calculation

LightSimtastic.py

@@ -339,6 +339,8 @@ class SiPMSimulation:
                     )
              )
         
+#         nAp = np.random.binomial(n=1, p=(ns.pAp*(1 - np.exp(-tAp/ns.taur))), size=nPG_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)
@@ -356,9 +358,6 @@ class SiPMSimulation:
         
         pXT_n = np.split(npXT, nPG_cumsum)[:-1]
         dXT_n = np.split(ndXT, nPG_cumsum)[:-1]
-        dXT_nt =ns.dist_ftdXT.rvs(size=sum(map(sum, dXT_n)))
-        dXT_nt = self.dXTtSplit(dXT_nt, dXT_n)  
- 
         Ap_n = np.split(nAp, nPG_cumsum)[:-1]
         Ap_nt = np.split(tAp, nPG_cumsum)[:-1]
         
@@ -421,21 +420,20 @@ class SiPMSimulation:
             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, elem_dXT_nt))
+                            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, elem_dXT_nt in zip(list_pg_d, list_pg_t, list_dXT_n, list_dXT_nt)
+                            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, list_dXT_nt in zip(pg_d, pg_t, dXT_n, dXT_nt)
+            for list_pg_d, list_pg_t, list_dXT_n in zip(pg_d, pg_t, dXT_n)
         ]
         
         
@@ -446,7 +444,7 @@ class SiPMSimulation:
             np.array(
                 list(chain.from_iterable(
                         [
-                            [elem_Ap_n]*elem_Ap_n
+                            [1-np.exp(-elem_Ap_t/ns.tslow)]*elem_Ap_n
                             if(elem_Ap_n > 0 and elem_pg_d >0)
                             else []
                             for elem_pg_d, elem_pg_t, elem_Ap_n, elem_Ap_t in zip(list_pg_d, list_pg_t, list_Ap_n, list_Ap_t)
@@ -493,7 +491,6 @@ class SiPMSimulation:
         
         pbar.update(1)
         pbar.set_description("{0}".format(self.pb_ag_text[6]))
-        pbar.close()
     
     
     def QsG(self, ns, t):
@@ -530,8 +527,6 @@ class SiPMSimulation:
     
     def QSiPM(self, ns):
         
-        
-        
         pbar = tqdm(range(len(self.pb_Q_text)-1))
         
         pbar.set_description("{0}".format(self.pb_Q_text[0]))
@@ -574,15 +569,13 @@ class SiPMSimulation:
                      
                      
 
-    def ISiPM(self, ns, ntim=2048, convolve=True, n_transients=20):
-        
-        
+    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(ns.n_sim*len(tim))
+        I = ns.dist_I.rvs(n_transients*len(tim))
         I = np.array(np.split(I, len(tim))).T
         Gt = ns.dist_Gt.pdf(tim)
         
@@ -594,13 +587,11 @@ class SiPMSimulation:
         ns.Sig1 = Sig1
 
         
-        I_temp = I[0:int(n_transients)]
         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(
@@ -622,16 +613,14 @@ class SiPMSimulation:
                     
                 ])
             
-        I_temp+=I_SiPM
+        I+=I_SiPM
         
-        I = np.apply_along_axis(signal.fftconvolve, 1, I_temp, Gt, 'same')
+        I = np.apply_along_axis(signal.fftconvolve, 1, I, Gt, 'same')
         
         return I, tim
     
    
 
-        
-
     def Simulate(self,
                  n_sim,
                  transients=False,
@@ -717,21 +706,29 @@ class SiPMSimulation:
             self.AllGeiger(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
-
             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, "ChargeSpectrum"] =  Qs
             self.pars.at[index, "TimeElapsed_ChargeSpectrum"] = end_time-start_time
 
             if(transients):
                 start_time = time.time()
-                I, t = self.ISiPM(ns)
+                I, t = self.ISiPM(ns, n_transients=n_transients)
                 end_time = time.time()
-                Trans = np.vstack([I,t])
+                
+            
+                Trans = {
+                    "I":I,
+                    "t":t
+                }
+                
                 self.pars.at[index, "Transients"] =  Trans
                 self.pars.at[index, "TimeElapsed_Transients"] = end_time-start_time
                
@@ -775,3 +772,4 @@ class SiPMSimulation:
                 formatted = (["{:.3f}".format(__elem) if __elem is not None else " " for __elem in _elem])
                 formatted = ','.join(formatted)
                 f.write(formatted + '\n')
+