diff --git a/benchmarks/dvcs/scripts/dvcs_ps_gen.cxx b/benchmarks/dvcs/analysis/dvcs_ps_gen.cxx
similarity index 91%
rename from benchmarks/dvcs/scripts/dvcs_ps_gen.cxx
rename to benchmarks/dvcs/analysis/dvcs_ps_gen.cxx
index c0642a7c05cd487a5668741ed00ca997d26c782f..1631d930fb096f99a3c121959a3c082bd1fa3ab0 100644
--- a/benchmarks/dvcs/scripts/dvcs_ps_gen.cxx
+++ b/benchmarks/dvcs/analysis/dvcs_ps_gen.cxx
@@ -1,3 +1,9 @@
+#include "TH2F.h"
+#include "TLorentzVector.h"
+#include "TGenPhaseSpace.h"
+
+#include <iostream>
+
 void dvcs_ps_gen() {
   double E_p = 100.0;
   double M_p = 0.938;
diff --git a/benchmarks/dvcs/analysis/dvcs_tests.cxx b/benchmarks/dvcs/analysis/dvcs_tests.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..43695195dc12a7b05058d4666d0a2d00b9292351
--- /dev/null
+++ b/benchmarks/dvcs/analysis/dvcs_tests.cxx
@@ -0,0 +1,54 @@
+#include <cmath>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "ROOT/RDataFrame.hxx"
+#include "Math/Vector4D.h"
+#include "TCanvas.h"
+
+#include <nlohmann/json.hpp>
+using json = nlohmann::json;
+
+R__LOAD_LIBRARY(libfmt.so)
+#include "fmt/core.h"
+#include "fmt/color.h"
+
+R__LOAD_LIBRARY(libeicd.so)
+R__LOAD_LIBRARY(libDD4pod.so)
+
+#include "eicd/InclusiveKinematicsCollection.h"
+#include "eicd/ReconstructedParticleCollection.h"
+
+void dvcs_tests(const char* fname = "rec_dvcs.root"){
+
+  fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green), "Running DVCS analysis...\n");
+
+  // Run this in multi-threaded mode if desired
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame df("events", fname);
+
+  auto df0 = df.Define("n_parts", "ReconstructedParticles.size()")
+               .Define("isQ2gt1", "InclusiveKinematicsTruth.Q2 > 1.0")
+               .Define("n_Q2gt1", "isQ2gt1.size()");
+
+  auto h_n_parts = df0.Histo1D({"h_n_parts", "; h_n_parts n", 10, 0, 10}, "n_parts");
+  auto h_Q2      = df0.Histo1D({"h_Q2", "; Q^{2} [GeV^{2}/c^{2}]", 100, 0, 30}, "InclusiveKinematicsTruth.Q2");
+  auto n_Q2gt1   = df0.Mean("n_Q2gt1");
+  auto n_parts   = df0.Mean("n_parts");
+
+  // ---------------------------
+  // Do evaluation
+
+  auto c = new TCanvas();
+  h_Q2->DrawCopy();
+  c->SaveAs("results/dvcs/Q2.png");
+  c->SaveAs("results/dvcs/Q2.pdf");
+  fmt::print("{} DVCS events Q2>1\n",*n_Q2gt1);
+  fmt::print("{} tracks per event\n",*n_parts);
+
+  c = new TCanvas();
+  h_n_parts->DrawCopy();
+  c->SaveAs("results/dvcs/n_parts.png");
+  c->SaveAs("results/dvcs/n_parts.pdf");
+}
diff --git a/benchmarks/dvcs/config.yml b/benchmarks/dvcs/config.yml
index ee77a0f3b6498bfc30da81457fbff4d7598a9d9e..9b00dfebafd6f75f413700ef4f9d7565c1bf3957 100644
--- a/benchmarks/dvcs/config.yml
+++ b/benchmarks/dvcs/config.yml
@@ -5,6 +5,7 @@ dvcs:process:
     - s3
   needs: ["common:detector"]
   script:
+    - compile_analyses.py dvcs
     - bash benchmarks/dvcs/dvcs.sh --all
 
 dvcs:results:
diff --git a/benchmarks/dvcs/dvcs.sh b/benchmarks/dvcs/dvcs.sh
index 7a47ac14b2c9522f08454dc4a84f4420bd52d050..c05c3486b292eb2aa2acb8d6809777aaf9cc50ba 100644
--- a/benchmarks/dvcs/dvcs.sh
+++ b/benchmarks/dvcs/dvcs.sh
@@ -145,7 +145,7 @@ if [[ -n "${DO_ANALYSIS}" || -n "${DO_ALL}" ]] ; then
   mkdir -p results/dvcs
 
   # here you can add as many scripts as you want.
-  root -b -q "benchmarks/dvcs/scripts/dvcs_tests.cxx(\"${JUGGLER_REC_FILE}\")"
+  root -b -q "benchmarks/dvcs/analysis/dvcs_tests.cxx+(\"${JUGGLER_REC_FILE}\")"
   if [[ "$?" -ne "0" ]] ; then
     echo "ERROR running root script"
     exit 1
diff --git a/benchmarks/dvcs/scripts/dvcs_tests.cxx b/benchmarks/dvcs/scripts/dvcs_tests.cxx
deleted file mode 100644
index 2ebf9f45f73b91baf65fa56eb9b111da7025e7dc..0000000000000000000000000000000000000000
--- a/benchmarks/dvcs/scripts/dvcs_tests.cxx
+++ /dev/null
@@ -1,145 +0,0 @@
-#include <cmath>
-#include <iostream>
-#include <string>
-#include <vector>
-
-#include "ROOT/RDataFrame.hxx"
-
-#include <nlohmann/json.hpp>
-using json = nlohmann::json;
-
-R__LOAD_LIBRARY(libfmt.so)
-#include "fmt/core.h"
-#include "fmt/color.h"
-
-R__LOAD_LIBRARY(libeicd.so)
-R__LOAD_LIBRARY(libDD4pod.so)
-
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "eicd/TrackParametersCollection.h"
-#include "eicd/ClusterCollection.h"
-#include "eicd/ReconstructedParticleCollection.h"
-#include "eicd/BasicParticleCollection.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;
-};
-
-
-auto 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].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y));
-  }
-  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].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
-    result.push_back(lv);
-  }
-  return result;
-};
-auto recfourvec = [](ROOT::VecOps::RVec<eic::ReconstructedParticleData> 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].p.x, in[i].p.y, in[i].p.z, 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;
-};
-
-
-void dvcs_tests(const char* fname = "rec_dvcs.root"){
-
-  fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green), "Running DVCS analysis...\n");
-
-  // Run this in multi-threaded mode if desired
-  ROOT::EnableImplicitMT();
-  ROOT::RDataFrame df("events", fname);
-
-  using ROOT::Math::PxPyPzMVector;
-  PxPyPzMVector p_ebeam = {0,0,-10, 0.000511};
-  PxPyPzMVector p_pbeam = {0,0,275,  0.938 };
-
-  auto eprime = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
-    for(const auto& p : in){
-      if(p.pdgID == 11 ) {
-        return PxPyPzMVector(p.ps.x,p.ps.y,p.ps.z,p.mass);
-      }
-    }
-    return PxPyPzMVector(0,0,0,0);
-  };
-  auto q_vec = [=](PxPyPzMVector const& p) {
-    return p_ebeam - p;
-  };
-
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
-                 .Define("thrownP", fourvec, {"thrownParticles"})
-                 .Define("recP", recfourvec, {"ReconstructedParticles"})
-                 .Define("NPart", "recP.size()")
-                 .Define("p_thrown", momentum, {"thrownP"})
-                 .Define("nTracks", "outputTrackParameters.size()")
-                 .Define("p_track", p_track, {"outputTrackParameters"})
-                 .Define("delta_p",delta_p, {"p_track", "p_thrown"})
-                 .Define("eprime", eprime, {"thrownParticles"})
-                 .Define("q",  q_vec, {"eprime"})
-                 .Define("Q2", "-1.0*(q.Dot(q))");
-
-  auto h_n_dummy = df0.Histo1D({"h_n_part", "; h_n_part n", 10, 0, 10}, "NPart");
-  auto h_Q2      = df0.Histo1D({"h_Q2", "; Q^{2} [GeV^{2}/c^{2}]", 100, 0, 30}, "Q2");
-  auto n_Q2      = df0.Filter("Q2>1").Count();
-  auto n_tracks  = df0.Mean("nTracks");
-
-  // ---------------------------
-  // Do evaluation
-
-  auto c = new TCanvas();
-  h_Q2->DrawCopy();
-  c->SaveAs("results/dvcs/Q2.png");
-  c->SaveAs("results/dvcs/Q2.pdf");
-  fmt::print("{} DVCS events\n",*n_Q2);
-  fmt::print("{} tracks per event\n",*n_tracks);
-
-  c = new TCanvas();
-  h_n_dummy->DrawCopy();
-  c->SaveAs("results/dvcs/n_dummy.png");
-  //c->SaveAs("results/dvcs/n_dummy.pdf");
-
-}
diff --git a/benchmarks/single/config.yml b/benchmarks/single/config.yml
index 35ed2a0d046972fd4ccf7a18b60d9a334c79d155..a84da3f88880307899bfe1f16df1668cfe30bb4f 100644
--- a/benchmarks/single/config.yml
+++ b/benchmarks/single/config.yml
@@ -3,4 +3,5 @@ single:process:
   timeout: 24 hours
   stage: process
   script:
+    - compile_analyses.py single
     - bash benchmarks/single/single.sh e-_1GeV_45to135deg
diff --git a/benchmarks/u_omega/scripts/demo.cxx b/benchmarks/u_omega/analysis/demo.cxx
similarity index 74%
rename from benchmarks/u_omega/scripts/demo.cxx
rename to benchmarks/u_omega/analysis/demo.cxx
index afdc0d9a110e447a0bac7dab499015cf0d9ce34a..ad63c8959a2eace4e18c283bd3f16ec7c76a536e 100644
--- a/benchmarks/u_omega/scripts/demo.cxx
+++ b/benchmarks/u_omega/analysis/demo.cxx
@@ -4,6 +4,8 @@
 #include <vector>
 
 #include "ROOT/RDataFrame.hxx"
+#include "Math/Vector4D.h"
+#include "TCanvas.h"
 
 #include <nlohmann/json.hpp>
 using json = nlohmann::json;
@@ -20,6 +22,7 @@ R__LOAD_LIBRARY(libDD4pod.so)
 #include "eicd/ClusterCollection.h"
 #include "eicd/ReconstructedParticleCollection.h"
 #include "eicd/BasicParticleCollection.h"
+#include "eicd/InclusiveKinematicsCollection.h"
 
 using ROOT::RDataFrame;
 using namespace ROOT::VecOps;
@@ -109,23 +112,15 @@ void demo(const char* fname = "rec_dvcs.root"){
     return p_ebeam - p;
   };
 
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
-                 .Define("thrownP", fourvec, {"thrownParticles"})
-                 .Define("recP", recfourvec, {"ReconstructedParticles"})
-                 .Define("NPart", "recP.size()")
-                 .Define("p_thrown", momentum, {"thrownP"})
-                 .Define("nTracks", "outputTrackParameters.size()")
-                 .Define("p_track", p_track, {"outputTrackParameters"})
-                 .Define("delta_p",delta_p, {"p_track", "p_thrown"})
-                 .Define("eprime", eprime, {"thrownParticles"})
-                 .Define("q",  q_vec, {"eprime"})
-                 .Define("Q2", "-1.0*(q.Dot(q))");
-
-  auto h_n_dummy = df0.Histo1D({"h_n_part", "; h_n_part n", 10, 0, 10}, "NPart");
-  auto h_Q2      = df0.Histo1D({"h_Q2", "; Q^{2} [GeV^{2}/c^{2}]", 100, 0, 30}, "Q2");
-  auto n_Q2      = df0.Filter("Q2>1").Count();
-  auto n_tracks  = df0.Mean("nTracks");
+  auto df0 = df.Define("n_parts", "ReconstructedParticles.size()")
+               .Define("isQ2gt1", "InclusiveKinematicsTruth.Q2 > 1.0")
+               .Define("n_Q2gt1", "isQ2gt1.size()");
+
+
+  auto h_n_parts = df0.Histo1D({"h_n_parts", "; h_n_parts n", 10, 0, 10}, "n_parts");
+  auto h_Q2      = df0.Histo1D({"h_Q2", "; Q^{2} [GeV^{2}/c^{2}]", 100, 0, 30}, "InclusiveKinematicsTruth.Q2");
+  auto n_Q2gt1   = df0.Mean("n_Q2gt1");
+  auto n_parts   = df0.Mean("n_parts");
 
   // ---------------------------
   // Do evaluation
@@ -134,12 +129,12 @@ void demo(const char* fname = "rec_dvcs.root"){
   h_Q2->DrawCopy();
   c->SaveAs("results/u_omega/Q2.png");
   c->SaveAs("results/u_omega/Q2.pdf");
-  fmt::print("{} u_omega events\n",*n_Q2);
-  fmt::print("{} tracks per event\n",*n_tracks);
+  fmt::print("{} u_omega events Q2>1\n",*n_Q2gt1);
+  fmt::print("{} tracks per event\n",*n_parts);
 
   c = new TCanvas();
-  h_n_dummy->DrawCopy();
-  c->SaveAs("results/u_omega/n_dummy.png");
-  //c->SaveAs("results/u_omega/n_dummy.pdf");
+  h_n_parts->DrawCopy();
+  c->SaveAs("results/u_omega/n_parts.png");
+  c->SaveAs("results/u_omega/n_parts.pdf");
 
 }
diff --git a/benchmarks/u_omega/config.yml b/benchmarks/u_omega/config.yml
index 9e4925b6d431097afec2363d590e7c2aeb26480a..40d5f4be7fc677d18f10492de31b484ace90d422 100644
--- a/benchmarks/u_omega/config.yml
+++ b/benchmarks/u_omega/config.yml
@@ -5,6 +5,7 @@ u_omega:process:
     - s3
   needs: ["common:detector"]
   script:
+    - compile_analyses.py u_omega
     - bash benchmarks/u_omega/u_omega.sh --all
 
 u_omega:results:
diff --git a/benchmarks/u_omega/u_omega.sh b/benchmarks/u_omega/u_omega.sh
index 6d051589e1851b77ce0485ee72e0a4929fc41d89..23bbf403af8899ab7cd9481db0a811d7967cbf1a 100644
--- a/benchmarks/u_omega/u_omega.sh
+++ b/benchmarks/u_omega/u_omega.sh
@@ -144,7 +144,7 @@ if [[ -n "${DO_ANALYSIS}" || -n "${DO_ALL}" ]] ; then
   mkdir -p "results/${FILE_NAME_TAG}"
 
   # here you can add as many scripts as you want.
-  root -b -q "benchmarks/${FILE_NAME_TAG}/scripts/demo.cxx(\"${JUGGLER_REC_FILE}\")"
+  root -b -q "benchmarks/${FILE_NAME_TAG}/analysis/demo.cxx+(\"${JUGGLER_REC_FILE}\")"
   if [[ "$?" -ne "0" ]] ; then
     echo "ERROR running root script"
     exit 1
diff --git a/options/reconstruction.py b/options/reconstruction.py
index 66b449dcee6216d84dbd4143fc70df29ea70973c..927a8ffccca7fb7d1e3a3ae4f4145e302a70996a 100644
--- a/options/reconstruction.py
+++ b/options/reconstruction.py
@@ -62,6 +62,7 @@ from Configurables import Jug__Fast__SmearedFarForwardParticles as FFSmearedPart
 from Configurables import Jug__Fast__MatchClusters as MatchClusters
 from Configurables import Jug__Fast__ClusterMerger as ClusterMerger
 from Configurables import Jug__Fast__TruthEnergyPositionClusterMerger as EnergyPositionClusterMerger
+from Configurables import Jug__Fast__InclusiveKinematicsTruth as InclusiveKinematicsTruth
 
 from Configurables import Jug__Digi__PhotoMultiplierDigi as PhotoMultiplierDigi
 from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
@@ -129,6 +130,13 @@ dummy = MC2DummyParticle("dummy",
         smearing=0)
 algorithms.append(dummy)
 
+# Truth level kinematics
+truth_incl_kin = InclusiveKinematicsTruth("truth_incl_kin",
+        inputMCParticles="mcparticles",
+        outputData="InclusiveKinematicsTruth"
+)
+algorithms.append(truth_incl_kin)
+
 # Crystal Endcap Ecal
 ce_ecal_daq = dict(
         dynamicRangeADC=5.*units.GeV,