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/MyFCN.h

@@ -0,0 +1,112 @@
+#ifndef WIRECELLCLUS_MYFCN_H
+#define WIRECELLCLUS_MYFCN_H
+
+#include "WireCellClus/PRGraph.h"
+#include "WireCellClus/Facade_Cluster.h"
+#include "WireCellUtil/Units.h"
+#include "WireCellClus/TrackFitting.h"
+
+
+#include <vector>
+#include <tuple>
+
+namespace WireCell::Clus::PR {
+
+    /// MyFCN: Vertex fitting class for pattern recognition
+    /// Performs vertex position fitting by minimizing distances from segments
+    class MyFCN {
+    public:
+        /// Constructor
+        /// @param vtx The vertex to fit
+        /// @param flag_vtx_constraint Whether to constrain the vertex position
+        /// @param vtx_constraint_range Range for vertex constraint
+        /// @param vertex_protect_dis Protection distance for vertex
+        /// @param vertex_protect_dis_short_track Protection distance for short tracks
+        /// @param fit_dis Fitting distance
+        MyFCN(VertexPtr vtx, 
+              bool flag_vtx_constraint = false, 
+              double vtx_constraint_range = 1*units::cm, 
+              double vertex_protect_dis = 1.5*units::cm, 
+              double vertex_protect_dis_short_track = 0.9*units::cm, 
+              double fit_dis = 6*units::cm);
+
+        ~MyFCN();
+
+        /// Update the fitting range parameters
+        void update_fit_range(double tmp_vertex_protect_dis = 1.5*units::cm, 
+                             double tmp_vertex_protect_dis_short_track = 0.9*units::cm, 
+                             double tmp_fit_dis = 6*units::cm);
+
+        /// Add a segment to the fitting
+        void AddSegment(SegmentPtr sg);
+
+        /// Fit the vertex position
+        /// @return pair of (success flag, fitted position)
+        std::pair<bool, Facade::geo_point_t> FitVertex();
+
+        /// Update information after fitting
+        /// @param fit_pos The fitted position
+        /// @param temp_cluster The cluster being processed
+        /// @param default_dis_cut Default distance cut
+        void UpdateInfo(Facade::geo_point_t fit_pos, 
+                       Facade::Cluster& temp_cluster,
+                       TrackFitting& track_fitter, IDetectorVolumes::pointer dv, 
+                       double default_dis_cut = 4.0*units::cm);
+
+        /// Get segment information at index i
+        /// @return pair of (segment pointer, index)
+        std::pair<SegmentPtr, int> get_seg_info(int i);
+
+        /// Get number of fittable tracks
+        int get_fittable_tracks();
+
+        /// Get vertex constraint flag
+        bool get_flag_vtx_constraint() { return flag_vtx_constraint; }
+
+        /// Set vertex constraint flag
+        void set_flag_vtx_constraint(bool val) { flag_vtx_constraint = val; }
+
+        /// Set vertex constraint range
+        void set_vtx_constraint_range(double val) { vtx_constraint_range = val; }
+
+        /// Get the vector of segments
+        std::vector<SegmentPtr>& get_segments() { return segments; }
+
+        /// Get the vector of point vectors
+        std::vector<std::vector<Facade::geo_point_t>>& get_vec_points() { return vec_points; }
+
+        /// Print points for debugging
+        void print_points();
+
+        /// Set enforce two track fit flag
+        void set_enforce_two_track_fit(bool val) { enforce_two_track_fit = val; }
+
+        /// Get enforce two track fit flag
+        bool get_enforce_two_track_fit() { return enforce_two_track_fit; }
+
+    private:
+        VertexPtr vtx;
+        bool enforce_two_track_fit;
+        bool flag_vtx_constraint;
+        double vtx_constraint_range;
+
+        double vertex_protect_dis;
+        double vertex_protect_dis_short_track;
+        double fit_dis;
+
+        std::vector<SegmentPtr> segments;
+        std::vector<std::vector<Facade::geo_point_t>> vec_points;
+
+        // PCA directions: tuple of (dir1, dir2, dir3)
+        std::vector<std::tuple<Facade::geo_point_t, Facade::geo_point_t, Facade::geo_point_t>> vec_PCA_dirs;
+
+        // PCA eigenvalues: tuple of (val1, val2, val3)
+        std::vector<std::tuple<double, double, double>> vec_PCA_vals;
+
+        // Centers for each segment
+        std::vector<Facade::geo_point_t> vec_centers;
+    };
+
+} // namespace WireCell::Clus::PR
+
+#endif // WIRECELLCLUS_MYFCN_H