diff --git a/scripts/muonscanfit.jl b/scripts/muonscanfit.jl
index 74b217f94d32031b47894b669a91fb7d1e6fb1bc..d6b309d30e0b4a6674747e045e6ee6222e7205eb 100644
--- a/scripts/muonscanfit.jl
+++ b/scripts/muonscanfit.jl
@@ -37,6 +37,8 @@ struct MuonScanfitResult
     dir_z::Float64
     t::Float64
     Q::Float64
+    S1::Float64
+    S2::Float64
     NDF::Int
     walltime::Float64
     mc_angular_error::Float64
@@ -93,6 +95,8 @@ function main()
                 muon.dir...,
                 muon.t,
                 muon.Q,
+                muon.S1,
+                muon.S2,
                 muon.NDF,
                 walltime,
                 mc_angular_error,
diff --git a/src/scanfit.jl b/src/scanfit.jl
index 3eed0ff90693e8b242f28745130b6b891975df34..0bf2612ef428b2c940d6527f79f5ba16beddb754 100644
--- a/src/scanfit.jl
+++ b/src/scanfit.jl
@@ -1,3 +1,6 @@
+const N_FITS_SPREAD = 10  # number of prefits to calculate the spread of the fit
+
+
 Base.@kwdef struct MuonScanfitParameters
     tmaxlocal::Float64 = 18.0  # [ns]
     roadwidth::Float64 = 200.0  # [m]
@@ -59,6 +62,8 @@ function (msf::MuonScanfit)(hits::Vector{T}) where T<:KM3io.AbstractHit
     isempty(candidates) && return candidates
     sort!(candidates, by=m->m.Q; rev=true)
 
+    S1 = spread(candidates[1:min(N_FITS_SPREAD, length(candidates))])
+
     # Second round on directed cones pointing towards the previous best directions
     # TODO: currently disabled until all the allocations are minimised
     # here, reusing a Vector{Direction} (attached to msf as buffer) might be a good idea.
@@ -76,9 +81,12 @@ function (msf::MuonScanfit)(hits::Vector{T}) where T<:KM3io.AbstractHit
         sort!(candidates, by=m->m.Q; rev=true)
     end
 
-    candidates[1:msf.params.nfits]
+    S2 = spread(candidates[1:min(N_FITS_SPREAD, length(candidates))])
+
+    [setproperties(c, (S1=S1, S2=S2)) for c in candidates[1:msf.params.nfits]]
 end
 
+
 """
 
 Performs the scanfit for each given direction and returns a
@@ -105,9 +113,14 @@ struct MuonScanfitCandidate
     pos::Position{Float64}
     dir::Direction{Float64}
     t::Float64
-    Q::Float64
+    Q::Float64  # quality of the fit
+    S1::Float64  # spread of the prefits, the smaller the better
+    S2::Float64  # spread of the last fits, the smaller the better
     NDF::Int
 end
+MuonScanfitCandidate(pos::Position{Float64}, dir::Direction{Float64}, t::Float64, Q::Float64, NDF::Int) = MuonScanfitCandidate(pos, dir, t, Q, π, π, NDF)
+MuonScanfitCandidate() = MuonScanfitCandidate(Position(0.0, 0.0, 0.0), Direction(0.0, 0.0, 1.0), NaN, NaN, NaN, NaN, 0)
+Base.angle(m1::T, m2::T) where T<:MuonScanfitCandidate = angle(m1.dir, m2.dir)
 
 """
 The quality of the fit, the larger the better, as used in e.g. Jpp.
@@ -384,7 +397,7 @@ function (s::XYTSolver)(hits::Vector{T}, dir::Direction{Float64}, α::Float64) w
     hits = s.hits_buffer  # just for convenience
     n_final_hits = length(hits)
 
-    n_final_hits <= s.est.NUMBER_OF_PARAMETERS && return MuonScanfitCandidate(Position(0, 0, 0), dir, 0, -Inf, 0)
+    n_final_hits <= s.est.NUMBER_OF_PARAMETERS && return MuonScanfitCandidate()
 
     NDF = n_final_hits - s.est.NUMBER_OF_PARAMETERS
     N = hitcount(hits)
@@ -394,7 +407,7 @@ function (s::XYTSolver)(hits::Vector{T}, dir::Direction{Float64}, α::Float64) w
         estimate!(s.est, hits)
     catch ex
         # isa(ex, SingularSVDException) && @warn "Singular SVD"
-        return MuonScanfitCandidate(Position(0, 0, 0), dir, 0, -Inf, 0)
+        return MuonScanfitCandidate()
     end
 
     # TODO: consider creating a "pos()" getter for everything
diff --git a/test/math_tests.jl b/test/math_tests.jl
index 8e3559480c47e943718b7c72ed5f3a073005845e..ba53b261acba0f0f6de672fd944f81c9cbf40b03 100644
--- a/test/math_tests.jl
+++ b/test/math_tests.jl
@@ -8,6 +8,9 @@ using Test
     @test π/2 ≈ angle([1,0,0], [0,0,1])
     @test π ≈ angle([1,0,0], [-1,0,0])
 
+    m1 = NeRCA.MuonScanfitCandidate(Position(1, 2, 3), Direction(1.0, 0.0, 0.0), 0, 0, 0, 0)
+    m2 = NeRCA.MuonScanfitCandidate(Position(1, 2, 3), Direction(0.0, 1.0, 0.0), 0, 0, 0, 0)
+    @test π/2 ≈ angle(m1, m2)
 end