Merge pull request #9 from dkoh0207/develop

Added new cut varaibles and heuristics for nue selection.

Author: francois-drielsma

spine/data/out/interaction.py

@@ -279,6 +279,14 @@ def from_particles(cls, particles):
 @inherit_docstring(RecoBase, InteractionBase)
 class RecoInteraction(InteractionBase, RecoBase):
     """Reconstructed interaction information."""
+    
+    leading_shower_vertex_distance: float = -1.
+    leading_shower_axial_spread: float = None
+    leading_shower_directional_spread: float = None
+    leading_shower_start_dedx: float = None
+    leading_shower_trunk_straightness: float = None
+    leading_shower_adjacent_bragg_pearsonr: float = None
+    leading_shower_split_angle: float = None
 
     # Attributes that must never be stored to file
     _skip_attrs = (
@@ -307,7 +315,145 @@ def __str__(self):
             Basic information about the interaction properties
         """
         return 'Reco' + super().__str__()
+    
+    
+    @property
+    def leading_shower(self):
+        """Leading shower of this interaction
 
+        Returns
+        -------
+        RecoParticle
+            Primary shower with the highest reco kinetic energy.
+        """
+        showers = [part for part in self.primary_particles if part.shape == 0]
+        if len(showers) == 0:
+            return None
+        return max(showers, key=lambda x: x.ke)
+    
+    @leading_shower.setter
+    def leading_shower(self, leading_shower):
+        pass
+    
+    @property
+    def leading_shower_start_dedx(self):
+        """Leading shower dE/dx
+
+        Returns
+        -------
+        float
+            dE/dx of the leading shower
+        """
+        leading_shower = self.leading_shower
+        if leading_shower is None:
+            return -1.
+        return leading_shower.start_dedx
+    
+    @leading_shower_start_dedx.setter
+    def leading_shower_start_dedx(self, leading_shower_start_dedx):
+        pass
+    
+    @property
+    def leading_shower_directional_spread(self):
+        """Distance-weighted mean cosine distance from mean direction 
+        of the leading shower.
+        
+        See ShowerSpreadProcessor under post/reco/shower.py for more details.
+
+        Returns
+        -------
+        float
+            Measure of shower directional spread
+        """
+        leading_shower = self.leading_shower
+        if leading_shower is None:
+            return -1.
+        return leading_shower.directional_spread
+    
+    @leading_shower_directional_spread.setter
+    def leading_shower_directional_spread(self, leading_shower_directional_spread):
+        pass
+    
+    @property
+    def leading_shower_axial_spread(self):
+        """Pearson R correlation coefficient between the distance from 
+        shower start to a shower point (x) and the perpendicular distance
+        from the shower point to the shower axis. 
+
+        Returns
+        -------
+        coeff
+            Measure of shower axial spread
+        """
+        leading_shower = self.leading_shower
+        if leading_shower is None:
+            return -np.inf
+        return leading_shower.axial_spread
+    
+    @leading_shower_axial_spread.setter
+    def leading_shower_axial_spread(self, leading_shower_axial_spread):
+        pass
+    
+    @property
+    def leading_shower_trunk_straightness(self):
+        """Measure of how straight the trunk of the leading shower is.
+
+        Returns
+        -------
+        float
+            Measure of shower trunk straightness (between 0 and 1)
+        """
+        leading_shower = self.leading_shower
+        if leading_shower is None:
+            return -1.
+        return leading_shower.trunk_straightness
+    
+    @leading_shower_trunk_straightness.setter
+    def leading_shower_trunk_straightness(self, leading_shower_trunk_straightness):
+        pass
+    
+    @property
+    def leading_shower_adjacent_bragg_pearsonr(self):
+        """Pearson R correlation coefficient between:
+        x: distance from shower start to a point on a track adjacent to the
+        shower. 
+        y: the local dedx value at that point.
+        
+        If this value is strongly negative (close to -1), it suggests that 
+        the shower is likely to be adjacent to a bragg peak. 
+
+        Returns
+        -------
+        float
+            Measure of how likely the leading shower is adjacent 
+            to a bragg peak (between -1 and 1)
+        """
+        leading_shower = self.leading_shower
+        if leading_shower is None:
+            return -np.inf
+        return leading_shower.adjacent_bragg_pearsonr
+    
+    @leading_shower_adjacent_bragg_pearsonr.setter
+    def leading_shower_adjacent_bragg_pearsonr(self, leading_shower_adjacent_bragg_pearsonr):
+        pass
+    
+    @property
+    def leading_shower_split_angle(self):
+        """Angle between the two branches of the leading shower.
+
+        Returns
+        -------
+        float
+            Angle between the two branches of the leading shower
+        """
+        leading_shower = self.leading_shower
+        if leading_shower is None:
+            return -1.
+        return leading_shower.split_angle
+    
+    @leading_shower_split_angle.setter
+    def leading_shower_split_angle(self, leading_shower_split_angle):
+        pass
 
 @dataclass(eq=False)
 @inherit_docstring(TruthBase, InteractionBase)

spine/data/out/particle.py

@@ -218,13 +218,30 @@ class RecoParticle(ParticleBase, RecoBase):
     shower_split_angle: float
         Estimate of the opening angle of the shower. If particle is not a
         shower, then this is set to -1. (units of degrees)
+    shower_spread: float
+        Estimate of the spread of the shower (roughly mean cosine 
+        from mean direction)
+    axial_pearsonr: float
+        Pearson correlation coefficient of the axial profile of the shower.
+    trunk_validity: float
+        Explained variance ratio of the trunk of the particle. 
+    shower_dedx: float
+        dE/dx of the shower in MeV/cm, using direction reconstruction. 
+    adjacent_bragg: float
+        A number that represents whether the shower has a adjacent track
+        with a Bragg peak.
     """
     pid_scores: np.ndarray = None
     primary_scores: np.ndarray = None
     ppn_ids: np.ndarray = None
     ppn_points: np.ndarray = None
     vertex_distance: float = -1.
-    shower_split_angle: float = -1.
+    split_angle: float = -1.
+    directional_spread: float = -1.
+    axial_spread: float = -np.inf
+    trunk_straightness: float = -1.
+    start_dedx: float = -1.
+    adjacent_bragg_pearsonr: float = -np.inf
 
     # Fixed-length attributes
     _fixed_length_attrs = (

spine/driver.py

@@ -448,7 +448,7 @@ def get_prefixes(file_paths, split_output):
             log_prefix += f'--{suffix}'
 
         # Truncate file names that are too long
-        max_length = 230
+        max_length = 150
         if len(log_prefix) > max_length:
             log_prefix = log_prefix[:max_length-3] + '---'
 

spine/post/reco/__init__.py

@@ -13,3 +13,4 @@
 from .kinematics import *
 from .label import *
 from .shower import *
+from .topology import *

spine/post/reco/geometry.py

@@ -26,7 +26,8 @@ class ContainmentProcessor(PostBase):
 
     def __init__(self, margin, cathode_margin=None, detector=None,
                  geometry_file=None, mode='module',
-                 allow_multi_module=False, min_particle_sizes=0,
+                 allow_multi_module=False, skip_showers=False, 
+                 min_particle_sizes=0,
                  obj_type=('particle', 'interaction'),
                  truth_point_mode='points', run_mode='both'):
         """Initialize the containment conditions.
@@ -63,6 +64,8 @@ def __init__(self, margin, cathode_margin=None, detector=None,
               Note that this does not guarantee containment within the detector.
         allow_multi_module : bool, default False
             Whether to allow particles/interactions to span multiple modules
+        skip_showers : bool, default False
+            Whether to skip checking containment for showers
         min_particle_sizes : Union[int, dict], default 0
             When checking interaction containment, ignore particles below the
             size (in voxel count) specified by this parameter. If specified
@@ -87,6 +90,7 @@ def __init__(self, margin, cathode_margin=None, detector=None,
 
         # Store parameters
         self.allow_multi_module = allow_multi_module
+        self.skip_showers = skip_showers
 
         # Store the particle size thresholds in a dictionary
         if np.isscalar(min_particle_sizes):
@@ -118,6 +122,12 @@ def process(self, data):
             for obj in data[k]:
                 # Make sure the particle coordinates are expressed in cm
                 self.check_units(obj)
+                
+                # Skip showers for containment check if requested
+                if self.skip_showers:
+                    if hasattr(obj, 'shape') and obj.shape == 0:
+                        obj.is_contained = True
+                        continue
 
                 # Get point coordinates
                 points = self.get_points(obj)