Add functions to calculate pvalue

core/CMakeLists.txt

@@ -57,6 +57,8 @@ traccc_add_library( traccc_core core TYPE SHARED
   "include/traccc/finding/finding_config.hpp"
   "include/traccc/finding/interaction_register.hpp"
   # Fitting algorithmic code
+  "include/traccc/fitting/detail/gamma_functions.hpp"
+  "src/fitting/detail/gamma_functions.cpp"
   "include/traccc/fitting/kalman_filter/gain_matrix_smoother.hpp"
   "include/traccc/fitting/kalman_filter/gain_matrix_updater.hpp"
   "include/traccc/fitting/kalman_filter/kalman_actor.hpp"

core/include/traccc/definitions/common.hpp

@@ -11,13 +11,19 @@
 #include "traccc/definitions/primitives.hpp"
 
 // Detray include(s).
+#include "detray/definitions/math.hpp"
 #include "detray/definitions/units.hpp"
 
 namespace traccc {
 
+namespace math_ns = detray::math_ns;
+
 template <typename scalar_t>
 using unit = detray::unit<scalar_t>;
 
+template <typename scalar_t>
+using constant = detray::constant<scalar_t>;
+
 // epsilon for float variables
 constexpr scalar float_epsilon = 1e-5;
 

core/include/traccc/edm/track_state.hpp

@@ -32,6 +32,9 @@ struct fitter_info {
 
     /// Chi square of fitted track
     scalar_type chi2{0};
+
+    /// pvalue
+    scalar_type pval{0};
 };
 
 /// Fitting result per measurement

core/include/traccc/fitting/detail/chi2_cdf.hpp

@@ -0,0 +1,30 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#include "traccc/definitions/qualifiers.hpp"
+#include "traccc/fitting/detail/gamma_functions.hpp"
+
+// Reference: ProbFuncMathCore.cxx of ROOT library
+namespace traccc::detail {
+
+// Funtions to calculate the upper incomplete gamma function from a given chi2
+// and ndf
+//
+// @param x chi square
+// @param r ndof
+// @return upper incomplete gamma function (pvalue)
+template <typename scalar_t>
+TRACCC_HOST_DEVICE scalar_t chisquared_cdf_c(const scalar_t x,
+                                             const scalar_t r) {
+    double retval =
+        igamc(0.5 * static_cast<double>(r), 0.5 * static_cast<double>(x));
+    return static_cast<scalar_t>(retval);
+}
+
+}  // namespace traccc::detail

core/include/traccc/fitting/detail/gamma_functions.hpp

@@ -0,0 +1,63 @@
+/** TRACCC library, part of the ACTS project (R&D line)
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+#include "traccc/definitions/qualifiers.hpp"
+
+// gamma and related functions from Cephes library
+// see:  http://www.netlib.org/cephes
+//
+// Copyright 1985, 1987, 2000 by Stephen L. Moshier
+namespace traccc::detail {
+
+// incomplete gamma function (complement integral)
+//  igamc(a,x)   =   1 - igam(a,x)
+//
+//                            inf.
+//                              -
+//                     1       | |  -t  a-1
+//               =   -----     |   e   t   dt.
+//                    -      | |
+//                   | (a)    -
+//                             x
+//
+//
+
+// In this implementation both arguments must be positive.
+// The integral is evaluated by either a power series or
+// continued fraction expansion, depending on the relative
+// values of a and x.
+TRACCC_HOST_DEVICE double igamc(const double a, const double x);
+
+// left tail of incomplete gamma function:
+//
+//          inf.      k
+//   a  -x   -       x
+//  x  e     >   ----------
+//           -     -
+//          k=0   | (a+k+1)
+//
+TRACCC_HOST_DEVICE double igam(const double a, const double x);
+
+TRACCC_HOST_DEVICE double lgam(const double x);
+
+/*
+ * calculates a value of a polynomial of the form:
+ * a[0]x^N+a[1]x^(N-1) + ... + a[N]
+ */
+TRACCC_HOST_DEVICE double Polynomialeval(const double x, const double* a,
+                                         const unsigned int N);
+
+/*
+ * calculates a value of a polynomial of the form:
+ * x^N+a[0]x^(N-1) + ... + a[N-1]
+ */
+TRACCC_HOST_DEVICE double Polynomial1eval(const double x, const double* a,
+                                          const unsigned int N);
+
+}  // namespace traccc::detail

core/include/traccc/fitting/fitting_config.hpp

@@ -16,7 +16,7 @@ namespace traccc {
 template <typename scalar_t>
 struct fitting_config {
 
-    std::size_t n_iterations = 1;
+    std::size_t max_iterations = 10;
     scalar_t pathlimit = std::numeric_limits<scalar_t>::max();
     scalar_t overstep_tolerance = -10 * detray::unit<scalar_t>::um;
     scalar_t step_constraint = std::numeric_limits<scalar_t>::max();

core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp

@@ -12,6 +12,7 @@
 #include "traccc/edm/track_candidate.hpp"
 #include "traccc/edm/track_parameters.hpp"
 #include "traccc/edm/track_state.hpp"
+#include "traccc/fitting/detail/chi2_cdf.hpp"
 #include "traccc/fitting/fitting_config.hpp"
 #include "traccc/fitting/kalman_filter/gain_matrix_smoother.hpp"
 #include "traccc/fitting/kalman_filter/kalman_actor.hpp"
@@ -127,25 +128,64 @@ class kalman_fitter {
         const seed_parameters_t& seed_params, state& fitter_state,
         vector_type<intersection_type>&& nav_candidates = {}) {
 
+        scalar_type pval = 0.f;
+
         // Run the kalman filtering for a given number of iterations
-        for (std::size_t i = 0; i < m_cfg.n_iterations; i++) {
+        for (std::size_t i_it = 0; i_it < m_cfg.max_iterations; i_it++) {
 
             // Reset the iterator of kalman actor
             fitter_state.m_fit_actor_state.reset();
 
-            if (i == 0) {
+            if (i_it == 0) {
                 filter(seed_params, fitter_state, std::move(nav_candidates));
+                pval = fitter_state.m_fit_info.pval;
             }
             // From the second iteration, seed parameter is the smoothed track
             // parameter at the first surface
-            else {
-                const auto& new_seed_params =
+            else if (i_it > 0) {
+
+                auto& smoothed =
+                    fitter_state.m_fit_actor_state.m_track_states[0].smoothed();
+
+                // Get intersection on surface
+                const detray::surface<detector_type> sf{
+                    m_detector, smoothed.surface_link()};
+                using cxt_t = typename detector_type::geometry_context;
+                const cxt_t ctx{};
+                const auto free_vec =
+                    sf.bound_to_free_vector(ctx, smoothed.vector());
+
+                intersection_type sfi;
+                sfi.sf_desc = m_detector.surface(smoothed.surface_link());
+                sf.template visit_mask<detray::intersection_update>(
+                    detray::detail::ray<transform3_type>(free_vec), sfi,
+                    m_detector.transform_store());
+
+                // Apply material interaction backwardly to track state
+                typename interactor::state interactor_state;
+                interactor_state.do_multiple_scattering = false;
+                interactor{}.update(
+                    smoothed, interactor_state,
+                    static_cast<int>(detray::navigation::direction::e_backward),
+                    sf, sfi.cos_incidence_angle);
+
+                // Make new seed parameter
+                auto new_seed_params =
                     fitter_state.m_fit_actor_state.m_track_states[0].smoothed();
 
                 filter(new_seed_params, fitter_state,
                        std::move(nav_candidates));
+
+                const auto dpval = pval - fitter_state.m_fit_info.pval;
+                pval = fitter_state.m_fit_info.pval;
+                if (std::abs(dpval) < 1e-2) {
+                    std::cout << "N iterations: " << i_it + 1 << std::endl;
+                    return;
+                }
             }
         }
+        std::cout << "N iterations (FULL): " << m_cfg.max_iterations
+                  << std::endl;
     }
 
     /// Run the kalman fitter for an iteration
@@ -226,6 +266,12 @@ class kalman_fitter {
 
     TRACCC_HOST_DEVICE
     void update_statistics(state& fitter_state) {
+
+        // Reset the ndf and chi2 of fitter info
+        fitter_state.m_fit_info.ndf = 0.f;
+        fitter_state.m_fit_info.chi2 = 0.f;
+        fitter_state.m_fit_info.pval = 0.f;
+
         auto& fit_info = fitter_state.m_fit_info;
         auto& track_states = fitter_state.m_fit_actor_state.m_track_states;
 
@@ -242,6 +288,9 @@ class kalman_fitter {
 
         // Subtract the NDoF with the degree of freedom of the bound track (=5)
         fit_info.ndf = fit_info.ndf - 5.f;
+
+        // p value
+        fit_info.pval = detail::chisquared_cdf_c(fit_info.chi2, fit_info.ndf);
     }
 
     private:

core/include/traccc/fitting/kalman_filter/statistics_updater.hpp

@@ -10,6 +10,7 @@
 // Project include(s).
 #include "traccc/definitions/qualifiers.hpp"
 #include "traccc/edm/track_state.hpp"
+#include "traccc/fitting/detail/chi2_cdf.hpp"
 
 namespace traccc {