take into account elastic scattering in parameter estimation

include/aare/Models.hpp

@@ -600,8 +600,17 @@ struct GaussianChargeSharing {
         return par[3] > 0.0;
     }
 
-    static std::array<double, npar> estimate_par(NDView<double, 1> x,
-                                                 NDView<double, 1> y) {
+    /**
+     * @brief Data-driven initial parameter estimates.
+     * @param x Independent variable values
+     * @param y Dependent variable values
+     * @param elastic_scattering If true, take into account elastic scattering
+     * effects
+     * @return Initial parameter estimates
+     */
+    static std::array<double, npar>
+    estimate_par(NDView<double, 1> x, NDView<double, 1> y,
+                 const bool elastic_scattering = true) {
         const ssize_t n = y.size();
 
         const auto max_it = std::max_element(y.begin(), y.end());
@@ -609,7 +618,9 @@ struct GaussianChargeSharing {
 
         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);
+        const ssize_t tail = std::min<ssize_t>(
+            std::max<ssize_t>(n / 10, 2),
+            n); // TODO: why hardcoded to 10 % - doesnt it depend on x range??
 
         double left_mean = 0.0;
         double right_mean = 0.0;
@@ -618,18 +629,22 @@ struct GaussianChargeSharing {
             left_mean += y[i];
         left_mean /= tail;
 
-        for (ssize_t i = n - tail; i < n; ++i)
-            right_mean += y[i];
-        right_mean /= tail;
+        if (elastic_scattering) {
+            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);
+        const double mu = x[i_max];
+        const double N =
+            std::max(*max_it - right_mean, 1e-9); // amplitude K_\alpha
 
         // Left plateau excess is roughly N*C.
-        const double C = (left_mean - right_mean) / N;
+        const double C =
+            (left_mean - right_mean) / N; // charge-sharing amplitude ratio
 
         double sigma = 0.1 * x_range;
         const double half = right_mean + 0.5 * ((*max_it) - right_mean);
@@ -811,9 +826,16 @@ struct GaussianChargeSharingKb {
         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);
+    /**
+     * @brief Data-driven initial parameter estimates.
+     * @param elastic_scattering If true, take into account elastic scattering
+     * effects
+     */
+    static std::array<double, npar>
+    estimate_par(NDView<double, 1> x, NDView<double, 1> y,
+                 const bool elastic_scattering = true) {
+        const auto base =
+            GaussianChargeSharing::estimate_par(x, y, elastic_scattering);
 
         const double p0 = base[0];
         const double p1 = base[1];

python/src/fit.hpp

@@ -21,10 +21,10 @@ fit_dispatch(const aare::FitModel<Model> &model,
 
 template <typename Model> void bind_fit_model(py::module &m, const char *name) {
     using FM = aare::FitModel<Model>;
-    py::class_<FM>(m, name)
-        .def(py::init<unsigned int, unsigned int, double, bool>(),
-             py::arg("strategy") = 0, py::arg("max_calls") = 100,
-             py::arg("tolerance") = 0.5, py::arg("compute_errors") = false)
+    auto cls = py::class_<FM>(m, name);
+    cls.def(py::init<unsigned int, unsigned int, double, bool>(),
+            py::arg("strategy") = 0, py::arg("max_calls") = 100,
+            py::arg("tolerance") = 0.5, py::arg("compute_errors") = false)
         .def("SetParLimits",
              py::overload_cast<unsigned int, double, double>(&FM::SetParLimits),
              py::arg("idx"), py::arg("lo"), py::arg("hi"))
@@ -101,6 +101,45 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
             )doc",
             py::arg("x"), py::arg("y"), py::arg("y_err") = py::none(),
             py::arg("n_threads") = 4);
+
+    if constexpr (std::is_same_v<Model, model::GaussianChargeSharing> or
+                  std::is_same_v<Model, model::GaussianChargeSharingKb>) {
+        cls.def(
+            "estimate_par",
+            [](const FM & /*self*/,
+               py::array_t<double, py::array::c_style | py::array::forcecast> x,
+               py::array_t<double, py::array::c_style | py::array::forcecast> y,
+               const bool elastic_scattering) {
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_1d(y);
+
+                return Model::estimate_par(x_view, y_view, elastic_scattering);
+            },
+            py::arg("x"), py::arg("y"), py::arg("elastic_scattering") = true,
+            R"doc(
+                 Estimate starting parameters from 1D data.
+                    Parameters
+                    ----------
+                    elastic_scattering : bool
+                        If true, take into account elastic scattering
+                        effects.
+                 )doc");
+    } else {
+        cls.def(
+            "estimate_par",
+            [](const FM & /*self*/,
+               py::array_t<double, py::array::c_style | py::array::forcecast> x,
+               py::array_t<double, py::array::c_style | py::array::forcecast>
+                   y) {
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_1d(y);
+
+                return Model::estimate_par(x_view, y_view);
+            },
+            py::arg("x"), py::arg("y"), R"doc(
+             Estimate starting parameters from 1D data.
+             )doc");
+    }
 }
 
 template <typename Model>