Filter in Fourier space does not behave like it

This is a follow-up to an answered question that I asked and that can be found here.

I have several points (x,y,z coordinates) in a 3D box with associated masses. I want to draw an histogram of the mass-density that is found in spheres of a given radius R. The idea is to compute a 3D histogram of my box (with binning much smaller than the radius), take its FFT, multiply by the filter (a ball in real space) and inverse FFT the result. From there, I just compute the 1D histogram of the values obtained in each 3D-bin.

Following the issue I had by using an analytic expression of the filter in Fourier space, I am now generating the ball in real space and take its FFT to obtain my filter. However, the histogram I get out of this method is really strange, where I would expect a Gaussian I am getting this:

My code is the following:

import numpy as np
import matplotlib.pyplot as plt
import random
from numba import njit


# 1. Generate a bunch of points with masses from 1 to 3 separated by a radius of 1 cm

size = 100

radius = 1
rangeX = (0, size)
rangeY = (0, size)
rangeZ = (0, size)
rangem = (1,3)
qty = 300000  # or however many points you want


deltas = set()
for x in range(-radius, radius+1):
    for y in range(-radius, radius+1):
        for z in range(-radius, radius+1):
            if x*x + y*y + z*z<= radius*radius:
                deltas.add((x,y,z))

X = []
Y = []
Z = []
M = []
excluded = set()
for i in range(qty):
    x = random.randrange(*rangeX)
    y = random.randrange(*rangeY)
    z = random.randrange(*rangeZ)
    m = random.uniform(*rangem)
    if (x,y,z) in excluded: continue
    X.append(x)
    Y.append(y)
    Z.append(z)
    M.append(1)
    excluded.update((x+dx, y+dy, z+dz) for (dx,dy,dz) in deltas)

#print("There is ",len(X)," points in the box")



# Compute the 3D histogram

a = np.vstack((X, Y, Z)).T

b = 200
R = 10

H, edges = np.histogramdd(a, weights=M, bins = b)  

Fh = np.fft.fftn(H, axes=(-3,-2, -1))

# Generate the filter in real space

Kreal = np.zeros((b,b,b))


X = edges[0]
Y = edges[1]
Z = edges[2]

mid = int(b/2)

s = (X.max()-X.min()+Y.max()-Y.min()+Z.max()-Z.min())/(3*b)

cst = 1/2 + (1/12 - (R/s)**2)*np.arctan((0.5*np.sqrt((R/s)**2-0.5))/(0.5-(R/s)**2)) + 1/3*np.sqrt((R/s)**2-0.5) + ((R/s)**2 - 1/12)*np.arctan(0.5/(np.sqrt((R/s)**2-0.5))) - 4/3*(R/s)**3*np.arctan(0.25/((R/s)*np.sqrt((R/s)**2-0.5)))


@njit(parallel=True)
def remp(Kreal):
    for i in range(b):
        for j in range(b):
            for k in range(b):
                a = cst - np.sqrt((X[i]-X[mid])**2 + (Y[j]-Y[mid])**2 + (Z[k]-Z[mid])**2)/s
                if a >= 0.1 and a < 0.2:
                    Kreal[i][j][k] = 0.1
                elif a >= 0.2 and a < 0.3:
                   Kreal[i][j][k] = 0.2 
                elif a >= 0.3 and a < 0.4:
                   Kreal[i][j][k] = 0.3
                elif a >= 0.4 and a < 0.5:
                   Kreal[i][j][k] = 0.4
                elif a >= 0.5 and a < 0.6:
                   Kreal[i][j][k] = 0.5
                elif a >= 0.6 and a < 0.7:
                   Kreal[i][j][k] = 0.6
                elif a >= 0.7 and a < 0.8:
                   Kreal[i][j][k] = 0.7
                elif a >= 0.8 and a < 0.9:
                   Kreal[i][j][k] = 0.8
                elif a >= 0.9 and a < 0.99:
                   Kreal[i][j][k] = 0.9
                elif a >= 0.99:
                   Kreal[i][j][k] = 1
    return Kreal


Kreal = remp(Kreal)

Kreal = np.fft.ifftshift(Kreal)

Kh = np.fft.fftn(Kreal, axes=(-3,-2, -1))

Gh = np.multiply(Fh, Kh)

Density = np.real(np.fft.ifftn(Gh,axes=(-3,-2, -1)))

# Generate the filter in fourier space using its analytic expression

kx = 2*np.pi*np.fft.fftfreq(len(edges[0][:-1]))*len(edges[0][:-1])/(np.amax(X)-np.amin(X))
ky = 2*np.pi*np.fft.fftfreq(len(edges[1][:-1]))*len(edges[1][:-1])/(np.amax(Y)-np.amin(Y))
kz = 2*np.pi*np.fft.fftfreq(len(edges[2][:-1]))*len(edges[2][:-1])/(np.amax(Z)-np.amin(Z))

kr = np.sqrt(kx[:,None,None]**2 + ky[None,:,None]**2 + kz[None,None,:]**2)
kr *= R
Kh = (np.sin(kr)-kr*np.cos(kr))*3/(kr)**3


Kh[0,0,0] = 1


Gh = np.multiply(Fh, Kh)

Density2 = np.real(np.fft.ifftn(Gh,axes=(-3,-2, -1)))


D = Density.flatten()
N = np.mean(D)

D2 = Density2.flatten()
N2 = np.mean(D2)


# I then compute the histogram I want

hist, bins = np.histogram(D/N, bins='auto', density=True)
bin_centers = (bins[1:]+bins[:-1])*0.5
plt.plot(bin_centers, hist,'.',label = "Defining the Filter in real space")


hist, bins = np.histogram(D2/N2, bins='auto', density=True)
bin_centers = (bins[1:]+bins[:-1])*0.5
plt.plot(bin_centers, hist,'.',label = "Using analytic expression")



plt.xlabel('Normalised Density')
plt.ylabel('Probability density')
plt.legend()
plt.show()

Do you understand why this happens ? Thank you very much for your help.

PS: the long list of if statements when I define the Filter in real space comes from how i'm drawing the sphere on the grid. I assign the value 1 to all the bins that are 100% within the sphere, and then the value decreases as the volume occupied by the sphere in the bin decreases. I checked that it gives me a sphere of the radius wanted. Details on the subject can be found here(part 2.5 and figure 8 for accuracy).

--EDIT--

The code only seems to behave like this when all the particle masses are identical

My problem comes from how I am generating my filter. In my code, the way I associate weight to voxels not entirely in the sphere is discontinuous: For example I give the weight 0.1 to a voxel whose volume ratio is between 0.1 et 0.2.

Thus what happens when all points have the same mass is: I have multiples of 1 in my grid that I multiply with a finite number of coefficients, thus there is a finite nummber of possible values that my grid can take, and thus some bins are empty or at least 'less full'. This is less likely to happen when the masses of my particle are more continuously distributed.

A fix is thus to appoint the right weight to the voxels.

def remp(Kreal):
for i in range(b):
    for j in range(b):
        for k in range(b):
            a = cst - np.sqrt((X[i]-X[mid])**2 + (Y[j]-Y[mid])**2 + (Z[k]-Z[mid])**2)/s
            if a >= 0.1 and a < 0.99:
                Kreal[i][j][k] = a
            elif a >= 0.99:
               Kreal[i][j][k] = 1