Double PMT ap probability for photons emitting double PE.

wfsim/core/afterpulse.py

@@ -187,9 +187,20 @@ def photon_afterpulse(signal_pulse, resource, config):
             # Assign each photon FRIST random uniform number rU0 from (0, 1] for timing
             rU0 = 1 - np.random.rand(len(signal_pulse._photon_timings))
 
+            # delaytime_cdf is intentionally not normalized to 1 but the probability of the AP 
+            prob_ap = delaytime_cdf[signal_pulse._photon_channels, -1]
+            if prob_ap.max() * config['pmt_ap_modifier'] > 0.5:
+                prob = prob_ap.max() * config['pmt_ap_modifier']
+                log.warning(f'PMT after pulse probability is {prob} larger than 0.5?')
+
+            # Scaling down (up) rU0 effectivly increase (decrease) the ap rate
+            rU0 /= config['pmt_ap_modifier']
+
+            # Double the probability for those photon emitting dpe
+            rU0[signal_pulse._photon_is_dpe] /= 2
+
             # Select those photons with U <= max of cdf of specific channel
-            cdf_max = delaytime_cdf[signal_pulse._photon_channels, -1]
-            sel_photon_id = np.where(rU0 <= cdf_max * config['pmt_ap_modifier'])[0]
+            sel_photon_id = np.where(rU0 <= prob_ap)[0]
             if len(sel_photon_id) == 0:
                 continue
             sel_photon_channel = signal_pulse._photon_channels[sel_photon_id]

wfsim/core/pulse.py

@@ -63,11 +63,19 @@ def __call__(self, *args):
         self._raw_area = self._raw_area_bottom = 0
         self._raw_area_trigger = self._raw_area_trigger_bottom = 0
 
+        # Assign DPE, and save it for compute after-pulses
+        p_dpe = self.config['p_double_pe_emision']
+        self._photon_is_dpe = np.random.binomial(n=1,
+                                                p=p_dpe,
+                                                size=len(self._photon_timings)).astype(np.bool_)
+
         counts_start = 0  # Secondary loop index for assigning channel
         for channel, counts in zip(*np.unique(self._photon_channels, return_counts=True)):
 
             # Use 'counts' amount of photon for this channel
             _channel_photon_timings = self._photon_timings[counts_start:counts_start+counts]
+            _channel_photon_is_dpe = self._photon_is_dpe[counts_start:counts_start+counts]
+
             counts_start += counts
             if channel in self.config['turned_off_pmts']:
                 continue
@@ -81,10 +89,8 @@ def __call__(self, *args):
                     * self.uniform_to_pe_arr(np.random.random(len(_channel_photon_timings)), channel)
 
                 # Add some double photoelectron emission by adding another sampled gain
-                n_double_pe = np.random.binomial(len(_channel_photon_timings),
-                                                 p=self.config['p_double_pe_emision'])
-
-                _channel_photon_gains[:n_double_pe] += self.config['gains'][channel] \
+                n_double_pe = _channel_photon_is_dpe.sum()
+                _channel_photon_gains[_channel_photon_is_dpe] += self.config['gains'][channel] \
                     * self.uniform_to_pe_arr(np.random.random(n_double_pe), channel)
 
             else: