Skip to content

Commit

Permalink
Refactor s2 correction (#704)
Browse files Browse the repository at this point in the history 
* refactor fields for correctedareas

* add corrected s2 with s2xy correction only

* fix dtype

* use infer_dtype

* remove redundant variables

* remove elife_corr from infer_dtype

* line break after operator

* use strax.to_str_tuple

* remove cs2_bottom_wo_elifecorr

* fix bug

Co-authored-by: Joran Angevaare <jorana@nikhef.nl>
  • Loading branch information
@JoranAngevaare
Jingqiang Ye and JoranAngevaare authored Oct 29, 2021
1 parent cf868e7 commit 1e9b592
Showing 1 changed file with 39 additions and 36 deletions.
75 changes: 39 additions & 36 deletions straxen/plugins/event_processing.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -635,44 +635,41 @@ class CorrectedAreas(strax.Plugin):
cs2_top and cs2_bottom are corrected by the corresponding maps,
and cs2 is the sum of the two.
"""
__version__ = '0.1.2'
__version__ = '0.1.3'

depends_on = ['event_basics', 'event_positions']
dtype = [('cs1', np.float32, 'Corrected S1 area [PE]'),
('cs2', np.float32, 'Corrected S2 area [PE]'),
('cs2_top', np.float32, 'Corrected S2 area in the top PMT array [PE]'),
('cs2_bottom', np.float32, 'Corrected S2 area in the bottom PMT array [PE]'),
('alt_cs1', np.float32, 'Corrected area of the alternate S1 [PE]'),
('alt_cs2', np.float32, 'Corrected area of the alternate S2 [PE]'),
('alt_cs2_top', np.float32, 'Corrected area of the alternate S2 in the top PMT array [PE]'),
('alt_cs2_bottom', np.float32, 'Corrected area of the alternate S2 in the bottom PMT array [PE]'),
('cs2_area_fraction_top', np.float32, 'Main cS2 fraction of area seen by the top PMT array'),
('alt_cs2_area_fraction_top', np.float32, 'Alternate cS2 fraction of area seen by the top PMT array'),
('elife_correction', np.float32, 'Main cS2 correction factor due to electron lifetime'),
('alt_elife_correction', np.float32, 'Alt cS2 correction factor due to electron lifetime'),

] + strax.time_fields

def infer_dtype(self):
dtype = []
dtype += strax.time_fields

for peak_type, peak_name in zip(['', 'alt_'], ['main', 'alternate']):
dtype += [(f'{peak_type}cs1', np.float32, f'Corrected area of {peak_name} S1 [PE]'),
(f'{peak_type}cs2_wo_elifecorr', np.float32,
f'Corrected area of {peak_name} S2 before elife correction (s2 xy correction only) [PE]'),
(f'{peak_type}cs2_area_fraction_top', np.float32,
f'Fraction of area seen by the top PMT array for corrected {peak_name} S2'),
(f'{peak_type}cs2_bottom', np.float32,
f'Corrected area of {peak_name} S2 in the bottom PMT array [PE]'),
(f'{peak_type}cs2', np.float32, f'Corrected area of {peak_name} S2 [PE]'),]

return dtype

def setup(self):
self.elife = get_correction_from_cmt(self.run_id, self.config['elife_conf'])

if isinstance(self.config['s1_xyz_correction_map'], str):
self.config['s1_xyz_correction_map'] = [self.config['s1_xyz_correction_map']]
if isinstance(self.config['s2_xy_correction_map'], str):
self.config['s2_xy_correction_map'] = [self.config['s2_xy_correction_map']]

self.s1_map = InterpolatingMap(
get_cmt_resource(self.run_id,
tuple(['suffix',
self.config['default_reconstruction_algorithm'],
*self.config['s1_xyz_correction_map']]),
*strax.to_str_tuple(self.config['s1_xyz_correction_map'])]),
fmt='text'))

self.s2_map = InterpolatingMap(
get_cmt_resource(self.run_id,
tuple(['suffix',
self.config['default_reconstruction_algorithm'],
*self.config['s2_xy_correction_map']]),
*strax.to_str_tuple(self.config['s2_xy_correction_map'])]),
fmt='json'))

def compute(self, events):
Expand All @@ -699,23 +696,29 @@ def compute(self, events):
s2_bottom_map_name = "map"

for peak_type in ["", "alt_"]:
# For electron lifetime corrections to the S2s,
# use lifetimes computed using the main S1.
el_string = peak_type + "s2_interaction_" if peak_type == "alt_" else peak_type
result[f"{peak_type}elife_correction"] = np.exp(events[f'{el_string}drift_time'] / self.elife)

# S2(x,y) corrections use the observed S2 positions
s2_positions = np.vstack([events[f'{peak_type}s2_x'], events[f'{peak_type}s2_y']]).T

result[f"{peak_type}cs2_top"] = (events[f'{peak_type}s2_area'] * events[f'{peak_type}s2_area_fraction_top']
* result[f"{peak_type}elife_correction"]
/ self.s2_map(s2_positions, map_name=s2_top_map_name))
result[f"{peak_type}cs2_bottom"] = (events[f'{peak_type}s2_area']
* (1 - events[f'{peak_type}s2_area_fraction_top'])
* result[f"{peak_type}elife_correction"]
/ self.s2_map(s2_positions, map_name=s2_bottom_map_name))
result[f"{peak_type}cs2"] = result[f"{peak_type}cs2_top"] + result[f"{peak_type}cs2_bottom"]
result[f"{peak_type}cs2_area_fraction_top"] = result[f"{peak_type}cs2_top"] / result[f"{peak_type}cs2"]
# corrected s2 with s2 xy map only, i.e. no elife correction
# this is for s2-only events which don't have drift time info
cs2_top_wo_elifecorr = (events[f'{peak_type}s2_area'] * events[f'{peak_type}s2_area_fraction_top'] /
self.s2_map(s2_positions, map_name=s2_top_map_name))
cs2_bottom_wo_elifecorr = (events[f'{peak_type}s2_area'] *
(1 - events[f'{peak_type}s2_area_fraction_top']) /
self.s2_map(s2_positions, map_name=s2_bottom_map_name))
result[f"{peak_type}cs2_wo_elifecorr"] = (cs2_top_wo_elifecorr + cs2_bottom_wo_elifecorr)

# cs2aft doesn't need elife correction as it cancels
result[f"{peak_type}cs2_area_fraction_top"] = cs2_top_wo_elifecorr / result[f"{peak_type}cs2_wo_elifecorr"]

# For electron lifetime corrections to the S2s,
# use drift time computed using the main S1.
el_string = peak_type + "s2_interaction_" if peak_type == "alt_" else peak_type
elife_correction = np.exp(events[f'{el_string}drift_time'] / self.elife)

cs2_top = cs2_top_wo_elifecorr * elife_correction
result[f"{peak_type}cs2_bottom"] = cs2_bottom_wo_elifecorr * elife_correction
result[f"{peak_type}cs2"] = cs2_top + result[f"{peak_type}cs2_bottom"]

return result

Expand Down

0 comments on commit 1e9b592

Please sign in to comment.