Skip to content

Commit

Permalink
<833> <Refactor xy algorithms>
Browse files Browse the repository at this point in the history 
#833

[author: gonzaponte]

This PR provides a refactor of the basic reconstruction algorithms. It
includes, among other things:

- the decoupling of corona and barycenter
- the addition (and check) of annotations to the xy reconstruction functions
- a generalization of penthesilea and dorothea to allow different reconstruction algorithms
- a refactor of the functions related to hit-building used in penthesilea

Closes #587.

[reviewer: paolafer]

This is a very nice, long-due PR cleaning up the xy reconstruction algorithms. Approved!
  • Loading branch information
@paolafer @carhc
paolafer authored and carhc committed Feb 19, 2024
2 parents 94c2f21 + 8ee3052 commit 86e39c4
Show file tree
Hide file tree
Showing 18 changed files with 399 additions and 281 deletions.
78 changes: 47 additions & 31 deletions invisible_cities/cities/components.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@
from .. reco import wfm_functions as wfm
from .. reco import paolina_functions as plf
from .. reco .xy_algorithms import corona
from .. reco .xy_algorithms import barycenter
from .. filters.s1s2_filter import S12Selector
from .. filters.s1s2_filter import pmap_filter
from .. database import load_db
Expand All @@ -77,6 +78,7 @@
from .. types .symbols import EventRange
from .. types .symbols import HitEnergy
from .. types .symbols import SiPMCharge
from .. types .symbols import XYReco


def city(city_function):
Expand Down Expand Up @@ -778,13 +780,17 @@ def peak_classifier(**params):
return partial(pmap_filter, selector)


def compute_xy_position(dbfile, run_number, **reco_params):
# `reco_params` is the set of parameters for the corona
# algorithm either for the full corona or for barycenter
datasipm = load_db.DataSiPM(dbfile, run_number)
def compute_xy_position(dbfile, run_number, algo, **reco_params):
if algo is XYReco.corona:
datasipm = load_db.DataSiPM(dbfile, run_number)
reco_params = dict(all_sipms = datasipm, **reco_params)
algorithm = corona
else:
algorithm = barycenter

def compute_xy_position(xys, qs):
return corona(xys, qs, datasipm, **reco_params)
return algorithm(xys, qs, **reco_params)

return compute_xy_position


Expand Down Expand Up @@ -861,35 +867,46 @@ def build_pointlike_event(pmap, selector_output, event_number, timestamp):
return build_pointlike_event


def hit_builder(dbfile, run_number, drift_v, reco,
rebin_slices, rebin_method,
global_reco_params, charge_type):
def get_s1_time(pmap, selector_output):
# in order to compute z one needs to define one S1
# for time reference. By default the filter will only
# take events with exactly one s1. Otherwise, the
# convention is to take the first peak in the S1 object
# as reference.
if np.any(selector_output.s1_peaks):
first_s1 = np.where(selector_output.s1_peaks)[0][0]
s1_t = pmap.s1s[first_s1].time_at_max_energy
else:
first_s2 = np.where(selector_output.s2_peaks)[0][0]
s1_t = pmap.s2s[first_s2].times[0]

return s1_t


def try_global_reco(reco, xys, qs):
try : cluster = reco(xys, qs)[0]
except XYRecoFail: return xy.empty()
else : return xy(cluster.X, cluster.Y)


def sipm_positions(dbfile, run_number):
datasipm = load_db.DataSiPM(dbfile, run_number)
sipm_xs = datasipm.X.values
sipm_ys = datasipm.Y.values
sipm_xys = np.stack((sipm_xs, sipm_ys), axis=1)
return sipm_xys

sipm_noise = NoiseSampler(dbfile, run_number).signal_to_noise

barycenter = partial(corona, all_sipms = datasipm, **global_reco_params)

def empty_cluster():
return Cluster(NN, xy(0,0), xy(0,0), 0)
def hit_builder(dbfile, run_number, drift_v,
rebin_slices, rebin_method,
global_reco, slice_reco,
charge_type):
sipm_xys = sipm_positions(dbfile, run_number)
sipm_noise = NoiseSampler(dbfile, run_number).signal_to_noise

def build_hits(pmap, selector_output, event_number, timestamp):
hitc = HitCollection(event_number, timestamp * 1e-3)

# in order to compute z one needs to define one S1
# for time reference. By default the filter will only
# take events with exactly one s1. Otherwise, the
# convention is to take the first peak in the S1 object
# as reference.
if np.any(selector_output.s1_peaks):
first_s1 = np.where(selector_output.s1_peaks)[0][0]
s1_t = pmap.s1s[first_s1].time_at_max_energy
else:
first_s2 = np.where(selector_output.s2_peaks)[0][0]
s1_t = pmap.s2s[first_s2].times[0]
s1_t = get_s1_time(pmap, selector_output)

# here hits are computed for each peak and each slice.
# In case of an exception, a hit is still created with a NN cluster.
Expand All @@ -904,26 +921,25 @@ def build_hits(pmap, selector_output, event_number, timestamp):
xys = sipm_xys[peak.sipms.ids]
qs = peak.sipm_charge_array(sipm_noise, charge_type,
single_point = True)
try : cluster = barycenter(xys, qs)[0]
except XYRecoFail: xy_peak = xy(NN, NN)
else : xy_peak = xy(cluster.X, cluster.Y)
xy_peak = try_global_reco(global_reco, xys, qs)

sipm_charge = peak.sipm_charge_array(sipm_noise ,
charge_type ,
single_point=False)

for slice_no, (t_slice, qs) in enumerate(zip(peak.times ,
sipm_charge)):
z_slice = (t_slice - s1_t) * units.ns * drift_v
e_slice = peak.pmts.sum_over_sensors[slice_no]
try:
xys = sipm_xys[peak.sipms.ids]
clusters = reco(xys, qs)
clusters = slice_reco(xys, qs)
es = hif.split_energy(e_slice, clusters)
for c, e in zip(clusters, es):
hit = Hit(peak_no, c, z_slice, e, xy_peak)
hit = Hit(peak_no, c, z_slice, e, xy_peak)
hitc.hits.append(hit)
except XYRecoFail:
hit = Hit(peak_no, empty_cluster(), z_slice,
hit = Hit(peak_no, Cluster.empty(), z_slice,
e_slice, xy_peak)
hitc.hits.append(hit)

Expand Down
4 changes: 3 additions & 1 deletion invisible_cities/cities/components_test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
from .. core import system_of_units as units
from .. types.symbols import WfType
from .. types.symbols import EventRange as ER
from .. types.symbols import XYReco

from . components import event_range
from . components import collect
Expand Down Expand Up @@ -98,7 +99,8 @@ def test_compute_xy_position_depends_on_actual_run_number():
'msipm': 9,
'consider_masked': True}
run_number = 6977
find_xy_pos = compute_xy_position('new', run_number, **reco_parameters)
find_xy_pos = compute_xy_position( 'new', run_number
, XYReco.corona, **reco_parameters)

xs_to_test = np.array([-65, -65, -55, -55, -55, -45, -45, -45])
ys_to_test = np.array([ 5, 25, 5, 15, 25, 5, 15, 25])
Expand Down
50 changes: 26 additions & 24 deletions invisible_cities/cities/dorothea.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@
from .. io.run_and_event_io import run_and_event_writer
from .. io. event_filter_io import event_filter_writer
from .. types.symbols import SiPMCharge
from .. types.symbols import XYReco

from .. dataflow import dataflow as fl
from .. dataflow.dataflow import push
Expand All @@ -64,29 +65,29 @@


@city
def dorothea( files_in : OneOrManyFiles
, file_out : str
, compression : str
, event_range : EventRangeType
, print_mod : int
, detector_db : str
, run_number : int
, drift_v : float
, s1_nmin : int, s1_nmax : int
, s1_emin : float, s1_emax : float
, s1_wmin : float, s1_wmax : float
, s1_hmin : float, s1_hmax : float
, s1_ethr : float
, s2_nmin : int, s2_nmax : int
, s2_emin : float, s2_emax : float
, s2_wmin : float, s2_wmax : float
, s2_hmin : float, s2_hmax : float
, s2_ethr : float
, s2_nsipmmin : int, s2_nsipmmax : int
, global_reco_params: dict
, sipm_charge_type : SiPMCharge
, include_mc : Optional[bool] = False
):
def dorothea( files_in : OneOrManyFiles
, file_out : str
, compression : str
, event_range : EventRangeType
, print_mod : int
, detector_db : str
, run_number : int
, drift_v : float
, s1_nmin : int, s1_nmax : int
, s1_emin : float, s1_emax : float
, s1_wmin : float, s1_wmax : float
, s1_hmin : float, s1_hmax : float
, s1_ethr : float
, s2_nmin : int, s2_nmax : int
, s2_emin : float, s2_emax : float
, s2_wmin : float, s2_wmax : float
, s2_hmin : float, s2_hmax : float
, s2_ethr : float
, s2_nsipmmin : int, s2_nsipmmax : int
, global_reco_algo : XYReco, global_reco_params: dict
, sipm_charge_type : SiPMCharge
, include_mc : Optional[bool] = False
):
# global_reco_params are qth, qlm, lm_radius, new_lm_radius, msipm
# qlm = 0 * pes every Cluster must contain at least one SiPM with charge >= qlm
# lm_radius = -1 * mm by default, use overall barycenter for KrCity
Expand All @@ -101,7 +102,8 @@ def dorothea( files_in : OneOrManyFiles
pmap_passed = fl.map(attrgetter("passed"), args="selector_output", out="pmap_passed")
pmap_select = fl.count_filter(bool, args="pmap_passed")

reco_algo = compute_xy_position(detector_db, run_number, **global_reco_params)
reco_algo = compute_xy_position( detector_db, run_number
, global_reco_algo, **global_reco_params)
build_pointlike_event = fl.map(build_pointlike_event_( detector_db, run_number, drift_v
, reco_algo, sipm_charge_type),
args = ("pmap", "selector_output", "event_number", "timestamp"),
Expand Down
56 changes: 29 additions & 27 deletions invisible_cities/cities/penthesilea.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,8 @@
from .. dataflow.dataflow import push
from .. dataflow.dataflow import pipe
from .. types.symbols import RebinMethod
from .. types.symbols import SiPMCharge
from .. types.symbols import SiPMCharge
from .. types.symbols import XYReco

from . components import city
from . components import copy_mc_info
Expand All @@ -54,26 +55,28 @@
from typing import Optional

@city
def penthesilea( files_in : OneOrManyFiles
, file_out : str
, compression : str
, event_range : EventRangeType
, print_mod : int
, detector_db : str
, run_number : int
, drift_v : float
, rebin : int
, s1_nmin : int, s1_nmax : int
, s1_emin : float, s1_emax : float
, s1_wmin : float, s1_wmax : float
, s1_hmin : float, s1_hmax : float
, s1_ethr : float
, s2_nmin : int, s2_nmax : int
, s2_emin : float, s2_emax : float
, s2_wmin : float, s2_wmax : float
, s2_hmin : float, s2_hmax : float
, s2_ethr : float
, s2_nsipmmin : int, s2_nsipmmax : int
def penthesilea( files_in : OneOrManyFiles
, file_out : str
, compression : str
, event_range : EventRangeType
, print_mod : int
, detector_db : str
, run_number : int
, drift_v : float
, rebin : int
, s1_nmin : int, s1_nmax : int
, s1_emin : float, s1_emax : float
, s1_wmin : float, s1_wmax : float
, s1_hmin : float, s1_hmax : float
, s1_ethr : float
, s2_nmin : int, s2_nmax : int
, s2_emin : float, s2_emax : float
, s2_wmin : float, s2_wmax : float
, s2_hmin : float, s2_hmax : float
, s2_ethr : float
, s2_nsipmmin : int, s2_nsipmmax : int
, slice_reco_algo : XYReco
, global_reco_algo : XYReco
, slice_reco_params : dict
, global_reco_params : dict
, rebin_method : RebinMethod
Expand All @@ -82,6 +85,8 @@ def penthesilea( files_in : OneOrManyFiles

# slice_reco_params are qth, qlm, lm_radius, new_lm_radius, msipm used for hits reconstruction
# global_reco_params are qth, qlm, lm_radius, new_lm_radius, msipm used for overall global (pointlike event) reconstruction
slice_reco = compute_xy_position(detector_db, run_number, slice_reco_algo, ** slice_reco_params)
global_reco = compute_xy_position(detector_db, run_number, global_reco_algo, **global_reco_params)


classify_peaks = df.map(peak_classifier(**locals()),
Expand All @@ -91,17 +96,14 @@ def penthesilea( files_in : OneOrManyFiles
pmap_passed = df.map(attrgetter("passed"), args="selector_output", out="pmap_passed")
pmap_select = df.count_filter(bool, args="pmap_passed")

reco_algo_slice = compute_xy_position(detector_db, run_number, **slice_reco_params)
build_hits = df.map(hit_builder(detector_db, run_number, drift_v,
reco_algo_slice, rebin,
rebin_method,
global_reco_params,
rebin, rebin_method,
global_reco, slice_reco,
sipm_charge_type),
args = ("pmap", "selector_output", "event_number", "timestamp"),
out = "hits" )
reco_algo_global = compute_xy_position(detector_db, run_number, **global_reco_params)
build_pointlike_event = df.map(build_pointlike_event_( detector_db, run_number, drift_v
, reco_algo_global, sipm_charge_type),
, global_reco, sipm_charge_type),
args = ("pmap", "selector_output", "event_number", "timestamp"),
out = "pointlike_event" )

Expand Down
9 changes: 1 addition & 8 deletions invisible_cities/config/barycenter.conf
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,3 @@
barycenter_params = dict(
Qthr = 1 * pes,
lm_radius = -1 * mm ,

# the previous parameter sets barycenter as reco algorithm.
# The following arguments are not necessary
#Qlm
#new_lm_radius
#msipm
Qthr = 1 * pes
)
1 change: 1 addition & 0 deletions invisible_cities/config/dorothea.conf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ s2_nsipmmax = 1000

include('$ICDIR/config/barycenter.conf')

global_reco_algo = barycenter
global_reco_params = barycenter_params
del barycenter_params

Expand Down
1 change: 1 addition & 0 deletions invisible_cities/config/dorothea_corona.conf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,5 +44,6 @@ s2_nsipmmax = 1000

include('$ICDIR/config/corona.conf')

global_reco_algo = corona
global_reco_params = corona_params
del corona_params
9 changes: 3 additions & 6 deletions invisible_cities/config/penthesilea.conf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ rebin_method = stride

sipm_charge_type = raw

slice_reco_algo = corona
slice_reco_params = dict(
Qthr = 2 * pes,
Qlm = 5 * pes,
Expand All @@ -40,9 +41,5 @@ slice_reco_params = dict(
msipm = 1 )


global_reco_params = dict(
Qthr = 1 * pes,
Qlm = 0 * pes,
lm_radius = -1 * mm ,
new_lm_radius = -1 * mm ,
msipm = 1 )
global_reco_algo = barycenter
global_reco_params = dict(Qthr = 1 * pes)
Loading

0 comments on commit 86e39c4

Please sign in to comment.