Merge branch 'main' into find-clusters-refactor

CMakeLists.txt

@@ -169,7 +169,11 @@ if(AARE_FETCH_MINUIT2)
   set(BUILD_TESTING
       OFF
       CACHE BOOL "")
+  set(MINUIT2_INSTALL
+      ON
+      CACHE BOOL "")
   FetchContent_MakeAvailable(Minuit2)
+
   set_property(TARGET Minuit2Math PROPERTY POSITION_INDEPENDENT_CODE ON)
   set_property(TARGET Minuit2 PROPERTY POSITION_INDEPENDENT_CODE ON)
 else()
@@ -238,17 +242,8 @@ if(AARE_FETCH_FMT)
   set(FMT_INSTALL
       ON
       CACHE BOOL "")
-  # set(FMT_CMAKE_DIR "")
   FetchContent_MakeAvailable(fmt)
   set_property(TARGET fmt PROPERTY POSITION_INDEPENDENT_CODE ON)
-  install(
-    TARGETS fmt
-    EXPORT ${project}-targets
-    LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
-    ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
-    RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
-    INCLUDES
-    DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
 else()
   find_package(fmt 6 REQUIRED)
 endif()
@@ -261,10 +256,6 @@ if(AARE_FETCH_JSON)
       ON
       CACHE BOOL "")
   FetchContent_MakeAvailable(json)
-  set(NLOHMANN_JSON_TARGET_NAME nlohmann_json)
-
-  install(TARGETS nlohmann_json EXPORT "${TARGETS_EXPORT_NAME}")
-  message(STATUS "target: ${NLOHMANN_JSON_TARGET_NAME}")
 else()
   find_package(nlohmann_json 3.11.3 REQUIRED)
 endif()

include/aare/Chi2.hpp

@@ -138,6 +138,12 @@ class Chi2Model1DGrad : public ROOT::Minuit2::FCNGradientBase {
 // ── Convenient aliases ──────────────────────────────────────────────
 
 using Chi2Gaussian = Chi2Model1DGrad<aare::model::Gaussian>;
+using Chi2GaussianErfcPlateau =
+    Chi2Model1DGrad<aare::model::GaussianErfcPlateau>;
+using Chi2GaussianChargeSharing =
+    Chi2Model1DGrad<aare::model::GaussianChargeSharing>;
+using Chi2GaussianChargeSharingKb =
+    Chi2Model1DGrad<aare::model::GaussianChargeSharingKb>;
 using Chi2RisingScurve = Chi2Model1DGrad<aare::model::RisingScurve>;
 using Chi2FallingScurve = Chi2Model1DGrad<aare::model::FallingScurve>;
 using Chi2Pol1 = Chi2Model1DGrad<aare::model::Pol1>;

include/aare/ClusterFile.hpp

@@ -383,7 +383,7 @@ ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
             "Unexpected error (not feof or ferror) when reading frame number");
     }
 
-    int32_t n_clusters; // Saved as 32bit integer in the cluster file
+    uint32_t n_clusters;
     if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
         throw std::runtime_error(LOCATION +
                                  "Could not read number of clusters");

include/aare/Models.hpp

@@ -340,6 +340,510 @@ struct Gaussian {
     }
 };
 
+// _____________________________________________________________________
+//
+// GaussianErfcPlateau
+// _____________________________________________________________________
+
+/**
+ * @brief Gaussian peak plus a falling error-function plateau.
+ *
+ * Equivalent to the ROOT-style expression:
+ *
+ *   [0] * exp(-0.5 * ((x - [2])/[3])^2)
+ * + [1] * (1 - erf((x - [2]) / ([3] * sqrt(2)))) / 2
+ *
+ * f(x) = A * exp(-(x - mu)^2 / (2*sigma^2))
+ *      + S * 0.5 * (1 - erf((x - mu) / (sqrt(2)*sigma)))
+ *
+ * Parameters:
+ *   par[0] = A      (Gaussian amplitude)
+ *   par[1] = S      (left plateau height)
+ *   par[2] = mu     (shared Gaussian center / step inflection point)
+ *   par[3] = sigma  (shared Gaussian width / step width, must be > 0)
+ *
+ * Analytic partial derivatives:
+ *   Let dx = x - mu
+ *       z  = dx / (sqrt(2)*sigma)
+ *       E  = exp(-z^2)
+ *       H  = 0.5 * (1 - erf(z))
+ *
+ *   df/dA     = E
+ *   df/dS     = H
+ *   df/dmu    = A*E*dx/sigma^2 + S*E/(sqrt(2*pi)*sigma)
+ *   df/dsigma = A*E*dx^2/sigma^3 + S*E*dx/(sqrt(2*pi)*sigma^2)
+ */
+struct GaussianErfcPlateau {
+    static constexpr std::size_t npar = 4;
+
+    static constexpr std::array<ParamInfo, npar> param_info = {{
+        {"A", -no_bound, no_bound},
+        {"S", -no_bound, no_bound},
+        {"mu", -no_bound, no_bound},
+        {"sigma", 1e-12, no_bound},
+    }};
+
+    static double eval(double x, const std::vector<double> &par) {
+        const double A = par[0];
+        const double S = par[1];
+        const double mu = par[2];
+        const double sig = par[3];
+
+        const double dx = x - mu;
+        const double z = dx * inv_sqrt2 / sig;
+
+        const double e = std::exp(-z * z);
+        const double step = 0.5 * (1.0 - fast_erf(z));
+
+        return A * e + S * step;
+    }
+
+    static void eval_and_grad(double x, const std::vector<double> &par,
+                              double &f, std::array<double, npar> &g) {
+        const double A = par[0];
+        const double S = par[1];
+        const double mu = par[2];
+        const double sig = par[3];
+
+        const double dx = x - mu;
+        const double z = dx * inv_sqrt2 / sig;
+
+        const double e = std::exp(-z * z);
+        const double step = 0.5 * (1.0 - fast_erf(z));
+
+        f = A * e + S * step;
+
+        const double inv_sig = 1.0 / sig;
+        const double inv_sig2 = inv_sig * inv_sig;
+        const double inv_sig3 = inv_sig2 * inv_sig;
+
+        g[0] = e;    // df/dA
+        g[1] = step; // df/dS
+
+        g[2] = A * e * dx * inv_sig2 + S * inv_sqrt_2pi * e * inv_sig; // df/dmu
+
+        g[3] = A * e * dx * dx * inv_sig3 +
+               S * inv_sqrt_2pi * e * dx * inv_sig2; // df/dsigma
+    }
+
+    /** @brief Reject degenerate or negative width. */
+    static bool is_valid(const std::vector<double> &par) {
+        return par[3] > 0.0;
+    }
+
+    /**
+     * @brief Data-driven initial parameter estimates.
+     *
+     * Assumes x is sorted ascending and the plateau is on the left.
+     *
+     * S      ~ average of first 10% of y values
+     * mu     ~ x at maximum y
+     * A      ~ y_max - 0.5*S
+     * sigma  ~ right-side half-width of the peak, fallback to 10% x-range
+     */
+    static std::array<double, npar> estimate_par(NDView<double, 1> x,
+                                                 NDView<double, 1> y) {
+        const ssize_t n = y.size();
+
+        const auto max_it = std::max_element(y.begin(), y.end());
+        const ssize_t i_max = std::distance(y.begin(), max_it);
+
+        const double y_max = *max_it;
+
+        const double x_range = std::max(x[n - 1] - x[0], 1e-9);
+
+        const ssize_t tail = std::min<ssize_t>(std::max<ssize_t>(n / 10, 2), n);
+
+        double left_mean = 0.0;
+        double right_mean = 0.0;
+
+        for (ssize_t i = 0; i < tail; ++i)
+            left_mean += y[i];
+        left_mean /= tail;
+
+        for (ssize_t i = n - tail; i < n; ++i)
+            right_mean += y[i];
+        right_mean /= tail;
+
+        const double S = left_mean;
+        const double mu = x[i_max];
+
+        // At x = mu, the erfc step contributes roughly S/2.
+        const double A = y_max - 0.5 * S;
+
+        // Estimate sigma from the right half-width of the peak.
+        // This avoids the left plateau biasing the FWHM estimate.
+        const double half = right_mean + 0.5 * (y_max - right_mean);
+
+        double sig = 0.1 * x_range;
+        for (ssize_t i = i_max; i < n; ++i) {
+            if (y[i] <= half) {
+                const double half_width = std::abs(x[i] - mu);
+                if (half_width > 1e-12) {
+                    sig = half_width / 1.1774100225154747; // sqrt(2*ln(2))
+                }
+                break;
+            }
+        }
+
+        sig = std::max(sig, 1e-6);
+
+        return {A, S, mu, sig};
+    }
+
+    static void compute_steps(const std::array<double, npar> &start,
+                              double x_range, double y_range,
+                              double /*slope_scale*/,
+                              std::array<double, npar> &steps) {
+        steps[0] = std::max(0.1 * std::abs(start[0]), 0.1 * y_range);
+        steps[1] = std::max(0.1 * std::abs(start[1]), 0.1 * y_range);
+        steps[2] = 0.05 * x_range;
+        steps[3] = 0.05 * x_range;
+    }
+};
+
+// _____________________________________________________________________
+//
+// GaussianChargeSharing
+// _____________________________________________________________________
+
+/**
+ * @brief Gaussian peak with erfc charge-sharing tail and linear pedestal.
+ *
+ * Equivalent to the old ROOT energyCalibration function:
+ *
+ *   gaussChargeSharing(x, par)
+ *
+ * for sign = +1:
+ *
+ *   f(x) = p0 - p1*x
+ *        + N * [ G(x; mu, sigma) + C * H(x; mu, sigma) ]
+ *
+ * where:
+ *
+ *   G = exp(-0.5 * ((x - mu)/sigma)^2)
+ *   H = 0.5 * erfc((x - mu)/(sqrt(2)*sigma))
+ *
+ * Parameters:
+ *   par[0] = p0     background pedestal
+ *   par[1] = p1     background slope, ROOT convention uses p0 - p1*x
+ *   par[2] = mu     peak position
+ *   par[3] = sigma  noise RMS / Gaussian width
+ *   par[4] = N      number of photons / global amplitude
+ *   par[5] = C      charge-sharing pedestal
+ */
+struct GaussianChargeSharing {
+    static constexpr std::size_t npar = 6;
+
+    static constexpr std::array<ParamInfo, npar> param_info = {{
+        {"p0", -no_bound, no_bound},
+        {"p1", -no_bound, no_bound},
+        {"mu", -no_bound, no_bound},
+        {"sigma", 1e-12, no_bound},
+        {"N", -no_bound, no_bound},
+        {"C", -no_bound, no_bound},
+    }};
+
+    static double eval(double x, const std::vector<double> &par) {
+        const double p0 = par[0];
+        const double p1 = par[1];
+        const double mu = par[2];
+        const double sig = par[3];
+        const double N = par[4];
+        const double C = par[5];
+
+        const double dx = x - mu;
+        const double u = dx / sig;
+
+        const double G = std::exp(-0.5 * u * u);
+        const double H = 0.5 * (1.0 - fast_erf(u * inv_sqrt2));
+
+        return p0 - p1 * x + N * (G + C * H);
+    }
+
+    static void eval_and_grad(double x, const std::vector<double> &par,
+                              double &f, std::array<double, npar> &g) {
+        const double p0 = par[0];
+        const double p1 = par[1];
+        const double mu = par[2];
+        const double sig = par[3];
+        const double N = par[4];
+        const double C = par[5];
+
+        const double dx = x - mu;
+        const double u = dx / sig;
+
+        const double G = std::exp(-0.5 * u * u);
+        const double H = 0.5 * (1.0 - fast_erf(u * inv_sqrt2));
+
+        f = p0 - p1 * x + N * (G + C * H);
+
+        const double inv_sig = 1.0 / sig;
+        const double inv_sig2 = inv_sig * inv_sig;
+        const double inv_sig3 = inv_sig2 * inv_sig;
+
+        const double dG_dmu = G * dx * inv_sig2;
+        const double dG_dsig = G * dx * dx * inv_sig3;
+
+        const double dH_dmu = inv_sqrt_2pi * G * inv_sig;
+        const double dH_dsig = inv_sqrt_2pi * G * dx * inv_sig2;
+
+        g[0] = 1.0;
+        g[1] = -x;
+        g[2] = N * (dG_dmu + C * dH_dmu);
+        g[3] = N * (dG_dsig + C * dH_dsig);
+        g[4] = G + C * H;
+        g[5] = N * H;
+    }
+
+    static bool is_valid(const std::vector<double> &par) {
+        return par[3] > 0.0;
+    }
+
+    static std::array<double, npar> estimate_par(NDView<double, 1> x,
+                                                 NDView<double, 1> y) {
+        const ssize_t n = y.size();
+
+        const auto max_it = std::max_element(y.begin(), y.end());
+        const ssize_t i_max = std::distance(y.begin(), max_it);
+
+        const double x_range = std::max(x[n - 1] - x[0], 1e-9);
+
+        const ssize_t tail = std::min<ssize_t>(std::max<ssize_t>(n / 10, 2), n);
+
+        double left_mean = 0.0;
+        double right_mean = 0.0;
+
+        for (ssize_t i = 0; i < tail; ++i)
+            left_mean += y[i];
+        left_mean /= tail;
+
+        for (ssize_t i = n - tail; i < n; ++i)
+            right_mean += y[i];
+        right_mean /= tail;
+
+        const double p0 = right_mean;
+        const double p1 = 0.0;
+        const double mu = x[i_max];
+
+        const double N = std::max(*max_it - right_mean, 1e-9);
+
+        // Left plateau excess is roughly N*C.
+        const double C = (left_mean - right_mean) / N;
+
+        double sigma = 0.1 * x_range;
+        const double half = right_mean + 0.5 * ((*max_it) - right_mean);
+
+        for (ssize_t i = i_max; i < n; ++i) {
+            if (y[i] <= half) {
+                const double half_width = std::abs(x[i] - mu);
+                if (half_width > 1e-12)
+                    sigma = half_width / 1.1774100225154747; // sqrt(2*ln(2))
+                break;
+            }
+        }
+
+        sigma = std::max(sigma, 1e-6);
+
+        return {p0, p1, mu, sigma, N, C};
+    }
+
+    static void compute_steps(const std::array<double, npar> &start,
+                              double x_range, double y_range,
+                              double slope_scale,
+                              std::array<double, npar> &steps) {
+        steps[0] = std::max(0.1 * std::abs(start[0]), 0.1 * y_range);
+        steps[1] = 0.1 * slope_scale;
+        steps[2] = 0.05 * x_range;
+        steps[3] = 0.05 * x_range;
+        steps[4] = std::max(0.1 * std::abs(start[4]), 0.1 * y_range);
+        steps[5] = std::max(0.1 * std::abs(start[5]), 0.01);
+    }
+};
+
+// _____________________________________________________________________
+//
+// GaussianChargeSharingKb
+// _____________________________________________________________________
+
+/**
+ * @brief Double-Gaussian Kα/Kβ spectrum with erfc charge-sharing tails.
+ *
+ * Equivalent to the old ROOT energyCalibration function:
+ *
+ *   gaussChargeSharingKb(x, par)
+ *
+ * for sign = +1:
+ *
+ *   mu_b = r * mu
+ *
+ *   f(x) = p0 - p1*x
+ *        + N * [
+ *              G(x; mu, sigma)
+ *            + C * H(x; mu, sigma)
+ *            + q * G(x; mu_b, sigma)
+ *            + q * C * H(x; mu_b, sigma)
+ *          ]
+ *
+ * Parameters:
+ *   par[0] = p0      background pedestal
+ *   par[1] = p1      background slope, ROOT convention uses p0 - p1*x
+ *   par[2] = mu      Kα peak position
+ *   par[3] = sigma   shared width / noise RMS
+ *   par[4] = N       global amplitude / number of photons
+ *   par[5] = C       charge-sharing pedestal
+ *   par[6] = r       Kβ mean ratio, mu_Kb = r * mu
+ *   par[7] = q       Kβ fraction
+ */
+struct GaussianChargeSharingKb {
+    static constexpr std::size_t npar = 8;
+
+    static constexpr std::array<ParamInfo, npar> param_info = {{
+        {"p0", -no_bound, no_bound},
+        {"p1", -no_bound, no_bound},
+        {"mu", -no_bound, no_bound},
+        {"sigma", 1e-12, no_bound},
+        {"N", -no_bound, no_bound},
+        {"C", -no_bound, no_bound},
+        {"kb_mean", -no_bound, no_bound},
+        {"kb_frac", -no_bound, no_bound},
+    }};
+
+    static double eval(double x, const std::vector<double> &par) {
+        const double p0 = par[0];
+        const double p1 = par[1];
+        const double mu = par[2];
+        const double sig = par[3];
+        const double N = par[4];
+        const double C = par[5];
+        const double r = par[6];
+        const double q = par[7];
+
+        const double mu_b = r * mu;
+
+        const double dx = x - mu;
+        const double dxb = x - mu_b;
+
+        const double u = dx / sig;
+        const double ub = dxb / sig;
+
+        const double G = std::exp(-0.5 * u * u);
+        const double Gb = std::exp(-0.5 * ub * ub);
+
+        const double H = 0.5 * (1.0 - fast_erf(u * inv_sqrt2));
+        const double Hb = 0.5 * (1.0 - fast_erf(ub * inv_sqrt2));
+
+        const double ka = G + C * H;
+        const double kb = Gb + C * Hb;
+
+        return p0 - p1 * x + N * (ka + q * kb);
+    }
+
+    static void eval_and_grad(double x, const std::vector<double> &par,
+                              double &f, std::array<double, npar> &g) {
+        const double p0 = par[0];
+        const double p1 = par[1];
+        const double mu = par[2];
+        const double sig = par[3];
+        const double N = par[4];
+        const double C = par[5];
+        const double r = par[6];
+        const double q = par[7];
+
+        const double mu_b = r * mu;
+
+        const double dx = x - mu;
+        const double dxb = x - mu_b;
+
+        const double u = dx / sig;
+        const double ub = dxb / sig;
+
+        const double G = std::exp(-0.5 * u * u);
+        const double Gb = std::exp(-0.5 * ub * ub);
+
+        const double H = 0.5 * (1.0 - fast_erf(u * inv_sqrt2));
+        const double Hb = 0.5 * (1.0 - fast_erf(ub * inv_sqrt2));
+
+        const double ka = G + C * H;
+        const double kb = Gb + C * Hb;
+
+        f = p0 - p1 * x + N * (ka + q * kb);
+
+        const double inv_sig = 1.0 / sig;
+        const double inv_sig2 = inv_sig * inv_sig;
+        const double inv_sig3 = inv_sig2 * inv_sig;
+
+        const double dG_dmu = G * dx * inv_sig2;
+        const double dG_dsig = G * dx * dx * inv_sig3;
+        const double dH_dmu = inv_sqrt_2pi * G * inv_sig;
+        const double dH_dsig = inv_sqrt_2pi * G * dx * inv_sig2;
+
+        const double dGb_dmu_b = Gb * dxb * inv_sig2;
+        const double dGb_dsig = Gb * dxb * dxb * inv_sig3;
+        const double dHb_dmu_b = inv_sqrt_2pi * Gb * inv_sig;
+        const double dHb_dsig = inv_sqrt_2pi * Gb * dxb * inv_sig2;
+
+        const double dka_dmu = dG_dmu + C * dH_dmu;
+        const double dka_dsig = dG_dsig + C * dH_dsig;
+
+        const double dkb_dmu_b = dGb_dmu_b + C * dHb_dmu_b;
+        const double dkb_dsig = dGb_dsig + C * dHb_dsig;
+
+        g[0] = 1.0;
+        g[1] = -x;
+
+        // mu_b = r * mu, so dmu_b/dmu = r
+        g[2] = N * (dka_dmu + q * r * dkb_dmu_b);
+
+        g[3] = N * (dka_dsig + q * dkb_dsig);
+
+        g[4] = ka + q * kb;
+
+        g[5] = N * (H + q * Hb);
+
+        // mu_b = r * mu, so dmu_b/dr = mu
+        g[6] = N * q * mu * dkb_dmu_b;
+
+        g[7] = N * kb;
+    }
+
+    static bool is_valid(const std::vector<double> &par) {
+        return par[3] > 0.0;
+    }
+
+    static std::array<double, npar> estimate_par(NDView<double, 1> x,
+                                                 NDView<double, 1> y) {
+        const auto base = GaussianChargeSharing::estimate_par(x, y);
+
+        const double p0 = base[0];
+        const double p1 = base[1];
+        const double mu = base[2];
+        const double sigma = base[3];
+        const double N = base[4];
+        const double C = base[5];
+
+        // Good generic starting values for Cu Kβ relative to Kα.
+        const double kb_mean = 1.10;
+        const double kb_frac = 0.10;
+
+        return {p0, p1, mu, sigma, N, C, kb_mean, kb_frac};
+    }
+
+    static void compute_steps(const std::array<double, npar> &start,
+                              double x_range, double y_range,
+                              double slope_scale,
+                              std::array<double, npar> &steps) {
+        steps[0] = std::max(0.1 * std::abs(start[0]), 0.1 * y_range);
+        steps[1] = 0.1 * slope_scale;
+        steps[2] = 0.05 * x_range;
+        steps[3] = 0.05 * x_range;
+        steps[4] = std::max(0.1 * std::abs(start[4]), 0.1 * y_range);
+        steps[5] = std::max(0.1 * std::abs(start[5]), 0.01);
+        steps[6] = std::max(0.01 * std::abs(start[6]), 1e-3);
+        steps[7] = std::max(0.1 * std::abs(start[7]), 0.01);
+    }
+};
+
 // _____________________________________________________________________
 //
 // RisingScurve

include/aare/VarClusterFinder.hpp

@@ -45,7 +45,11 @@ template <typename T> class VarClusterFinder {
         {-1, -1, 0, 1}}; // col ### 8-neighbour by scaning from top to bottom
     const std::array<int, 8> di_{{0, 0, -1, 1, -1, 1, -1, 1}}; // row
     const std::array<int, 8> dj_{{-1, 1, 0, 0, 1, -1, -1, 1}}; // col
-    std::map<int, int> child; // heirachy: key: child; val: parent
+    std::map<int, int> child;   // heirachy: key: child; val: parent
+    int numberOfNeighbours = 8; // 4 or 8
+    bool empty_surroundingPixels =
+        true; // whether to set peripheral pixels to 0, to avoid potential
+              // influence for pedestal updating
     std::unordered_map<int, Hit> h_size;
     std::vector<Hit> hits;
     // std::vector<std::vector<int16_t>> row
@@ -67,6 +71,16 @@ template <typename T> class VarClusterFinder {
     void set_peripheralThresholdFactor(int factor) {
         peripheralThresholdFactor_ = factor;
     }
+    void set_numberOfNeighbours(int num) {
+        if (num == 4 || num == 8)
+            numberOfNeighbours = num;
+        else
+            throw std::runtime_error(
+                LOCATION + "number of neighbours should be either 4 or 8");
+    }
+    void set_empty_surroundingPixels(bool empty) {
+        empty_surroundingPixels = empty;
+    }
     void find_clusters(NDView<T, 2> img);
     void find_clusters_X(NDView<T, 2> img);
     void rec_FillHit(int clusterIndex, int i, int j);
@@ -201,7 +215,7 @@ void VarClusterFinder<T>::rec_FillHit(int clusterIndex, int i, int j) {
     }
     original_(i, j) = 0;
 
-    for (int k = 0; k < 8; ++k) { // 8 for 8-neighbour
+    for (int k = 0; k < numberOfNeighbours; ++k) { //
         const auto row = i + di_[k];
         const auto col = j + dj_[k];
         if (row >= 0 && col >= 0 && row < shape_[0] && col < shape_[1]) {
@@ -218,9 +232,11 @@ void VarClusterFinder<T>::rec_FillHit(int clusterIndex, int i, int j) {
                 //     h_size[clusterIndex].enes[h_size[clusterIndex].size] =
                 //     original_(row, col);
                 // }// ? weather to include peripheral pixels
-                original_(row, col) =
-                    0; // remove peripheral pixels, to avoid potential influence
-                       // for pedestal updating
+                if (empty_surroundingPixels) {
+                    original_(row, col) =
+                        0; // remove peripheral pixels, to avoid potential
+                           // influence for pedestal updating
+                }
             }
         }
     }

python/aare/__init__.py

@@ -17,7 +17,7 @@ from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, Clu
 from .ClusterVector import ClusterVector
 from .Cluster import Cluster
 
-from ._aare import Gaussian, RisingScurve, FallingScurve, Pol1, Pol2
+from ._aare import Gaussian, RisingScurve, FallingScurve, Pol1, Pol2, GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb
 from ._aare import fit
 from ._aare import fit_gaus, fit_pol1, fit_scurve, fit_scurve2
 from ._aare import Interpolator

python/src/fit.hpp

@@ -679,6 +679,11 @@ void define_fit_bindings(py::module &m) {
 
     // ── Bind model classes ──────────────────────────────────────────
     bind_fit_model<aare::model::Gaussian>(m, "Gaussian");
+    bind_fit_model<aare::model::GaussianErfcPlateau>(m, "GaussianErfcPlateau");
+    bind_fit_model<aare::model::GaussianChargeSharing>(m,
+                                                       "GaussianChargeSharing");
+    bind_fit_model<aare::model::GaussianChargeSharingKb>(
+        m, "GaussianChargeSharingKb");
     bind_fit_model<aare::model::RisingScurve>(m, "RisingScurve");
     bind_fit_model<aare::model::FallingScurve>(m, "FallingScurve");
     bind_fit_model<aare::model::Pol1>(m, "Pol1");
@@ -716,6 +721,38 @@ void define_fit_bindings(py::module &m) {
                     mdl, x, y, y_err_obj, n_threads);
             }
 
+            // ── GaussianErfcPlateau ───────
+            if (py::isinstance<aare::FitModel<GaussianErfcPlateau>>(
+                    model_obj)) {
+                const auto &mdl =
+                    model_obj
+                        .cast<const aare::FitModel<GaussianErfcPlateau> &>();
+                return fit_dispatch<GaussianErfcPlateau,
+                                    Chi2GaussianErfcPlateau>(
+                    mdl, x, y, y_err_obj, n_threads);
+            }
+
+            // ── GaussianChargeSharing ───────
+            if (py::isinstance<aare::FitModel<GaussianChargeSharing>>(
+                    model_obj)) {
+                const auto &mdl =
+                    model_obj
+                        .cast<const aare::FitModel<GaussianChargeSharing> &>();
+                return fit_dispatch<GaussianChargeSharing,
+                                    Chi2GaussianChargeSharing>(
+                    mdl, x, y, y_err_obj, n_threads);
+            }
+
+            // ── GaussianChargeSharingKb ───────
+            if (py::isinstance<aare::FitModel<GaussianChargeSharingKb>>(
+                    model_obj)) {
+                const auto &mdl = model_obj.cast<
+                    const aare::FitModel<GaussianChargeSharingKb> &>();
+                return fit_dispatch<GaussianChargeSharingKb,
+                                    Chi2GaussianChargeSharingKb>(
+                    mdl, x, y, y_err_obj, n_threads);
+            }
+
             // ── Rising Scurve ───────
             if (py::isinstance<aare::FitModel<RisingScurve>>(model_obj)) {
                 const auto &mdl =

python/src/var_cluster.hpp

@@ -34,6 +34,10 @@ void define_var_cluster_finder_bindings(py::module &m) {
                  auto noise_map_span = make_view_2d(noise_map);
                  self.set_noiseMap(noise_map_span);
              })
+        .def("set_numberOfNeighbours",
+             &VarClusterFinder<double>::set_numberOfNeighbours)
+        .def("set_empty_surroundingPixels",
+             &VarClusterFinder<double>::set_empty_surroundingPixels)
         .def("set_peripheralThresholdFactor",
              &VarClusterFinder<double>::set_peripheralThresholdFactor)
         .def("find_clusters",