diff --git a/src/theory.jl b/src/theory.jl
index 6eb577e..5988772 100755
--- a/src/theory.jl
+++ b/src/theory.jl
@@ -25,7 +25,8 @@ function Theory(cosmology::Cosmology,
             dz = get(Nuisances, string(name, "_", "dz"), 0.0)
             tracer = NumberCountsTracer(cosmology, zs .+ dz, nz;
                                         b=b, res=res_gc, 
-                                        nz_interpolation=int_gc)
+                                        nz_interpolation=int_gc,
+                                        smooth=true)
         elseif t_type == "galaxy_shear"
             zs_mean, nz_mean = files[string("nz_", name)]
             m = get(Nuisances, string(name, "_", "m"), 0.0)
@@ -57,7 +58,9 @@ function Theory(cosmology::Cosmology,
         ls = files[string("ls_", name1, "_", name2)]
         tracer1 = tracers[name1]
         tracer2 = tracers[name2]
-        cls[idx[i]+1:idx[i+1]] = angularCℓs(cosmology, tracer1, tracer2, ls)
+        lss = [[i for i in l-3:l+3] for l in ls]
+        clss = [mean(angularCℓs(cosmology, tracer1, tracer2, _lls)) for _lls in lss]
+        cls[idx[i]+1:idx[i+1]] = clss
     end
     
     return cls
diff --git a/src/tracers.jl b/src/tracers.jl
index ca3acb3..fadf547 100755
--- a/src/tracers.jl
+++ b/src/tracers.jl
@@ -28,7 +28,7 @@ Returns:
 - `NumberCountsTracer::NumberCountsTracer` : Number counts tracer structure.
 """
 NumberCountsTracer(cosmo::Cosmology, z_n, nz; 
-    b=1.0, res=1000, nz_interpolation="linear") = begin
+    b=1.0, res=1000, nz_interpolation="linear", smooth=false) = begin
     z_w, nz_w = nz_interpolate(z_n, nz, res;
         mode=nz_interpolation)
     nz_norm = integrate(z_w, nz_w, SimpsonEven())
@@ -37,6 +37,9 @@ NumberCountsTracer(cosmo::Cosmology, z_n, nz;
     hz = Hmpc(cosmo, z_w)
     
     w_arr = @. (nz_w*hz/nz_norm)
+    if smooth
+        w_arr = smooth_w_neighbors(w_arr)
+    end
     wint = linear_interpolation(chi, b .* w_arr, extrapolation_bc=Line())
     F::Function = ℓ -> 1
     NumberCountsTracer(wint, F)
@@ -175,3 +178,21 @@ function nz_interpolate(z, nz, res; mode="linear")
         return z, nz
     end
 end
+
+function smooth_w_neighbors(arr::Vector{T}, k::Int = 5) where T
+    N = length(arr)
+    neighbors = Vector{}(undef, N)  # Create an array to hold neighbor arrays
+
+    half_k = div(k, 2)  # This is 2 if k=5, so we look at two neighbors on each side
+
+    for i in 1:N
+        # Define the range around the current element, clamping to avoid out-of-bounds
+        start_idx = max(1, i - half_k)
+        end_idx = min(N, i + half_k)
+        
+        # Collect neighbors and assign to the current position
+        neighbors[i] = mean(arr[start_idx:end_idx])
+    end
+
+    return neighbors
+end
diff --git a/test/test_output.npz b/test/test_output.npz
index 845d1b6..91df7d1 100644
Binary files a/test/test_output.npz and b/test/test_output.npz differ