Calculate halo mass from a density profile

import h5py
import numpy as np

# constants and load data
h      = 0.7
rho_c  =  2.77536627e11 * 0.7**2;
Omegam = 0.3
rho_h  = Omegam * rho_c * 200 # at z=3
f      = h5py.File('fdm_profiles_p8640_b10000_m7e-23.hdf5', 'r')
allrho    = f['rho_ms_kpc3']
r      = f['r_bins_kpc'][:]

# change units
allrho   = allrho[:].T * 1e9  * h**2 # in unit M_solar/Mpc**3
rmid     = (r[1:] + r[:-1])/2 / h  / 1e3 # in unit Mpc
dr       = (r[1:] - r[:-1]) / h / 1e3 # in unit Mpc

#filter
select_simon_halo = 0
rho_halo = allrho[:,select_simon_halo]

nanmask  = np.isfinite(rho_halo)
rho_halo = rho_halo[nanmask]
rmid     = rmid[nanmask]
dr       = dr[nanmask]

count = f['count'][select_simon_halo][nanmask]
cell_len_Mpc = 10 / h / 8640
bin_masses = rho_halo * count * cell_len_Mpc**3
sphere_masses = np.cumsum(bin_masses)

mask     = np.where(sphere_masses / (4/3 * np.pi * rmid**3) <= rho_h)[0][0]
rho_halo = rho_halo[0:mask]
rmid     = rmid[0:mask]
dr       = dr[0:mask]

#size     = mask.shape[0]
size     = mask
rho_halo = rho_halo.reshape(size,1)
rmid     = rmid.reshape(size,1)
dr       = dr.reshape(size,1)

Mh        = sphere_masses[mask]
print ('Mh = %e'%(Mh))

import h5py
import numpy
import astropy.cosmology

h = 0.7
omega_m = 0.3
omega_l = 0.7

c = astropy.cosmology.FlatLambdaCDM(100 * h, omega_l)
rho0 = c.critical_density0.to('M_sun/kpc**3').value * omega_m / h**2
with h5py.File('fdm_profiles_p8640_b10000_m7e-23.hdf5', 'r') as f:
	p = f['rho_ms_kpc3'][...]
	p[numpy.isnan(p)] = 0
	r = f['r_bins_kpc'][...]

for halo in [0,1,36,62,75]:
	spherical = numpy.cumsum(p[halo] * 4/3 * numpy.pi * (r[1:]**3 - r[:-1]**3)) / (
		4/3 * numpy.pi * r[1:]**3
	)
	overdense_mask = spherical > 200 * rho0
	overdense_bins, = numpy.nonzero(overdense_mask)
	contiguous = numpy.arange(len(overdense_bins))
	overdense_bins = overdense_bins[overdense_bins == contiguous]
	i = overdense_bins[-1]
	mass = spherical[i] * 4/3 * numpy.pi * r[i + 1]**3
	print(f'{halo}: {mass / h:e}')