From bd35bd992cd6f5565e6027fb648212287c3679a6 Mon Sep 17 00:00:00 2001 From: JanWeldert Date: Thu, 23 Oct 2025 08:09:18 -0500 Subject: [PATCH 1/5] Allow relative pid improvement and remove pid stage --- pisa/stages/README.md | 3 +- pisa/stages/pid/README.md | 3 - pisa/stages/pid/__init__.py | 0 pisa/stages/pid/shift_scale_pid.py | 116 ----------------------------- pisa/stages/reco/resolutions.py | 19 ++++- 5 files changed, 16 insertions(+), 125 deletions(-) delete mode 100644 pisa/stages/pid/README.md delete mode 100644 pisa/stages/pid/__init__.py delete mode 100644 pisa/stages/pid/shift_scale_pid.py diff --git a/pisa/stages/README.md b/pisa/stages/README.md index 28450dcda..0e433c0c0 100644 --- a/pisa/stages/README.md +++ b/pisa/stages/README.md @@ -12,8 +12,7 @@ Directories are PISA stages, and within each directory can be found the services * `flux/` - All stages relating to the atmospheric neutrino flux. * `likelihood/` - A stage that pre-computes some quantities needed for the "generalized likelihood" * `osc/` - All stages relating to neutrino oscillations. -* `pid/` - All stages relating to particle identification. -* `reco/` - All stages relating to applying reconstruction kernels. +* `reco/` - All stages relating to applying reconstruction kernels (including pid). * `utils/` - All "utility" stages (not representing physics effects). * `xsec/` - All stages relating to cross sections. * `GLOBALS.md` - File that describes globally available variables within PISA together with the services that implement them. diff --git a/pisa/stages/pid/README.md b/pisa/stages/pid/README.md deleted file mode 100644 index c48cd8ec1..000000000 --- a/pisa/stages/pid/README.md +++ /dev/null @@ -1,3 +0,0 @@ -# Stage: Particle ID - -The purpose of this stage is to change the event classification, classifying the reconstructed events the track and cascade channels. diff --git a/pisa/stages/pid/__init__.py b/pisa/stages/pid/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/pisa/stages/pid/shift_scale_pid.py b/pisa/stages/pid/shift_scale_pid.py deleted file mode 100644 index 0e2e1dedc..000000000 --- a/pisa/stages/pid/shift_scale_pid.py +++ /dev/null @@ -1,116 +0,0 @@ -""" -The purpose of this stage is shift and/or scale the pid values. -""" - -from __future__ import absolute_import, print_function, division - -from numba import guvectorize -import numpy as np - -from pisa import FTYPE, TARGET -from pisa.core.param import Param, ParamSet -from pisa.core.stage import Stage -from pisa.utils import vectorizer - -__all__ = ['shift_scale_pid', 'calculate_pid_function', 'init_test'] - -__author__ = 'L. Fischer' - - -class shift_scale_pid(Stage): # pylint: disable=invalid-name - """ - Shift/scale pid. - - Parameters - ---------- - params - bias : float - shift pid values by given bias - scale : float - scale pid values by given scale factor - """ - - def __init__(self, - **std_kwargs, - ): - - # register expected parameters - expected_params = ('bias', 'scale',) - - expected_container_keys = ('pid', ) - - # init base class - super().__init__( - expected_params=expected_params, - expected_container_keys=expected_container_keys, - **std_kwargs, - ) - - assert self.calc_mode == 'events' - - def setup_function(self): - """Setup the stage""" - - # set the correct data mode - self.data.representation = self.calc_mode - for container in self.data: - container['calculated_pid'] = np.empty((container.size), dtype=FTYPE) - container['original_pid'] = np.empty((container.size), dtype=FTYPE) - vectorizer.assign(vals=container['pid'], out=container['original_pid']) - - def compute_function(self): - """Perform computation""" - - # bias/scale have no units. - bias = self.params.bias.m_as('dimensionless') - scale = self.params.scale.m_as('dimensionless') - - for container in self.data: - calculate_pid_function(bias, - scale, - container['original_pid'], - out=container['calculated_pid']) - container.mark_changed('calculated_pid') - - def apply_function(self): - for container in self.data: - # set the pid value to the calculated one - vectorizer.assign(vals=container['calculated_pid'], out=container['pid']) - -signatures = [ - '(f4[:], f4[:], f4[:], f4[:])', - '(f8[:], f8[:], f8[:], f8[:])' -] - -layout = '(),(),()->()' - - -@guvectorize(signatures, layout, target=TARGET) -def calculate_pid_function(bias_value, scale_factor, pid, out): - """This function selects a pid cut by shifting the pid variable so - the default cut at 1.0 is at the desired cut position. - - Parameters - ---------- - bias_value : scalar - shift pid values by this bias - scale_factor : scalar - scale pid values with this factor - pid : scalar - pid variable - out : scalar - shifted pid values - - """ - - out[0] = (scale_factor[0] * pid[0]) + bias_value[0] - - -def init_test(**param_kwargs): - """Instantiation example""" - param_set = ParamSet([ - Param(name='bias', value=0.0, **param_kwargs), - Param(name='scale', value=1.0, **param_kwargs) - ]) - - return shift_scale_pid(calc_mode='events', params=param_set) diff --git a/pisa/stages/reco/resolutions.py b/pisa/stages/reco/resolutions.py index e8fbad45d..996f42190 100644 --- a/pisa/stages/reco/resolutions.py +++ b/pisa/stages/reco/resolutions.py @@ -4,6 +4,7 @@ """ from __future__ import absolute_import, print_function, division +import numpy as np from pisa.core.param import Param, ParamSet from pisa.core.stage import Stage @@ -28,7 +29,8 @@ class resolutions(Stage): # pylint: disable=invalid-name coszen_improvement : quantity (dimensionless) scale the reco error down by this fraction pid_improvement : quantity (dimensionless) - applies a shift to the classification parameter + applies a shift to the classification parameter [if relative_pid=False] + scales the pid error down by this fraction [if relative_pid=True] Expected container keys are .. :: @@ -41,8 +43,10 @@ class resolutions(Stage): # pylint: disable=invalid-name """ def __init__( self, + relative_pid=False, **std_kwargs ): + expected_params = ( 'energy_improvement', 'coszen_improvement', @@ -64,6 +68,7 @@ def __init__( **std_kwargs, ) + self.relative_pid = relative_pid def setup_function(self): @@ -76,14 +81,20 @@ def setup_function(self): logging.info('Changing coszen resolutions') container['reco_coszen'] += (container['true_coszen'] - container['reco_coszen']) * self.params.coszen_improvement.m_as('dimensionless') + container['reco_coszen'] = np.clip(container['reco_coszen'], -1, 1) container.mark_changed('reco_coszen') - # TODO: make sure coszen is within -1/1 ? logging.info('Changing PID resolutions') if container.name in ['numu_cc', 'numubar_cc']: - container['pid'] += self.params.pid_improvement.m_as('dimensionless') + if self.relative_pid: + container['pid'] += (1 - container['pid']) * self.params.pid_improvement.m_as('dimensionless') + else: + container['pid'] += self.params.pid_improvement.m_as('dimensionless') else: - container['pid'] -= self.params.pid_improvement.m_as('dimensionless') + if self.relative_pid: + container['pid'] += (0 - container['pid']) * self.params.pid_improvement.m_as('dimensionless') + else: + container['pid'] -= self.params.pid_improvement.m_as('dimensionless') container.mark_changed('pid') From fe68d4d1c78870cc20426938d8724facdb3f610d Mon Sep 17 00:00:00 2001 From: Jan Weldert <32642322+JanWeldert@users.noreply.github.com> Date: Mon, 27 Oct 2025 17:08:20 +0100 Subject: [PATCH 2/5] Remove pid stage from globals --- pisa/stages/GLOBALS.md | 1 - 1 file changed, 1 deletion(-) diff --git a/pisa/stages/GLOBALS.md b/pisa/stages/GLOBALS.md index fc2bb0cf2..a55b02916 100644 --- a/pisa/stages/GLOBALS.md +++ b/pisa/stages/GLOBALS.md @@ -66,7 +66,6 @@ Also note that where a service implements `FTYPE` and relies on C extension code | `osc.nusquids` | :heavy_minus_sign: | :heavy_check_mark: | :heavy_check_mark: | | `osc.prob3` | :heavy_check_mark: | :heavy_check_mark: | :heavy_check_mark: | | `osc.two_nu_osc` | :heavy_check_mark: | :heavy_check_mark: | :heavy_check_mark: | -| `pid.shift_scale_pid` | :heavy_check_mark: | :heavy_check_mark: | :heavy_check_mark: | | `reco.resolutions` | :heavy_minus_sign: | :heavy_minus_sign: | :heavy_minus_sign: | | `reco.simple_param` | :heavy_minus_sign: | :heavy_minus_sign: | :heavy_check_mark: | | `utils.add_indices` | :heavy_minus_sign: | :heavy_minus_sign: | :heavy_minus_sign: | From aed4234c5db302c635662d1af49d152208b77907 Mon Sep 17 00:00:00 2001 From: JanWeldert Date: Tue, 16 Dec 2025 10:57:10 -0600 Subject: [PATCH 3/5] Clip number of events after Gaussian smearing --- pisa/core/map.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pisa/core/map.py b/pisa/core/map.py index 5140f550c..94afb21fd 100755 --- a/pisa/core/map.py +++ b/pisa/core/map.py @@ -1217,6 +1217,7 @@ def fluctuate(self, method, random_state=None, jumpahead=None): loc=orig_hist[valid_mask], scale=sigma[valid_mask], random_state=random_state ) + gauss = np.clip(gauss, a_min=0, a_max=None) hist_vals = np.empty_like(orig_hist, dtype=np.float64) hist_vals[valid_mask] = poisson.rvs( From 986089edc3fe82a33bb93e1c54ca0fda481ac8f0 Mon Sep 17 00:00:00 2001 From: JanWeldert Date: Tue, 16 Dec 2025 11:02:06 -0600 Subject: [PATCH 4/5] Do not create linked container if none of the containers exist --- pisa/core/container.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/pisa/core/container.py b/pisa/core/container.py index 7ede83f1f..a1cd0217f 100644 --- a/pisa/core/container.py +++ b/pisa/core/container.py @@ -125,10 +125,12 @@ def link_containers(self, key, names): % (set(names) - set(self.names)) ) containers = [self.__getitem__(name) for name in link_names] - logging.trace('Linking containers %s into %s'%(link_names, key)) - new_container = VirtualContainer(key, containers) - self.linked_containers.append(new_container) - + if len(containers) > 0: + logging.trace('Linking containers %s into %s'%(link_names, key)) + new_container = VirtualContainer(key, containers) + self.linked_containers.append(new_container) + else: + logging.warning("None of the containers exist, skipping %s"%(key)) def unlink_containers(self): """Unlink all container""" From 4fd7f0c59a06b0fe1be12253336b3b3b219629ad Mon Sep 17 00:00:00 2001 From: JanWeldert Date: Wed, 17 Dec 2025 02:40:14 -0600 Subject: [PATCH 5/5] Make adaptive KDE an option --- pisa/stages/utils/kde.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pisa/stages/utils/kde.py b/pisa/stages/utils/kde.py index a09bd947e..669c82634 100644 --- a/pisa/stages/utils/kde.py +++ b/pisa/stages/utils/kde.py @@ -55,6 +55,7 @@ def __init__( coszen_name="reco_coszen", oversample=10, coszen_reflection=0.25, + adaptive=True, alpha=0.1, stack_pid=True, stash_hists=False, @@ -70,6 +71,7 @@ def __init__( self.oversample = int(oversample) self.coszen_reflection = float(coszen_reflection) self.alpha = float(alpha) + self.adaptive = adaptive self.stack_pid = stack_pid self.stash_hists = stash_hists self.stash_valid = False @@ -185,6 +187,7 @@ def apply_function(self): bw_method=self.bw_method, coszen_name=self.coszen_name, coszen_reflection=self.coszen_reflection, + adaptive=self.adaptive, alpha=self.alpha, oversample=self.oversample, use_cuda=False,