Skip to content

Commit

Permalink
Merge branch 'develop' into feature/ehinkle_geo_resource_charge_coord…
Browse files Browse the repository at this point in the history
…_fix
  • Loading branch information
edhinkle authored Dec 15, 2023
2 parents 9393d6d + 3a03882 commit 67676c1
Show file tree
Hide file tree
Showing 4 changed files with 36 additions and 10 deletions.
7 changes: 3 additions & 4 deletions src/proto_nd_flow/reco/charge/calib_hit_merger.py
Original file line number Diff line number Diff line change
Expand Up @@ -241,7 +241,6 @@ def merge_hits(self,hits, weights, seg_fracs, dt_cut, sum_fields=None, weighted_

# calculate segment contributions for each merged hit
if has_mc_truth:
tmp_bt = np.full(shape=new_hits.shape+(2,self.max_contrib_segments),fill_value=0.)
back_track = np.full(shape=new_hits.shape,fill_value=0.,dtype=self.hit_frac_dtype)
# loop over hits
for hit_it, hit in np.ndenumerate(new_hits):
Expand All @@ -250,12 +249,12 @@ def merge_hits(self,hits, weights, seg_fracs, dt_cut, sum_fields=None, weighted_
# renormalize the fractional contributions given the charge weighted average
norm = np.sum(np.multiply(hit_contr[0],hit_contr[1]))
if norm == 0.: norm = 1.
tmp_bt[hit_it][0] = np.multiply(hit_contr[0],hit_contr[1])/norm # fractional contributions
tmp_bt[hit_it][1] = hit_contr[2] # segment_ids
tmp_bt_0 = np.multiply(hit_contr[0],hit_contr[1])/norm # fractional contributions
tmp_bt_1 = hit_contr[2] # segment_ids

# merge unique track contributions
track_dict = defaultdict(lambda:0)
for track in zip(tmp_bt[hit_it][0],tmp_bt[hit_it][1]):
for track in zip(tmp_bt_0,tmp_bt_1):
track_dict[track[1]] += track[0]
track_dict = dict(track_dict)
bt_unique_segs = np.array(list(track_dict.keys()))
Expand Down
6 changes: 6 additions & 0 deletions src/proto_nd_flow/reco/charge/calib_prompt_hits.py
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,8 @@ class CalibHitBuilder(H5FlowStage):
z f8, pixel z location [cm]
t_drift u8, drift time [ticks???]
ts_pps f8, PPS packet timestamp [ns]
io_group u8, io group ID (PACMAN number)
io_channel u8, io channel ID (related to PACMAN number & PACMAN UART Number)
Q f8, hit charge [ke-]
E f8, hit energy [MeV]
Expand All @@ -82,6 +84,8 @@ class CalibHitBuilder(H5FlowStage):
('z', 'f8'),
('t_drift', 'f8'),
('ts_pps', 'u8'),
('io_group', 'u8'),
('io_channel', 'u8'),
('Q', 'f8'),
('E', 'f8')
])
Expand Down Expand Up @@ -212,6 +216,8 @@ def run(self, source_name, source_slice, cache):
calib_hits_arr['z'] = zy[:,0]
calib_hits_arr['ts_pps'] = raw_hits_arr['ts_pps']
calib_hits_arr['t_drift'] = drift_t
calib_hits_arr['io_group'] = packets_arr['io_group']
calib_hits_arr['io_channel'] = packets_arr['io_channel']
calib_hits_arr['Q'] = self.charge_from_dataword(packets_arr['dataword'],vref,vcm,ped)
calib_hits_arr['E'] = self.charge_from_dataword(packets_arr['dataword'],vref,vcm,ped) * 23.6e-3 # hardcoding W_ion and not accounting for finite electron lifetime

Expand Down
31 changes: 26 additions & 5 deletions src/proto_nd_flow/resources/geometry.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
import numpy as np
import numpy.ma as ma
import logging
import warnings
import yaml

from h5flow.core import H5FlowResource
Expand All @@ -22,6 +23,8 @@ class Geometry(H5FlowResource):
Parameters:
- ``path``: ``str``, path to stored geometry data within file
- ``network_agnostic``: ``bool``, optional, ignore the (io_channel % 4), useful if the io_channel mapping changes run-to-run (``True == ignore``)
- ``n_io_channels_per_tile``: ``int``, optional, only used with ``network_agnostic == True``, sets the number of io channels to have same geometry info
- ``crs_geometry_file``: ``str``, path to yaml file describing charge readout system geometry
- ``det_geometry_file``: ``str``, path to yaml file describing overall detector geometry
- ``lrs_geometry_file``: ``str``, path to yaml file describing light readout system
Expand Down Expand Up @@ -72,9 +75,11 @@ class Geometry(H5FlowResource):
lrs_geometry_file: 'data/proto_nd_flow/light_module_desc-0.0.0.yaml'
'''
class_version = '0.1.0'
class_version = '0.2.0'

default_path = 'geometry_info'
default_network_agnostic = False
default_n_io_channels_per_tile = 4
default_det_geometry_file = '-'
default_crs_geometry_file = '-'
default_lrs_geometry_file = '-'
Expand All @@ -86,6 +91,8 @@ def __init__(self, **params):
super(Geometry, self).__init__(**params)

self.path = params.get('path', self.default_path)
self.network_agnostic = params.get('network_agnostic', self.default_network_agnostic)
self.n_io_channels_per_tile = params.get('n_io_channels_per_tile', self.default_n_io_channels_per_tile)
self.crs_geometry_file = params.get('crs_geometry_file', self.default_crs_geometry_file)
self.det_geometry_file = params.get('det_geometry_file', self.default_crs_geometry_file)
self.lrs_geometry_file = params.get('lrs_geometry_file', self.default_lrs_geometry_file)
Expand Down Expand Up @@ -119,6 +126,8 @@ def init(self, source_name):
max_drift_distance=self.max_drift_distance,
module_RO_bounds=self.module_RO_bounds,
pixel_pitch=self.pixel_pitch,
network_agnostic=self.network_agnostic,
n_io_channels_per_tile=self.n_io_channels_per_tile
)
write_lut(self.data_manager, self.path, self.anode_drift_coordinate, 'anode_drift_coordinate')
write_lut(self.data_manager, self.path, self.drift_dir, 'drift_dir')
Expand Down Expand Up @@ -150,7 +159,7 @@ def init(self, source_name):
lut_size = (self.anode_drift_coordinate.nbytes + self.drift_dir.nbytes
+ self.pixel_coordinates_2D.nbytes + self.tile_id.nbytes
+ self.tpc_id.nbytes + self.det_id.nbytes
+ self.det_bounds.nbytes) * 4
+ self.det_bounds.nbytes)

if self.rank == 0:
logging.info(f'Geometry LUT(s) size: {lut_size/1024/1024:0.02f}MB')
Expand Down Expand Up @@ -558,7 +567,7 @@ def _load_charge_geometry(self):

tiles = np.arange(1,len(geometry_yaml['tile_chip_to_io'])*len(det_geometry_yaml['module_to_io_groups'])+1)
io_groups = [
geometry_yaml['tile_chip_to_io'][tile][chip] // 1000 * (mod-1)*2
geometry_yaml['tile_chip_to_io'][tile][chip] // 1000 + (mod-1)*2
for tile in geometry_yaml['tile_chip_to_io']
for chip in geometry_yaml['tile_chip_to_io'][tile]
for mod in module_to_io_groups
Expand All @@ -582,7 +591,7 @@ def _load_charge_geometry(self):

pixel_coordinates_2D_min_max = [(min(v), max(v)) for v in (io_groups, io_channels, chip_ids, channel_ids)]
self._pixel_coordinates_2D = LUT('f4', *pixel_coordinates_2D_min_max, shape=(2,))
self._pixel_coordinates_2D.default = 0.
self._pixel_coordinates_2D.default = np.nan

tile_min_max = [(min(v), len(module_to_io_groups)*max(v)) for v in (io_groups, io_channels)]
self._tile_id = LUT('i4', *tile_min_max)
Expand Down Expand Up @@ -612,14 +621,26 @@ def _load_charge_geometry(self):
io_channel = io_group_io_channel % 1000
self._tile_id[([io_group], [io_channel])] = tile+(module_id-1)*len(tile_chip_to_io)

if self.network_agnostic == True:
# if we don't care about the network configuration, then we
# can just loop over every N io channels and add them to the LUT
start_io_channel = ((io_channel-1)//self.n_io_channels_per_tile)*self.n_io_channels_per_tile + 1
for io_channel in range(start_io_channel, start_io_channel+self.n_io_channels_per_tile):
self._tile_id[([io_group], [io_channel])] = tile

for chip_channel in chip_channel_to_position:
chip = chip_channel // 1000
channel = chip_channel % 1000

try:
io_group_io_channel = tile_chip_to_io[tile][chip]
except KeyError:
continue
if self.network_agnostic == True:
warnings.warn('Encountered an out-of-network chip, but because you enabled ``network_agnostic``, we will carry on with assumptions about the io group and io channel')
# using the info about the first chip on the tile for all others
io_group_io_channel = list(geometry_yaml['tile_chip_to_io'][tile].values())[0]
else:
continue

io_group = io_group_io_channel // 1000 + (module_id-1)*len(det_geometry_yaml['module_to_io_groups'][module_id])
io_channel = io_group_io_channel % 1000
Expand Down
2 changes: 1 addition & 1 deletion src/proto_nd_flow/resources/run_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,7 @@ def _lookup_mc_info(self):
# mc info has already exists, return
return

if self.input_filename[-3:] == '.h5' or '.h5' in self.input_filename:
if any(self.input_filename.endswith(ext) for ext in ['.h5', '.hdf5']):
if H5FLOW_MPI:
with h5py.File(self.input_filename, 'r', driver='mpio', comm=self.comm) as f:
is_mc = 'mc_packets_assn' in f
Expand Down

0 comments on commit 67676c1

Please sign in to comment.