Double PMT ap probability for photons emitting double PE.

wfsim/core/afterpulse.py

@@ -189,6 +189,13 @@ def photon_afterpulse(signal_pulse, resource, config):
 
             # Select those photons with U <= max of cdf of specific channel
             cdf_max = delaytime_cdf[signal_pulse._photon_channels, -1]
+            if cdf_max.max() * config['pmt_ap_modifier'] > 0.5:
+                prob = cdf_max.max() * config['pmt_ap_modifier']
+                log.warning(f'PMT after pulse probability is {prob} larger than 0.5?')
+
+            # Double the probability for those photon emitting dpe
+            cdf_max[signal_pulse._photon_emit_dpe] *= 2
+
             sel_photon_id = np.where(rU0 <= cdf_max * config['pmt_ap_modifier'])[0]
             if len(sel_photon_id) == 0:
                 continue

wfsim/core/pulse.py

@@ -62,6 +62,7 @@ def __call__(self, *args):
         self._n_pe_trigger = self._n_pe_trigger_bottom = 0
         self._raw_area = self._raw_area_bottom = 0
         self._raw_area_trigger = self._raw_area_trigger_bottom = 0
+        self._photon_emit_dpe = np.zeros(len(self._photon_timings), dtype=np.bool_)  # For PMT afterpulse
 
         counts_start = 0  # Secondary loop index for assigning channel
         for channel, counts in zip(*np.unique(self._photon_channels, return_counts=True)):
@@ -87,6 +88,8 @@ def __call__(self, *args):
                 _channel_photon_gains[:n_double_pe] += self.config['gains'][channel] \
                     * self.uniform_to_pe_arr(np.random.random(n_double_pe), channel)
 
+                self._photon_emit_dpe[counts_start: counts_start+n_double_pe] = True
+
             else:
                 n_double_pe = 0
                 _channel_photon_gains = np.array(self._photon_gains[self._photon_channels == channel])