centralmembraneanalysis/CentralMembraneStripeMatching/GlobalFieldFitter.h

0001 #ifndef TPCCALIB_GLOBALFIELDFITTER_H
0002 #define TPCCALIB_GLOBALFIELDFITTER_H
0003
0004 #include <array>
0005 #include <vector>
0006
0007 // A field observation is std::array<double, 6>.
0008 // It represents one measured/reference stripe displacement used in the fit.
0009 inline constexpr int observation_phi = 0;
0010 inline constexpr int observation_r = 1;
0011 inline constexpr int observation_delta_r = 2;
0012 inline constexpr int observation_r_delta_phi = 3;
0013 inline constexpr int observation_sigma_delta_r = 4;
0014 inline constexpr int observation_sigma_r_delta_phi = 5;
0015
0016 // A fitted control point is std::array<double, 7>.
0017 // The last entry is stored as 0.0 or 1.0 because this version avoids custom
0018 // record objects and keeps everything in basic STL containers.
0019 inline constexpr int control_point_phi = 0;
0020 inline constexpr int control_point_r = 1;
0021 inline constexpr int control_point_delta_r = 2;
0022 inline constexpr int control_point_r_delta_phi = 3;
0023 inline constexpr int control_point_delta_phi = 4;
0024 inline constexpr int control_point_support = 5;
0025 inline constexpr int control_point_active = 6;
0026
0027 class GlobalFieldFitter {
0028 public:
0029   // observations: each entry is indexed with observation_* above.
0030   // control_r_positions: optional radial grid from lamination gap centers.
0031   explicit GlobalFieldFitter(const std::vector<std::array<double, 6>> &observations, const std::vector<double> &control_r_positions = {});
0032
0033   // Fits DeltaR and R*DeltaPhi on one shared phi/R control grid.
0034   bool fit();
0035
0036   // True only after fit() succeeds and coefficient arrays are usable.
0037   bool is_valid() const;
0038
0039   // Evaluate the fitted distortion components at any phi/R point.
0040   double evaluate_delta_r(double phi, double r) const;
0041   double evaluate_r_delta_phi(double phi, double r) const;
0042   double evaluate_delta_phi(double phi, double r) const;
0043
0044   // Returns robust global fallback uncertainties from the fitted residuals.
0045   double predictive_sigma_delta_r(double phi, double r) const;
0046   double predictive_sigma_r_delta_phi(double phi, double r) const;
0047
0048   // Returns std::array<double, 7> points indexed with control_point_*.
0049   std::vector<std::array<double, 7>> control_points() const;
0050
0051 private:
0052   // Build the rectangular phi/R control grid before solving.
0053   void initialize_grid();
0054
0055   // Converts a 2D control-grid coordinate into a flat vector index.
0056   int control_index(int phi_index, int r_index) const;
0057
0058   // Finds the four surrounding controls and bilinear weights for one point.
0059   void accumulate_bilinear_weights(double phi, double r, std::vector<int> &indices, std::vector<double> &weights) const;
0060
0061   // Adds finite-difference smoothing terms to the normal equations.
0062   void add_smoothness_penalty(std::vector<std::vector<double>> &normal_matrix, const std::vector<int> &full_to_fit_index) const;
0063
0064   // Plain Gaussian elimination solver for the small dense normal equations.
0065   bool solve_linear_system(std::vector<std::vector<double>> matrix, std::vector<double> rhs, std::vector<double> &solution) const;
0066
0067   // Evaluates one fitted component from a coefficient vector.
0068   double evaluate_component(const std::vector<double> &coefficients, double phi, double r) const;
0069
0070   // Counts how much observation weight reaches each control point.
0071   void compute_control_support();
0072
0073   // Chooses active controls when sparse interpolation is enabled.
0074   void select_active_controls(std::vector<int> &full_to_fit_index);
0075
0076   // Fills inactive controls from neighboring active controls for evaluation.
0077   void build_dense_evaluation_grid();
0078
0079   // Control-grid coordinate helpers.
0080   double control_phi(int phi_index) const;
0081   double control_r(int r_index) const;
0082
0083   // Interpolates a value at an inactive control from nearby active controls.
0084   double interpolate_active_value(const std::vector<double> &coefficients, double phi, double r) const;
0085
0086   // Robust residual scale used when writing uncertainty diagnostics.
0087   double estimate_global_residual_sigma(bool fit_delta_r) const;
0088
0089   // Raw fitted observations indexed by observation_*.
0090   std::vector<std::array<double, 6>> m_observations;
0091
0092   // Compact fitted coefficient vectors. These may contain only active controls.
0093   std::vector<double> m_coefficients_delta_r;
0094   std::vector<double> m_coefficients_r_delta_phi;
0095
0096   // Dense coefficient vectors used for fast evaluation everywhere on the grid.
0097   std::vector<double> m_evaluation_coefficients_delta_r;
0098   std::vector<double> m_evaluation_coefficients_r_delta_phi;
0099
0100   // Per-control support and active flags.
0101   std::vector<double> m_controlSupport;
0102   std::vector<bool> m_activeControl;
0103
0104   // Control-grid dimensions and numeric ranges.
0105   int m_nControlPhi = 0;
0106   int m_nControlR = 0;
0107   double m_phiMin = 0.0;
0108   double m_phiMax = 0.0;
0109   double m_phiStep = 0.0;
0110   double m_rMin = 0.0;
0111   double m_rMax = 0.0;
0112   double m_rStep = 0.0;
0113
0114   // Requested radial controls from the caller and final sorted usable controls.
0115   std::vector<double> m_requestedControlRPositions;
0116   std::vector<double> m_controlRPositions;
0117
0118   // Set true only after both DeltaR and R*DeltaPhi solves succeed.
0119   bool m_isValid = false;
0120 };
0121
0122 #endif