diff --git a/garlic/fig/nbody/test.R b/garlic/fig/nbody/test.R
index 0d0b619..21cb225 100644
--- a/garlic/fig/nbody/test.R
+++ b/garlic/fig/nbody/test.R
@@ -1,37 +1,53 @@
 library(ggplot2)
 library(dplyr)
 library(scales)
+library(jsonlite)
 
-# Load the dataset
-df=read.table("data.csv", col.names=c("blocksize", "time"))
+args=commandArgs(trailingOnly=TRUE)
 
-bs_unique = unique(df$blocksize)
-nbs=length(bs_unique)
+# Read the timetable from args[1]
+input_file = "input.json"
+if (length(args)>0) { input_file = args[1] }
+
+# Load the dataset in NDJSON format
+dataset = jsonlite::stream_in(file(input_file)) %>%
+	jsonlite::flatten()
+
+particles = unique(dataset$config.particles)
+
+# We only need the nblocks and time
+df = select(dataset, config.nblocks, config.hw.cpusPerSocket, time) %>%
+	rename(nblocks=config.nblocks,
+		cpusPerSocket=config.hw.cpusPerSocket)
+
+df = df %>% mutate(blocksPerCpu = nblocks / cpusPerSocket)
+df$nblocks = as.factor(df$nblocks)
+df$blocksPerCpuFactor = as.factor(df$blocksPerCpu)
 
 # Normalize the time by the median
-D=group_by(df, blocksize) %>%
-	mutate(tnorm = time / median(time) - 1) # %>%
-#	mutate(bad = (abs(tnorm) >= 0.01)) %>%
-#	mutate(color = ifelse(bad,"red","black"))
+D=group_by(df, nblocks) %>%
+	mutate(tnorm = time / median(time) - 1)
 
-D$bad = cut(abs(D$tnorm), breaks=c(-Inf, 0.01, +Inf), labels=c("good", "bad"))
+bs_unique = unique(df$nblocks)
+nbs=length(bs_unique)
 
 print(D)
 
-#ppi=300
-#h=5
-#w=5
-#png("box.png", width=w*ppi, height=h*ppi, res=ppi)
+ppi=300
+h=5
+w=5
+
+png("box.png", width=w*ppi, height=h*ppi, res=ppi)
 #
 #
 #
-# Create the plot with the normalized time vs blocksize
-p = ggplot(D, aes(x=blocksize, y=tnorm)) +
+# Create the plot with the normalized time vs nblocks
+p = ggplot(data=D, aes(x=blocksPerCpuFactor, y=tnorm)) +
 
 	# Labels
-	labs(x="Block size", y="Normalized time",
-              title="Nbody normalized time",
-              subtitle="@expResult@") +
+	labs(x="Blocks/CPU", y="Normalized time",
+              title=sprintf("Nbody normalized time. Particles=%d", particles), 
+              subtitle=input_file) +
 
 	# Center the title
 	#theme(plot.title = element_text(hjust = 0.5)) +
@@ -39,44 +55,49 @@ p = ggplot(D, aes(x=blocksize, y=tnorm)) +
 	# Black and white mode (useful for printing)
 	#theme_bw() +
 
-	# Draw boxplots
-	geom_boxplot(aes(group=blocksize)) +
+	# Add the maximum allowed error lines
+	geom_hline(yintercept=c(-0.01, 0.01),
+		linetype="dashed", color="red") +
 
-	# Use log2 scale in x
-	scale_x_continuous(trans=log2_trans(),
-			   breaks=bs_unique) +
+	# Draw boxplots
+	geom_boxplot() +
 
 	scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
 
-	# Add the maximum allowed error lines
-	geom_hline(yintercept=c(-0.01, 0.01),
-		linetype="dashed", color="red")
+	theme_bw() +
+
+	theme(plot.subtitle=element_text(size=8)) +
+
+	theme(legend.position = c(0.85, 0.85)) #+
+
+	# Place each variant group in one separate plot
+	#facet_wrap(~jemalloc)
 
 
 
 # Render the plot
 print(p)
-#
-## Save the png image
-#dev.off()
-#
-#png("scatter.png", width=w*ppi, height=h*ppi, res=ppi)
 
-## Create the plot with the normalized time vs blocksize
-#p = ggplot(D, aes(x=blocksize, y=time, color=bad)) +
+## Save the png image
+dev.off()
 #
-#	labs(x="Blocksize", y="Time (s)",
-#              title="Nbody granularity",
-#              subtitle="@expResult@") +
+png("scatter.png", width=w*ppi, height=h*ppi, res=ppi)
 #
-#	geom_point(shape=21, size=1.5) +
-#	scale_color_manual(values=c("black", "red")) +
-#	scale_x_continuous(trans=log2_trans(),
-#			   breaks=bs_unique) +
-#	scale_y_continuous(trans=log2_trans())
-#
-## Render the plot
-#print(p)
+## Create the plot with the normalized time vs nblocks
+p = ggplot(D, aes(x=blocksPerCpuFactor, y=time)) +
+
+	labs(x="Blocks/CPU", y="Time (s)",
+              title=sprintf("Nbody granularity. Particles=%d", particles), 
+              subtitle=input_file) +
+	theme_bw() +
+	theme(plot.subtitle=element_text(size=8)) +
+	theme(legend.position = c(0.5, 0.88)) +
+
+	geom_point(shape=21, size=3) +
+	scale_y_continuous(trans=log2_trans())
+
+# Render the plot
+print(p)
 
 # Save the png image
-#dev.off()
+dev.off()
diff --git a/overlay.nix b/overlay.nix
index 75cd33f..43abed0 100644
--- a/overlay.nix
+++ b/overlay.nix
@@ -335,10 +335,11 @@ let
       };
 
       # Datasets used in the figures
-      ds = with self.bsc.garlic; {
-        nbody = {
-          jemalloc = with exp.nbody; pp.merge [ baseline jemalloc ];
-          freeCpu = with exp.nbody; pp.merge [ baseline freeCpu ];
+      ds = with self.bsc.garlic; with pp; {
+        nbody = with exp.nbody; {
+          test      = merge [ baseline ];
+          jemalloc  = merge [ baseline jemalloc ];
+          freeCpu   = merge [ baseline freeCpu ];
         };
       };
 
@@ -346,6 +347,11 @@ let
       fig = with self.bsc.garlic; {
         nbody = {
 
+          test = pp.rPlot {
+            script = ./garlic/fig/nbody/test.R;
+            dataset = ds.nbody.test;
+          };
+
           jemalloc = pp.rPlot {
             script = ./garlic/fig/nbody/jemalloc.R;
             dataset = ds.nbody.jemalloc;