diff --git a/.clang-format b/.clang-format
index 05b10dc8ce7802fa782322aca13d6d67d14c6cf8..c43c6165c7c67310628bd74e3b38b2c0a5a072a7 100644
--- a/.clang-format
+++ b/.clang-format
@@ -44,7 +44,7 @@ BreakConstructorInitializersBeforeComma: false
 BreakConstructorInitializers: BeforeColon
 BreakAfterJavaFieldAnnotations: false
 BreakStringLiterals: true
-ColumnLimit:     120
+ColumnLimit:     100
 CommentPragmas:  '^ IWYU pragma:'
 CompactNamespaces: false
 ConstructorInitializerAllOnOneLineOrOnePerLine: false
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 4f8190234dfb414ebcf22cf9b68be15bcad2345c..d522499b09cd15557bea954c55c5d1c56848c87d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -33,9 +33,9 @@ detector:
     - ./util/build_detector.sh
 
 include:
-  - local: 'dis/config.yml'
-  - local: 'dvmp/config.yml'
-  - local: 'dvcs/config.yml'
+  - local: 'benchmarks/dis/config.yml'
+  - local: 'benchmarks/dvmp/config.yml'
+  - local: 'benchmarks/dvcs/config.yml'
 
 summary:
   stage: finish
diff --git a/.rootlogon.C b/.rootlogon.C
index 0be78dd3599f186dcabb06c3c8d0c20177a8f509..9c11c3a48d08428a9b1bab910a66c30b50245309 100644
--- a/.rootlogon.C
+++ b/.rootlogon.C
@@ -1,10 +1,16 @@
 {
   // Ensure fmt is loaded
   R__LOAD_LIBRARY(libfmt);
+  //
+  // top-level include-dir
+  gROOT->ProcessLine(".include include");
 
   // setup a local build directory so we don't polute our source code with
-  // ROOT dictionaries etc.
-  gSystem->SetBuildDir("/tmp/root_build");
+  // ROOT dictionaries etc. if desired
+  const char* build_dir = gSystem->Getenv("ROOT_BUILD_DIR");
+  if (build_dir) {
+    gSystem->SetBuildDir(build_dir);
+  }
 
   // style definition based off the ATLAS style
   TStyle* s = gStyle;
@@ -37,7 +43,7 @@
 
   // use large fonts
   // Int_t font=72; // Helvetica italics
-  Int_t font = 43; // Helvetica
+  Int_t    font  = 43; // Helvetica
   Double_t tsize = 26;
   s->SetTextFont(font);
 
diff --git a/accelerator b/accelerator
new file mode 160000
index 0000000000000000000000000000000000000000..f3ff428e3b926a41e95beaa984d8dc05cec37cc7
--- /dev/null
+++ b/accelerator
@@ -0,0 +1 @@
+Subproject commit f3ff428e3b926a41e95beaa984d8dc05cec37cc7
diff --git a/benchmarks.json b/benchmarks/benchmarks.json
similarity index 100%
rename from benchmarks.json
rename to benchmarks/benchmarks.json
diff --git a/dis/.gitignore b/benchmarks/dis/.gitignore
similarity index 100%
rename from dis/.gitignore
rename to benchmarks/dis/.gitignore
diff --git a/dis/README.md b/benchmarks/dis/README.md
similarity index 100%
rename from dis/README.md
rename to benchmarks/dis/README.md
diff --git a/benchmarks/dis/analysis/dis.h b/benchmarks/dis/analysis/dis.h
new file mode 100644
index 0000000000000000000000000000000000000000..0e88a01a008726cd6c6f2301fa07aef9d0da03ce
--- /dev/null
+++ b/benchmarks/dis/analysis/dis.h
@@ -0,0 +1,26 @@
+#ifndef DVMP_H
+#define DVMP_H
+
+#include <util.h>
+
+#include <algorithm>
+#include <cmath>
+#include <exception>
+#include <fmt/core.h>
+#include <limits>
+#include <string>
+#include <vector>
+
+#include <Math/Vector4D.h>
+
+// Additional utility functions for DVMP benchmarks. Where useful, these can be
+// promoted to the top-level util library
+namespace util {
+
+  // ADD EXTRA DIS UTILTIY FUNCTIONS HERE
+
+  //=========================================================================================================
+
+} // namespace util
+
+#endif
diff --git a/benchmarks/dis/analysis/dis_electrons.cxx b/benchmarks/dis/analysis/dis_electrons.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..02dc7b0dac383580c671aa98591f0a64cca3671b
--- /dev/null
+++ b/benchmarks/dis/analysis/dis_electrons.cxx
@@ -0,0 +1,116 @@
+#include "dis.h"
+#include "plot.h"
+
+#include <benchmark.h>
+#include <mt.h>
+#include <util.h>
+
+#include "ROOT/RDataFrame.hxx"
+#include <cmath>
+#include <fmt/color.h>
+#include <fmt/core.h>
+#include <fstream>
+#include <iostream>
+#include <nlohmann/json.hpp>
+#include <string>
+#include <vector>
+
+int dis_electrons(const std::string& config_name)
+{
+  // read our configuration
+  std::ifstream  config_file{config_name};
+  nlohmann::json config;
+  config_file >> config;
+
+  const std::string rec_file      = config["rec_file"];
+  const std::string detector      = config["detector"];
+  std::string       output_prefix = config["output_prefix"];
+  const std::string test_tag      = config["test_tag"];
+
+  fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
+             "Running DIS electron analysis...\n");
+  fmt::print(" - Detector package: {}\n", detector);
+  fmt::print(" - input file: {}\n", rec_file);
+  fmt::print(" - output prefix: {}\n", output_prefix);
+  fmt::print(" - test tag: {}\n", test_tag);
+
+  // create our test definition
+  // test_tag
+  eic::util::Test dis_Q2_resolution{
+      {{"name", fmt::format("{}_Q2_resolution", test_tag)},
+       {"title", "DIS Q2 resolution"},
+       {"description",
+        fmt::format("DIS Q2 resolution with {}, estimated using a Gaussian fit.", detector)},
+       {"quantity", "resolution (in %)"},
+       {"target", "0.1"}}};
+
+  // Run this in multi-threaded mode if desired
+  ROOT::EnableImplicitMT(kNumThreads);
+
+  const double electron_mass = util::get_pdg_mass("electron");
+
+  // Ensure our output prefix always ends on a dot, a slash or a dash
+  // Necessary when generating output plots
+  if (output_prefix.back() != '.' && output_prefix.back() != '/' && output_prefix.back() != '-') {
+    output_prefix += "-";
+  }
+
+  ROOT::RDataFrame d("events", rec_file);
+
+  // utility lambda functions to bind the reconstructed particle type
+  // (as we have no PID yet)
+  auto momenta_from_tracking =
+      [electron_mass](const std::vector<eic::TrackParametersData>& tracks) {
+        return util::momenta_from_tracking(tracks, electron_mass);
+      };
+
+  auto d0 = d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"})
+                .Define("N", "p_rec.size()")
+                .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
+                .Define("mom_sim", util::mom, {"p_sim"})
+                .Define("mom_rec", util::mom, {"p_rec"});
+
+  auto h_mom_sim = d0.Histo1D({"h_mom_sim", "; GeV; counts", 100, 0, 50}, "mom_sim");
+  auto h_mom_rec = d0.Histo1D({"h_mom_rec", "; GeV; counts", 100, 0, 50}, "mom_rec");
+
+  auto c = new TCanvas();
+
+  // Plot our histograms.
+  // TODO: to start I'm explicitly plotting the histograms, but want to
+  // factorize out the plotting code moving forward.
+  {
+    TCanvas c{"canvas", "canvas", 1200, 1200};
+    c.cd();
+    // gPad->SetLogx(false);
+    gPad->SetLogy(true);
+    auto& h1 = *h_mom_sim;
+    auto& h2 = *h_mom_rec;
+    // histogram style
+    h1.SetLineColor(plot::kMpBlue);
+    h1.SetLineWidth(2);
+    h2.SetLineColor(plot::kMpOrange);
+    h2.SetLineWidth(2);
+    // axes
+    h1.GetXaxis()->CenterTitle();
+    h1.GetYaxis()->CenterTitle();
+    // draw everything
+    h1.DrawClone("hist");
+    h2.DrawClone("hist same");
+    // FIXME hardcoded beam configuration
+    plot::draw_label(18, 275, detector);
+    TText* tptr1;
+    auto   t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+    t1->SetFillColorAlpha(kWhite, 0);
+    t1->SetTextFont(43);
+    t1->SetTextSize(25);
+    tptr1 = t1->AddText("simulated");
+    tptr1->SetTextColor(plot::kMpBlue);
+    tptr1 = t1->AddText("reconstructed");
+    tptr1->SetTextColor(plot::kMpOrange);
+    t1->Draw();
+
+    c.Print(fmt::format("{}momentum.png", output_prefix).c_str());
+  }
+
+  return 0;
+}
diff --git a/benchmarks/dis/benchmark.json b/benchmarks/dis/benchmark.json
new file mode 100644
index 0000000000000000000000000000000000000000..1754f81aaa83a4c2c448c00a67c242a41b005b1b
--- /dev/null
+++ b/benchmarks/dis/benchmark.json
@@ -0,0 +1,6 @@
+{
+  "name": "DIS/SIDIS",
+  "title": "DIS/SIDIS Benchmarks",
+  "description": "Benchmark for (Semi-inclusive) DIS",
+  "target": "0.8"
+}
diff --git a/benchmarks/dis/config.yml b/benchmarks/dis/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..6331fd56f4bde5000a7108c53c1296ccf7df4831
--- /dev/null
+++ b/benchmarks/dis/config.yml
@@ -0,0 +1,50 @@
+dis:generate:
+  stage: initialize
+  needs: []
+  timeout: 1 hours
+  cache:
+    key:
+      files:
+        - benchmarks/dis/generator/pythia_dis.cxx
+      prefix: "$CI_COMMIT_REF_SLUG"
+    paths:
+      - input/dis
+  artifacts:
+    paths:
+      - input
+  script:
+    - bash benchmarks/dis/gen.sh --config barrel --ebeam 18 --pbeam 275
+
+dis:process:
+  stage: process
+  needs: ["detector", "dis:generate"]
+  timeout: 1 hour
+  script:
+    - source options/env.sh
+    - ./util/compile_analyses.py dis
+    - ./benchmarks/dis/dis.sh --config barrel --ebeam 18 --pbeam 275
+  artifacts:
+    paths:
+      - results
+  retry:
+    max: 2
+    when:
+      - runner_system_failure
+  cache:
+    key:
+      files:
+        - .rootlogon.C
+        - util/compile_analyses.py
+      prefix: "$CI_COMMIT_REF_SLUG"
+    paths:
+      - .local/root_build
+
+dis:results:
+  stage: collect
+  needs: ["dis:process"]
+  script:
+    - ./util/collect_tests.py dis
+  artifacts:
+    paths:
+      - results/dis.json
+      - results/dis
diff --git a/benchmarks/dis/dis.sh b/benchmarks/dis/dis.sh
new file mode 100755
index 0000000000000000000000000000000000000000..77fb85a577d6726e413166a153f2893b81aab1ef
--- /dev/null
+++ b/benchmarks/dis/dis.sh
@@ -0,0 +1,141 @@
+#!/bin/bash
+
+## =============================================================================
+## Run the DIS benchmarks in 5 steps:
+## 1. Parse the command line and setup environment
+## 2. Detector simulation through npsim
+## 3. Digitization and reconstruction through Juggler
+## 4. Root-based Physics analyses
+## 5. Finalize
+## =============================================================================
+
+## make sure we launch this script from the project root directory
+PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/../..
+pushd ${PROJECT_ROOT}
+
+echo "Running the DIS benchmarks"
+
+## =============================================================================
+## Step 1: Setup the environment variables
+##
+## First parse the command line flags.
+## This sets the following environment variables:
+## - CONFIG:   The specific generator configuration
+## - EBEAM:    The electron beam energy
+## - PBEAM:    The ion beam energy
+source util/parse_cmd.sh $@
+
+## To run the reconstruction, we need the following global variables:
+## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon)
+## - JUGGLER_DETECTOR:       the detector package we want to use for this benchmark
+## - DETECTOR_PATH:          full path to the detector definitions
+##
+## You can ready options/env.sh for more in-depth explanations of the variables
+## and how they can be controlled.
+source options/env.sh
+
+## We also need the following benchmark-specific variables:
+##
+## - BENCHMARK_TAG: Unique identified for this benchmark process.
+## - BEAM_TAG:      Identifier for the chosen beam configuration
+## - INPUT_PATH:    Path for generator-level input to the benchmarks
+## - TMP_PATH:      Path for temporary data (not exported as artifacts)
+## - RESULTS_PATH:  Path for benchmark output figures and files
+##
+## You can read dvmp/env.sh for more in-depth explanations of the variables.
+source benchmarks/dis/env.sh
+
+## Get a unique file names based on the configuration options
+GEN_FILE=${INPUT_PATH}/gen-${CONFIG}_${JUGGLER_N_EVENTS}.hepmc
+
+SIM_FILE=${TMP_PATH}/sim-${CONFIG}.root
+SIM_LOG=${TMP_PATH}/sim-${CONFIG}.log
+
+
+REC_FILE=${TMP_PATH}/rec-${CONFIG}.root
+REC_LOG=${TMP_PATH}/sim-${CONFIG}.log
+
+PLOT_TAG=${CONFIG}
+
+## =============================================================================
+## Step 2: Run the simulation
+echo "Running Geant4 simulation"
+npsim --runType batch \
+      --part.minimalKineticEnergy 1000*GeV  \
+      -v WARNING \
+      --numberOfEvents ${JUGGLER_N_EVENTS} \
+      --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+      --inputFiles ${GEN_FILE} \
+      --outputFile ${SIM_FILE}
+if [ "$?" -ne "0" ] ; then
+  echo "ERROR running npsim"
+  exit 1
+fi
+
+## =============================================================================
+## Step 3: Run digitization & reconstruction
+echo "Running the digitization and reconstruction"
+## FIXME Need to figure out how to pass file name to juggler from the commandline
+## the tracker_reconstruction.py options file uses the following environment
+## variables:
+## - JUGGLER_SIM_FILE:    input detector simulation
+## - JUGGLER_REC_FILE:    output reconstructed data
+## - JUGGLER_DETECTOR_PATH: Location of the detector geometry
+## - JUGGLER_N_EVENTS:    number of events to process (part of global environment)
+## - JUGGLER_DETECTOR:    detector package (part of global environment)
+export JUGGLER_SIM_FILE=${SIM_FILE}
+export JUGGLER_REC_FILE=${REC_FILE}
+export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+  gaudirun.py options/tracker_reconstruction.py \
+  2>&1 > ${REC_LOG}
+## on-error, first retry running juggler again as there is still a random
+## crash we need to address FIXME
+if [ "$?" -ne "0" ] ; then
+  echo "Juggler crashed, retrying..."
+  xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+    gaudirun.py options/tracker_reconstruction.py \
+    2>&1 > ${REC_LOG}
+  if [ "$?" -ne "0" ] ; then
+    echo "ERROR running juggler, both attempts failed"
+    exit 1
+  fi
+fi
+
+## =============================================================================
+## Step 4: Analysis
+## write a temporary configuration file for the analysis script
+echo "Running analysis"
+CONFIG="${TMP_PATH}/${PLOT_TAG}.json"
+cat << EOF > ${CONFIG}
+{
+  "rec_file": "${REC_FILE}",
+  "detector": "${JUGGLER_DETECTOR}",
+  "output_prefix": "${RESULTS_PATH}/${PLOT_TAG}",
+  "test_tag": "${BEAM_TAG}"
+}
+EOF
+#cat ${CONFIG}
+root -b -q "benchmarks/dis/analysis/dis_electrons.cxx+(\"${CONFIG}\")"
+#root -b -q "benchmarks/dis/analysis/dis_electrons.cxx(\"${CONFIG}\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running rec_dis_electron script"
+  exit 1
+fi
+
+## =============================================================================
+## Step 5: finalize
+echo "Finalizing DIS benchmark"
+
+## Move over reconsturction artifacts as long as we don't have
+## too many events
+if [ "${JUGGLER_N_EVENTS}" -lt "500" ] ; then 
+  cp ${REC_FILE} ${RESULTS_PATH}
+fi
+
+## Always move over log files to the results path
+cp ${REC_LOG} ${RESULTS_PATH}
+
+## =============================================================================
+## All done!
+echo "DIS benchmarks complete"
diff --git a/benchmarks/dis/env.sh b/benchmarks/dis/env.sh
new file mode 100644
index 0000000000000000000000000000000000000000..1a5b153f74f55b5a764019826ad061c7a6d196ba
--- /dev/null
+++ b/benchmarks/dis/env.sh
@@ -0,0 +1,50 @@
+#!/bin/bash
+
+## =============================================================================
+## Local configuration variables for this particular benchmark 
+## It defines the following additional variables: 
+##
+##  - BENCHMARK_TAG:          Tag to identify this particular benchmark
+##  - BEAM_TAG                Tag to identify the beam configuration
+##  - INPUT_PATH:             Path for generator-level input to the benchmarks
+##  - TMP_PATH:               Path for temporary data (not exported as artifacts)
+##  - RESULTS_PATH:           Path for benchmark output figures and files
+##
+## This script assumes that EBEAM and PBEAM are set as part of the
+## calling script (usually as command line argument).
+##
+## =============================================================================
+
+## Tag for the local benchmark. Should be the same as the directory name for
+## this particular benchmark set (for clarity). 
+## This tag is used for the output artifacts directory (results/${JUGGLER_TAG}) 
+## and a tag in some of the output files.
+export BENCHMARK_TAG="dis"
+echo "Setting up the local environment for the ${BENCHMARK_TAG^^} benchmarks"
+
+## Extra beam tag to identify the desired beam configuration
+export BEAM_TAG="${EBEAM}on${PBEAM}"
+
+## Data path for input data (generator-level hepmc file)
+INPUT_PATH="input/${BENCHMARK_TAG}/${BEAM_TAG}"
+mkdir -p ${INPUT_PATH}
+export INPUT_PATH=`realpath ${INPUT_PATH}`
+echo "INPUT_PATH:             ${INPUT_PATH}"
+
+
+## Data path for temporary data (not exported as artifacts)
+TMP_PATH=${LOCAL_PREFIX}/tmp/${BEAM_TAG}
+mkdir -p ${TMP_PATH}
+export TMP_PATH=`realpath ${TMP_PATH}`
+echo "TMP_PATH:               ${TMP_PATH}"
+
+## Data path for benchmark output (plots and reconstructed files
+## if not too big).
+RESULTS_PATH="results/${BENCHMARK_TAG}/${BEAM_TAG}"
+mkdir -p ${RESULTS_PATH}
+export RESULTS_PATH=`realpath ${RESULTS_PATH}`
+echo "RESULTS_PATH:           ${RESULTS_PATH}"
+
+## =============================================================================
+## That's all!
+echo "Local environment setup complete."
diff --git a/benchmarks/dis/gen.sh b/benchmarks/dis/gen.sh
new file mode 100755
index 0000000000000000000000000000000000000000..4a42081ae6f437d542de677613eef22f4c6c9516
--- /dev/null
+++ b/benchmarks/dis/gen.sh
@@ -0,0 +1,96 @@
+#!/bin/bash
+
+## =============================================================================
+## Standin for a proper pythia generation process, similar to how we
+## generate events for DVMP
+## Runs in 5 steps:
+##   1. Parse the command line and setup the environment
+##   2. Check if we can load the requested file from the cache
+##   3. Build generator exe 
+##   4. Run the actual generator
+##   5. Finalize
+## =============================================================================
+## =============================================================================
+
+## make sure we launch this script from the project root directory
+PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/../..
+pushd ${PROJECT_ROOT}
+
+## =============================================================================
+## Step 1: Setup the environment variables
+## First parse the command line flags.
+## This sets the following environment variables:
+## - CONFIG:   The specific generator configuration --> not currenlty used FIXME
+## - EBEAM:    The electron beam energy --> not currently used FIXME
+## - PBEAM:    The ion beam energy --> not currently used FIXME
+source util/parse_cmd.sh $@
+
+## To run the generator, we need the following global variables:
+##
+## - LOCAL_PREFIX:      Place to cache local packages and data
+## - JUGGLER_N_EVENTS:  Number of events to process
+## - JUGGLER_RNG_SEED:  Random seed for event generation.
+##
+## You can read options/env.sh for more in-depth explanations of the variables
+## and how they can be controlled.
+source options/env.sh
+
+## We also need the following benchmark-specific variables:
+##
+## - BENCHMARK_TAG: Unique identified for this benchmark process.
+## - INPUT_PATH:    Path for generator-level input to the benchmarks
+## - TMP_PATH:      Path for temporary data (not exported as artifacts)
+##
+## You can read dvmp/env.sh for more in-depth explanations of the variables.
+source benchmarks/dis/env.sh
+
+## Get a unique file name prefix based on the configuration options
+GEN_TAG=gen-${CONFIG}_${JUGGLER_N_EVENTS} ## Generic file prefix
+
+## =============================================================================
+## Step 2: Check if we really need to run, or can use the cache.
+if [ -f "${INPUT_PATH}/${GEN_TAG}.hepmc" ]; then
+  echo "Found cached generator output for $GEN_TAG, no need to rerun"
+  exit 0
+fi
+
+echo "Generator output for $GEN_TAG not found in cache, need to run generator"
+
+## =============================================================================
+## Step 3: Build generator exe 
+##         TODO: need to configurability to the generator exe 
+
+echo "Compiling   benchmarks/dis/generator/pythia_dis.cxx ..."
+g++ benchmarks/dis/generator/pythia_dis.cxx -o pythia_dis  \
+   -I/usr/local/include  -Iinclude \
+   -O2 -std=c++11 -pedantic -W -Wall -Wshadow -fPIC  \
+   -L/usr/local/lib -Wl,-rpath,/usr/local/lib -lpythia8 -ldl \
+   -L/usr/local/lib -Wl,-rpath,/usr/local/lib -lHepMC3
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR compiling pythia"
+  exit 1
+fi
+echo "done"
+
+## =============================================================================
+## Step 4: Run the event generator
+echo "Running the generator"
+./pythia_dis ${TMP_PATH}/${GEN_TAG}.hepmc
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running pythia"
+  exit 1
+fi
+
+
+## =============================================================================
+## Step 5: Finally, move relevant output into the artifacts directory and clean up
+## =============================================================================
+echo "Moving generator output into ${INPUT_PATH}"
+mv ${TMP_PATH}/${GEN_TAG}.hepmc ${INPUT_PATH}/${GEN_TAG}.hepmc
+## this step only matters for local execution
+echo "Cleaning up"
+## does nothing
+
+## =============================================================================
+## All done!
+echo "$BENCHMARK_TAG event generation complete"
diff --git a/dis/generator/gen_central_electrons.cxx b/benchmarks/dis/generator/gen_central_electrons.cxx
similarity index 100%
rename from dis/generator/gen_central_electrons.cxx
rename to benchmarks/dis/generator/gen_central_electrons.cxx
diff --git a/dis/src/pythia_dis.cc b/benchmarks/dis/generator/pythia_dis.cxx
similarity index 94%
rename from dis/src/pythia_dis.cc
rename to benchmarks/dis/generator/pythia_dis.cxx
index ff819a7e3243373539aace192d0990003ca3368d..a48679683baf6ae692825275e713f79b8317f685 100644
--- a/dis/src/pythia_dis.cc
+++ b/benchmarks/dis/generator/pythia_dis.cxx
@@ -1,7 +1,7 @@
 #include "Pythia8/Pythia.h"
 #include "Pythia8Plugins/HepMC3.h"
-#include <unistd.h>
 #include <cassert>
+#include <unistd.h>
 
 using namespace Pythia8;
 
@@ -13,7 +13,7 @@ using std::string;
 enum class mode { none, help, list, part };
 
 struct settings {
-  double      E_electron   = 10.0;  // GeV
+  double      E_electron   = 18.0;  // GeV
   double      E_ion        = 275.0; // GeV
   std::string outfile      = "dis.hepmc";
   int         ion_PID      = 2212;
@@ -22,7 +22,7 @@ struct settings {
   bool        success      = false;
   double      Q2_min       = 4.0;
   int         N_events     = 1000;
-  mode    selected       = mode::none;
+  mode        selected     = mode::none;
 };
 
 template <typename T>
@@ -63,16 +63,15 @@ void print_man_page(T cli, const char* argv0)
 }
 //______________________________________________________________________________
 
-settings cmdline_settings(int argc, char* argv[]) {
+settings cmdline_settings(int argc, char* argv[])
+{
   settings s;
-  auto     lastOpt =
-      " options:" % (option("-h", "--help").set(s.selected, mode::help) % "show help",
-                     value("file", s.outfile).if_missing([] {
-                       std::cout << "You need to provide an output filename argument!\n";
-                     }) % "output file");
+  auto     lastOpt = " options:" % (option("-h", "--help").set(s.selected, mode::help) % "show help",
+                                value("file", s.outfile).if_missing([] {
+                                  std::cout << "You need to provide an output filename argument!\n";
+                                }) % "output file");
 
-  auto cli =
-      (command("help").set(s.selected, mode::help) | lastOpt);
+  auto cli = (command("help").set(s.selected, mode::help) | lastOpt);
 
   assert(cli.flags_are_prefix_free());
 
@@ -95,7 +94,8 @@ settings cmdline_settings(int argc, char* argv[]) {
 }
 //______________________________________________________________________________
 
-int main(int argc, char* argv[]) {
+int main(int argc, char* argv[])
+{
 
   settings s = cmdline_settings(argc, argv);
   if (!s.success) {
diff --git a/dvcs/config.yml b/benchmarks/dvcs/config.yml
similarity index 80%
rename from dvcs/config.yml
rename to benchmarks/dvcs/config.yml
index 0f63429f92eaf871c30b8bb9189e9a058cef87ee..d2c1a44710489b1e0e2af7ba590e70e2545b5a7a 100644
--- a/dvcs/config.yml
+++ b/benchmarks/dvcs/config.yml
@@ -3,7 +3,7 @@ dvcs:process:
   timeout: 1 hour
   needs: ["detector"]
   script:
-    - bash dvcs/dvcs.sh
+    - bash benchmarks/dvcs/dvcs.sh
   artifacts:
     paths:
       - results
@@ -18,5 +18,3 @@ dvcs:results:
     paths:
       - results
         #  reports:
-        #    junit: ["results/dvcs/dvcs_report.xml"]
-
diff --git a/dvcs/dvcs.sh b/benchmarks/dvcs/dvcs.sh
similarity index 87%
rename from dvcs/dvcs.sh
rename to benchmarks/dvcs/dvcs.sh
index 436406f26997f0f265f342f5222b00c1486493f1..35bbcafeddfe3e2e218ba286f209464cf8a297d7 100644
--- a/dvcs/dvcs.sh
+++ b/benchmarks/dvcs/dvcs.sh
@@ -15,9 +15,9 @@ echo "JUGGLER_FILE_NAME_TAG = ${JUGGLER_FILE_NAME_TAG}"
 ## - JUGGLER_DETECTOR:       the detector package we want to use for this benchmark
 ## - DETECTOR_PATH:          full path to the detector definitions
 ##
-## You can ready config/env.sh for more in-depth explanations of the variables
+## You can ready options/env.sh for more in-depth explanations of the variables
 ## and how they can be controlled.
-source config/env.sh
+source options/env.sh
 
 
 curl -o test_proton_dvcs_eic.hepmc "https://eicweb.phy.anl.gov/api/v4/projects/345/jobs/artifacts/master/raw/data/test_proton_dvcs_eic.hepmc?job=compile"
@@ -51,7 +51,7 @@ fi
 
 mkdir -p results/dvcs
 
-root -b -q "dvcs/scripts/dvcs_tests.cxx(\"${JUGGLER_REC_FILE}\")"
+root -b -q "benchmarks/dvcs/scripts/dvcs_tests.cxx(\"${JUGGLER_REC_FILE}\")"
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running root script"
   exit 1
@@ -62,10 +62,6 @@ if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
 cp ${JUGGLER_REC_FILE} results/dvcs/.
 fi
 
-# Collect the results
-#cp dvcs/report.xml results/dvcs/.
-#cp dvcs/report2.xml results/dvcs/.
-
 
 
 
diff --git a/dvcs/scripts/dvcs_tests.cxx b/benchmarks/dvcs/scripts/dvcs_tests.cxx
similarity index 100%
rename from dvcs/scripts/dvcs_tests.cxx
rename to benchmarks/dvcs/scripts/dvcs_tests.cxx
diff --git a/dvcs/scripts/merge_results.py b/benchmarks/dvcs/scripts/merge_results.py
similarity index 100%
rename from dvcs/scripts/merge_results.py
rename to benchmarks/dvcs/scripts/merge_results.py
diff --git a/benchmarks/dvmp/analysis/dvmp.h b/benchmarks/dvmp/analysis/dvmp.h
new file mode 100644
index 0000000000000000000000000000000000000000..371b454f8d92f4c0cd5b17edd4a8f5e1ecc77e77
--- /dev/null
+++ b/benchmarks/dvmp/analysis/dvmp.h
@@ -0,0 +1,56 @@
+#ifndef DVMP_H
+#define DVMP_H
+
+#include <util.h>
+
+#include <algorithm>
+#include <cmath>
+#include <exception>
+#include <fmt/core.h>
+#include <limits>
+#include <string>
+#include <vector>
+
+#include <Math/Vector4D.h>
+
+// Additional utility functions for DVMP benchmarks. Where useful, these can be
+// promoted to the top-level util library
+namespace util {
+
+  // Calculate the 4-vector sum of a given pair of particles
+  inline ROOT::Math::PxPyPzMVector
+  get_sum(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>& particle_pair)
+  {
+    return (particle_pair.first + particle_pair.second);
+  }
+
+  //========================================================================================================
+  // for structure functions
+
+  struct inv_quant { // add more when needed
+    double nu, Q2, x;
+  };
+
+  // for simu
+  inline inv_quant calc_inv_quant_simu(const std::vector<ROOT::Math::PxPyPzMVector>& parts)
+  {
+    ROOT::Math::PxPyPzMVector q(parts[0] - parts[2]);
+    ROOT::Math::PxPyPzMVector P(parts[3]);
+
+    double    nu         = q.Dot(P) / P.mass();
+    double    Q2         = -q.Dot(q);
+    inv_quant quantities = {nu, Q2, Q2 / 2. / P.mass() / nu};
+    return quantities;
+  }
+
+  inline double get_nu_simu(inv_quant quantities) { return quantities.nu / 1000.; }
+  inline double get_Q2_simu(inv_quant quantities) { return quantities.Q2; }
+  inline double get_x_simu(inv_quant quantities) { return quantities.x; }
+
+  // for tracking, add later
+
+  //=========================================================================================================
+
+} // namespace util
+
+#endif
diff --git a/dvmp/analysis/vm_invar.cxx b/benchmarks/dvmp/analysis/vm_invar.cxx
similarity index 54%
rename from dvmp/analysis/vm_invar.cxx
rename to benchmarks/dvmp/analysis/vm_invar.cxx
index f8723733fc38bc811fc14cdae22498d3debf68b2..03349825902a7e1210251fab0b4d8b467ddfae9d 100644
--- a/dvmp/analysis/vm_invar.cxx
+++ b/benchmarks/dvmp/analysis/vm_invar.cxx
@@ -1,7 +1,9 @@
-#include "benchmark.hh"
-#include "mt.h"
+#include "dvmp.h"
 #include "plot.h"
-#include "util.h"
+
+#include <benchmark.h>
+#include <mt.h>
+#include <util.h>
 
 #include <ROOT/RDataFrame.hxx>
 #include <cmath>
@@ -17,36 +19,35 @@
 // a desired vector meson (e.g. jpsi) and a desired decay particle (e.g. muon)
 // Output figures are written to our output prefix (which includes the output
 // file prefix), and labeled with our detector name.
-// TODO: I think it would be better to pass small json configuration file to
-//       the test, instead of this ever-expanding list of function arguments.
 // FIXME: MC does not trace back into particle history. Need to fix that
-int vm_invar(const std::string& config_name) {
+int vm_invar(const std::string& config_name)
+{
   // read our configuration
-  std::ifstream config_file{config_name};
+  std::ifstream  config_file{config_name};
   nlohmann::json config;
   config_file >> config;
 
-  const std::string rec_file = config["rec_file"];
-  const std::string vm_name = config["vm_name"];
-  const std::string decay_name = config["decay"];
-  const std::string detector = config["detector"];
-  std::string output_prefix = config["output_prefix"];
-  const std::string test_tag = config["test_tag"];
+  const std::string rec_file      = config["rec_file"];
+  const std::string vm_name       = config["vm_name"];
+  const std::string decay_name    = config["decay"];
+  const std::string detector      = config["detector"];
+  std::string       output_prefix = config["output_prefix"];
+  const std::string test_tag      = config["test_tag"];
 
   fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
              "Running VM invariant mass analysis...\n");
   fmt::print(" - Vector meson: {}\n", vm_name);
   fmt::print(" - Decay particle: {}\n", decay_name);
   fmt::print(" - Detector package: {}\n", detector);
+  fmt::print(" - input file: {}\n", rec_file);
   fmt::print(" - output prefix: {}\n", output_prefix);
 
   // create our test definition
   // test_tag
-  eic::util::Test vm_mass_resolution_test{
-      {{"name",
-        fmt::format("{}_{}_{}_mass_resolution", test_tag, vm_name, decay_name)},
+  eic::util::Test Q2_resolution_test{
+      {{"name", fmt::format("{}_{}_{}_Q2_resolution", test_tag, vm_name, decay_name)},
        {"title",
-        fmt::format("{} -> {} Invariant Mass Resolution", vm_name, decay_name)},
+        fmt::format("Q^2 Resolution for {} -> {} events with {}", vm_name, decay_name, detector)},
        {"description", "Invariant Mass Resolution calculated from raw "
                        "tracking data using a Gaussian fit."},
        {"quantity", "resolution"},
@@ -56,12 +57,11 @@ int vm_invar(const std::string& config_name) {
   ROOT::EnableImplicitMT(kNumThreads);
 
   // The particles we are looking for. E.g. J/psi decaying into e+e-
-  const double vm_mass = util::get_pdg_mass(vm_name);
+  const double vm_mass    = util::get_pdg_mass(vm_name);
   const double decay_mass = util::get_pdg_mass(decay_name);
 
   // Ensure our output prefix always ends on a dot, a slash or a dash
-  if (output_prefix.back() != '.' && output_prefix.back() != '/' &&
-      output_prefix.back() != '-') {
+  if (output_prefix.back() != '.' && output_prefix.back() != '/' && output_prefix.back() != '-') {
     output_prefix += "-";
   }
 
@@ -70,74 +70,68 @@ int vm_invar(const std::string& config_name) {
 
   // utility lambda functions to bind the vector meson and decay particle
   // types
-  auto momenta_from_tracking =
-      [decay_mass](const std::vector<eic::TrackParametersData>& tracks) {
-        return util::momenta_from_tracking(tracks, decay_mass);
-      };
+  auto momenta_from_tracking = [decay_mass](const std::vector<eic::TrackParametersData>& tracks) {
+    return util::momenta_from_tracking(tracks, decay_mass);
+  };
 
   //====================================================================
-    
+
   // Define analysis flow
-  auto d_im =
-      d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"})
-          .Define("N", "p_rec.size()")
-          .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
-          //================================================================
-          .Define("invariant_quantities", util::calc_inv_quant_simu, {"p_sim"})
-          .Define("nu_sim" , util::get_nu_simu, {"invariant_quantities"})
-          .Define("Q2_sim" , util::get_Q2_simu, {"invariant_quantities"})
-          .Define("x_sim" ,  util::get_x_simu, {"invariant_quantities"});
-          //================================================================
+  auto d_im = d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"})
+                  .Define("N", "p_rec.size()")
+                  .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
+                  //================================================================
+                  .Define("invariant_quantities", util::calc_inv_quant_simu, {"p_sim"})
+                  .Define("nu_sim", util::get_nu_simu, {"invariant_quantities"})
+                  .Define("Q2_sim", util::get_Q2_simu, {"invariant_quantities"})
+                  .Define("x_sim", util::get_x_simu, {"invariant_quantities"});
+  //================================================================
 
   // Define output histograms
-  
-  auto h_nu_sim = d_im.Histo1D(
-      {"h_nu_sim", ";#nu/1000;#", 100, 0., 2.}, "nu_sim");
-  auto h_Q2_sim = d_im.Histo1D(
-      {"h_Q2_sim", ";Q^{2};#", 100, 0., 15.}, "Q2_sim");
-  auto h_x_sim = d_im.Histo1D(
-      {"h_x_sim", ";x;#", 100, 0., 0.1}, "x_sim");
 
+  auto h_nu_sim = d_im.Histo1D({"h_nu_sim", ";#nu/1000;#", 100, 0., 2.}, "nu_sim");
+  auto h_Q2_sim = d_im.Histo1D({"h_Q2_sim", ";Q^{2};#", 100, 0., 15.}, "Q2_sim");
+  auto h_x_sim  = d_im.Histo1D({"h_x_sim", ";x;#", 100, 0., 0.1}, "x_sim");
 
   // Plot our histograms.
   // TODO: to start I'm explicitly plotting the histograms, but want to
   // factorize out the plotting code moving forward.
   {
-    
+
     // Print canvas to output file
-    
+
     TCanvas c{"canvas2", "canvas2", 1800, 600};
     c.Divide(3, 1, 0.0001, 0.0001);
-    //pad 1 nu
+    // pad 1 nu
     c.cd(1);
-    //gPad->SetLogx(false);
-    //gPad->SetLogy(false);
+    // gPad->SetLogx(false);
+    // gPad->SetLogy(false);
     auto& hnu = *h_nu_sim;
     // histogram style
     hnu.SetLineColor(plot::kMpBlue);
     hnu.SetLineWidth(2);
     // axes
     hnu.GetXaxis()->CenterTitle();
-    //hnu.GetXaxis()->SetTitle("#times1000");
+    // hnu.GetXaxis()->SetTitle("#times1000");
     // draw everything
     hnu.DrawClone("hist");
     // FIXME hardcoded beam configuration
-    plot::draw_label(10, 100, detector, vm_name, "#nu");
+    plot::draw_label(10, 100, detector);
     TText* tptr21;
-    auto t21 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+    auto   t21 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
     t21->SetFillColorAlpha(kWhite, 0);
     t21->SetTextFont(43);
     t21->SetTextSize(25);
     tptr21 = t21->AddText("simulated");
     tptr21->SetTextColor(plot::kMpBlue);
-    //tptr1 = t1->AddText("reconstructed");
-    //tptr1->SetTextColor(plot::kMpOrange);
+    // tptr1 = t1->AddText("reconstructed");
+    // tptr1->SetTextColor(plot::kMpOrange);
     t21->Draw();
-    
-    //pad 2 Q2
+
+    // pad 2 Q2
     c.cd(2);
-    //gPad->SetLogx(false);
-    //gPad->SetLogy(false);
+    // gPad->SetLogx(false);
+    // gPad->SetLogy(false);
     auto& hQ2 = *h_Q2_sim;
     // histogram style
     hQ2.SetLineColor(plot::kMpBlue);
@@ -147,22 +141,22 @@ int vm_invar(const std::string& config_name) {
     // draw everything
     hQ2.DrawClone("hist");
     // FIXME hardcoded beam configuration
-    plot::draw_label(10, 100, detector, vm_name, "Q^{2}");
+    plot::draw_label(10, 100, detector);
     TText* tptr22;
-    auto t22 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+    auto   t22 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
     t22->SetFillColorAlpha(kWhite, 0);
     t22->SetTextFont(43);
     t22->SetTextSize(25);
     tptr22 = t22->AddText("simulated");
     tptr22->SetTextColor(plot::kMpBlue);
-    //tptr1 = t1->AddText("reconstructed");
-    //tptr1->SetTextColor(plot::kMpOrange);
+    // tptr1 = t1->AddText("reconstructed");
+    // tptr1->SetTextColor(plot::kMpOrange);
     t22->Draw();
-    
-    //pad 1 nu
+
+    // pad 1 nu
     c.cd(3);
-    //gPad->SetLogx(false);
-    //gPad->SetLogy(false);
+    // gPad->SetLogx(false);
+    // gPad->SetLogy(false);
     auto& hx = *h_x_sim;
     // histogram style
     hx.SetLineColor(plot::kMpBlue);
@@ -172,28 +166,27 @@ int vm_invar(const std::string& config_name) {
     // draw everything
     hx.DrawClone("hist");
     // FIXME hardcoded beam configuration
-    plot::draw_label(10, 100, detector, vm_name, "x");
+    plot::draw_label(10, 100, detector);
     TText* tptr23;
-    auto t23 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+    auto   t23 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
     t23->SetFillColorAlpha(kWhite, 0);
     t23->SetTextFont(43);
     t23->SetTextSize(25);
     tptr23 = t23->AddText("simulated");
     tptr23->SetTextColor(plot::kMpBlue);
-    //tptr1 = t1->AddText("reconstructed");
-    //tptr1->SetTextColor(plot::kMpOrange);
+    // tptr1 = t1->AddText("reconstructed");
+    // tptr1->SetTextColor(plot::kMpOrange);
     t23->Draw();
-    
+
     c.Print(fmt::format("{}InvariantQuantities.png", output_prefix).c_str());
   }
 
-  // TODO we're not actually doing an IM fit yet, so for now just return an
+  // TODO we're not actually getting the resolutions yet
   // error for the test result
-  vm_mass_resolution_test.error(-1);
+  Q2_resolution_test.error(-1);
 
   // write out our test data
-  eic::util::write_test(vm_mass_resolution_test,
-                           fmt::format("{}vm_invar.json", output_prefix));
+  eic::util::write_test(Q2_resolution_test, fmt::format("{}vm_invar.json", output_prefix));
 
   // That's all!
   return 0;
diff --git a/dvmp/analysis/vm_mass.cxx b/benchmarks/dvmp/analysis/vm_mass.cxx
similarity index 56%
rename from dvmp/analysis/vm_mass.cxx
rename to benchmarks/dvmp/analysis/vm_mass.cxx
index fb736909c0e957805f8a3bec3ef91358d4f54b76..82abba6f77d8343df259e89d823e0b7d14419534 100644
--- a/dvmp/analysis/vm_mass.cxx
+++ b/benchmarks/dvmp/analysis/vm_mass.cxx
@@ -1,7 +1,9 @@
-#include "benchmark.hh"
-#include "mt.h"
+#include "dvmp.h"
 #include "plot.h"
-#include "util.h"
+
+#include <benchmark.h>
+#include <mt.h>
+#include <util.h>
 
 #include <ROOT/RDataFrame.hxx>
 #include <cmath>
@@ -20,33 +22,34 @@
 // TODO: I think it would be better to pass small json configuration file to
 //       the test, instead of this ever-expanding list of function arguments.
 // FIXME: MC does not trace back into particle history. Need to fix that
-int vm_mass(const std::string& config_name) {
+int vm_mass(const std::string& config_name)
+{
   // read our configuration
-  std::ifstream config_file{config_name};
+  std::ifstream  config_file{config_name};
   nlohmann::json config;
   config_file >> config;
 
-  const std::string rec_file = config["rec_file"];
-  const std::string vm_name = config["vm_name"];
-  const std::string decay_name = config["decay"];
-  const std::string detector = config["detector"];
-  std::string output_prefix = config["output_prefix"];
-  const std::string test_tag = config["test_tag"];
+  const std::string rec_file      = config["rec_file"];
+  const std::string vm_name       = config["vm_name"];
+  const std::string decay_name    = config["decay"];
+  const std::string detector      = config["detector"];
+  std::string       output_prefix = config["output_prefix"];
+  const std::string test_tag      = config["test_tag"];
 
   fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
              "Running VM invariant mass analysis...\n");
   fmt::print(" - Vector meson: {}\n", vm_name);
   fmt::print(" - Decay particle: {}\n", decay_name);
   fmt::print(" - Detector package: {}\n", detector);
+  fmt::print(" - input file: {}\n", rec_file);
   fmt::print(" - output prefix: {}\n", output_prefix);
 
   // create our test definition
   // test_tag
-  eic::util::Test vm_mass_resolution_test{
-      {{"name",
-        fmt::format("{}_{}_{}_mass_resolution", test_tag, vm_name, decay_name)},
-       {"title",
-        fmt::format("{} -> {} Invariant Mass Resolution", vm_name, decay_name)},
+  eic::util::Test mass_resolution_test{
+      {{"name", fmt::format("{}_{}_{}_mass_resolution", test_tag, vm_name, decay_name)},
+       {"title", fmt::format("{} Invariant Mass Resolution for {} -> {} with {}", vm_name, vm_name,
+                             decay_name, detector)},
        {"description", "Invariant Mass Resolution calculated from raw "
                        "tracking data using a Gaussian fit."},
        {"quantity", "resolution"},
@@ -56,12 +59,11 @@ int vm_mass(const std::string& config_name) {
   ROOT::EnableImplicitMT(kNumThreads);
 
   // The particles we are looking for. E.g. J/psi decaying into e+e-
-  const double vm_mass = util::get_pdg_mass(vm_name);
+  const double vm_mass    = util::get_pdg_mass(vm_name);
   const double decay_mass = util::get_pdg_mass(decay_name);
 
   // Ensure our output prefix always ends on a dot, a slash or a dash
-  if (output_prefix.back() != '.' && output_prefix.back() != '/' &&
-      output_prefix.back() != '-') {
+  if (output_prefix.back() != '.' && output_prefix.back() != '/' && output_prefix.back() != '-') {
     output_prefix += "-";
   }
 
@@ -70,56 +72,46 @@ int vm_mass(const std::string& config_name) {
 
   // utility lambda functions to bind the vector meson and decay particle
   // types
-  auto momenta_from_tracking =
-      [decay_mass](const std::vector<eic::TrackParametersData>& tracks) {
-        return util::momenta_from_tracking(tracks, decay_mass);
-      };
-  auto find_decay_pair =
-      [vm_mass](const std::vector<ROOT::Math::PxPyPzMVector>& parts) {
-        return util::find_decay_pair(parts, vm_mass);
-      };
-  
-      
-    //util::PrintGeant4(mcparticles2);
+  auto momenta_from_tracking = [decay_mass](const std::vector<eic::TrackParametersData>& tracks) {
+    return util::momenta_from_tracking(tracks, decay_mass);
+  };
+  auto find_decay_pair = [vm_mass](const std::vector<ROOT::Math::PxPyPzMVector>& parts) {
+    return util::find_decay_pair(parts, vm_mass);
+  };
+
+  // util::PrintGeant4(mcparticles2);
   // Define analysis flow
-  auto d_im =
-      d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"})
-          .Define("N", "p_rec.size()")
-          .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
-          .Define("decay_pair_rec", find_decay_pair, {"p_rec"})
-          .Define("decay_pair_sim", find_decay_pair, {"p_sim"})
-          .Define("mass_rec", util::get_im, {"decay_pair_rec"})
-          .Define("mass_sim", util::get_im, {"decay_pair_sim"})
-          .Define("pt_rec", util::get_pt, {"decay_pair_rec"})
-          .Define("pt_sim", util::get_pt, {"decay_pair_sim"})
-          .Define("phi_rec" , util::get_phi, {"decay_pair_rec"})
-          .Define("phi_sim" , util::get_phi, {"decay_pair_sim"})
-          .Define("rapidity_rec" , util::get_y, {"decay_pair_rec"})
-          .Define("rapidity_sim" , util::get_y, {"decay_pair_sim"});
-          
+  auto d_im = d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"})
+                  .Define("N", "p_rec.size()")
+                  .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
+                  .Define("decay_pair_rec", find_decay_pair, {"p_rec"})
+                  .Define("decay_pair_sim", find_decay_pair, {"p_sim"})
+                  .Define("p_vm_rec", "decay_pair_rec.first + decay_pair_rec.second")
+                  .Define("p_vm_sim", "decay_pair_sim.first + decay_pair_sim.second")
+                  //.Define("p_vm_sim", util::get_sum, {"decay_pair_sim"})
+                  .Define("mass_rec", "p_vm_rec.M()")
+                  .Define("mass_sim", "p_vm_sim.M()")
+                  .Define("pt_rec", "p_vm_rec.pt()")
+                  .Define("pt_sim", "p_vm_sim.pt()")
+                  .Define("phi_rec", "p_vm_rec.phi()")
+                  .Define("phi_sim", "p_vm_sim.phi()")
+                  .Define("eta_rec", "p_vm_rec.eta()")
+                  .Define("eta_sim", "p_vm_sim.eta()");
 
   // Define output histograms
-  auto h_im_rec = d_im.Histo1D(
-      {"h_im_rec", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_rec");
-  auto h_im_sim = d_im.Histo1D(
-      {"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_sim");
-      
-  auto h_pt_rec = d_im.Histo1D(
-      {"h_pt_rec", ";p_{T} (GeV/c);#", 400, 0., 40.}, "pt_rec");
-  auto h_pt_sim = d_im.Histo1D(
-      {"h_pt_sim", ";p_{T} (GeV/c);#", 400, 0., 40.}, "pt_sim"); 
-      
-  auto h_phi_rec = d_im.Histo1D(
-      {"h_phi_rec", ";#phi_{ll'};#", 90, -M_PI, M_PI}, "phi_rec");
-  auto h_phi_sim = d_im.Histo1D(
-      {"h_phi_sim", ";#phi_{ll'};#", 90, -M_PI, M_PI}, "phi_sim");
-      
-  auto h_y_rec = d_im.Histo1D(
-      {"h_y_rec", ";y_{ll'};#", 1000, -5., 5.}, "rapidity_rec");
-  auto h_y_sim = d_im.Histo1D(
-      {"h_y_sim", ";y_{ll'};#", 1000, -5., 5.}, "rapidity_sim");
+  auto h_im_rec =
+      d_im.Histo1D({"h_im_rec", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_rec");
+  auto h_im_sim =
+      d_im.Histo1D({"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_sim");
+
+  auto h_pt_rec = d_im.Histo1D({"h_pt_rec", ";p_{T} (GeV/c);#", 400, 0., 40.}, "pt_rec");
+  auto h_pt_sim = d_im.Histo1D({"h_pt_sim", ";p_{T} (GeV/c);#", 400, 0., 40.}, "pt_sim");
 
+  auto h_phi_rec = d_im.Histo1D({"h_phi_rec", ";#phi_{ll'};#", 90, -M_PI, M_PI}, "phi_rec");
+  auto h_phi_sim = d_im.Histo1D({"h_phi_sim", ";#phi_{ll'};#", 90, -M_PI, M_PI}, "phi_sim");
 
+  auto h_eta_rec = d_im.Histo1D({"h_eta_rec", ";#eta_{ll'};#", 1000, -5., 5.}, "eta_rec");
+  auto h_eta_sim = d_im.Histo1D({"h_eta_sim", ";#eta_{ll'};#", 1000, -5., 5.}, "eta_sim");
 
   // Plot our histograms.
   // TODO: to start I'm explicitly plotting the histograms, but want to
@@ -127,10 +119,10 @@ int vm_mass(const std::string& config_name) {
   {
     TCanvas c{"canvas", "canvas", 1200, 1200};
     c.Divide(2, 2, 0.0001, 0.0001);
-    //pad 1 mass
+    // pad 1 mass
     c.cd(1);
-    //gPad->SetLogx(false);
-    //gPad->SetLogy(false);
+    // gPad->SetLogx(false);
+    // gPad->SetLogy(false);
     auto& h11 = *h_im_sim;
     auto& h12 = *h_im_rec;
     // histogram style
@@ -145,9 +137,9 @@ int vm_mass(const std::string& config_name) {
     h11.DrawClone("hist");
     h12.DrawClone("hist same");
     // FIXME hardcoded beam configuration
-    plot::draw_label(10, 100, detector, vm_name, "Invariant mass");
+    plot::draw_label(10, 100, detector);
     TText* tptr1;
-    auto t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+    auto   t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
     t1->SetFillColorAlpha(kWhite, 0);
     t1->SetTextFont(43);
     t1->SetTextSize(25);
@@ -156,11 +148,11 @@ int vm_mass(const std::string& config_name) {
     tptr1 = t1->AddText("reconstructed");
     tptr1->SetTextColor(plot::kMpOrange);
     t1->Draw();
-    
-    //pad 2 pt
+
+    // pad 2 pt
     c.cd(2);
-    //gPad->SetLogx(false);
-    //gPad->SetLogy(false);
+    // gPad->SetLogx(false);
+    // gPad->SetLogy(false);
     auto& h21 = *h_pt_sim;
     auto& h22 = *h_pt_rec;
     // histogram style
@@ -175,9 +167,9 @@ int vm_mass(const std::string& config_name) {
     h21.DrawClone("hist");
     h22.DrawClone("hist same");
     // FIXME hardcoded beam configuration
-    plot::draw_label(10, 100, detector, vm_name, "Transverse Momentum");
+    plot::draw_label(10, 100, detector);
     TText* tptr2;
-    auto t2 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+    auto   t2 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
     t2->SetFillColorAlpha(kWhite, 0);
     t2->SetTextFont(43);
     t2->SetTextSize(25);
@@ -186,11 +178,11 @@ int vm_mass(const std::string& config_name) {
     tptr2 = t2->AddText("reconstructed");
     tptr2->SetTextColor(plot::kMpOrange);
     t2->Draw();
-    
-    //pad 3 phi
+
+    // pad 3 phi
     c.cd(3);
-    //gPad->SetLogx(false);
-    //gPad->SetLogy(false);
+    // gPad->SetLogx(false);
+    // gPad->SetLogy(false);
     auto& h31 = *h_phi_sim;
     auto& h32 = *h_phi_rec;
     // histogram style
@@ -205,9 +197,9 @@ int vm_mass(const std::string& config_name) {
     h31.DrawClone("hist");
     h32.DrawClone("hist same");
     // FIXME hardcoded beam configuration
-    plot::draw_label(10, 100, detector, vm_name, "#phi");
+    plot::draw_label(10, 100, detector);
     TText* tptr3;
-    auto t3 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+    auto   t3 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
     t3->SetFillColorAlpha(kWhite, 0);
     t3->SetTextFont(43);
     t3->SetTextSize(25);
@@ -216,13 +208,13 @@ int vm_mass(const std::string& config_name) {
     tptr3 = t3->AddText("reconstructed");
     tptr3->SetTextColor(plot::kMpOrange);
     t3->Draw();
-    
-    //pad 4 rapidity
+
+    // pad 4 rapidity
     c.cd(4);
-    //gPad->SetLogx(false);
-    //gPad->SetLogy(false);
-    auto& h41 = *h_y_sim;
-    auto& h42 = *h_y_rec;
+    // gPad->SetLogx(false);
+    // gPad->SetLogy(false);
+    auto& h41 = *h_eta_sim;
+    auto& h42 = *h_eta_rec;
     // histogram style
     h41.SetLineColor(plot::kMpBlue);
     h41.SetLineWidth(2);
@@ -235,9 +227,9 @@ int vm_mass(const std::string& config_name) {
     h41.DrawClone("hist");
     h42.DrawClone("hist same");
     // FIXME hardcoded beam configuration
-    plot::draw_label(10, 100, detector, vm_name, "Rapidity");
+    plot::draw_label(10, 100, detector);
     TText* tptr4;
-    auto t4 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+    auto   t4 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
     t4->SetFillColorAlpha(kWhite, 0);
     t4->SetTextFont(43);
     t4->SetTextSize(25);
@@ -248,16 +240,14 @@ int vm_mass(const std::string& config_name) {
     t4->Draw();
 
     c.Print(fmt::format("{}vm_mass_pt_phi_rapidity.png", output_prefix).c_str());
-
   }
 
   // TODO we're not actually doing an IM fit yet, so for now just return an
   // error for the test result
-  vm_mass_resolution_test.error(-1);
+  mass_resolution_test.error(-1);
 
   // write out our test data
-  eic::util::write_test(vm_mass_resolution_test,
-                           fmt::format("{}vm_mass.json", output_prefix));
+  eic::util::write_test(mass_resolution_test, fmt::format("{}_mass.json", output_prefix));
 
   // That's all!
   return 0;
diff --git a/dvmp/benchmark.json b/benchmarks/dvmp/benchmark.json
similarity index 100%
rename from dvmp/benchmark.json
rename to benchmarks/dvmp/benchmark.json
diff --git a/dvmp/config.yml b/benchmarks/dvmp/config.yml
similarity index 62%
rename from dvmp/config.yml
rename to benchmarks/dvmp/config.yml
index 5a17fddd9ab77ae7225b22f4299ab49cdf404ab5..941bafd42481774bee08e8d6017051b5777926ce 100644
--- a/dvmp/config.yml
+++ b/benchmarks/dvmp/config.yml
@@ -6,8 +6,8 @@ dvmp:generate:
   cache:
     key:
       files:
-        - dvmp/generator/jpsi_central.json
-        - dvmp/scripts/jpsi_central-generate.sh
+        - benchmarks/dvmp/generator/jpsi_central.json
+        - benchmarks/dvmp/scripts/jpsi_central-generate.sh
       prefix: "$CI_COMMIT_REF_SLUG"
     paths:
       - input/dvmp
@@ -15,7 +15,7 @@ dvmp:generate:
     paths:
       - input
   script:
-    - ./util/run_many.py ./dvmp/gen.sh 
+    - ./util/run_many.py ./benchmarks/dvmp/gen.sh 
           -c jpsi_barrel 
           -e 10x100 
           --decay muon --decay electron
@@ -26,7 +26,9 @@ dvmp:process:
   needs: ["detector", "dvmp:generate"]
   timeout: 1 hour
   script:
-    - ./util/run_many.py ./dvmp/dvmp.sh 
+    - source options/env.sh
+    - ./util/compile_analyses.py dvmp
+    - ./util/run_many.py ./benchmarks/dvmp/dvmp.sh 
           -c jpsi_barrel 
           -e 10x100 
           --decay muon --decay electron
@@ -35,6 +37,18 @@ dvmp:process:
   artifacts:
     paths:
       - results
+  retry:
+    max: 2
+    when:
+      - runner_system_failure
+  cache:
+    key:
+      files:
+        - .rootlogon.C
+        - util/compile_analyses.py
+      prefix: "$CI_COMMIT_REF_SLUG"
+    paths:
+      - .local/root_build
 
 dvmp:results:
   stage: collect
diff --git a/dvmp/dvmp.sh b/benchmarks/dvmp/dvmp.sh
similarity index 91%
rename from dvmp/dvmp.sh
rename to benchmarks/dvmp/dvmp.sh
index 0d263a96ef1e2c61d5fe8506bedc44f22de13f4f..1746a8600c1a3e20a8968b43db51884d2311e2fc 100755
--- a/dvmp/dvmp.sh
+++ b/benchmarks/dvmp/dvmp.sh
@@ -10,7 +10,7 @@
 ## =============================================================================
 
 ## make sure we launch this script from the project root directory
-PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/..
+PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/../..
 pushd ${PROJECT_ROOT}
 
 echo "Running the DVMP benchmarks"
@@ -35,9 +35,9 @@ source util/parse_cmd.sh $@
 ## - JUGGLER_DETECTOR:       the detector package we want to use for this benchmark
 ## - DETECTOR_PATH:          full path to the detector definitions
 ##
-## You can ready config/env.sh for more in-depth explanations of the variables
+## You can ready options/env.sh for more in-depth explanations of the variables
 ## and how they can be controlled.
-source config/env.sh
+source options/env.sh
 
 ## We also need the following benchmark-specific variables:
 ##
@@ -48,7 +48,7 @@ source config/env.sh
 ## - RESULTS_PATH:  Path for benchmark output figures and files
 ##
 ## You can read dvmp/env.sh for more in-depth explanations of the variables.
-source dvmp/env.sh
+source benchmarks/dvmp/env.sh
 
 ## Get a unique file names based on the configuration options
 GEN_FILE=${INPUT_PATH}/gen-${CONFIG}_${DECAY}_${JUGGLER_N_EVENTS}.hepmc
@@ -107,7 +107,6 @@ if [ "$?" -ne "0" ] ; then
     exit 1
   fi
 fi
-#ls -l
 
 ## =============================================================================
 ## Step 4: Analysis
@@ -127,10 +126,14 @@ EOF
 #cat ${CONFIG}
 
 ## run the analysis script with this configuration
-root -b -q "dvmp/analysis/vm_mass.cxx(\"${CONFIG}\")"
-root -b -q "dvmp/analysis/vm_invar.cxx(\"${CONFIG}\")"
+root -b -q "benchmarks/dvmp/analysis/vm_mass.cxx+(\"${CONFIG}\")"
 if [ "$?" -ne "0" ] ; then
-  echo "ERROR running root script"
+  echo "ERROR running vm_mass script"
+  exit 1
+fi
+root -b -q "benchmarks/dvmp/analysis/vm_invar.cxx+(\"${CONFIG}\")"
+if [ "$?" -ne "0" ] ; then
+  echo "ERROR running vm_invar script"
   exit 1
 fi
 
@@ -145,11 +148,11 @@ if [ "${JUGGLER_N_EVENTS}" -lt "500" ] ; then
 fi
 
 ## Always move over log files to the results path
-mv ${SIM_LOG} ${REC_LOG} ${RESULTS_PATH}
+mv ${REC_LOG} ${RESULTS_PATH}
 
 
 ## cleanup output files
-rm -f ${REC_FILE} ${SIM_FILE}
+#rm -f ${REC_FILE} ${SIM_FILE} ## --> not needed for CI
 
 ## =============================================================================
 ## All done!
diff --git a/dvmp/env.sh b/benchmarks/dvmp/env.sh
similarity index 100%
rename from dvmp/env.sh
rename to benchmarks/dvmp/env.sh
diff --git a/dvmp/gen.sh b/benchmarks/dvmp/gen.sh
similarity index 91%
rename from dvmp/gen.sh
rename to benchmarks/dvmp/gen.sh
index 0c9544b7b414f5ee96a40690c46ba13f1695fbe1..8bef0a116f14014869d134c38720e7717ad5e49d 100755
--- a/dvmp/gen.sh
+++ b/benchmarks/dvmp/gen.sh
@@ -11,9 +11,10 @@
 ## =============================================================================
 
 ## make sure we launch this script from the project root directory
-PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/..
+PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/../..
 pushd ${PROJECT_ROOT}
 
+
 ## =============================================================================
 ## Step 1: Setup the environment variables
 ##
@@ -32,9 +33,9 @@ source util/parse_cmd.sh $@
 ## - JUGGLER_N_EVENTS:  Number of events to process
 ## - JUGGLER_RNG_SEED:  Random seed for event generation.
 ##
-## You can read config/env.sh for more in-depth explanations of the variables
+## You can read options/env.sh for more in-depth explanations of the variables
 ## and how they can be controlled.
-source config/env.sh
+source options/env.sh
 
 ## We also need the following benchmark-specific variables:
 ##
@@ -43,7 +44,7 @@ source config/env.sh
 ## - TMP_PATH:      Path for temporary data (not exported as artifacts)
 ##
 ## You can read dvmp/env.sh for more in-depth explanations of the variables.
-source dvmp/env.sh
+source benchmarks/dvmp/env.sh
 
 ## Get a unique file name prefix based on the configuration options
 GEN_TAG=gen-${CONFIG}_${DECAY}_${JUGGLER_N_EVENTS} ## Generic file prefix
@@ -72,7 +73,7 @@ elif [ $DECAY = "muon" ]; then
 fi
 
 ## generate the config file for this generator setup
-CONFIG_IN="${BENCHMARK_TAG}/generator/${CONFIG}.json.in"
+CONFIG_IN="benchmarks/${BENCHMARK_TAG}/generator/${CONFIG}.json.in"
 echo "Creating generator configuration file ${GEN_TAG}.json"
 if [ ! -f ${CONFIG_IN} ]; then
   echo "ERROR: cannot find master config file ${CONFIG_IN}"
@@ -106,4 +107,6 @@ done
 echo "Cleaning up"
 rm ${TMP_PATH}/${GEN_TAG}.json
 
+## =============================================================================
 ## All done!
+echo "$BENCHMARK_TAG event generation complete"
diff --git a/dvmp/generator/jpsi_barrel.json.in b/benchmarks/dvmp/generator/jpsi_barrel.json.in
similarity index 100%
rename from dvmp/generator/jpsi_barrel.json.in
rename to benchmarks/dvmp/generator/jpsi_barrel.json.in
diff --git a/dvmp/generator/jpsi_central.json.in b/benchmarks/dvmp/generator/jpsi_central.json.in
similarity index 100%
rename from dvmp/generator/jpsi_central.json.in
rename to benchmarks/dvmp/generator/jpsi_central.json.in
diff --git a/dis/analysis/rec_dis_electrons.cxx b/dis/analysis/rec_dis_electrons.cxx
deleted file mode 100644
index b7227a1ddb33c2564fd160e2124da774a9711161..0000000000000000000000000000000000000000
--- a/dis/analysis/rec_dis_electrons.cxx
+++ /dev/null
@@ -1,113 +0,0 @@
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "eicd/TrackParametersCollection.h"
-#include "eicd/ClusterCollection.h"
-#include "eicd/ClusterData.h"
-
-using ROOT::RDataFrame;
-using namespace ROOT::VecOps;
-
-auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
-  std::vector<double> result;
-  for (size_t i = 0; i < in.size(); ++i) {
-    result.push_back(std::abs(1.0/(in[i].qOverP)));
-  }
-  return result;
-};
-
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
-  std::vector<float> result;
-  for (size_t i = 0; i < in.size(); ++i) {
-    result.push_back(std::sqrt(in[i].psx * in[i].psx + in[i].psy * in[i].psy));
-  }
-  return result;
-}
-
-auto momentum = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
-  std::vector<double> result;
-  for (size_t i = 0; i < in.size(); ++i) {
-   result.push_back(in[i].E());
-  }
-  return result;
-};
-auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
-  std::vector<double> result;
-  for (size_t i = 0; i < in.size(); ++i) {
-   result.push_back(in[i].Theta()*180/M_PI);
-  }
-  return result;
-};
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
-  std::vector<ROOT::Math::PxPyPzMVector> result;
-  ROOT::Math::PxPyPzMVector lv;
-  for (size_t i = 0; i < in.size(); ++i) {
-    lv.SetCoordinates(in[i].psx, in[i].psy, in[i].psz, in[i].mass);
-    result.push_back(lv);
-  }
-  return result;
-};
-
-auto delta_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
-  std::vector<double> res;
-  for (const auto& p1 : thrown) {
-    for (const auto& p2 : tracks) {
-      res.push_back(p1 - p2);
-    }
-  }
-  return res;
-};
-
-int rec_dis_electrons(const char* fname = "topside/rec_central_electrons.root")
-{
-
-  ROOT::EnableImplicitMT();
-  ROOT::RDataFrame df("events", fname);
-
-  auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles2[isThrown]")
-                 .Define("thrownP", fourvec, {"thrownParticles"})
-                 .Define("p_thrown", momentum, {"thrownP"})
-                 .Define("nTracks", "outputTrackParameters.size()")
-                 .Define("p_track", p_track, {"outputTrackParameters"})
-                 .Define("p_track1", p_track, {"outputTrackParameters1"})
-                 .Define("p_track2", p_track, {"outputTrackParameters2"})
-                 .Define("delta_p",delta_p, {"p_track", "p_thrown"})
-                 .Define("delta_p1",delta_p, {"p_track1", "p_thrown"})
-                 .Define("delta_p2",delta_p, {"p_track2", "p_thrown"});
-
-  auto h_nTracks = df0.Histo1D({"h_nTracks", "; N tracks ", 10, 0, 10}, "nTracks");
-  auto h_pTracks = df0.Histo1D({"h_pTracks", "; GeV/c ", 100, 0, 10}, "p_track");
-
-  auto h_delta_p  = df0.Histo1D({"h_delta_p", "; GeV/c ",  100, -10, 10}, "delta_p");
-  auto h_delta_p1 = df0.Histo1D({"h_delta_p1", "; GeV/c ", 100, -10, 10}, "delta_p1");
-  auto h_delta_p2 = df0.Histo1D({"h_delta_p2", "; GeV/c ", 100, -10, 10}, "delta_p2");
-
-  auto c = new TCanvas();
-
-  h_nTracks->DrawCopy();
-  c->SaveAs("results/dis/rec_central_electrons_nTracks.png");
-  c->SaveAs("results/dis/rec_central_electrons_nTracks.pdf");
-
-  h_pTracks->DrawCopy();
-  c->SaveAs("results/dis/rec_central_electrons_pTracks.png");
-  c->SaveAs("results/dis/rec_central_electrons_pTracks.pdf");
-
-  THStack * hs = new THStack("hs_delta_p","; GeV/c "); 
-  TH1D* h1 = (TH1D*) h_delta_p->Clone();
-  hs->Add(h1);
-  h1 = (TH1D*) h_delta_p1->Clone();
-  h1->SetLineColor(2);
-  hs->Add(h1);
-  h1 = (TH1D*) h_delta_p2->Clone();
-  h1->SetLineColor(4);
-  h1->SetFillStyle(3001);
-  h1->SetFillColor(4);
-  hs->Add(h1);
-  hs->Draw("nostack");
-  c->SaveAs("results/dis/rec_central_electrons_delta_p.png");
-  c->SaveAs("results/dis/rec_central_electrons_delta_p.pdf");
-
-  return 0;
-}
diff --git a/dis/config.yml b/dis/config.yml
deleted file mode 100644
index 7db7df61e51aa673c989eb022888eaf2f9de2e35..0000000000000000000000000000000000000000
--- a/dis/config.yml
+++ /dev/null
@@ -1,26 +0,0 @@
-dis:generate:
-  stage: initialize
-  image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:latest
-  needs: []
-  timeout: 1 hours
-  artifacts:
-    paths:
-      - results
-  script:
-    - bash dis/gen.sh
-
-dis:process:
-  stage: process
-  needs: ["detector", "dis:generate"]
-  timeout: 1 hour
-  artifacts:
-    paths:
-      - results
-  script:
-    - echo "DIS benchmarks"
-
-dis:results:
-  stage: collect
-  needs: ["dis:process"]
-  script:
-    - echo "All DIS benchmarks successful"
diff --git a/dis/dis.sh b/dis/dis.sh
deleted file mode 100644
index 7c68e130ef94d42141529c1140bd0006f2eecd21..0000000000000000000000000000000000000000
--- a/dis/dis.sh
+++ /dev/null
@@ -1,100 +0,0 @@
-#!/bin/bash
-
-## =============================================================================
-## Run the DVMP benchmarks in 5 steps:
-## 1. Build/install detector package
-## 2. Detector simulation through npsim
-## 3. Digitization and reconstruction through Juggler
-## 4. Root-based Physics analyses
-## 5. Finalize
-## =============================================================================
-
-echo "Running the DIS benchmarks"
-
-## make sure we launch this script from the project root directory
-PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/..
-pushd ${PROJECT_ROOT}
-
-## =============================================================================
-## Load the environment variables. To build the detector we need the following
-## variables:
-##
-## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon)
-## - JUGGLER_DETECTOR:       the detector package we want to use for this benchmark
-## - DETECTOR_PATH:          full path to the detector definitions
-##
-## You can ready config/env.sh for more in-depth explanations of the variables
-## and how they can be controlled.
-source config/env.sh
-
-## Extra environment variables for DVMP:
-## file tag for these tests
-JUGGLER_FILE_NAME_TAG="dis"
-# TODO use the input file name, as we will be generating a lot of these
-# in the future...
-# FIXME Generator file hardcoded for now
-## note: these variables need to be exported to be accessible from
-##       the juggler options.py. We should really work on a dedicated
-##       juggler launcher to get rid of these "magic" variables. FIXME
-export JUGGLER_GEN_FILE="results/dis/${JUGGLER_FILE_NAME_TAG}.hepmc"
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
-export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
-
-
-## =============================================================================
-## Step 1: Build/install the desired detector package
-## TODO remove this
-#bash util/build_detector.sh
-
-
-## =============================================================================
-## Step 2: Run the simulation
-echo "Running Geant4 simulation"
-npsim --runType batch \
-      --part.minimalKineticEnergy 1000*GeV  \
-      -v WARNING \
-      --numberOfEvents ${JUGGLER_N_EVENTS} \
-      --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
-      --inputFiles ${JUGGLER_GEN_FILE} \
-      --outputFile ${JUGGLER_SIM_FILE}
-if [ "$?" -ne "0" ] ; then
-  echo "ERROR running npsim"
-  exit 1
-fi
-
-## =============================================================================
-## Step 3: Run digitization & reconstruction
-echo "Running the digitization and reconstruction"
-# FIXME Need to figure out how to pass file name to juggler from the commandline
-xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
-  gaudirun.py options/tracker_reconstruction.py
-if [ "$?" -ne "0" ] ; then
-  echo "ERROR running juggler"
-  exit 1
-fi
-ls -l
-
-## =============================================================================
-## Step 4: Analysis
-root -b -q "dis/analysis/rec_dis_electrons.cxx(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")"
-if [[ "$?" -ne "0" ]] ; then
-  echo "ERROR running root script"
-  exit 1
-fi
-
-## =============================================================================
-## Step 5: finalize
-echo "Finalizing ${JUGGLER_FILE_NAME_TAG} benchmark"
-
-## Copy over reconsturction artifacts as long as we don't have
-## too many events
-if [ "${JUGGLER_N_EVENTS}" -lt "500" ] ; then 
-  cp ${JUGGLER_REC_FILE} results/dis/.
-fi
-
-## cleanup output files
-rm ${JUGGLER_REC_FILE} ${JUGGLER_SIM_FILE}
-
-## =============================================================================
-## All done!
-echo "${JUGGLER_FILE_NAME_TAG} benchmarks complete"
diff --git a/dis/gen.sh b/dis/gen.sh
deleted file mode 100644
index e56e4edd49380efd16dc897b1589c6a56bac9909..0000000000000000000000000000000000000000
--- a/dis/gen.sh
+++ /dev/null
@@ -1,65 +0,0 @@
-#!/bin/bash
-
-## =============================================================================
-## Standin for a proper pythia generation process, similar to how we
-## generate events for DVMP
-## =============================================================================
-
-## TODO: use JUGGLER_FILE_NAME_TAG instead of explicitly refering to dis
-
-echo "Running the DIS benchmarks"
-
-## make sure we launch this script from the project root directory
-PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/..
-pushd ${PROJECT_ROOT}
-
-## =============================================================================
-## Load the environment variables. To build the detector we need the following
-## variables:
-##
-## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon)
-## - JUGGLER_DETECTOR:       the detector package we want to use for this benchmark
-## - DETECTOR_PATH:          full path to the detector definitions
-##
-## You can ready config/env.sh for more in-depth explanations of the variables
-## and how they can be controlled.
-source config/env.sh
-
-## Setup local environment
-export DATA_PATH=results/dis
-
-## Extra environment variables for DVMP:
-## file tag for these tests
-JUGGLER_FILE_NAME_TAG="dis"
-
-## =============================================================================
-## Step 1: Dummy event generator
-## TODO better file name that encodes the actual configuration we're running
-echo "Compiling   dis/src/pythia_dis.cc ..."
-g++ dis/src/pythia_dis.cc -o pythia_dis  \
-   -I/usr/local/include  -Iinclude \
-   -O2 -std=c++11 -pedantic -W -Wall -Wshadow -fPIC  \
-   -L/usr/local/lib -Wl,-rpath,/usr/local/lib -lpythia8 -ldl \
-   -L/usr/local/lib -Wl,-rpath,/usr/local/lib -lHepMC3
-if [[ "$?" -ne "0" ]] ; then
-  echo "ERROR compiling pythia"
-  exit 1
-fi
-echo "done"
-pwd
-ls -lrth
-
-./pythia_dis dis.hepmc
-if [[ "$?" -ne "0" ]] ; then
-  echo "ERROR running pythia"
-  exit 1
-fi
-
-## =============================================================================
-## Step 2: finalize
-echo "Moving event generator output into ${DATA_PATH}"
-#mv .local/${JUGGLER_FILE_NAME_TAG}.hepmc ${DATA_PATH}/${JUGGLER_FILE_NAME_TAG}.hepmc
-
-## =============================================================================
-## All done!
-echo "dis event generation complete"
diff --git a/dis/generator/placeholder b/dis/generator/placeholder
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/dis/util/placeholder b/dis/util/placeholder
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/dvcs/tests/report.xml b/dvcs/tests/report.xml
deleted file mode 100644
index 2ba9c6e484d8b0ce7c0c7164759f22f0d7d82d60..0000000000000000000000000000000000000000
--- a/dvcs/tests/report.xml
+++ /dev/null
@@ -1,10 +0,0 @@
-<testsuites>
-  <testsuite name="DVCS" tests="3" failures="1">
-    <testcase classname="foo1" name="ASuccessfulTest" time="10.0"/>
-    <testcase classname="foo2" name="AnotherSuccessfulTest"  time="20.0"/>
-    <testcase classname="foo3" name="AFailingTest"  time="30.0">
-      <failure type="NotEnoughFoo"> details about failure </failure>
-    </testcase>
-  </testsuite>
-</testsuites>
-
diff --git a/dvcs/tests/report2.xml b/dvcs/tests/report2.xml
deleted file mode 100644
index f437ef44d7b01cb33b69259fead289a20fc82072..0000000000000000000000000000000000000000
--- a/dvcs/tests/report2.xml
+++ /dev/null
@@ -1,12 +0,0 @@
-<testsuite tests="3">
-    <testcase time="1.00" classname="bar1" name="ASuccessfulTest2"/>
-    <testcase time="2.00" classname="bar2" name="AnotherSuccessfulTest2">
-      <!--
-      <system-out>[[ATTACHMENT|/absolute/path/to/some/file]]</system-out>
-      -->
-    </testcase>
-    <testcase time="3.00" classname="bar3" name="AFailingTest2">
-        <failure type="NotEnoughFoo"> MORE details about failure </failure>
-    </testcase>
-</testsuite>
-
diff --git a/dvmp/.gitignore b/dvmp/.gitignore
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/dvmp/analysis/benchmark.hh b/dvmp/analysis/benchmark.hh
deleted file mode 120000
index a96772424fe616aa97770a31a305b3ae74f9911b..0000000000000000000000000000000000000000
--- a/dvmp/analysis/benchmark.hh
+++ /dev/null
@@ -1 +0,0 @@
-../../util/benchmark.hh
\ No newline at end of file
diff --git a/dvmp/analysis/exception.hh b/dvmp/analysis/exception.hh
deleted file mode 120000
index e50e23b6af1bfcf85181374407bd3e80d9d65696..0000000000000000000000000000000000000000
--- a/dvmp/analysis/exception.hh
+++ /dev/null
@@ -1 +0,0 @@
-../../util/exception.hh
\ No newline at end of file
diff --git a/dvmp/analysis/plot.h b/dvmp/analysis/plot.h
deleted file mode 100644
index 69ebba341af239541496a410c66f78eb33e8adac..0000000000000000000000000000000000000000
--- a/dvmp/analysis/plot.h
+++ /dev/null
@@ -1,40 +0,0 @@
-#ifndef PLOT_H
-#define PLOT_H
-
-#include <TColor.h>
-#include <fmt/core.h>
-#include <vector>
-
-namespace plot {
-
-const int kMpBlue = TColor::GetColor(0x1f, 0x77, 0xb4);
-const int kMpOrange = TColor::GetColor(0xff, 0x7f, 0x0e);
-const int kMpGreen = TColor::GetColor(0x2c, 0xa0, 0x2c);
-const int kMpRed = TColor::GetColor(0xd6, 0x27, 0x28);
-const int kMpPurple = TColor::GetColor(0x94, 0x67, 0xbd);
-const int kMpBrown = TColor::GetColor(0x8c, 0x56, 0x4b);
-const int kMpPink = TColor::GetColor(0xe3, 0x77, 0xc2);
-const int kMpGrey = TColor::GetColor(0x7f, 0x7f, 0x7f);
-const int kMpMoss = TColor::GetColor(0xbc, 0xbd, 0x22);
-const int kMpCyan = TColor::GetColor(0x17, 0xbe, 0xcf);
-
-const std::vector<int> kPalette = {kMpBlue,   kMpOrange, kMpGreen, kMpRed,
-                                   kMpPurple, kMpBrown,  kMpPink,  kMpGrey,
-                                   kMpMoss,   kMpCyan};
-
-void draw_label(int ebeam, int pbeam, const std::string_view detector,
-                std::string_view vm, std::string_view what) {
-  auto t = new TPaveText(.15, 0.800, .7, .925, "NB NDC");
-  t->SetFillColorAlpha(kWhite, 0.4);
-  t->SetTextFont(43);
-  t->SetTextSize(25);
-  t->AddText(fmt::format("#bf{{{} }}SIMULATION", detector).c_str());
-  t->AddText(fmt::format("{} for {} DVMP.", what, vm).c_str());
-  t->AddText(fmt::format("{} GeV on {} GeV", ebeam, pbeam, what).c_str());
-  t->SetTextAlign(12);
-  t->Draw();
-}
-
-} // namespace plot
-
-#endif
diff --git a/dvmp/analysis/util.h b/dvmp/analysis/util.h
deleted file mode 100644
index dad7048e016d7c5d03cca9eb37481472ea4e5a74..0000000000000000000000000000000000000000
--- a/dvmp/analysis/util.h
+++ /dev/null
@@ -1,198 +0,0 @@
-#ifndef UTIL_H
-#define UTIL_H
-
-#include <algorithm>
-#include <cmath>
-#include <exception>
-#include <fmt/core.h>
-#include <limits>
-#include <string>
-#include <vector>
-
-#include <Math/Vector4D.h>
-
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "eicd/TrackParametersCollection.h"
-
-namespace util {
-
-// Exception definition for unknown particle errors
-// FIXME: A utility exception base class should be included in the analysis
-//        utility library, so we can skip most of this boilerplate
-class unknown_particle_error : public std::exception {
-public:
-  unknown_particle_error(std::string_view particle) : m_particle{particle} {}
-  virtual const char* what() const throw() {
-    return fmt::format("Unknown particle type: {}", m_particle).c_str();
-  }
-  virtual const char* type() const throw() { return "unknown_particle_error"; }
-
-private:
-  const std::string m_particle;
-};
-
-// Simple function to return the appropriate PDG mass for the particles
-// we care about for this process.
-// FIXME: consider something more robust (maybe based on hepPDT) to the
-//        analysis utility library
-inline double get_pdg_mass(std::string_view part) {
-  if (part == "electron") {
-    return 0.0005109989461;
-  } else if (part == "muon") {
-    return .1056583745;
-  } else if (part == "jpsi") {
-    return 3.0969;
-  } else if (part == "upsilon") {
-    return 9.49630;
-  } else {
-    throw unknown_particle_error{part};
-  }
-}
-
-// Get a vector of 4-momenta from raw tracking info, using an externally
-// provided particle mass assumption.
-// FIXME: should be part of utility library
-inline auto
-momenta_from_tracking(const std::vector<eic::TrackParametersData>& tracks,
-                      const double mass) {
-  std::vector<ROOT::Math::PxPyPzMVector> momenta{tracks.size()};
-  // transform our raw tracker info into proper 4-momenta
-  std::transform(tracks.begin(), tracks.end(), momenta.begin(),
-                 [mass](const auto& track) {
-                   // make sure we don't divide by zero
-                   if (fabs(track.qOverP) < 1e-9) {
-                     return ROOT::Math::PxPyPzMVector{};
-                   }
-                   const double p = fabs(1. / track.qOverP);
-                   const double px = p * cos(track.phi) * sin(track.theta);
-                   const double py = p * sin(track.phi) * sin(track.theta);
-                   const double pz = p * cos(track.theta);
-                   return ROOT::Math::PxPyPzMVector{px, py, pz, mass};
-                 });
-  return momenta;
-}
-
-// Get a vector of 4-momenta from the simulation data.
-// FIXME: should be part of utility library
-// TODO: Add PID selector (maybe using ranges?)
-inline auto
-momenta_from_simulation(const std::vector<dd4pod::Geant4ParticleData>& parts) {
-  std::vector<ROOT::Math::PxPyPzMVector> momenta{parts.size()};
-  // transform our simulation particle data into 4-momenta
-  std::transform(parts.begin(), parts.end(), momenta.begin(),
-                 [](const auto& part) {
-                   return ROOT::Math::PxPyPzMVector{part.psx, part.psy,
-                                                    part.psz, part.mass};
-                 });
-  return momenta;
-}
-
-// Find the decay pair candidates from a vector of particles (parts),
-// with invariant mass closest to a desired value (pdg_mass)
-// FIXME: not sure if this belongs here, or in the utility library. Probably the
-//        utility library
-inline std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>
-find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts,
-                const double pdg_mass) {
-  int first = -1;
-  int second = -1;
-  double best_mass = -1;
-
-  // go through all particle combinatorics, calculate the invariant mass
-  // for each combination, and remember which combination is the closest
-  // to the desired pdg_mass
-  for (int i = 0; i < parts.size(); ++i) {
-    for (int j = i + 1; j < parts.size(); ++j) {
-      const double new_mass{(parts[i] + parts[j]).mass()};
-      if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
-        first = i;
-        second = j;
-        best_mass = new_mass;
-      }
-    }
-  }
-  if (first < 0) {
-    return {{}, {}};
-  }
-  return {parts[first], parts[second]};
-}
-
-// Calculate the invariant mass of a given pair of particles
-inline double
-get_im(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
-           particle_pair) {
-  return (particle_pair.first + particle_pair.second).mass();
-}
-
-// Calculate the transverse momentum of a given pair of particles
-inline double
-get_pt(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
-           particle_pair) {
-  double px_pair = (particle_pair.first + particle_pair.second).px();
-  double py_pair = (particle_pair.first + particle_pair.second).py();  
-  return sqrt(px_pair*px_pair + py_pair*py_pair);
-}
-
-// Calculate the azimuthal angle of a given pair of particles
-inline double
-get_phi(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
-           particle_pair) {
-  double px_pair = (particle_pair.first + particle_pair.second).px();
-  double py_pair = (particle_pair.first + particle_pair.second).py();
-  double phi_pair = std::atan2(py_pair,px_pair);
-  //if(py_pair <= 0.) phi_pair = - phi_pair;
-  return phi_pair;
-}
-
-// Calculate the rapidity of a given pair of particles
-inline double
-get_y(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
-           particle_pair) {
-  double px_pair = (particle_pair.first + particle_pair.second).px();
-  double py_pair = (particle_pair.first + particle_pair.second).py();
-  double pz_pair = (particle_pair.first + particle_pair.second).pz();
-  double mass_pair = (particle_pair.first + particle_pair.second).mass();
-  double energy_pair = sqrt(mass_pair*mass_pair + px_pair*px_pair + py_pair*py_pair + pz_pair*pz_pair);
-  return 0.5*log((energy_pair + pz_pair)/(energy_pair - pz_pair));
-}
-
-//========================================================================================================
-//for structure functions
-
-struct inv_quant{    //add more when needed
-    double nu, Q2, x;
-};
-
-//for simu
-inline inv_quant calc_inv_quant_simu(const std::vector<ROOT::Math::PxPyPzMVector>& parts){
-  ROOT::Math::PxPyPzMVector q(parts[0] - parts[2]);
-  ROOT::Math::PxPyPzMVector P(parts[3]);
-  
-  double nu = q.Dot(P) / P.mass();
-  double Q2 = - q.Dot(q);  
-  inv_quant quantities = {nu, Q2, Q2/2./P.mass()/nu};
-  return quantities;
-}
-
-inline double get_nu_simu(inv_quant quantities) {
-  return quantities.nu/1000.;
-}
-inline double get_Q2_simu(inv_quant quantities) {
-  return quantities.Q2;
-}
-inline double get_x_simu(inv_quant quantities) {
-  return quantities.x;
-}
-
-//for tracking, add later
-
-//=========================================================================================================
-
-
-
-
-
-
-} // namespace util
-
-#endif
diff --git a/include/benchmark.h b/include/benchmark.h
new file mode 100644
index 0000000000000000000000000000000000000000..0634be0a1494739c4e7337f99e35ffbac0be7f7f
--- /dev/null
+++ b/include/benchmark.h
@@ -0,0 +1,124 @@
+#ifndef BENCHMARK_H
+#define BENCHMARK_H
+
+#include "exception.h"
+#include <fmt/core.h>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <nlohmann/json.hpp>
+#include <string>
+
+// Bookkeeping of test data to store data of one or more tests in a json file to
+// facilitate future accounting.
+//
+// Usage Example 1 (single test):
+// ==============================
+// 1. define our test
+//      eic::util::Test test1{
+//        {{"name", "example_test"},
+//        {"title", "Example Test"},
+//        {"description", "This is an example of a test definition"},
+//        {"quantity", "efficiency"},
+//        {"target", "1"}}};
+// 2. set pass/fail/error status and return value (in this case .99)
+//      test1.pass(0.99)
+// 3. write our test data to a json file
+//      eic::util::write_test(test1, "test1.json");
+//
+// Usage Example 2 (multiple tests):
+// =================================
+// 1. define our tests
+//      eic::util::Test test1{
+//        {{"name", "example_test"},
+//        {"title", "Example Test"},
+//        {"description", "This is an example of a test definition"},
+//        {"quantity", "efficiency"},
+//        {"target", "1"}}};
+//      eic::util::Test test2{
+//        {{"name", "another_test"},
+//        {"title", "Another example Test"},
+//        {"description", "This is a second example of a test definition"},
+//        {"quantity", "resolution"},
+//        {"target", "3."}}};
+// 2. set pass/fail/error status and return value (in this case .99)
+//      test1.fail(10)
+// 3. write our test data to a json file
+//      eic::util::write_test({test1, test2}, "test.json");
+
+// Namespace for utility scripts, FIXME this should be part of an independent
+// library
+namespace eic::util {
+
+  struct TestDefinitionError : Exception {
+    TestDefinitionError(std::string_view msg) : Exception(msg, "test_definition_error") {}
+  };
+
+  // Wrapper for our test data json, with three methods to set the status
+  // after test completion (to pass, fail or error). The default value
+  // is error.
+  // The following fields should be defined in the definitions json
+  // for the test to make sense:
+  //  - name: unique identifier for this test
+  //  - title: Slightly more verbose identifier for this test
+  //  - description: Concise description of what is tested
+  //  - quantity: What quantity is tested? Unites of value/target
+  //  - target: Target value of <quantity> that we want to reach
+  //  - value: Actual value of <quantity>
+  //  - weight: Weight for this test (this is defaulted to 1.0 if not specified)
+  //  - result: pass/fail/error
+  struct Test {
+    // note: round braces for the json constructor, as else it will pick the wrong
+    //       initializer-list constructur (it will put everything in an array)
+    Test(const std::map<std::string, std::string>& definition) : json(definition)
+    {
+      // std::cout << json.dump() << std::endl;
+      // initialize with error (as we don't have a value yet)
+      error();
+      // Check that all required fields are present
+      for (const auto& field :
+           {"name", "title", "description", "quantity", "target", "value", "result"}) {
+        if (json.find(field) == json.end()) {
+          throw TestDefinitionError{
+              fmt::format("Error in test definition: field '{}' missing", field)};
+        }
+      }
+      // Default "weight" to 1 if not set
+      if (json.find("weight") == json.end()) {
+        json["weight"] = 1.0;
+      }
+    }
+    // Set this test to pass/fail/error
+    void pass(double value) { update_result("pass", value); }
+    void fail(double value) { update_result("fail", value); }
+    void error(double value = 0) { update_result("error", value); }
+
+    nlohmann::json json;
+
+  private:
+    void update_result(std::string_view status, double value)
+    {
+      json["result"] = status;
+      json["value"]  = value;
+    }
+  };
+
+  void write_test(const std::vector<Test>& data, const std::string& fname)
+  {
+    nlohmann::json test;
+    for (auto& entry : data) {
+      test["tests"].push_back(entry.json);
+    }
+    std::cout << fmt::format("Writing test data to {}\n", fname);
+    std::ofstream output_file(fname);
+    output_file << std::setw(4) << test << "\n";
+  }
+  void write_test(const Test& data, const std::string& fname)
+  {
+    std::vector<Test> vtd{data};
+    write_test(vtd, fname);
+  }
+
+} // namespace eic::util
+
+#endif
diff --git a/include/exception.h b/include/exception.h
new file mode 100644
index 0000000000000000000000000000000000000000..b630bd008a0d8855211eca5ee0af60650aec1314
--- /dev/null
+++ b/include/exception.h
@@ -0,0 +1,22 @@
+#ifndef UTIL_EXCEPTION_H
+#define UTIL_EXCEPTION_H
+
+#include <exception>
+#include <string>
+
+namespace eic::util {
+  class Exception : public std::exception {
+  public:
+    Exception(std::string_view msg, std::string_view type = "exception") : msg_{msg}, type_{type} {}
+
+    virtual const char* what() const throw() { return msg_.c_str(); }
+    virtual const char* type() const throw() { return type_.c_str(); }
+    virtual ~Exception() throw() {}
+
+  private:
+    std::string msg_;
+    std::string type_;
+  };
+} // namespace eic::util
+
+#endif
diff --git a/dvmp/analysis/mt.h b/include/mt.h
similarity index 100%
rename from dvmp/analysis/mt.h
rename to include/mt.h
diff --git a/include/plot.h b/include/plot.h
new file mode 100644
index 0000000000000000000000000000000000000000..c198616325d54d73ae5e4c08f5ba7113ae77f817
--- /dev/null
+++ b/include/plot.h
@@ -0,0 +1,42 @@
+#ifndef PLOT_H
+#define PLOT_H
+
+#include <TCanvas.h>
+#include <TColor.h>
+#include <TPad.h>
+#include <TPaveText.h>
+#include <TStyle.h>
+#include <fmt/core.h>
+#include <vector>
+
+namespace plot {
+
+  const int kMpBlue   = TColor::GetColor(0x1f, 0x77, 0xb4);
+  const int kMpOrange = TColor::GetColor(0xff, 0x7f, 0x0e);
+  const int kMpGreen  = TColor::GetColor(0x2c, 0xa0, 0x2c);
+  const int kMpRed    = TColor::GetColor(0xd6, 0x27, 0x28);
+  const int kMpPurple = TColor::GetColor(0x94, 0x67, 0xbd);
+  const int kMpBrown  = TColor::GetColor(0x8c, 0x56, 0x4b);
+  const int kMpPink   = TColor::GetColor(0xe3, 0x77, 0xc2);
+  const int kMpGrey   = TColor::GetColor(0x7f, 0x7f, 0x7f);
+  const int kMpMoss   = TColor::GetColor(0xbc, 0xbd, 0x22);
+  const int kMpCyan   = TColor::GetColor(0x17, 0xbe, 0xcf);
+
+  const std::vector<int> kPalette = {kMpBlue,  kMpOrange, kMpGreen, kMpRed,  kMpPurple,
+                                     kMpBrown, kMpPink,   kMpGrey,  kMpMoss, kMpCyan};
+
+  void draw_label(int ebeam, int pbeam, const std::string_view detector)
+  {
+    auto t = new TPaveText(.15, 0.800, .7, .925, "NB NDC");
+    t->SetFillColorAlpha(kWhite, 0.4);
+    t->SetTextFont(43);
+    t->SetTextSize(25);
+    t->AddText(fmt::format("#bf{{{} }}SIMULATION", detector).c_str());
+    t->AddText(fmt::format("{} GeV on {} GeV", ebeam, pbeam).c_str());
+    t->SetTextAlign(12);
+    t->Draw();
+  }
+
+} // namespace plot
+
+#endif
diff --git a/include/util.h b/include/util.h
new file mode 100644
index 0000000000000000000000000000000000000000..56fb12893787e1d1fff00bd01653ad23b149f48e
--- /dev/null
+++ b/include/util.h
@@ -0,0 +1,159 @@
+#ifndef UTIL_H
+#define UTIL_H
+
+// TODO: should probably be moved to a global benchmark utility library
+
+#include <algorithm>
+#include <cmath>
+#include <exception>
+#include <fmt/core.h>
+#include <limits>
+#include <string>
+#include <vector>
+
+#include <Math/Vector4D.h>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/TrackParametersCollection.h"
+
+namespace util {
+
+  // Exception definition for unknown particle errors
+  // FIXME: A utility exception base class should be included in the analysis
+  //        utility library, so we can skip most of this boilerplate
+  class unknown_particle_error : public std::exception {
+  public:
+    unknown_particle_error(std::string_view particle) : m_particle{particle} {}
+    virtual const char* what() const throw()
+    {
+      return fmt::format("Unknown particle type: {}", m_particle).c_str();
+    }
+    virtual const char* type() const throw() { return "unknown_particle_error"; }
+
+  private:
+    const std::string m_particle;
+  };
+
+  // Simple function to return the appropriate PDG mass for the particles
+  // we care about for this process.
+  // FIXME: consider something more robust (maybe based on hepPDT) to the
+  //        analysis utility library
+  inline double get_pdg_mass(std::string_view part)
+  {
+    if (part == "electron") {
+      return 0.0005109989461;
+    } else if (part == "muon") {
+      return .1056583745;
+    } else if (part == "jpsi") {
+      return 3.0969;
+    } else if (part == "upsilon") {
+      return 9.49630;
+    } else {
+      throw unknown_particle_error{part};
+    }
+  }
+
+  // Get a vector of 4-momenta from raw tracking info, using an externally
+  // provided particle mass assumption.
+  inline auto momenta_from_tracking(const std::vector<eic::TrackParametersData>& tracks,
+                                    const double                                 mass)
+  {
+    std::vector<ROOT::Math::PxPyPzMVector> momenta{tracks.size()};
+    // transform our raw tracker info into proper 4-momenta
+    std::transform(tracks.begin(), tracks.end(), momenta.begin(), [mass](const auto& track) {
+      // make sure we don't divide by zero
+      if (fabs(track.qOverP) < 1e-9) {
+        return ROOT::Math::PxPyPzMVector{};
+      }
+      const double p  = fabs(1. / track.qOverP);
+      const double px = p * cos(track.phi) * sin(track.theta);
+      const double py = p * sin(track.phi) * sin(track.theta);
+      const double pz = p * cos(track.theta);
+      return ROOT::Math::PxPyPzMVector{px, py, pz, mass};
+    });
+    return momenta;
+  }
+
+  // Get a vector of 4-momenta from the simulation data.
+  // TODO: Add PID selector (maybe using ranges?)
+  inline auto momenta_from_simulation(const std::vector<dd4pod::Geant4ParticleData>& parts)
+  {
+    std::vector<ROOT::Math::PxPyPzMVector> momenta{parts.size()};
+    // transform our simulation particle data into 4-momenta
+    std::transform(parts.begin(), parts.end(), momenta.begin(), [](const auto& part) {
+      return ROOT::Math::PxPyPzMVector{part.psx, part.psy, part.psz, part.mass};
+    });
+    return momenta;
+  }
+
+  // Find the decay pair candidates from a vector of particles (parts),
+  // with invariant mass closest to a desired value (pdg_mass)
+  inline std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>
+  find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass)
+  {
+    int    first     = -1;
+    int    second    = -1;
+    double best_mass = -1;
+
+    // go through all particle combinatorics, calculate the invariant mass
+    // for each combination, and remember which combination is the closest
+    // to the desired pdg_mass
+    for (size_t i = 0; i < parts.size(); ++i) {
+      for (size_t j = i + 1; j < parts.size(); ++j) {
+        const double new_mass{(parts[i] + parts[j]).mass()};
+        if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
+          first     = i;
+          second    = j;
+          best_mass = new_mass;
+        }
+      }
+    }
+    if (first < 0) {
+      return {{}, {}};
+    }
+    return {parts[first], parts[second]};
+  }
+
+  // Calculate the magnitude of the momentum of a vector of 4-vectors
+  inline auto mom(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
+  {
+    std::vector<double> P(momenta.size());
+    // transform our raw tracker info into proper 4-momenta
+    std::transform(momenta.begin(), momenta.end(), P.begin(),
+                   [](const auto& mom) { return mom.P(); });
+    return P;
+  }
+  // Calculate the transverse momentum of a vector of 4-vectors
+  inline auto pt(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
+  {
+    std::vector<double> pt(momenta.size());
+    // transform our raw tracker info into proper 4-momenta
+    std::transform(momenta.begin(), momenta.end(), pt.begin(),
+                   [](const auto& mom) { return mom.pt(); });
+    return pt;
+  }
+
+  // Calculate the azimuthal angle phi of a vector of 4-vectors
+  inline auto phi(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
+  {
+    std::vector<double> phi(momenta.size());
+    // transform our raw tracker info into proper 4-momenta
+    std::transform(momenta.begin(), momenta.end(), phi.begin(),
+                   [](const auto& mom) { return mom.phi(); });
+    return phi;
+  }
+  // Calculate the pseudo-rapidity of a vector of particles
+  inline auto eta(const std::vector<ROOT::Math::PxPyPzMVector>& momenta)
+  {
+    std::vector<double> eta(momenta.size());
+    // transform our raw tracker info into proper 4-momenta
+    std::transform(momenta.begin(), momenta.end(), eta.begin(),
+                   [](const auto& mom) { return mom.eta(); });
+    return eta;
+  }
+
+  //=========================================================================================================
+
+} // namespace util
+
+#endif
diff --git a/config/env.sh b/options/env.sh
similarity index 96%
rename from config/env.sh
rename to options/env.sh
index 5abe41947886239dabd08153fe89e1353e1d26fa..542ee7dec4fe387a7e13069e16ecb53972e95fa7 100755
--- a/config/env.sh
+++ b/options/env.sh
@@ -87,6 +87,10 @@ echo "DETECTOR_PREFIX:        ${DETECTOR_PREFIX}"
 export DETECTOR_PATH="${DETECTOR_PREFIX}/${JUGGLER_DETECTOR}"
 echo "DETECTOR_PATH:          ${DETECTOR_PATH}"
 
+## build dir for ROOT to put its binaries etc.
+export ROOT_BUILD_DIR=$LOCAL_PREFIX/root_build
+echo "ROOT_BUILD_DIR:         ${ROOT_BUILD_DIR}"
+
 ## =============================================================================
 ## Setup PATH and LD_LIBRARY_PATH to include our prefixes
 echo "Adding JUGGLER_INSTALL_PREFIX and LOCAL_PREFIX to PATH and LD_LIBRARY_PATH"
diff --git a/pythia_dis b/pythia_dis
new file mode 100755
index 0000000000000000000000000000000000000000..91a83f279189632193eaa87181405aaad17b8378
Binary files /dev/null and b/pythia_dis differ
diff --git a/tools/start_dev_shell.sh b/tools/dev-shell
similarity index 100%
rename from tools/start_dev_shell.sh
rename to tools/dev-shell
diff --git a/util/benchmark.hh b/util/benchmark.hh
deleted file mode 100644
index f8abec063205d43beba254d1af3a919c387d2acd..0000000000000000000000000000000000000000
--- a/util/benchmark.hh
+++ /dev/null
@@ -1,117 +0,0 @@
-#ifndef BENCHMARK_LOADED
-#define BENCHMARK_LOADED
-
-#include "exception.hh"
-#include <fmt/core.h>
-#include <fstream>
-#include <iostream>
-#include <nlohmann/json.hpp>
-#include <string>
-
-// Bookkeeping of test data to store data of one or more tests in a json file to
-// facilitate future accounting.
-//
-// Usage Example 1 (single test):
-// ==============================
-// 1. define our test
-//      eic::util::Test test1{
-//        {{"name", "example_test"},
-//        {"title", "Example Test"},
-//        {"description", "This is an example of a test definition"},
-//        {"quantity", "efficiency"},
-//        {"target", "1"}}};
-// 2. set pass/fail/error status and return value (in this case .99)
-//      test1.pass(0.99)
-// 3. write our test data to a json file
-//      eic::util::write_test(test1, "test1.json");
-//
-// Usage Example 2 (multiple tests):
-// =================================
-// 1. define our tests
-//      eic::util::Test test1{
-//        {{"name", "example_test"},
-//        {"title", "Example Test"},
-//        {"description", "This is an example of a test definition"},
-//        {"quantity", "efficiency"},
-//        {"target", "1"}}};
-//      eic::util::Test test2{
-//        {{"name", "another_test"},
-//        {"title", "Another example Test"},
-//        {"description", "This is a second example of a test definition"},
-//        {"quantity", "resolution"},
-//        {"target", "3."}}};
-// 2. set pass/fail/error status and return value (in this case .99)
-//      test1.fail(10)
-// 3. write our test data to a json file
-//      eic::util::write_test({test1, test2}, "test.json");
-
-// Namespace for utility scripts, FIXME this should be part of an independent
-// library
-namespace eic::util {
-
-struct TestDefinitionError : Exception {
-  TestDefinitionError(std::string_view msg)
-      : Exception(msg, "test_definition_error") {}
-};
-
-// Wrapper for our test data json, with three methods to set the status
-// after test completion (to pass, fail or error). The default value
-// is error.
-// The following fields should be defined in the definitions json
-// for the test to make sense:
-//  - name: unique identifier for this test
-//  - title: Slightly more verbose identifier for this test
-//  - description: Concise description of what is tested
-//  - quantity: What quantity is tested? Unites of value/target
-//  - target: Target value of <quantity> that we want to reach
-//  - value: Actual value of <quantity>
-//  - weight: Weight for this test (this is defaulted to 1.0 if not specified)
-//  - result: pass/fail/error
-struct Test {
-  Test(nlohmann::json definition) : json{std::move(definition)} {
-    // initialize with error (as we don't have a value yet)
-    error();
-    // Check that all required fields are present
-    for (const auto& field : {"name", "title", "description", "quantity",
-                              "target", "value", "result"}) {
-      if (json.find(field) == json.end()) {
-        throw TestDefinitionError{
-            fmt::format("Error in test definition: field '{}' missing", field)};
-      }
-    }
-    // Default "weight" to 1 if not set
-    if (json.find("weight") == json.end()) {
-      json["weight"] = 1.0;
-    }
-  }
-  // Set this test to pass/fail/error
-  void pass(double value) { update_result("pass", value); }
-  void fail(double value) { update_result("fail", value); }
-  void error(double value = 0) { update_result("error", value); }
-
-  nlohmann::json json;
-
-private:
-  void update_result(std::string_view status, double value) {
-    json["result"] = status;
-    json["value"] = value;
-  }
-};
-
-void write_test(const std::vector<Test>& data, const std::string& fname) {
-  nlohmann::json test;
-  for (auto& entry : data) {
-    test["tests"].push_back(entry.json);
-  }
-  std::cout << fmt::format("Writing test data to {}\n", fname);
-  std::ofstream output_file(fname);
-  output_file << std::setw(4) << test << "\n";
-}
-void write_test(const Test& data, const std::string& fname) {
-  std::vector<Test> vtd{data};
-  write_test(vtd, fname);
-}
-
-} // namespace eic::util
-
-#endif
diff --git a/util/build_detector.sh b/util/build_detector.sh
index 03641f48151b95011700f80ffd1c0814417b47aa..02313585d56e0ca28af161bc19a7f7b3ae5f53c7 100755
--- a/util/build_detector.sh
+++ b/util/build_detector.sh
@@ -18,9 +18,9 @@ pushd ${PROJECT_ROOT}
 ## - DETECTOR_PATH:    full path for the detector definitions
 ##                     this is the same as ${DETECTOR_PREFIX}/${JUGGLER_DETECTOR}
 ##
-## You can read config/env.sh for more in-depth explanations of the variables
+## You can read options/env.sh for more in-depth explanations of the variables
 ## and how they can be controlled.
-source config/env.sh
+source options/env.sh
 
 ## =============================================================================
 ## Step 1: download/update the detector definitions (if needed)
diff --git a/util/collect_benchmarks.py b/util/collect_benchmarks.py
index 8f9fd675e5d824580788a13690596515d333b988..0af7e9a12b37eb7616be3400d4046f0088bcf223 100755
--- a/util/collect_benchmarks.py
+++ b/util/collect_benchmarks.py
@@ -9,7 +9,7 @@ directory.
 """
 
 ## Our master definition file, the benchmark project directory
-MASTER_FILE=r'benchmarks.json'
+MASTER_FILE=r'benchmarks/benchmarks.json'
 
 ## Our results directory
 RESULTS_PATH=r'results'
diff --git a/util/collect_tests.py b/util/collect_tests.py
index c56d8c8b2ffd6a1f7ef525c88f66ed6614bb8857..4d860ca79d9f996204a5ca9dc447fa10ed8ec4f4 100755
--- a/util/collect_tests.py
+++ b/util/collect_tests.py
@@ -15,7 +15,7 @@ files to identify them as benchmark components.
 """
 
 ## Our benchmark definition file, stored in the benchmark root directory
-BENCHMARK_FILE=r'{}/benchmark.json'
+BENCHMARK_FILE=r'benchmarks/{}/benchmark.json'
 
 ## Our benchmark results directory
 RESULTS_PATH=r'results/{}'
diff --git a/util/compile_analyses.py b/util/compile_analyses.py
new file mode 100755
index 0000000000000000000000000000000000000000..153f2ea2f61429b6864a21b3ad625e6b53373ba3
--- /dev/null
+++ b/util/compile_analyses.py
@@ -0,0 +1,119 @@
+#!/usr/bin/env python3
+
+"""
+Compile all root analysis scripts under
+benchmarks/<BENCHMARK>/analysis/*.cxx
+
+Doing this step here rather than during the main benchmark script has
+multiple advantages:
+    1. Get feedback on syntax errors early on, without wasting compute resources
+    2. Avoid race conditions for large benchmarks run in parallel
+    3. Make it easier to properly handle the root build directory, as
+       this has to exist prior to our attempt to compile, else all will
+       fail (this is probably an old bug in root...)
+
+Analysis scripts are expected to have extension 'cxx' and be located in the analysis
+subdirectory
+"""
+
+## Our analysis path and file extension for glob
+ANALYSIS_PATH=r'benchmarks/{}/analysis'
+ANALYSIS_EXT = r'cxx'
+
+import argparse
+import os
+from pathlib import Path
+
+## Exceptions for this module
+class Error(Exception):
+    '''Base class for exceptions in this module.'''
+    pass
+
+class PathNotFoundError(Exception):
+    '''Path does not exist.
+
+    Attributes:
+        path: the path name
+        message: error message
+    '''
+    def __init__(self, path):
+        self.file = file
+        self.message = 'No such directory: {}'.format(file)
+class NoAnalysesFoundError(Exception):
+    '''Did not find any analysis scripts to complile
+
+    Attributes:
+        path: the analysis path
+        message: error message
+    '''
+    def __init__(self, path):
+        self.file = file
+        self.message = 'No analysis found (extension \'{}\' in path: {}'.format(file,
+                ANALYSIS_EXT)
+
+class CompilationError(Exception):
+    '''Raised when we failed to compile an analysis script
+
+    Attributes:
+        file: analysis file name
+        path: analysis path
+        message: error message
+    '''
+    def __init__(self, file):
+        self.file = file
+        self.message = "Analysis '{}' failed to compile".format(file)
+
+parser = argparse.ArgumentParser()
+parser.add_argument(
+        'benchmark',
+        help='A benchmarks for which to compile the analysis scripts.')
+
+def compile_analyses(benchmark):
+    '''Compile all analysis scripts for a benchmark.'''
+    print("Compiling all analyis scripts for '{}'".format(benchmark))
+
+    ## Ensure our build directory exists
+    _init_build_dir(benchmark)
+
+    ## Get a list of all analysis scripts
+    _compile_all(benchmark)
+
+    ## All done!
+    print('All analyses for', benchmark, 'compiled successfully')
+
+def _init_build_dir(benchmark):
+    '''Initialize our ROOT build directory (if using one).'''
+    print(' --> Initializing ROOT build directory ...')
+    build_prefix = os.getenv('ROOT_BUILD_DIR')
+    if build_prefix is None:
+        print('    --> ROOT_BUILD_DIR not set, no action needed.')
+        return
+    ## deduce the root build directory
+    pwd = os.getenv('PWD')
+    build_dir = '{}/{}/{}'.format(build_prefix, pwd, ANALYSIS_PATH.format(benchmark))
+    print("    --> Ensuring directory '{}' exists".format(build_dir))
+    os.system('mkdir -p {}'.format(build_dir))
+
+def _compile_all(benchmark):
+    '''Compile all analysis for this benchmark.'''
+    print(' --> Compiling analysis scripts')
+    anadir = Path(ANALYSIS_PATH.format(benchmark))
+    if not anadir.exists():
+        raise PathNotFoundError(anadir)
+    ana_list = []
+    for file in anadir.glob('*.{}'.format(ANALYSIS_EXT)):
+        ana_list.append(file)
+        print('    --> Compiling:', file, flush=True)
+        err = os.system(_compile_cmd(file))
+        if err:
+            raise CompilationError(file)
+    if len(ana_list) == 0:
+        raise NoAnalysesFoundError(anadir)
+
+def _compile_cmd(file):
+    '''Return a one-line shell command to compile an analysis script.'''
+    return r'bash -c "root -q -b -e \".L {}+\""'.format(file)
+
+if __name__ == "__main__":
+    args = parser.parse_args()
+    compile_analyses(args.benchmark)
diff --git a/util/exception.hh b/util/exception.hh
deleted file mode 100644
index d97ef27cd2f1e3ed4237d489af4cf0b742d88143..0000000000000000000000000000000000000000
--- a/util/exception.hh
+++ /dev/null
@@ -1,23 +0,0 @@
-#ifndef UTIL_EXCEPTION
-#define UTIL_EXCEPTION
-
-#include <exception>
-#include <string>
-
-namespace eic::util {
-class Exception : public std::exception {
-public:
-  Exception(std::string_view msg, std::string_view type = "exception")
-      : msg_{msg}, type_{type} {}
-
-  virtual const char* what() const throw() { return msg_.c_str(); }
-  virtual const char* type() const throw() { return type_.c_str(); }
-  virtual ~Exception() throw() {}
-
-private:
-  std::string msg_;
-  std::string type_;
-};
-} // namespace eic::util
-
-#endif
diff --git a/util/parse_cmd.sh b/util/parse_cmd.sh
index 9bf2f55319c4cf13d8dd73165580e148e94186dd..04028d8958d03df22f7b2d964a75acf2d9b47317 100755
--- a/util/parse_cmd.sh
+++ b/util/parse_cmd.sh
@@ -11,9 +11,11 @@
 ##   - REQUIRE_DECAY:     require the --decay flag to be set
 ## =============================================================================
 
+## Commented out because this should be taken care of by the 
+## calling script to not enforce a fixed directory structure.
 ## make sure we launch this script from the project root directory
-PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/..
-pushd ${PROJECT_ROOT}
+#PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/..
+#pushd ${PROJECT_ROOT}
 
 ## =============================================================================
 ## Step 1: Process the command line arguments