From e80adb070a8c68fb2b38fb9e39a52fdefa73d07c Mon Sep 17 00:00:00 2001 From: tianyu zhu Date: Fri, 10 Dec 2021 05:22:21 -0500 Subject: [PATCH] SPE spec bug fix, adding trigger truth (#283) * bug fix * new truth fields * delete * Update wfsim/core/pulse.py Co-authored-by: Joran Angevaare * add more * sneaky bug fix * change in other files * add comments from Tianyu Co-authored-by: Joran Angevaare Co-authored-by: Peter Gaemers --- wfsim/core/pulse.py | 83 +++++++++++++++++++++++++++++++++------- wfsim/core/rawdata.py | 17 ++------ wfsim/strax_interface.py | 20 +++++++--- 3 files changed, 88 insertions(+), 32 deletions(-) diff --git a/wfsim/core/pulse.py b/wfsim/core/pulse.py index 09ecd236..df1a37ec 100644 --- a/wfsim/core/pulse.py +++ b/wfsim/core/pulse.py @@ -29,7 +29,11 @@ def __init__(self, config): self.init_pmt_current_templates() self.init_spe_scaling_factor_distributions() self.config['turned_off_pmts'] = np.arange(len(config['gains']))[np.array(config['gains']) == 0] - + self.current_max = np.max(np.array(self._pmt_current_templates), axis=1) + self.current_2_adc = self.config['pmt_circuit_load_resistor'] \ + * self.config['external_amplification'] \ + / (self.config['digitizer_voltage_range'] / 2 ** (self.config['digitizer_bits'])) + self.clear_pulse_cache() def __call__(self, *args): @@ -52,7 +56,12 @@ def __call__(self, *args): len(self._photon_timings)).astype(np.int64) dt = self.config.get('sample_duration', 10) # Getting dt from the lib just once - self._n_double_pe = self._n_double_pe_bot = 0 # For truth aft output + self._n_photon = self._n_photon_bottom = 0 # For truth output + self._n_pe = self._n_pe_bottom = 0 + self._n_photon_trigger = self._n_photon_trigger_bottom = 0 + 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 counts_start = 0 # Secondary loop index for assigning channel for channel, counts in zip(*np.unique(self._photon_channels, return_counts=True)): @@ -69,24 +78,27 @@ def __call__(self, *args): if '_photon_gains' not in self.__dict__: _channel_photon_gains = self.config['gains'][channel] \ - * self.uniform_to_pe_arr(np.random.random(len(_channel_photon_timings))) + * 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']) - self._n_double_pe += n_double_pe - if channel in self.config['channels_bottom']: - self._n_double_pe_bot += n_double_pe - - if self.config['detector'] == 'XENON1T': - _channel_photon_gains[:n_double_pe] += self.config['gains'][channel] \ - * self.uniform_to_pe_arr(np.random.random(n_double_pe), channel) - else: - _channel_photon_gains[:n_double_pe] += self.config['gains'][channel] \ - * self.uniform_to_pe_arr(np.random.random(n_double_pe)) + + _channel_photon_gains[:n_double_pe] += self.config['gains'][channel] \ + * self.uniform_to_pe_arr(np.random.random(n_double_pe), channel) + else: + n_double_pe = 0 _channel_photon_gains = np.array(self._photon_gains[self._photon_channels == channel]) + # += truth per channel to internal truth fields + self.add_truth( + _channel_photon_timings, + _channel_photon_gains, + dt, + channel, + n_double_pe,) + # Build a simulated waveform, length depends on min and max of photon timings min_timing, max_timing = np.min( _channel_photon_timings), np.max(_channel_photon_timings) @@ -199,6 +211,51 @@ def uniform_to_pe_arr(self, p, channel=0): indices = (p * 2000).astype(np.int64) + 1 return self.__uniform_to_pe_arr[channel, indices] + def add_truth(self, + photon_timings, + photon_gains, + dt, + channel, + n_double_pe, + ): + """Add required information to the fields used for the truth information""" + if channel in self.config['turned_off_pmts']: + return + + if str(channel) in self.config.get('special_thresholds', {}): + threshold = self.config['special_thresholds'][str(channel)] - 0.5 + else: + threshold = self.config['zle_threshold'] - 0.5 + # Figure out if we were above threshold + # - Current max is the highest value in the SPE pulse model given a + # remainder [0-9] ns (offset from 10 ns sample time). + # The SPE pulse model sums up to 0.1 (pe/ns), and the peak is ~0.03 (pe/ns) + # - Multiply this by the sampled gain of each photon, we get (electron/ns) + # - The current_2_adc take into account n_electron -> current -> voltage + # -> amplification -> digitization, converting the electron/ns to adc value. + remainder = (photon_timings % dt).astype(int) + max_amplitude_adc = photon_gains * self.current_max[remainder] * self.current_2_adc + above_threshold = max_amplitude_adc > threshold + trigger_photon = np.sum(above_threshold) + trigger_dpe = np.sum(above_threshold[:n_double_pe]) + raw_area = np.sum(photon_gains) / self.config['gains'][channel] + trigger_raw_area = np.sum(photon_gains[above_threshold]) / self.config['gains'][channel] + + self._n_photon += len(photon_timings) + self._n_pe += len(photon_timings) + n_double_pe + self._n_photon_trigger += trigger_photon + self._n_pe_trigger += trigger_photon + trigger_dpe + self._raw_area += raw_area + self._raw_area_trigger += trigger_raw_area + + if channel in self.config['channels_bottom']: + self._n_photon_bottom += len(photon_timings) + self._n_pe_bottom += len(photon_timings) + n_double_pe + self._n_photon_trigger_bottom += trigger_photon + self._n_pe_trigger_bottom += trigger_photon + trigger_dpe + self._raw_area_bottom += raw_area + self._raw_area_trigger_bottom += trigger_raw_area + def clear_pulse_cache(self): self._pulses = [] diff --git a/wfsim/core/rawdata.py b/wfsim/core/rawdata.py index 178e207a..9a05cace 100644 --- a/wfsim/core/rawdata.py +++ b/wfsim/core/rawdata.py @@ -342,19 +342,10 @@ def get_truth(self, instruction, truth_buffer): (self.config['samples_before_pulse_center'] + self.config['samples_after_pulse_center'] + 1) \ * self.config['sample_duration'] - channels = getattr(pulse, '_photon_channels', []) - - n_dpe = getattr(pulse, '_n_double_pe', 0) - n_dpe_bot = getattr(pulse, '_n_double_pe_bot', 0) - - tb['n_photon'] -= np.sum(np.isin(channels, getattr(pulse, 'turned_off_pmts', []))) - tb['n_pe'] += tb['n_photon']+n_dpe - # this turned_off guy, check how this works with a config['turned_off_guys'] - channels_bottom = list( - set(self.config['channels_bottom']).difference(getattr(pulse, 'turned_off_pmts', []))) - tb['n_photon_bottom'] = ( - np.sum(np.isin(channels, channels_bottom))) - tb['n_pe_bottom'] = tb['n_photon_bottom'] + n_dpe_bot + # Copy single valued fields directly from pulse class + for field in ['n_pe', 'n_pe_trigger', 'n_photon', 'n_photon_trigger', 'raw_area', 'raw_area_trigger']: + for suffix in ['', '_bottom']: + tb[field+suffix] = getattr(pulse, '_' + field + suffix, 0) # Summarize the instruction cluster in one row of the truth file for field in instruction.dtype.names: diff --git a/wfsim/strax_interface.py b/wfsim/strax_interface.py index cb281686..148ffb9f 100644 --- a/wfsim/strax_interface.py +++ b/wfsim/strax_interface.py @@ -43,11 +43,19 @@ truth_extra_dtype = [ (('End time of the interaction [ns]', 'endtime'), np.int64), - (('Number of simulated electrons', 'n_electron'), np.float64), - (('Number of detected photons', 'n_photon'), np.float64), - (('Number of detected photoelectrons(n_photons+dpe)', 'n_pe'), np.float64), - (('number of photons detected in bottom array', 'n_photon_bottom'), np.float64), - (('number of photoelectrons detected in bottom array', 'n_pe_bottom'), np.float64), + (('Number of simulated electrons', 'n_electron'), np.int32), + (('Number of photons reaching PMT', 'n_photon'), np.int32), + (('Number of photons + dpe passing', 'n_pe'), np.int32), + (('Number of photons passing trigger', 'n_photon_trigger'), np.int32), + (('Number of photons + dpe passing trigger', 'n_pe_trigger'), np.int32), + (('Raw area in pe', 'raw_area'), np.float64), + (('Raw area in pe passing trigger', 'raw_area_trigger'), np.float64), + (('Number of photons reaching PMT (bottom)', 'n_photon_bottom'), np.int32), + (('Number of photons + dpe passing (bottom)', 'n_pe_bottom'), np.int32), + (('Number of photons passing trigger (bottom)', 'n_photon_trigger_bottom'), np.int32), + (('Number of photons + dpe passing trigger (bottom)', 'n_pe_trigger_bottom'), np.int32), + (('Raw area in pe (bottom)', 'raw_area_bottom'), np.float64), + (('Raw area in pe passing trigger (bottom)', 'raw_area_trigger_bottom'), np.float64), (('Arrival time of the first photon [ns]', 't_first_photon'), np.float64), (('Arrival time of the last photon [ns]', 't_last_photon'), np.float64), (('Mean time of the photons [ns]', 't_mean_photon'), np.float64), @@ -522,7 +530,7 @@ def set_config(self,): # Update some values stored in CMT if self.config['fax_config_override_from_cmt'] is not None: for fax_field, cmt_option in self.config['fax_config_override_from_cmt'].items(): - if (fax_field in ['fdc_3d', 's1_light_yield_map'] + if (fax_field in ['fdc_3d', 's1_lce_correction_map'] and self.config.get('default_reconstruction_algorithm', False)): cmt_option = tuple(['suffix', self.config['default_reconstruction_algorithm'],