Tracking dev progress

aux/src/ClusterHelpers.cxx

@@ -319,6 +319,8 @@ std::map<std::string, size_t> Aux::count(const cluster_graph_t& cgraph, bool nod
 }
 
 
+
+
 std::unordered_map<int, std::set<cluster_vertex_t> > Aux::blob_clusters(
     const cluster_graph_t& cg)
 {
@@ -368,4 +370,4 @@ std::unordered_map<int, std::set<cluster_vertex_t> > Aux::blob_clusters(
     }
 
     return groups;
-}
+}

cfg/pgrapher/common/clus.jsonnet

@@ -94,6 +94,17 @@ local wc = import "wirecell.jsonnet";
             } + dv_cfg + pcts_cfg
         },
 
+        tagger_check_neutrino(name="", trackfitting_config_file="", particle_dataset="", recombination_model="") :: {
+            type: "TaggerCheckNeutrino",
+            name: prefix + name,
+            data: {
+                grouping: "live",           // Which grouping to process
+                trackfitting_config_file: trackfitting_config_file, 
+                particle_dataset: particle_dataset,
+                recombination_model: recombination_model,
+            } + dv_cfg + pcts_cfg
+        },
+
         pointed(name="", groupings=["live"]) :: {
             type: "ClusteringPointed",
             name: prefix+name,

clus/inc/WireCellClus/Facade_Cluster.h

@@ -280,6 +280,30 @@ namespace WireCell::Clus::Facade {
             double charge_cut = 4000.0,
             bool disable_dead_mix_cell = true) const;
 
+        /// Get segment IDs for all points (computed from graph analysis)
+        /// @return Vector of segment IDs, -1 for unassigned points
+        std::vector<int> segment_ids() const;
+
+        /// Get shower flags for all points (computed from graph analysis)
+        /// @return Vector of flags, 1=shower, 0=track
+        std::vector<int> shower_flags() const;
+
+        /// Get segment ID for a specific point
+        /// @param point_index Global point index in cluster
+        /// @return Segment ID or -1 if unassigned
+        int segment_id(size_t point_index) const;
+
+        /// Get shower flag for a specific point
+        /// @param point_index Global point index in cluster
+        /// @return 1 if shower, 0 if track
+        int shower_flag(size_t point_index) const;
+
+        /// Invalidate cached segment data (IDs and shower flags)
+        /// Call this before re-computing segment information
+        void invalidate_segment_data() {
+            cache().invalidate_segment_data();
+        }
+
 
 
         // Return vector is size 3 holding vectors of size npoints providing k-d tree coordinate points.

clus/inc/WireCellClus/Facade_ClusterCache.h

@@ -74,11 +74,23 @@ namespace WireCell::Clus::Facade {
         // Set of point indices excluded during graph operations (equivalent to prototype's excluded_points)
         std::set<size_t> excluded_points;
 
+        // Segment IDs by point index (computed from graph analysis)
+        std::vector<int> point_segment_ids;
+
+        // Shower flags by point index (computed from graph analysis) 
+        std::vector<int> point_shower_flags;
+
         // Steiner point cloud k-d tree cache
         mutable std::unique_ptr<KDTree::MultiQuery> steiner_kd;
         mutable decltype(std::declval<KDTree::MultiQuery>().get<double>(std::vector<std::string>{})) steiner_query3d;
         mutable std::string cached_steiner_pc_name;
         mutable bool steiner_kd_built{false};
+
+        // Invalidate segment-related cached data
+        void invalidate_segment_data() {
+            point_segment_ids.clear();
+            point_shower_flags.clear();
+        }
     };
 
 }

clus/inc/WireCellClus/Facade_Util.h

@@ -145,6 +145,9 @@ namespace WireCell::Clus::Facade {
         results_type get_closest_index(const geo_point_t& p, const size_t N) const;
         /// @return index, geo_point_t
         std::pair<size_t, geo_point_t> get_closest_wcpoint(const geo_point_t& p) const;
+        double get_closest_dis(const geo_point_t& p) const;
+
+        std::vector<std::pair<size_t, geo_point_t>> get_closest_wcpoints_radius(const geo_point_t& p, const double radius) const;
 
         /// @param p_test1 is the point to start from
         /// @param dir is the direction to search along

clus/inc/WireCellClus/NeutrinoPatternBase.h

@@ -0,0 +1,110 @@
+#include "WireCellClus/PRGraph.h"
+#include "WireCellClus/Facade_Cluster.h"
+#include "WireCellClus/TrackFitting.h"
+
+namespace WireCell::Clus::PR {
+    class PatternAlgorithms{
+        public:
+        std::set<VertexPtr> find_cluster_vertices(Graph& graph, const Facade::Cluster& cluster);
+        std::set<SegmentPtr> find_cluster_segments(Graph& graph, const Facade::Cluster& cluster);
+        bool clean_up_graph(Graph& graph, const Facade::Cluster& cluster);
+
+        SegmentPtr init_first_segment(Graph& graph, Facade::Cluster& cluster, Facade::Cluster* main_cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv, bool flag_back_search = true);
+
+        // find the shortest path using steiner graph
+        std::vector<Facade::geo_point_t> do_rough_path(const Facade::Cluster& cluster,Facade::geo_point_t& first_point, Facade::geo_point_t& last_point);
+        std::vector<Facade::geo_point_t> do_rough_path_reg_pc(const Facade::Cluster& cluster, Facade::geo_point_t& first_point, Facade::geo_point_t& last_point,  std::string graph_name = "relaxed_pid");
+        // create a segment given a path
+        SegmentPtr create_segment_for_cluster(WireCell::Clus::Facade::Cluster& cluster, IDetectorVolumes::pointer dv, const std::vector<Facade::geo_point_t>& path_points, int dir = 0);
+        // create a segment given two vertices, null, if failed
+        SegmentPtr create_segment_from_vertices(Graph& graph, Facade::Cluster& cluster, VertexPtr v1, VertexPtr v2, IDetectorVolumes::pointer dv);
+        // replace a segment and vertex with another segment and vertex, assuming the original vertex only connect to this segment
+        bool replace_segment_and_vertex(Graph& graph, SegmentPtr& seg, VertexPtr& vtx, std::list<Facade::geo_point_t>& path_point_list, Facade::geo_point_t& break_point, IDetectorVolumes::pointer dv);
+        bool replace_segment_and_vertex(Graph& graph, SegmentPtr& seg, VertexPtr old_vertex, VertexPtr new_vertex, IDetectorVolumes::pointer dv);
+        bool break_segment_into_two(Graph& graph, VertexPtr vtx1, SegmentPtr seg, VertexPtr vtx2, std::list<Facade::geo_point_t>& path_point_list1, Facade::geo_point_t& break_point, std::list<Facade::geo_point_t>& path_point_list2, IDetectorVolumes::pointer dv);
+
+
+        // return the point and its index in the steiner tree as a pair
+        std::pair<Facade::geo_point_t, size_t> proto_extend_point(const Facade::Cluster& cluster, Facade::geo_point_t& p, Facade::geo_vector_t& dir, Facade::geo_vector_t& dir_other, bool flag_continue);
+        // return Steiner Graph path in wcps_list1 and wcps_list2
+        bool proto_break_tracks(const Facade::Cluster& cluster, const Facade::geo_point_t& first_wcp, const Facade::geo_point_t& curr_wcp, const Facade::geo_point_t& last_wcp, std::list<Facade::geo_point_t>& wcps_list1, std::list<Facade::geo_point_t>& wcps_list2, bool flag_pass_check);
+        // breaking segments ...
+        bool break_segments(Graph& graph, TrackFitting& track_fitter, IDetectorVolumes::pointer dv, std::vector<SegmentPtr>& remaining_segments, float dis_cut = 0);
+        // merge two segments to one
+        bool merge_two_segments_into_one(Graph& graph, SegmentPtr& seg1, VertexPtr& vtx, SegmentPtr& seg2, IDetectorVolumes::pointer dv);
+        // merge vertex into another
+        bool merge_vertex_into_another(Graph& graph, VertexPtr& vtx_from, VertexPtr& vtx_to, IDetectorVolumes::pointer dv);
+
+        // get direction with  distance cut ... 
+        Facade::geo_vector_t vertex_get_dir(VertexPtr& vertex, Graph& graph, double dis_cut = 5*units::cm);
+        Facade::geo_vector_t vertex_segment_get_dir(VertexPtr& vertex, SegmentPtr& segment, Graph& graph, double dis_cut = 2*units::cm);
+
+
+        // Structure examination
+        void examine_structure(Graph& graph, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv); // call examine_structure_1 and examine_structure_2      
+        bool examine_structure_1(Graph& graph, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool examine_structure_2(Graph& graph, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+
+        bool examine_structure_3(Graph& graph, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool examine_structure_4(VertexPtr vertex, bool flag_final_vertex, Graph& graph, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+
+        // identify other segments giving the graph ...
+        void find_other_segments(Graph& graph, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv , bool flag_break_track =true, double search_range = 1.5*units::cm, double scaling_2d = 0.8);
+
+        // examine segment
+        void examine_segment(Graph& graph, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool crawl_segment(Graph& graph, Facade::Cluster& cluster, SegmentPtr seg, VertexPtr vertex, TrackFitting& track_fitter, IDetectorVolumes::pointer dv );
+        void examine_partial_identical_segments(Graph& graph, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+
+        //examine vertices
+        void examine_vertices(Graph& graph, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool examine_vertices_1(Graph&graph, Facade::Cluster&cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool examine_vertices_1p(Graph&graph, VertexPtr v1, VertexPtr v2, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool examine_vertices_2(Graph&graph, Facade::Cluster&cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool examine_vertices_4(Graph&graph, Facade::Cluster&cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool examine_vertices_4p(Graph&graph, VertexPtr v1, VertexPtr v2, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        Facade::geo_point_t get_local_extension(Facade::Cluster& cluster, Facade::geo_point_t& wcp);
+        void examine_vertices_3(Graph& graph, Facade::Cluster& main_cluster, std::pair<VertexPtr, VertexPtr> main_cluster_initial_pair_vertices, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+
+        // master pattern recognition function
+        bool find_proto_vertex(Graph& graph, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv, bool flag_break_track = true, int nrounds_find_other_tracks = 2, bool flag_back_search = true);
+
+        void init_point_segment(Graph& graph, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+
+        // examine the structure of the patterns ... 
+        bool examine_structure_final_1(Graph& graph, VertexPtr main_vertex, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool examine_structure_final_1p(Graph& graph, VertexPtr main_vertex, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool examine_structure_final_2(Graph& graph, VertexPtr main_vertex, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool examine_structure_final_3(Graph& graph, VertexPtr main_vertex, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+        bool examine_structure_final(Graph& graph, VertexPtr main_vertex, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv);
+
+        // EM shower related
+        void clustering_points(Graph& graph, Facade::Cluster& cluster, const IDetectorVolumes::pointer& dv, const std::string& cloud_name = "associated_points", double search_range = 1.2*units::cm, double scaling_2d = 0.7);
+        void separate_track_shower(Graph&graph, Facade::Cluster& cluster);
+        // Direction
+        void determine_direction(Graph& graph, Facade::Cluster& cluster, const Clus::ParticleDataSet::pointer& particle_data, const IRecombinationModel::pointer& recomb_model);
+        std::pair<int, double> calculate_num_daughter_showers(Graph& graph, VertexPtr vertex, SegmentPtr segment, bool flag_count_shower = true); 
+        void examine_good_tracks(Graph& graph, Facade::Cluster& cluster, const Clus::ParticleDataSet::pointer& particle_data);
+        // about fix maps
+        void fix_maps_multiple_tracks_in(Graph& graph, Facade::Cluster& cluster);
+        void fix_maps_shower_in_track_out(Graph& graph, Facade::Cluster& cluster);
+        void improve_maps_one_in(Graph& graph, Facade::Cluster& cluster, const Clus::ParticleDataSet::pointer& particle_data, const IRecombinationModel::pointer& recomb_model, bool flag_strong_check = true);
+        void improve_maps_shower_in_track_out(Graph& graph, Facade::Cluster& cluster, const Clus::ParticleDataSet::pointer& particle_data, const IRecombinationModel::pointer& recomb_model, bool flag_strong_check = true);
+        void improve_maps_no_dir_tracks(Graph& graph, Facade::Cluster& cluster, const Clus::ParticleDataSet::pointer& particle_data, const IRecombinationModel::pointer& recomb_model);
+        void improve_maps_multiple_tracks_in(Graph& graph, Facade::Cluster& cluster, const Clus::ParticleDataSet::pointer& particle_data, const IRecombinationModel::pointer& recomb_model);
+        void judge_no_dir_tracks_close_to_showers(Graph& graph, Facade::Cluster& cluster, const Clus::ParticleDataSet::pointer& particle_data, IDetectorVolumes::pointer dv);
+        bool examine_maps(Graph&graph, Facade::Cluster& cluster);
+        void examine_all_showers(Graph& graph, Facade::Cluster& cluster, const Clus::ParticleDataSet::pointer& particle_data);
+        void shower_determining_in_main_cluster(Graph& graph, Facade::Cluster& cluster, const Clus::ParticleDataSet::pointer& particle_data, const IRecombinationModel::pointer& recomb_model, IDetectorVolumes::pointer dv);
+
+        // vertex related functions 
+        bool search_for_vertex_activities(Graph& graph, VertexPtr vertex, std::set<SegmentPtr>& segments_set, Facade::Cluster& cluster, TrackFitting& track_fitter, IDetectorVolumes::pointer dv, double search_range = 1.5*units::cm);
+
+        // global information transfer
+        void transfer_info_from_segment_to_cluster(Graph& graph, Facade::Cluster& cluster, const std::string& cloud_name = "associated_points");
+
+        // print information
+        void print_segs_info(Graph& graph, Facade::Cluster& cluster, VertexPtr vertex= nullptr);
+
+    };
+}

clus/inc/WireCellClus/PRCommon.h

@@ -95,18 +95,18 @@ namespace WireCell::Clus::PR {
     // multi APA/face detectors?
     struct WCPoint {
         WireCell::Point point;   // 3D point 
-        int uvw[3] = {-1,-1,-1}; // wire indices
-        int index{-1};           // point index in some container
+        // int uvw[3] = {-1,-1,-1}; // wire indices
+        // int index{-1};           // point index in some container
 
         // FIXME: WCP had this, does WCT need it?
         // blob* b;
 
 
         // Return true if the point information has been filled.
-        bool valid() const {
-            if (index < 0) return false;
-            return true;
-        }
+        // bool valid() const {
+        //     if (index < 0) return false;
+        //     return true;
+        // }
     };
     using WCPointVector = std::vector<WCPoint>;
 

clus/inc/WireCellClus/PRGraph.h

@@ -115,9 +115,27 @@ namespace WireCell::Clus::PR {
     /// The two Vertex objects are those associated with the source/target nodes
     /// of the segment's edge.  The pair is ordered.  The first Vertex is the
     /// one with a "wcpoint" closest to the segment's initial "wcpoint".
-    std::pair<VertexPtr, VertexPtr> find_endpoints(Graph& graph, SegmentPtr seg);
+    std::pair<VertexPtr, VertexPtr> find_vertices(Graph& graph, SegmentPtr seg);
 
-
+    /// Return the other vertex connected to a segment, given one known vertex.
+    ///
+    /// Returns nullptr if:
+    /// - The segment edge is not in the graph
+    /// - The given vertex is not connected to the segment
+    ///
+    /// This is useful when you have one vertex of a segment and need to find
+    /// the vertex at the other end.
+    VertexPtr find_other_vertex(Graph& graph, SegmentPtr seg, VertexPtr vertex);
+
+    /// Find the segment connecting two vertices.
+    ///
+    /// Returns nullptr if:
+    /// - Either vertex is not in the graph
+    /// - No edge exists between the two vertices
+    ///
+    /// This function searches for an edge between the two given vertices and
+    /// returns the associated segment if found.
+    SegmentPtr find_segment(Graph& graph, VertexPtr vtx1, VertexPtr vtx2);
 
 };
 #endif

clus/inc/WireCellClus/PRSegment.h

@@ -67,6 +67,8 @@ namespace WireCell::Clus::PR {
         // Getters
         const std::shared_ptr<Aux::ParticleInfo>& particle_info() const { return m_particle_info; }
         std::shared_ptr<Aux::ParticleInfo>& particle_info() { return m_particle_info; }
+        double particle_score() const { return m_particle_score; }
+        void particle_score(double score) { m_particle_score = score; }
 
         // Chainable setter
         Segment& particle_info(std::shared_ptr<Aux::ParticleInfo> pinfo) {
@@ -115,10 +117,30 @@ namespace WireCell::Clus::PR {
         bool fit_flag_skip(int i);
         void fit_flag_skip(int i, bool flag);
 
-        void set_fit_associate_vec(std::vector<WireCell::Point >& tmp_fit_pt_vec, std::vector<int>& tmp_fit_index, std::vector<bool>& tmp_fit_skip, const IDetectorVolumes::pointer& dv,const std::string& cloud_name="fit");
+        void set_fit_associate_vec(std::vector<PR::Fit >& tmp_fit_vec, const IDetectorVolumes::pointer& dv,const std::string& cloud_name="fit");
 
+        // Global indices management for point clouds
+        void set_global_indices(const std::string& cloud_name, std::vector<size_t> indices) {
+            m_pc_global_indices[cloud_name] = std::move(indices);
+        }
+
+        const std::vector<size_t>& global_indices(const std::string& cloud_name) const {
+            static const std::vector<size_t> empty_vec;
+            auto it = m_pc_global_indices.find(cloud_name);
+            return (it != m_pc_global_indices.end()) ? it->second : empty_vec;
+        }
+
+        bool has_global_indices(const std::string& cloud_name) const {
+            return m_pc_global_indices.find(cloud_name) != m_pc_global_indices.end();
+        }
+
+        // Add public accessor:
+        int id() const { return m_id; }
+        void set_id(int id) { m_id = id; }
 
     private:
+        // Add to PRSegment.h private section:
+        int m_id{-1};
 
         std::vector<WCPoint> m_wcpts;
         std::vector<Fit> m_fits;
@@ -127,6 +149,10 @@ namespace WireCell::Clus::PR {
         bool m_dir_weak{false};
 
         std::shared_ptr<Aux::ParticleInfo> m_particle_info{nullptr};
+        double m_particle_score{100};
+
+        // Mapping from local DPC index to global cluster point index
+        std::map<std::string, std::vector<size_t>> m_pc_global_indices;