Muonscanfit

src/scanfit.jl

@@ -2,6 +2,7 @@ Base.@kwdef struct MuonScanfitParameters
     tmaxlocal::Float64 = 18.0  # [ns]
     roadwidth::Float64 = 200.0  # [m]
     nmaxhits::Int = 50  # maximum number of hits to use
+    number_of_fits_to_keep::Int = 1
 end
 
 struct MuonScanfit
@@ -14,6 +15,9 @@ function Base.show(io::IO, m::MuonScanfit)
    print(io, "$(typeof(m)) in $(length(m.directions)) directions")
 end
 
+"""
+Performs a Muon track fit for a given set of hits (usually snapshot hits).
+"""
 function (msf::MuonScanfit)(hits::Vector{T}) where T<:KM3io.AbstractHit
     l1builder = L1Builder(L1BuilderParameters(msf.params.tmaxlocal, false))
     rhits = l1builder(HitR1, msf.detector, hits)
@@ -24,9 +28,15 @@ function (msf::MuonScanfit)(hits::Vector{T}) where T<:KM3io.AbstractHit
     clique = Clique(Match3B(msf.params.roadwidth, msf.params.tmaxlocal))
     clusterize!(rhits, clique)
 
+    # candidates_pool = let candidates_dict = Dict{Int, Vector{MuonScanFitResult}}()
+    #     n_candidates_per_thread = Int(ceil((length(msf.directions) / Threads.nthreads())))
+    #     for thread_id ∈ 1:Threads.nthreads()
+    #         candidates_dict[thread_id] = sizehint!(MuonScanFitResult[], n_candidates_per_thread)
+    #     end
+    #     candidates_dict
+    # end
     candidates = MuonScanFitResult[]
 
-    # TODO threading is bugged
     # Threads.@threads for dir ∈ msf.directions
     for dir ∈ msf.directions
         est = Line1ZEstimator(Line1Z(Position(0, 0, 0), 0))
@@ -50,12 +60,10 @@ function (msf::MuonScanfit)(hits::Vector{T}) where T<:KM3io.AbstractHit
 
         clusterize!(rhits_copy, clique1D)
 
-        length(rhits_copy) < est.NUMBER_OF_PARAMETERS && continue
         NDF = length(rhits_copy) - est.NUMBER_OF_PARAMETERS
         N = hitcount(rhits_copy)
-        if length(rhits_copy) <= est.NUMBER_OF_PARAMETERS
-            continue
-        end
+
+        length(rhits_copy) <= est.NUMBER_OF_PARAMETERS && continue
 
         sort!(rhits_copy)
 
@@ -63,6 +71,7 @@ function (msf::MuonScanfit)(hits::Vector{T}) where T<:KM3io.AbstractHit
             estimate!(est, rhits_copy)
         catch ex
             # if isa(ex, SingularSVDException)
+            # @warn "Singular SVD"
             continue
         end
 
@@ -70,28 +79,32 @@ function (msf::MuonScanfit)(hits::Vector{T}) where T<:KM3io.AbstractHit
         # TODO: pass alpha and sigma, like V.set(*this, data.begin(), __end1, gridAngle_deg, sigma_ns);  // JMatrixNZ
         V = covmatrix(est.model.pos, rhits_copy)
         Y = timeresvec(est.model, rhits_copy)
-
         V⁻¹ = inv(V)
-
         χ² = transpose(Y) * V⁻¹ * Y
-
-
         fit_pos = R \ est.model.pos
 
-        push!(candidates, MuonScanFitResult(fit_pos, dir, quality(χ², N, NDF), NDF))
+        # push!(candidates_pool[Threads.threadid()], MuonScanFitResult(fit_pos, dir, est.model.t, quality(χ², N, NDF), NDF))
+        push!(candidates, MuonScanFitResult(fit_pos, dir, est.model.t, quality(χ², N, NDF), NDF))
     end
-    sort!(candidates, by=m->m.Q; rev=true)
-    candidates
+    # candidates = vcat(values(candidates_pool)...)
+    isempty(candidates) && return candidates
+
+    # sort!(candidates, by=m->m.Q; rev=true)
+    candidates[1:msf.params.number_of_fits_to_keep]
 end
 
 struct MuonScanFitResult
     pos::Position
     dir::Direction
+    t::Float64
     Q::Float64
     NDF::Int
 end
 
-quality(χ², N, NDF) = N / 4.0 * χ² / NDF
+"""
+The quality of the fit, the larger the better, as used in e.g. Jpp.
+"""
+quality(χ², N, NDF) = N  -  0.25 * χ² / NDF
 
 abstract type EstimatorModel end
 
@@ -211,8 +224,8 @@ function estimate!(est::Line1ZEstimator, hits)
 
     est.model = Line1Z(
         Position(
-            est.V[1, 1] * y₀ + est.V[1, 2] * y₁ + est.V[1, 3] * y₂,
-            est.V[2, 1] * y₀ + est.V[2, 2] * y₁ + est.V[2, 3] * y₂,
+            pos.x + est.V[1, 1] * y₀ + est.V[1, 2] * y₁ + est.V[1, 3] * y₂,
+            pos.y + est.V[2, 1] * y₀ + est.V[2, 2] * y₁ + est.V[2, 3] * y₂,
             posz(lz)
         ),
         (est.V[3, 1] * y₀ + est.V[3, 2] * y₁ + est.V[3, 3] * y₂) * KM3io.Constants.KAPPA_WATER * KM3io.Constants.C_INVERSE + t₀