diff --git a/plots/densityprofile.png b/plots/densityprofile.png
index ce7f7065bba88fbdb463434a86f3899c6d2f7701..d663e09822e03f5a2e037d61d4a382c6e766cb10 100644
Binary files a/plots/densityprofile.png and b/plots/densityprofile.png differ
diff --git a/plots/densityprofile.py b/plots/densityprofile.py
index 3c6560f80cc5230c3038a9f7cdeafe896dd58984..aecddd63022511d6761696fdbadab439bac9ea47 100644
--- a/plots/densityprofile.py
+++ b/plots/densityprofile.py
@@ -29,9 +29,19 @@ rhoHernquist = M / (2*np.pi) * a / (bins[1:]-dr/2) / (bins[1:]-dr/2 + a)**3
 print("Number of bins:",nbins)
 print("Bin size:",dr)
 
-for i in range(len(hist)):
-    if hist[i] == 0:
-        hist[i] = np.nan
+# for i in range(len(hist)):
+#     if hist[i] == 0:
+#         hist[i] = np.nan
+
+halfmassradius = 0
+tmpmass = 0
+rsorted = np.sort(r)
+for i in range(len(rsorted)):
+    tmpmass += m[i]
+    if tmpmass > M/2:
+        halfmassradius = rsorted[i]
+        break
+print("Half mass radius:",halfmassradius)
 
 # plot poissonian error bars
 plt.errorbar(bins[1:], rhoHernquist, yerr=np.sqrt(rhoHernquist), color='grey', fmt='.', label='Hernquist with Poisson error', markersize=0)