diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 9c744b9a393ac4af34ae04d99d8aef509ba120f9..956e7ddfc08701b43031229535276ede0c29d86f 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -14,6 +14,7 @@ default:
       - .local/bin
       - .local/include
       - .local/share
+      - .local/root_build
       - results
       - config
       - juggler.env
@@ -23,6 +24,7 @@ default:
 stages:
   - config
   - initialize
+  - compile
   - run
   - benchmarks1
   - benchmarks2
@@ -52,6 +54,12 @@ common:detector:
     - mkdir -p config
     - print_env.sh
 
+.compile_benchmark:
+  needs:
+    - ["common:detector"]
+  before_script:
+    - source .local/bin/env.sh
+
 .rec_benchmark:
   needs:
     - ["common:detector"]
diff --git a/benchmarks/clustering/full_cal_clusters.sh b/benchmarks/clustering/full_cal_clusters.sh
index 301373594a0ecb611433d156f0a0b89c2dd6bebf..9d7f8864a36a309ae7379e72e49299013ca873f0 100644
--- a/benchmarks/clustering/full_cal_clusters.sh
+++ b/benchmarks/clustering/full_cal_clusters.sh
@@ -79,7 +79,7 @@ compact_path=${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml
 export JUGGLER_FILE_NAME_TAG="${nametag}"
 export JUGGLER_GEN_FILE="gen_${JUGGLER_FILE_NAME_TAG}.hepmc"
 
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
 export JUGGLER_DIGI_FILE="digi_${JUGGLER_FILE_NAME_TAG}.root"
 export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
 
@@ -100,9 +100,10 @@ fi
 ls -lh ${JUGGLER_GEN_FILE}
 
 # Run geant4 simulations
-npsim --runType batch \
+ddsim --runType batch \
       -v WARNING \
       --part.minimalKineticEnergy "1*TeV" \
+      --filter.tracker edep0 \
       --numberOfEvents ${JUGGLER_N_EVENTS} \
       --compactFile ${compact_path} \
       --inputFiles ${JUGGLER_GEN_FILE} \
diff --git a/benchmarks/clustering/options/full_cal_reco.py b/benchmarks/clustering/options/full_cal_reco.py
index 24ef7bbe35411f49b9fe23463d00e2cfbf1f7da6..1f7834536b34fec435f4ad6825c56a7b8af2cb43 100644
--- a/benchmarks/clustering/options/full_cal_reco.py
+++ b/benchmarks/clustering/options/full_cal_reco.py
@@ -53,7 +53,7 @@ from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
 
 # branches needed from simulation root file
 sim_coll = [
-    "mcparticles",
+    "MCParticles",
     "EcalEndcapNHits",
     "EcalEndcapPHits",
     "EcalBarrelHits",
@@ -87,6 +87,7 @@ ce_ecal_reco = CalHitReco("ce_ecal_reco",
         thresholdFactor=4,          # 4 sigma cut on pedestal sigma
         readoutClass="EcalEndcapNHits",
         sectorField="sector",
+        samplingFraction=0.998,      # this accounts for a small fraction of leakage
         **ce_ecal_daq)
 
 ce_ecal_cl = IslandCluster("ce_ecal_cl",
@@ -100,11 +101,8 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl",
         dimScaledLocalDistXY=[1.8, 1.8])          # dimension scaled dist is good for hybrid sectors with different module size
 
 ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
-        inputHitCollection=ce_ecal_cl.inputHitCollection,
         inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
         outputClusterCollection="EcalEndcapNClusters",
-        mcHits="EcalEndcapNHits",
-        samplingFraction=0.998,      # this accounts for a small fraction of leakage
         logWeightBase=4.6)
 
 # Endcap Sampling Ecal
@@ -123,6 +121,7 @@ ci_ecal_reco = CalHitReco("ci_ecal_reco",
         inputHitCollection=ci_ecal_digi.outputHitCollection,
         outputHitCollection="EcalEndcapPHitsReco",
         thresholdFactor=5.0,
+        samplingFraction=ci_ecal_sf,
         **ci_ecal_daq)
 
 # merge hits in different layer (projection to local x-y plane)
@@ -145,12 +144,9 @@ ci_ecal_cl = IslandCluster("ci_ecal_cl",
         localDistXY=[10*mm, 10*mm])
 
 ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
-        inputHitCollection=ci_ecal_cl.inputHitCollection,
         inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection,
         outputClusterCollection="EcalEndcapPClusters",
-        mcHits="EcalEndcapPHits",
-        logWeightBase=6.2,
-        samplingFraction=ci_ecal_sf)
+        logWeightBase=6.2)
 
 # Central Barrel Ecal (Imaging Cal.)
 cb_ecal_daq = dict(
@@ -172,6 +168,7 @@ cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
         readoutClass="EcalBarrelHits",  # readout class
         layerField="layer",             # field to get layer id
         sectorField="module",           # field to get sector id
+        samplingFraction=cb_ecal_sf,
         **cb_ecal_daq)
 
 cb_ecal_cl = ImagingCluster("cb_ecal_cl",
@@ -183,12 +180,10 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl",
         sectorDist=3.*cm)                       # different sector
 
 cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
-        samplingFraction=cb_ecal_sf,
-        inputHitCollection=cb_ecal_cl.inputHitCollection,
-        inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection,
+        inputProtoClusters=cb_ecal_cl.outputProtoClusterCollection,
         mcHits="EcalBarrelHits",
-        outputClusterCollection="EcalBarrelImagingClusters",
-        outputLayerCollection="EcalBarrelImagingLayers")
+        outputClusters="EcalBarrelImagingClusters",
+        outputLayers="EcalBarrelImagingLayers")
 
 #Central ECAL SciFi
 # use the same daq_setting for digi/reco pair
@@ -211,6 +206,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
         layerField="layer",
         sectorField="module",
         localDetFields=["system", "module"], # use local coordinates in each module (stave)
+        samplingFraction=scifi_barrel_sf,
         **scfi_barrel_daq)
 
 # merge hits in different layer (projection to local x-y plane)
@@ -231,12 +227,9 @@ scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
         localDistXZ=[30*mm, 30*mm])
 
 scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
-       inputHitCollection=scfi_barrel_cl.inputHitCollection,
        inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection,
        outputClusterCollection="EcalBarrelScFiClusters",
-       mcHits="EcalBarrelScFiHits",
-       logWeightBase=6.2,
-       samplingFraction= scifi_barrel_sf)
+       logWeightBase=6.2)
 
 
 # Central Barrel Hcal
@@ -258,6 +251,7 @@ cb_hcal_reco = CalHitReco("cb_hcal_reco",
         readoutClass="HcalBarrelHits",
         layerField="layer",
         sectorField="module",
+        samplingFraction=cb_hcal_sf,
         **cb_hcal_daq)
 
 cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
@@ -275,12 +269,9 @@ cb_hcal_cl = IslandCluster("cb_hcal_cl",
         localDistXY=[15.*cm, 15.*cm])
 
 cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
-        inputHitCollection=cb_hcal_cl.inputHitCollection,
         inputProtoClusterCollection=cb_hcal_cl.outputProtoClusterCollection,
         outputClusterCollection="HcalBarrelClusters",
-        mcHits="HcalBarrelHits",
-        logWeightBase=6.2,
-        samplingFraction=cb_hcal_sf)
+        logWeightBase=6.2)
 
 
 # Hcal Hadron Endcap
@@ -298,6 +289,7 @@ ci_hcal_reco = CalHitReco("ci_hcal_reco",
         inputHitCollection=ci_hcal_digi.outputHitCollection,
         outputHitCollection="HcalEndcapPHitsReco",
         thresholdFactor=5.0,
+        samplingFraction=ci_hcal_sf,
         **ci_hcal_daq)
 
 ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
@@ -315,12 +307,9 @@ ci_hcal_cl = IslandCluster("ci_hcal_cl",
         localDistXY=[15.*cm, 15.*cm])
 
 ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
-        inputHitCollection=ci_hcal_cl.inputHitCollection,
         inputProtoClusterCollection=ci_hcal_cl.outputProtoClusterCollection,
         outputClusterCollection="HcalEndcapPClusters",
-        mcHits="HcalEndcapPHits",
-        logWeightBase=6.2,
-        samplingFraction=ci_hcal_sf)
+        logWeightBase=6.2)
 
 # Hcal Electron Endcap
 ce_hcal_daq = dict(
@@ -338,6 +327,7 @@ ce_hcal_reco = CalHitReco("ce_hcal_reco",
         inputHitCollection=ce_hcal_digi.outputHitCollection,
         outputHitCollection="HcalEndcapNHitsReco",
         thresholdFactor=5.0,
+        samplingFraction=ce_hcal_sf,
         **ce_hcal_daq)
 
 ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
@@ -355,16 +345,13 @@ ce_hcal_cl = IslandCluster("ce_hcal_cl",
         localDistXY=[15.*cm, 15.*cm])
 
 ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
-        inputHitCollection=ce_hcal_cl.inputHitCollection,
         inputProtoClusterCollection=ce_hcal_cl.outputProtoClusterCollection,
         outputClusterCollection="HcalEndcapNClusters",
-        mcHits="HcalEndcapNHits",
-        logWeightBase=6.2,
-        samplingFraction=ce_hcal_sf)
+        logWeightBase=6.2)
 
 
 podout.outputCommands = ['drop *',
-        'keep mcparticles',
+        'keep MCParticles',
         'keep *Digi',
         'keep *Reco*',
         'keep *Cluster*',
diff --git a/benchmarks/clustering/scripts/cluster_plots.py b/benchmarks/clustering/scripts/cluster_plots.py
index 567c51515e21b2b0d60c07b13f27d251aeec9af9..c7b2fd46ef21773b584a560e15d5a3580cfb47b0 100644
--- a/benchmarks/clustering/scripts/cluster_plots.py
+++ b/benchmarks/clustering/scripts/cluster_plots.py
@@ -48,12 +48,12 @@ def flatten_collection(rdf, collection, cols=None):
     return dfp
 
 
-def thrown_particles_figure(rdf, save, mcbranch="mcparticles"):
+def thrown_particles_figure(rdf, save, mcbranch="MCParticles"):
     # define truth particle info
-    dft = flatten_collection(rdf, mcbranch, ['genStatus', 'pdgID', 'ps.x', 'ps.y', 'ps.z', 'mass'])
+    dft = flatten_collection(rdf, mcbranch, ['generatorStatus', 'PDG', 'momentum.x', 'momentum.y', 'momentum.z', 'mass'])
     dft.rename(columns={c: c.replace(mcbranch + '.', '') for c in dft.columns}, inplace=True)
     # select thrown particles
-    dft = dft[dft['genStatus'] == 1]
+    dft = dft[dft['generatorStatus'] == 1]
 
     # figure
     fig, axs = plt.subplots(2, 2, figsize=(16, 12), dpi=120)
@@ -62,7 +62,7 @@ def thrown_particles_figure(rdf, save, mcbranch="mcparticles"):
     get_pname = np.vectorize(lambda pid: pdgbase.GetParticle(int(pid)).GetName())
 
     # enumerate particle names
-    dft.loc[:, 'pname'] = get_pname(dft['pdgID'].values)
+    dft.loc[:, 'pname'] = get_pname(dft['PDG'].values)
     penum = {pname: i for i, pname in enumerate(dft['pname'].unique())}
     dft.loc[:, 'pname_id'] = dft['pname'].map(penum)
 
@@ -77,7 +77,7 @@ def thrown_particles_figure(rdf, save, mcbranch="mcparticles"):
 
     # calculate kinematics
     get_4vecs = np.vectorize(lambda x, y, z, m: ROOT.Math.PxPyPzMVector(x, y, z, m))
-    fourvecs = get_4vecs(*dft[['ps.x', 'ps.y', 'ps.z', 'mass']].values.T)
+    fourvecs = get_4vecs(*dft[['momentum.x', 'momentum.y', 'momentum.z', 'mass']].values.T)
 
     dft.loc[:, 'p'] = [v.P() for v in fourvecs]
     dft.loc[:, 'eta'] = [v.Eta() for v in fourvecs]
@@ -109,7 +109,7 @@ if __name__ == '__main__':
     parser.add_argument('-o', dest='outdir', default='.', help='Output directory.')
     parser.add_argument('-m', '--macros', dest='macros', default='rootlogon.C',
             help='Macro files to be loaded by root, separated by \",\".')
-    parser.add_argument('--mc-branch', dest='mc', default='mcparticles',
+    parser.add_argument('--mc-branch', dest='mc', default='MCParticles',
             help='Branch name of MC particles truth info.')
     parser.add_argument('--collections', dest='coll', default='',
             help='Collection name of clusters to plot')
diff --git a/benchmarks/clustering/scripts/deprecated/barrel_clusters.cxx b/benchmarks/clustering/scripts/deprecated/barrel_clusters.cxx
index 5720c7ea3fa19f07154a41e325bd508c368e645e..efc8291491b87c0748f940de386a82815a8fceeb 100644
--- a/benchmarks/clustering/scripts/deprecated/barrel_clusters.cxx
+++ b/benchmarks/clustering/scripts/deprecated/barrel_clusters.cxx
@@ -1,10 +1,8 @@
 #include <iostream>
 #include "ROOT/RDataFrame.hxx"
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
-// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
-//  export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
 
 using ROOT::RDataFrame;
 using namespace ROOT::VecOps;
@@ -23,34 +21,34 @@ auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
   }
   return result;
 };
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto fourvec = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv);
   }
   return result;
 };
-auto pt = [](std::vector<dd4pod::Geant4ParticleData> const& in) {
+auto pt = [](std::vector<edm4hep::MCParticleData> 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));
+    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
   }
   return result;
 };
-auto pid = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto pid = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> const& in) {
   std::vector<int> result;
   for (size_t i = 0; i < in.size(); ++i) {
-    result.push_back(in[i].pdgID);
+    result.push_back(in[i].PDG);
   }
   return result;
 };
-auto eta = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto eta = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> const& in) {
   std::vector<float> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv.Eta());
   }
   return result;
@@ -91,8 +89,8 @@ int barrel_clusters(const char* in_fname = "topside/rec_barrel_clusters.root")
   ROOT::EnableImplicitMT();
   ROOT::RDataFrame df("events", in_fname);
 
-  auto d0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                .Define("thrownParticles", "mcparticles[isThrown]")
+  auto d0 = df.Define("isThrown", "MCParticles.generatorStatus == 1")
+                .Define("thrownParticles", "MCParticles[isThrown]")
                 .Define("thrownP", fourvec, {"thrownParticles"})
                 .Define("thrownEta", eta, {"thrownParticles"})
                 .Define("thrownTheta", theta, {"thrownP"})
@@ -159,6 +157,6 @@ int barrel_clusters(const char* in_fname = "topside/rec_barrel_clusters.root")
 
   //c->SaveAs("results/crystal_cal_electrons_Ecluster.png");
   //std::string outfilename = "rdf_test.root";
-  //df2.Snapshot("events", outfilename, {"MCParticles_pt", "mcparticles"});
+  //df2.Snapshot("events", outfilename, {"MCParticles_pt", "MCParticles"});
   return 0;
 }
diff --git a/benchmarks/ecal/options/barrel.py b/benchmarks/ecal/options/barrel.py
index 1919fbf045a3012995163a88ff54f4405f27729b..88fcea482956bf5dbd20457fbd90807a8e004dbd 100644
--- a/benchmarks/ecal/options/barrel.py
+++ b/benchmarks/ecal/options/barrel.py
@@ -38,7 +38,7 @@ from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
 
 # branches needed from simulation root file
 sim_coll = [
-    "mcparticles",
+    "MCParticles",
     "EcalBarrelHits",
 ]
 
@@ -63,6 +63,7 @@ cb_ecal_digi = CalHitDigi("cb_ecal_digi",
 cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
         inputHitCollection=cb_ecal_digi.outputHitCollection,
         outputHitCollection="EcalBarrelImagingHitsReco",
+        samplingFraction=cb_ecal_sf,
         thresholdFactor=3,  # about 20 keV
         readoutClass="EcalBarrelHits",  # readout class
         layerField="layer",             # field to get layer id
@@ -78,15 +79,13 @@ cb_ecal_cl = ImagingCluster("cb_ecal_cl",
         sectorDist=3.*cm)                       # different sector
 
 cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
-        samplingFraction=cb_ecal_sf,
-        inputHitCollection=cb_ecal_cl.inputHitCollection,
-        inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection,
-        outputClusterCollection="EcalBarrelImagingClusters",
-        outputLayerCollection="EcalBarrelImagingLayers",
+        inputProtoClusters=cb_ecal_cl.outputProtoClusterCollection,
+        outputClusters="EcalBarrelImagingClusters",
+        outputLayers="EcalBarrelImagingLayers",
         mcHits="EcalBarrelHits")
 
 podout.outputCommands = ['drop *',
-        'keep mcparticles',
+        'keep MCParticles',
         'keep *HitsReco',
         'keep *HitsDigi',
         'keep *Cluster*']
diff --git a/benchmarks/ecal/options/endcap_e.py b/benchmarks/ecal/options/endcap_e.py
index 9eb12ffe8f34848c95a17a8e6bce8354ca014ab8..1b89fd08e630a71d90b8751b4e0b599ce3d60499 100644
--- a/benchmarks/ecal/options/endcap_e.py
+++ b/benchmarks/ecal/options/endcap_e.py
@@ -31,7 +31,7 @@ from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
 
 # branches needed from simulation root file
 sim_coll = [
-    "mcparticles",
+    "MCParticles",
     "EcalEndcapNHits",
 ]
 
@@ -56,6 +56,7 @@ ce_ecal_digi = CalHitDigi("ce_ecal_digi",
 ce_ecal_reco = CalHitReco("ce_ecal_reco",
         inputHitCollection=ce_ecal_digi.outputHitCollection,
         outputHitCollection="EcalEndcapNHitsReco",
+        samplingFraction=0.998,      # this accounts for a small fraction of leakage
         thresholdFactor=4,          # 4 sigma cut on pedestal sigma
         readoutClass="EcalEndcapNHits",
         sectorField="sector",
@@ -72,15 +73,12 @@ ce_ecal_cl = IslandCluster("ce_ecal_cl",
         dimScaledLocalDistXY=[1.8, 1.8])          # hybrid calorimeter with different dimensions, using a scaled dist
 
 ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
-        inputHitCollection=ce_ecal_cl.inputHitCollection,
         inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
         outputClusterCollection="EcalEndcapNClusters",
-        mcHits="EcalEndcapNHits",
-        samplingFraction=0.998,      # this accounts for a small fraction of leakage
         logWeightBase=4.6)
 
 podout.outputCommands = ['drop *',
-        'keep mcparticles',
+        'keep MCParticles',
         'keep *HitsReco',
         'keep *HitsDigi',
         'keep *Cluster*']
diff --git a/benchmarks/ecal/options/endcap_i.py b/benchmarks/ecal/options/endcap_i.py
index aaa00dc7145735456441e98dc29e0a0b762b67c1..e8a091b789b8daadd44ca1ec90844ead4e63e863 100644
--- a/benchmarks/ecal/options/endcap_i.py
+++ b/benchmarks/ecal/options/endcap_i.py
@@ -34,7 +34,7 @@ from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
 
 # branches needed from simulation root file
 sim_coll = [
-    "mcparticles",
+    "MCParticles",
     "EcalEndcapPHits",
 ]
 
@@ -58,6 +58,7 @@ ci_ecal_digi = CalHitDigi("ci_ecal_digi",
 ci_ecal_reco = CalHitReco("ci_ecal_reco",
         inputHitCollection=ci_ecal_digi.outputHitCollection,
         outputHitCollection="EcalEndcapPHitsReco",
+        samplingFraction=ci_ecal_sf,
         thresholdFactor=5.0,
         **ci_ecal_daq)
 
@@ -81,15 +82,12 @@ ci_ecal_cl = IslandCluster("ci_ecal_cl",
         localDistXY=[10*mm, 10*mm])
 
 ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
-        inputHitCollection=ci_ecal_cl.inputHitCollection,
         inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection,
         outputClusterCollection="EcalEndcapPClusters",
-        mcHits="EcalEndcapPHits",
-        logWeightBase=6.2,
-        samplingFraction=ci_ecal_sf)
+        logWeightBase=6.2)
 
 podout.outputCommands = ['drop *',
-        'keep mcparticles',
+        'keep MCParticles',
         'keep *HitsReco*',
         'keep *HitsDigi',
         'keep *Cluster*']
diff --git a/benchmarks/ecal/run_emcal_benchmarks.py b/benchmarks/ecal/run_emcal_benchmarks.py
index f5ef2801d212aac680257b731fa1d8c653b261f1..e44c455d9627f23bef0e52e3b87760e1023295b5 100755
--- a/benchmarks/ecal/run_emcal_benchmarks.py
+++ b/benchmarks/ecal/run_emcal_benchmarks.py
@@ -58,7 +58,7 @@ if args.run_type not in default_type:
 opt_scripts, ana_scripts, collections = default_type[args.run_type]
 
 gen_file = 'gen_{nametag}_{pmin}_{pmax}.hepmc'.format(**kwargs)
-sim_file = 'sim_{nametag}_{pmin}_{pmax}.root'.format(**kwargs)
+sim_file = 'sim_{nametag}_{pmin}_{pmax}.edm4hep.root'.format(**kwargs)
 rec_file = 'rec_{nametag}_{pmin}_{pmax}.root'.format(**kwargs)
 
 procs = [p.strip() for p in args.process.split(',')]
@@ -74,9 +74,10 @@ if 'sim' in procs:
             '--particles', args.particles]
     subprocess.run(gen_cmd)
     # simulation
-    sim_cmd = ['npsim',
+    sim_cmd = ['ddsim',
             '--part.minimalKineticEnergy', '1*TeV',
             '--numberOfEvents', '{}'.format(args.nev),
+            '--filter.tracker', 'edep0',
             '--runType', 'batch',
             '--inputFiles', gen_file,
             '--outputFile', sim_file,
diff --git a/benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx b/benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx
index 125135b03197f4496465bccbb13919004e913f60..5ac4ecaae55679140d568447c9151990df001725 100644
--- a/benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx
+++ b/benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx
@@ -8,8 +8,8 @@
 #include <iostream>
 #include <vector>
 
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "dd4pod/CalorimeterHitCollection.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
 
 #include "TCanvas.h"
 #include "TStyle.h"
@@ -30,7 +30,7 @@ R__LOAD_LIBRARY(libDD4pod.so)
 
 using ROOT::RDataFrame;
 
-void analysis_zdc_neutrons(const char* input_fname = "sim_zdc_uniform_neutron.root")
+void analysis_zdc_neutrons(const char* input_fname = "sim_zdc_uniform_neutron.edm4hep.root")
 {
   // Setting for graphs
   TStyle* kStyle = new TStyle("kStyle","Kim's Style");
@@ -67,38 +67,39 @@ void analysis_zdc_neutrons(const char* input_fname = "sim_zdc_uniform_neutron.ro
   ROOT::RDataFrame d0("events", input_fname);
 
   // Thrown Energy [GeV]
-  auto Ethr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+  auto Ethr = [](const std::vector<edm4hep::MCParticleData> & input) {
     auto p = input[2];
-    auto energy = TMath::Sqrt(p.ps.x * p.ps.x + p.ps.y * p.ps.y + p.ps.z * p.ps.z + p.mass * p.mass);
+    auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass);
     return energy;
   };
 
   // Theta [mrad]
-  auto Thetathr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+  auto Thetathr = [](const std::vector<edm4hep::MCParticleData> & input) {
     auto p = input[2];
-    auto theta = p.ps.theta();
+    auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
     return theta*1000.0;
   };
 
   // Phi [rad]
-  auto Phithr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+  auto Phithr = [](const std::vector<edm4hep::MCParticleData> & input) {
     auto p = input[2];
-    auto phi = p.ps.phi();
+    auto phi = std::atan2(p.momentum.y, p.momentum.x);
     return phi;
   };
   
   // Eta
-  auto Etathr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+  auto Etathr = [](const std::vector<edm4hep::MCParticleData> & input) {
     auto p = input[2];
-    auto eta = p.ps.eta();
+    auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
+    auto eta = -std::log(std::tan(theta/2));
     return eta;
   };  
 
   // Define variables
-  auto d1 = d0.Define("Ethr",      Ethr,     {"mcparticles"})
-              .Define("Thetathr",  Thetathr, {"mcparticles"})
-              .Define("Phithr",    Phithr,   {"mcparticles"})
-	      .Define("Etathr",    Etathr,   {"mcparticles"})
+  auto d1 = d0.Define("Ethr",      Ethr,     {"MCParticles"})
+              .Define("Thetathr",  Thetathr, {"MCParticles"})
+              .Define("Phithr",    Phithr,   {"MCParticles"})
+	      .Define("Etathr",    Etathr,   {"MCParticles"})
 	      ;
 
   // Define Histograms
diff --git a/benchmarks/far_forward/analysis/analysis_zdc_photons.cxx b/benchmarks/far_forward/analysis/analysis_zdc_photons.cxx
index e0cf310279ba9f98e9eb0413d18bef1ff15d7fb3..4f2ae18e765c98ceec0b34ff2732c9ef43c3bfe1 100644
--- a/benchmarks/far_forward/analysis/analysis_zdc_photons.cxx
+++ b/benchmarks/far_forward/analysis/analysis_zdc_photons.cxx
@@ -8,8 +8,8 @@
 #include <iostream>
 #include <vector>
 
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "dd4pod/CalorimeterHitCollection.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/CalorimeterHitCollection.h"
 
 #include "TCanvas.h"
 #include "TStyle.h"
@@ -30,7 +30,7 @@ R__LOAD_LIBRARY(libDD4pod.so)
 
 using ROOT::RDataFrame;
 
-void analysis_zdc_photons(const char* input_fname = "sim_zdc_uniform_photon.root")
+void analysis_zdc_photons(const char* input_fname = "sim_zdc_uniform_photon.edm4hep.root")
 {
   // Setting for graphs
   TStyle* kStyle = new TStyle("kStyle","Kim's Style");
@@ -67,38 +67,39 @@ void analysis_zdc_photons(const char* input_fname = "sim_zdc_uniform_photon.root
   ROOT::RDataFrame d0("events", input_fname);
 
   // Thrown Energy [GeV]
-  auto Ethr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+  auto Ethr = [](const std::vector<edm4hep::MCParticleData> & input) {
     auto p = input[2];
-    auto energy = TMath::Sqrt(p.ps.x * p.ps.x + p.ps.y * p.ps.y + p.ps.z * p.ps.z + p.mass * p.mass);
+    auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass);
     return energy;
   };
 
   // Theta [mrad]
-  auto Thetathr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+  auto Thetathr = [](const std::vector<edm4hep::MCParticleData> & input) {
     auto p = input[2];
-    auto theta = p.ps.theta();
+    auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
     return theta*1000.0;
   };
 
   // Phi [rad]
-  auto Phithr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+  auto Phithr = [](const std::vector<edm4hep::MCParticleData> & input) {
     auto p = input[2];
-    auto phi = p.ps.phi();
+    auto phi = std::atan2(p.momentum.y, p.momentum.x);
     return phi;
   };
   
   // Eta
-  auto Etathr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+  auto Etathr = [](const std::vector<edm4hep::MCParticleData> & input) {
     auto p = input[2];
-    auto eta = p.ps.eta();
+    auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
+    auto eta = -std::log(std::tan(theta/2));
     return eta;
   };  
 
   // Define variables
-  auto d1 = d0.Define("Ethr",      Ethr,     {"mcparticles"})
-              .Define("Thetathr",  Thetathr, {"mcparticles"})
-              .Define("Phithr",    Phithr,   {"mcparticles"})
-	      .Define("Etathr",    Etathr,   {"mcparticles"})
+  auto d1 = d0.Define("Ethr",      Ethr,     {"MCParticles"})
+              .Define("Thetathr",  Thetathr, {"MCParticles"})
+              .Define("Phithr",    Phithr,   {"MCParticles"})
+	      .Define("Etathr",    Etathr,   {"MCParticles"})
 	      ;
 
   // Define Histograms
diff --git a/benchmarks/far_forward/analysis/hits_far_forward_protons.cxx b/benchmarks/far_forward/analysis/hits_far_forward_protons.cxx
index 860a1b6c5a4774cce25de8a6668bfef29ef5ed14..976ad546ed7c15843926edfa64f23258f16618ec 100644
--- a/benchmarks/far_forward/analysis/hits_far_forward_protons.cxx
+++ b/benchmarks/far_forward/analysis/hits_far_forward_protons.cxx
@@ -10,9 +10,9 @@
 
 R__LOAD_LIBRARY(libeicd.so)
 R__LOAD_LIBRARY(libDD4pod.so)
-#include "dd4pod/TrackerHitCollection.h"
-#include "dd4pod/TrackerHitData.h"
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/SimTrackerHitData.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/TrackParametersCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
@@ -30,10 +30,10 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 };
 
 
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+std::vector<float> pt (std::vector<edm4hep::MCParticleData> 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));
+    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
   }
   return result;
 }
@@ -52,11 +52,11 @@ auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
   }
   return result;
 };
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto fourvec = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv);
   }
   return result;
@@ -82,14 +82,14 @@ auto delta_p_over_p = [](const std::vector<double>& tracks, const std::vector<do
   return res;
 };
 
-int hits_far_forward_protons(const char* fname = "sim_far_forward_protons.root")
+int hits_far_forward_protons(const char* fname = "sim_far_forward_protons.edm4hep.root")
 {
 
   ROOT::EnableImplicitMT();
   ROOT::RDataFrame df("events", fname);
 
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
+  auto df0 = df.Define("isThrown", "MCParticles.generatorStatus == 1")
+                 .Define("thrownParticles", "MCParticles[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
@@ -105,8 +105,8 @@ int hits_far_forward_protons(const char* fname = "sim_far_forward_protons.root")
                  //.Define("delta_p_over_p1",delta_p_over_p, {"p_track1", "p_thrown"})
                  //.Define("delta_p_over_p2",delta_p_over_p, {"p_track2", "p_thrown"})
                  //.Define("N_VtxBarrelHits",[](std::vector<eic::TrackerHitData> hits) { return hits.size();},{"VertexBarrelRecHits"})
-                 .Define("N_BarrelHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelHits"})
-                 .Define("N_EndcapHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapHits"})
+                 .Define("N_BarrelHits", [](std::vector<edm4hep::SimTrackerHitData> hits) { return hits.size();}, {"TrackerBarrelHits"})
+                 .Define("N_EndcapHits", [](std::vector<edm4hep::SimTrackerHitData> hits) { return hits.size();}, {"TrackerEndcapHits1"})
                  ;
 
   auto hBarrel_x_vs_y = df0.Histo2D({"hBarrel_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "TrackerBarrelHits.position.x", "TrackerBarrelHits.position.y");
diff --git a/benchmarks/far_forward/analysis/print_far_forward_protons.cxx b/benchmarks/far_forward/analysis/print_far_forward_protons.cxx
index 0fbbf60f500cc2cbdd1bb35e45ff4b0297ea9dc0..841c0bf89b3659d50d5341e8cf1c6f3e86ee05a8 100644
--- a/benchmarks/far_forward/analysis/print_far_forward_protons.cxx
+++ b/benchmarks/far_forward/analysis/print_far_forward_protons.cxx
@@ -8,7 +8,7 @@
 
 R__LOAD_LIBRARY(libeicd.so)
 R__LOAD_LIBRARY(libDD4pod.so)
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/TrackParametersCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
diff --git a/benchmarks/far_forward/analysis/rec_far_forward_protons.cxx b/benchmarks/far_forward/analysis/rec_far_forward_protons.cxx
index 9e0ae6b8c8dfdcd7e4cfde904fd5ea21c9219fd6..8a84c11cb1f6be7c0f97484c2d8add61652cca91 100644
--- a/benchmarks/far_forward/analysis/rec_far_forward_protons.cxx
+++ b/benchmarks/far_forward/analysis/rec_far_forward_protons.cxx
@@ -10,7 +10,7 @@
 
 R__LOAD_LIBRARY(libeicd.so)
 R__LOAD_LIBRARY(libDD4pod.so)
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/TrackParametersCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
@@ -28,10 +28,10 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 };
 
 
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+std::vector<float> pt (std::vector<edm4hep::MCParticleData> 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));
+    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
   }
   return result;
 }
@@ -50,11 +50,11 @@ auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
   }
   return result;
 };
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto fourvec = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv);
   }
   return result;
@@ -86,8 +86,8 @@ int rec_far_forward_protons(const char* fname = "topside/rec_far_forward_protons
   ROOT::EnableImplicitMT();
   ROOT::RDataFrame df("events", fname);
 
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
+  auto df0 = df.Define("isThrown", "MCParticles.generatorStatus == 1")
+                 .Define("thrownParticles", "MCParticles[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
diff --git a/benchmarks/far_forward/config.yml b/benchmarks/far_forward/config.yml
index 81dda4710e6edf8ad4d4d8d02e3e1f920db6c4b2..b965eaf677bbe126a2ee17d8eff6a5c320141a02 100644
--- a/benchmarks/far_forward/config.yml
+++ b/benchmarks/far_forward/config.yml
@@ -1,23 +1,26 @@
+far_forward:compile:
+  stage: compile
+  extends: .compile_benchmark
+  script:
+    - compile_analyses.py far_forward
+
 B0_far_forward_protons:
   extends: .rec_benchmark
   stage: run
-  timeout: 24 hours
+  needs: ["far_forward:compile"]
   script:
-    - compile_analyses.py far_forward
     - bash benchmarks/far_forward/far_forward_protons.sh
 
 ZDC_far_forward_neutrons:
   extends: .rec_benchmark
   stage: run
-  timeout: 24 hours
+  needs: ["far_forward:compile"]
   script:
-    - compile_analyses.py far_forward
     - bash benchmarks/far_forward/run_zdc_neutrons.sh
 
 ZDC_far_forward_photons:
   extends: .rec_benchmark
   stage: run
-  timeout: 24 hours
+  needs: ["far_forward:compile"]
   script:
-    - compile_analyses.py far_forward
     - bash benchmarks/far_forward/run_zdc_photons.sh
diff --git a/benchmarks/far_forward/far_forward_protons.sh b/benchmarks/far_forward/far_forward_protons.sh
index f91eec96e73deccb298293fd9dd956bd88e3218b..8b6bf299f720646736c2fc1514a1251964b2922b 100644
--- a/benchmarks/far_forward/far_forward_protons.sh
+++ b/benchmarks/far_forward/far_forward_protons.sh
@@ -59,7 +59,7 @@ fi
 export JUGGLER_FILE_NAME_TAG="far_forward_protons"
 export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
 
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
 export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
 
 echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
@@ -78,8 +78,9 @@ then
 
   echo "Running geant4 simulation"
   ## run geant4 simulations
-  npsim --runType batch \
+  ddsim --runType batch \
     --part.minimalKineticEnergy 1000*GeV  \
+    --filter.tracker edep0 \
     -v WARNING \
     --numberOfEvents ${JUGGLER_N_EVENTS} \
     --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
diff --git a/benchmarks/far_forward/options/far_forward_reconstruction.py b/benchmarks/far_forward/options/far_forward_reconstruction.py
index ad009bca6650f9958deb5f0dce7afd1be0f120f6..7aef3b7dccc184e5d8d9e429fa607b75e4bf56c8 100644
--- a/benchmarks/far_forward/options/far_forward_reconstruction.py
+++ b/benchmarks/far_forward/options/far_forward_reconstruction.py
@@ -16,7 +16,9 @@ geo_service = GeoSvc("GeoSvc", detectors=[compact_path], OutputLevel=WARNING)
 podioevent = EICDataSvc("EventDataSvc", inputs=input_sims, OutputLevel=WARNING)
 
 from Configurables import PodioInput
-from Configurables import Jug__Digi__UFSDTrackerDigi as TrackerDigi
+
+from Configurables import Jug__Digi__SimTrackerHitsCollector as SimTrackerHitsCollector
+from Configurables import Jug__Digi__SiliconTrackerDigi as TrackerDigi
 from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerReco
 
 from Configurables import Jug__Reco__TrackerSourceLinker as TrackerSourceLinker
@@ -30,12 +32,24 @@ from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrack
 from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
 
 sim_colls = [
-    "mcparticles",
-    "B0TrackerHits",
-    "ForwardRomanPotHits",
-    "ForwardOffMTrackerHits"
+    "MCParticles",
+    "B0TrackerHits"
+]
+
+ffi_romanpot_collections = [
+    "ForwardRomanPotHits1",
+    "ForwardRomanPotHits2"
 ]
 
+ffi_offmtracker_collections = [
+    "ForwardOffMTrackerHits1",
+    "ForwardOffMTrackerHits2",
+    "ForwardOffMTrackerHits3",
+    "ForwardOffMTrackerHits4"
+]
+
+sim_colls += ffi_romanpot_collections + ffi_offmtracker_collections
+
 # list of algorithms
 algorithms = []
 
@@ -44,8 +58,12 @@ podin = PodioInput("PodioReader", collections=sim_colls)
 algorithms.append(podin)
 
 ## Roman pots
+ffi_romanpot_coll = SimTrackerHitsCollector("ffi_romanpot_coll",
+        inputSimTrackerHits = ffi_romanpot_collections,
+        outputSimTrackerHits = "ForwardRomanPotAllHits")
+algorithms.append(ffi_romanpot_coll)
 ffi_romanpot_digi = TrackerDigi("ffi_romanpot_digi",
-        inputHitCollection = "ForwardRomanPotHits",
+        inputHitCollection = ffi_romanpot_coll.outputSimTrackerHits,
         outputHitCollection = "ForwardRomanPotRawHits",
         timeResolution = 8)
 algorithms.append(ffi_romanpot_digi)
@@ -61,8 +79,12 @@ ffi_romanpot_parts = FarForwardParticles("ffi_romanpot_parts",
 algorithms.append(ffi_romanpot_parts)
 
 ## Off momentum tracker
+ffi_offmtracker_coll = SimTrackerHitsCollector("ffi_offmtracker_coll",
+        inputSimTrackerHits = ffi_offmtracker_collections,
+        outputSimTrackerHits = "ForwardOffMTrackerAllHits")
+algorithms.append(ffi_offmtracker_coll)
 ffi_offmtracker_digi = TrackerDigi("ffi_offmtracker_digi",
-        inputHitCollection = "ForwardOffMTrackerHits",
+        inputHitCollection = ffi_offmtracker_coll.outputSimTrackerHits,
         outputHitCollection = "ForwardOffMTrackerRawHits",
         timeResolution = 8)
 algorithms.append(ffi_offmtracker_digi)
@@ -98,7 +120,7 @@ algorithms.append(sourcelinker)
 
 ## Track param init
 truth_trk_init = TrackParamTruthInit("truth_trk_init",
-        inputMCParticles="mcparticles",
+        inputMCParticles="MCParticles",
         outputInitialTrackParameters="InitTrackParams",
         OutputLevel=INFO)
 algorithms.append(truth_trk_init)
@@ -127,7 +149,7 @@ podout.outputCommands = [
     "drop trajectories",
     "drop outputSourceLinks",
     "drop outputInitialTrackParameters",
-    "keep mcparticles",
+    "keep MCParticles",
 ]
 algorithms.append(podout)
 
diff --git a/benchmarks/far_forward/options/zdc_reconstruction.py b/benchmarks/far_forward/options/zdc_reconstruction.py
index e95c2dd9d91189855cde0c82592f35413ec84f62..695d5fb3e1629c3f1a3cf5ba1c2e7335535a543a 100644
--- a/benchmarks/far_forward/options/zdc_reconstruction.py
+++ b/benchmarks/far_forward/options/zdc_reconstruction.py
@@ -14,7 +14,7 @@ from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
 kwargs = dict()
 
 # input arguments through environment variables
-kwargs['input'] = os.environ.get('JUGGLER_SIM_FILE', '../atheta/zdc_photons.root')
+kwargs['input'] = os.environ.get('JUGGLER_SIM_FILE', '../atheta/zdc_photons.edm4hep.root')
 kwargs['output'] = os.environ.get('JUGGLER_REC_FILE', 'rec_zdc_photons.root')
 kwargs['compact'] = os.environ.get('DETECTOR_COMPACT_PATH', '../atheta/athena.xml')
 kwargs['nev'] = int(os.environ.get('JUGGLER_N_EVENTS', 100))
@@ -30,7 +30,7 @@ geo_service = GeoSvc('GeoSvc', detectors=kwargs['compact'].split(','), OutputLev
 podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
 
 sim_colls = [
-    'mcparticles',
+    'MCParticles',
     'ZDCEcalHits',
     'ZDCHcalHits'
 ]
@@ -70,10 +70,8 @@ WSciFi_cl = IslandCluster('WSciFi_cl',
         splitCluster=True)
 
 WSciFi_clreco = RecoCoG('WSciFi_clreco',
-        inputHitCollection=WSciFi_cl.inputHitCollection,
         inputProtoClusterCollection = WSciFi_cl.outputProtoClusterCollection,
         outputClusterCollection='ZDCEcalClusters',
-        mcHits="ZDCEcalHits",
         logWeightBase=3.6)
 
 PbSci_cl = IslandCluster('PbSci_cl',
@@ -85,10 +83,8 @@ PbSci_cl = IslandCluster('PbSci_cl',
         splitCluster=False)
 
 PbSci_clreco = RecoCoG('PbSci_clreco',
-        inputHitCollection=PbSci_cl.inputHitCollection,
         inputProtoClusterCollection = PbSci_cl.outputProtoClusterCollection,
         outputClusterCollection='ZDCHcalClusters',
-        mcHits="ZDCHcalHits",
         logWeightBase=3.6,
         samplingFraction=kwargs['PbSci_sf'])
 
diff --git a/benchmarks/far_forward/run_zdc_neutrons.sh b/benchmarks/far_forward/run_zdc_neutrons.sh
index c4072140c09e7994c91d65cea1d9f6439ad2ff0f..9c54b656ae42d5bf6570071781338721f54086e8 100644
--- a/benchmarks/far_forward/run_zdc_neutrons.sh
+++ b/benchmarks/far_forward/run_zdc_neutrons.sh
@@ -72,7 +72,7 @@ export DETECTOR_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml
 export JUGGLER_FILE_NAME_TAG="zdc_neutrons"
 export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
 
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
 export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
 
 echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
@@ -88,8 +88,9 @@ then
   fi
 
   echo "Running Geant4 simulation"
-  npsim --runType batch \
+  ddsim --runType batch \
     --part.minimalKineticEnergy 0.5*MeV  \
+    --filter.tracker edep0 \
     -v WARNING \
     --numberOfEvents ${JUGGLER_N_EVENTS} \
     --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
diff --git a/benchmarks/far_forward/run_zdc_photons.sh b/benchmarks/far_forward/run_zdc_photons.sh
index 0baf77ec0c752c76f6e8d63ce467b65b8fb4c34a..ef22e2ab5d1324de7aee4450c8506a23753e1cc1 100644
--- a/benchmarks/far_forward/run_zdc_photons.sh
+++ b/benchmarks/far_forward/run_zdc_photons.sh
@@ -72,7 +72,7 @@ export DETECTOR_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml
 export JUGGLER_FILE_NAME_TAG="zdc_photons"
 export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
 
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
 export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
 
 echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
@@ -88,8 +88,9 @@ then
   fi
 
   echo "Running Geant4 simulation"
-  npsim --runType batch \
+  ddsim --runType batch \
     --part.minimalKineticEnergy 0.5*MeV  \
+    --filter.tracker edep0 \
     -v WARNING \
     --numberOfEvents ${JUGGLER_N_EVENTS} \
     --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
diff --git a/benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx b/benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx
index b07f5728b38685f19831b33fde5c01f905d3e2f1..288c59c5cacfbe5fdbf2f18a8cdb98a3e2b0f603 100644
--- a/benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx
+++ b/benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx
@@ -8,8 +8,8 @@
 #include <algorithm>
 #include <string>
 
-#include "dd4pod/CalorimeterHitCollection.h"
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/MCParticleCollection.h"
 
 #include "eicd/CalorimeterHitCollection.h"
 #include "eicd/CalorimeterHitData.h"
@@ -104,34 +104,34 @@ void emcal_barrel_pion_rejection_analysis(
   // MCParticles Functions/////////////////////////////////////////////////////////////////////////////////////
 
   // Filter function to get electrons
-  auto is_electron = [](std::vector<dd4pod::Geant4ParticleData> const& input){
-    return (input[2].pdgID == 11 ? true : false);
+  auto is_electron = [](std::vector<edm4hep::MCParticleData> const& input){
+    return (input[2].PDG == 11 ? true : false);
   };
 
   // Filter function to get just negative pions
-  auto is_piMinus = [](std::vector<dd4pod::Geant4ParticleData> const& input){
-    return (input[2].pdgID == -211 ? true : false);
+  auto is_piMinus = [](std::vector<edm4hep::MCParticleData> const& input){
+    return (input[2].PDG == -211 ? true : false);
   };
   
   // Thrown Energy [GeV]
-  auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-    return TMath::Sqrt(input[2].ps.x*input[2].ps.x + input[2].ps.y*input[2].ps.y + input[2].ps.z*input[2].ps.z + input[2].mass*input[2].mass);
+  auto Ethr = [](std::vector<edm4hep::MCParticleData> const& input) {
+    return TMath::Sqrt(input[2].momentum.x*input[2].momentum.x + input[2].momentum.y*input[2].momentum.y + input[2].momentum.z*input[2].momentum.z + input[2].mass*input[2].mass);
   };
 
   // Thrown Momentum [GeV]
-  auto Pthr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-    return TMath::Sqrt(input[2].ps.x*input[2].ps.x + input[2].ps.y*input[2].ps.y + input[2].ps.z*input[2].ps.z);
+  auto Pthr = [](std::vector<edm4hep::MCParticleData> const& input) {
+    return TMath::Sqrt(input[2].momentum.x*input[2].momentum.x + input[2].momentum.y*input[2].momentum.y + input[2].momentum.z*input[2].momentum.z);
   };
 
   // Thrown Eta 
-  auto Eta = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-    double E  = TMath::Sqrt(input[2].ps.x*input[2].ps.x + input[2].ps.y*input[2].ps.y + input[2].ps.z*input[2].ps.z + input[2].mass*input[2].mass);
-    return 0.5*TMath::Log((E + input[2].ps.z) / (E - input[2].ps.z));
+  auto Eta = [](std::vector<edm4hep::MCParticleData> const& input) {
+    double E  = TMath::Sqrt(input[2].momentum.x*input[2].momentum.x + input[2].momentum.y*input[2].momentum.y + input[2].momentum.z*input[2].momentum.z + input[2].mass*input[2].mass);
+    return 0.5*TMath::Log((E + input[2].momentum.z) / (E - input[2].momentum.z));
   };
 
   // Thrown pT [GeV]
-  auto pT = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-    return TMath::Sqrt(input[2].ps.x*input[2].ps.x + input[2].ps.y*input[2].ps.y);
+  auto pT = [](std::vector<edm4hep::MCParticleData> const& input) {
+    return TMath::Sqrt(input[2].momentum.x*input[2].momentum.x + input[2].momentum.y*input[2].momentum.y);
   };
 
   // RecoEcalBarrelImagingHits Functions/////////////////////////////////////////////////////////////////////////////////////
@@ -475,10 +475,10 @@ void emcal_barrel_pion_rejection_analysis(
 
   // Define variables
   auto d1 = d0
-              .Define("Ethr",                         Ethr,                           {"mcparticles"})
-              .Define("Pthr",                         Pthr,                           {"mcparticles"})
-              .Define("Eta",                          Eta,                            {"mcparticles"})
-              .Define("pT",                           pT,                             {"mcparticles"})
+              .Define("Ethr",                         Ethr,                           {"MCParticles"})
+              .Define("Pthr",                         Pthr,                           {"MCParticles"})
+              .Define("Eta",                          Eta,                            {"MCParticles"})
+              .Define("pT",                           pT,                             {"MCParticles"})
 
               .Define("ERecoImg",                     ERecoImg,                       {"RecoEcalBarrelImagingHits"})
               .Define("nhitsImgLayer",                nhitsImgLayer,                  {"RecoEcalBarrelImagingHits"})
@@ -552,8 +552,8 @@ void emcal_barrel_pion_rejection_analysis(
                                       "ERecoScFiLayer3OverP", "ERecoScFiLayer4OverP", "ERecoScFiLayer5OverP","ERecoScFiLayer6OverP"};
   std::vector<std::vector<double>> colRange = {{0,23},{0,19},{0,570},{-300,300},{0,280},{0.065,0.45},{-0.5,1},{0,1.5},{0,1.5},{0,2.5e-2},{0,2.5e-2},{0,2.5e-2},{0,2.5e-2},};
 
-  auto d_temp_ele = d1.Filter(is_electron, {"mcparticles"});
-  auto d_temp_pim = d1.Filter(is_piMinus,  {"mcparticles"});
+  auto d_temp_ele = d1.Filter(is_electron, {"MCParticles"});
+  auto d_temp_pim = d1.Filter(is_piMinus,  {"MCParticles"});
 
   int rangeCount = 0;
   for (auto && name : colList){
@@ -792,12 +792,12 @@ void emcal_barrel_pion_rejection_analysis(
   double effPim[7];
   double rejRatios[7];
 
-  auto he_uncut = d1.Filter(is_electron, {"mcparticles"}).Histo1D({"he_uncut",  "P", 7, &binEdges[0]}, "Pthr");
+  auto he_uncut = d1.Filter(is_electron, {"MCParticles"}).Histo1D({"he_uncut",  "P", 7, &binEdges[0]}, "Pthr");
   auto he_cut   = d_ele.Histo1D({"he_cut", "P", 7, &binEdges[0]}, "Pthr");
-  auto hp_uncut = d1.Filter(is_piMinus, {"mcparticles"}).Histo1D({"hp_uncut",  "P", 7, &binEdges[0]}, "Pthr");
+  auto hp_uncut = d1.Filter(is_piMinus, {"MCParticles"}).Histo1D({"hp_uncut",  "P", 7, &binEdges[0]}, "Pthr");
   auto hp_cut   = d_pim.Histo1D({"hp_cut", "P", 7, &binEdges[0]}, "Pthr");
 
-  double ele_eff_total = (double)*d_ele.Count() / (double)*d1.Filter(is_electron, {"mcparticles"}).Count();
+  double ele_eff_total = (double)*d_ele.Count() / (double)*d1.Filter(is_electron, {"MCParticles"}).Count();
 
   he_cut->Divide(he_uncut.GetPtr());//Effienciency 
   hp_uncut->Divide(hp_cut.GetPtr());//Rejection power
diff --git a/benchmarks/imaging_ecal/config.yml b/benchmarks/imaging_ecal/config.yml
index 152a518db99bf290a28f9fe39b195cc87dd8b3df..b269fb0d4c75ae01d28e5559a435d128b812dc40 100644
--- a/benchmarks/imaging_ecal/config.yml
+++ b/benchmarks/imaging_ecal/config.yml
@@ -4,6 +4,7 @@ imaging_ecal_electrons:
   stage: run
   script:
     - bash benchmarks/imaging_ecal/run_emcal_barrel.sh -t emcal_barrel_electrons -p "electron" -n 100
+  allow_failure: true
 
 imaging_ecal_photons:
   extends: .rec_benchmark
@@ -11,6 +12,7 @@ imaging_ecal_photons:
   stage: run
   script:
     - bash benchmarks/imaging_ecal/run_emcal_barrel.sh -t emcal_barrel_photons -p "photon" -n 100
+  allow_failure: true
 
 imaging_ecal_pions:
   extends: .rec_benchmark
@@ -18,6 +20,7 @@ imaging_ecal_pions:
   stage: run
   script:
     - bash benchmarks/imaging_ecal/run_emcal_barrel.sh -t emcal_barrel_pions -p "pion-" -n 100
+  allow_failure: true
 
 imaging_ecal_pion0:
   extends: .rec_benchmark
@@ -25,25 +28,27 @@ imaging_ecal_pion0:
   stage: run
   script:
     - bash benchmarks/imaging_ecal/run_imcal_pion0.sh -t imcal_barrel_pion0 -p "pion0" -n 100
+  allow_failure: true
 
-## SJJ removed as these are causing the CI to fail
-#imaging_ecal_pion_rejection:
-#  extends: .rec_benchmark
-#  timeout: 48 hours
-#  stage: run
-#  script:
-#    - compile_analyses.py imaging_ecal
-#    - bash benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh -t emcal_barrel_pion_rej_electron -p "electron" -n 100
-#    - bash benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh -t emcal_barrel_pion_rej_piminus -p "pion-" -n 100
+imaging_ecal_pion_rejection:
+  extends: .rec_benchmark
+  timeout: 48 hours
+  stage: run
+  script:
+    - compile_analyses.py imaging_ecal
+    - bash benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh -t emcal_barrel_pion_rej_electron -p "electron" -n 100
+    - bash benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh -t emcal_barrel_pion_rej_piminus -p "pion-" -n 100
+  allow_failure: true
 
-#imaging_ecal_pion_rejection:bench:
-#  extends: .rec_benchmark
-#  timeout: 1 hours
-#  stage: benchmarks2
-#  needs:
-#    - ["imaging_ecal_pion_rejection"]
-#  script:
-#    - ls -lhtR
-#    - compile_analyses.py imaging_ecal
-#    - root -b -q benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx+
-#    #- bash run_pion_rej_analysis.sh
+imaging_ecal_pion_rejection:bench:
+  extends: .rec_benchmark
+  timeout: 1 hours
+  stage: benchmarks2
+  needs:
+    - ["imaging_ecal_pion_rejection"]
+  script:
+    - ls -lhtR
+    - compile_analyses.py imaging_ecal
+    - root -b -q benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx+
+    #- bash run_pion_rej_analysis.sh
+  allow_failure: true
diff --git a/benchmarks/imaging_ecal/options/hybrid_cluster.py b/benchmarks/imaging_ecal/options/hybrid_cluster.py
index f4d3388fbf72d76f9903dd68922d800abab89420..1e6405c8c7e58538a930b5b2fc245e7a4f155d47 100644
--- a/benchmarks/imaging_ecal/options/hybrid_cluster.py
+++ b/benchmarks/imaging_ecal/options/hybrid_cluster.py
@@ -44,7 +44,7 @@ geo_service = GeoSvc('GeoSvc', detectors=kwargs['compact'].split(','), OutputLev
 podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
 
 sim_colls = [
-    'mcparticles',
+    'MCParticles',
     'EcalBarrelHits',
     'EcalBarrelScFiHits',
 ]
@@ -72,6 +72,7 @@ imcalreco = ImagingPixelReco('imcal_reco',
         readoutClass='EcalBarrelHits',
         layerField='layer',
         sectorField='module',
+        samplingFraction=kwargs['img_sf'],
         **imcaldaq)
 
 imcalcluster = ImagingTopoCluster('imcal_cluster',
@@ -84,12 +85,10 @@ imcalcluster = ImagingTopoCluster('imcal_cluster',
         sectorDist=3.*cm)
 clusterreco = ImagingClusterReco('imcal_clreco',
         #OutputLevel=DEBUG,
-        inputHitCollection=imcalcluster.inputHitCollection,
-        inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
-        outputLayerCollection='EcalBarrelImagingClustersLayers',
-        outputClusterCollection='EcalBarrelImagingClusters',
-        mcHits="EcalBarrelHits",
-        samplingFraction=kwargs['img_sf'])
+        inputProtoClusters=imcalcluster.outputProtoClusterCollection,
+        outputLayers='EcalBarrelImagingClustersLayers',
+        outputClusters='EcalBarrelImagingClusters',
+        mcHits="EcalBarrelHits")
 
 # scfi layers
 scfi_barrel_daq = dict(
@@ -111,6 +110,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
         layerField="layer",
         sectorField="module",
         localDetFields=["system", "module"], # use local coordinates in each module (stave)
+        samplingFraction=kwargs['scfi_sf'],
         **scfi_barrel_daq)
 
 # merge hits in different layer (projection to local x-y plane)
@@ -132,17 +132,14 @@ scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
         sectorDist=5.*cm)
 
 scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
-       inputHitCollection=scfi_barrel_cl.inputHitCollection,
        inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection,
        outputClusterCollection="EcalBarrelScFiClusters",
-        mcHits="EcalBarrelScFiHits",
-       logWeightBase=6.2,
-       samplingFraction=kwargs['scfi_sf'])
+       logWeightBase=6.2)
 
 # TODO: merge two types of clusters
 
 podout.outputCommands = ['drop *',
-        'keep mcparticles',
+        'keep MCParticles',
         'keep *Reco*',
         'keep *Digi*',
         'keep *Cluster*',
diff --git a/benchmarks/imaging_ecal/options/imaging_2dcluster.py b/benchmarks/imaging_ecal/options/imaging_2dcluster.py
index b4ac4a6d8c336f4b92c8b81aab7aac2f590699d1..70b66255a05c39414c0d256c26c6c0591c56130a 100644
--- a/benchmarks/imaging_ecal/options/imaging_2dcluster.py
+++ b/benchmarks/imaging_ecal/options/imaging_2dcluster.py
@@ -20,7 +20,7 @@ kwargs = dict()
 kwargs['img_sf'] = float(calib_data['sampling_fraction_img'])
 
 # input arguments through environment variables
-kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.root')
+kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.edm4hep.root')
 kwargs['output'] = os.environ.get('CB_EMCAL_REC_FILE', 'barrel_pion0_5GeV_cluster.root')
 kwargs['compact'] = os.environ.get('CB_EMCAL_COMPACT_PATH', '../topside/test.xml')
 kwargs['nev'] = int(os.environ.get('CB_EMCAL_NUMEV', 100))
@@ -34,7 +34,7 @@ print(kwargs)
 geo_service = GeoSvc('GeoSvc', detectors=kwargs['compact'].split(','), OutputLevel=INFO)
 podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
 
-podin = PodioInput('PodioReader', collections=['mcparticles', 'EcalBarrelHits'], OutputLevel=DEBUG)
+podin = PodioInput('PodioReader', collections=['MCParticles', 'EcalBarrelHits'], OutputLevel=DEBUG)
 podout = PodioOutput('out', filename=kwargs['output'])
 
 
@@ -76,7 +76,6 @@ imcal_barrel_cl = IslandCluster('imcal_barrel_cl',
 
 imcal_barrel_clreco = RecoCoG('imcal_barrel_clreco',
         # OutputLevel=DEBUG,
-        inputHitCollection=imcal_barrel_cl.inputHitCollection,
         inputProtoClusterCollection = imcal_barrel_cl.outputProtoClusterCollection,
         outputClusterCollection='EcalBarrelImagingClusters',
         mcHits="EcalBarrelHits",
diff --git a/benchmarks/imaging_ecal/options/imaging_topocluster.py b/benchmarks/imaging_ecal/options/imaging_topocluster.py
index fe243a1876a053d184ee0f3685a2d29bacff59ee..414e65d2c83fafb33637b55983c86eca92779587 100644
--- a/benchmarks/imaging_ecal/options/imaging_topocluster.py
+++ b/benchmarks/imaging_ecal/options/imaging_topocluster.py
@@ -21,7 +21,7 @@ kwargs = dict()
 kwargs['sf'] = float(calib_data['sampling_fraction_img'])
 
 # input arguments through environment variables
-kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.root')
+kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.edm4hep.root')
 kwargs['output'] = os.environ.get('CB_EMCAL_REC_FILE', 'barrel_pion0_5GeV_cluster.root')
 kwargs['compact'] = os.environ.get('CB_EMCAL_COMPACT_PATH', '../topside/test.xml')
 kwargs['nev'] = int(os.environ.get('CB_EMCAL_NUMEV', 100))
@@ -36,7 +36,7 @@ geo_service = GeoSvc("GeoSvc", detectors=kwargs['compact'].split(','), OutputLev
 podioevent = EICDataSvc("EventDataSvc", inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
 out = PodioOutput("out", filename=kwargs['output'])
 
-podioinput = PodioInput("PodioReader", collections=["mcparticles", "EcalBarrelHits"], OutputLevel=DEBUG)
+podioinput = PodioInput("PodioReader", collections=["MCParticles", "EcalBarrelHits"], OutputLevel=DEBUG)
 
 # use the same daq_setting for digi/reco pair
 daq_setting = dict(
@@ -70,7 +70,6 @@ imcalcluster = ImagingTopoCluster("imcal_cluster",
         minClusterNhits=5)
 clusterreco = ImagingClusterReco("imcal_clreco",
         # OutputLevel=DEBUG,
-        inputHitCollection=imcalcluster.inputHitCollection,
         inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
         outputLayerCollection="EcalBarrelImagingClustersLayers",
         outputClusterCollection="EcalBarrelImagingClusters",
diff --git a/benchmarks/imaging_ecal/options/scfi_cluster.py b/benchmarks/imaging_ecal/options/scfi_cluster.py
index 741d124a2480188aeef0ef7769e9648b0c98349f..b7f52ac80ce2e55c92adbd589f47e90cd989c107 100644
--- a/benchmarks/imaging_ecal/options/scfi_cluster.py
+++ b/benchmarks/imaging_ecal/options/scfi_cluster.py
@@ -20,7 +20,7 @@ kwargs = dict()
 kwargs['sf'] = float(calib_data['sampling_fraction_scfi'])
 
 # input arguments through environment variables
-kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.root')
+kwargs['input'] = os.environ.get('CB_EMCAL_SIM_FILE', '../topside/barrel_pion0_5GeV.edm4hep.root')
 kwargs['output'] = os.environ.get('CB_EMCAL_REC_FILE', 'barrel_pion0_5GeV_cluster.root')
 kwargs['compact'] = os.environ.get('CB_EMCAL_COMPACT_PATH', '../topside/test.xml')
 kwargs['nev'] = int(os.environ.get('CB_EMCAL_NUMEV', 100))
@@ -36,7 +36,7 @@ sf = float(os.environ.get('CB_EMCAL_SAMP_FRAC', '1.0'))
 geo_service = GeoSvc("GeoSvc", detectors=kwargs['compact'].split(','), OutputLevel=INFO)
 podioevent = EICDataSvc("EventDataSvc", inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
 
-podin = PodioInput("PodioReader", collections=["mcparticles", "EcalBarrelScFiHits"], OutputLevel=DEBUG)
+podin = PodioInput("PodioReader", collections=["MCParticles", "EcalBarrelScFiHits"], OutputLevel=DEBUG)
 podout = PodioOutput("out", filename=kwargs['output'])
 
 # use the same daq_setting for digi/reco pair
@@ -79,7 +79,6 @@ scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
         localDistXZ=[30*mm, 30*mm])
 
 scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
-       inputHitCollection=scfi_barrel_cl.inputHitCollection,
        inputProtoClusterCollection=scfi_barrel.outputProtoClusterCollection,
        outputClusterCollection="EcalBarrelScFiClusters",
        mcHits="EcalBarrelScFiHits",
diff --git a/benchmarks/imaging_ecal/run_emcal_barrel.sh b/benchmarks/imaging_ecal/run_emcal_barrel.sh
index d54844fbd2109ca6631699e838d408986de057bc..d81539546702af961a46f6410081aa7e433a8d2f 100644
--- a/benchmarks/imaging_ecal/run_emcal_barrel.sh
+++ b/benchmarks/imaging_ecal/run_emcal_barrel.sh
@@ -68,7 +68,7 @@ fi
 export CB_EMCAL_NAME_TAG="${nametag}"
 export CB_EMCAL_GEN_FILE="${CB_EMCAL_NAME_TAG}.hepmc"
 
-export CB_EMCAL_SIM_FILE="sim_${CB_EMCAL_NAME_TAG}.root"
+export CB_EMCAL_SIM_FILE="sim_${CB_EMCAL_NAME_TAG}.edm4hep.root"
 export CB_EMCAL_REC_FILE="rec_${CB_EMCAL_NAME_TAG}.root"
 
 echo "CB_EMCAL_NUMEV = ${CB_EMCAL_NUMEV}"
@@ -85,9 +85,10 @@ fi
 ls -lh ${CB_EMCAL_GEN_FILE}
 
 # Run geant4 simulations
-npsim --runType batch \
+ddsim --runType batch \
       -v WARNING \
       --part.minimalKineticEnergy "0.5*MeV" \
+      --filter.tracker edep0 \
       --numberOfEvents ${CB_EMCAL_NUMEV} \
       --compactFile ${CB_EMCAL_COMPACT_PATH} \
       --inputFiles ${CB_EMCAL_GEN_FILE} \
diff --git a/benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh b/benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh
index a43aaede64d338bc9a506c8a9b72c4ff7c4bd48c..a3e21fb7ed54d932718d893002e1734cc163ab64 100755
--- a/benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh
+++ b/benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh
@@ -68,7 +68,7 @@ fi
 export CB_EMCAL_NAME_TAG="${nametag}"
 export CB_EMCAL_GEN_FILE="${CB_EMCAL_NAME_TAG}.hepmc"
 
-export CB_EMCAL_SIM_FILE="sim_${CB_EMCAL_NAME_TAG}.root"
+export CB_EMCAL_SIM_FILE="sim_${CB_EMCAL_NAME_TAG}.edm4hep.root"
 export CB_EMCAL_REC_FILE="rec_${CB_EMCAL_NAME_TAG}.root"
 
 echo "CB_EMCAL_NUMEV = ${CB_EMCAL_NUMEV}"
@@ -85,9 +85,10 @@ fi
 ls -lh ${CB_EMCAL_GEN_FILE}
 
 # Run geant4 simulations
-npsim --runType batch \
+ddsim --runType batch \
       -v WARNING \
       --part.minimalKineticEnergy "0.5*MeV" \
+      --filter.tracker edep0 \
       --numberOfEvents ${CB_EMCAL_NUMEV} \
       --compactFile ${CB_EMCAL_COMPACT_PATH} \
       --inputFiles ${CB_EMCAL_GEN_FILE} \
diff --git a/benchmarks/imaging_ecal/run_imcal_pion0.sh b/benchmarks/imaging_ecal/run_imcal_pion0.sh
index 64cf9d07222d5ed138eb43edf8e864737e81365f..2866daea64d7f9b19b3ccab52cdcda4ff48f229a 100644
--- a/benchmarks/imaging_ecal/run_imcal_pion0.sh
+++ b/benchmarks/imaging_ecal/run_imcal_pion0.sh
@@ -64,7 +64,7 @@ fi
 export CB_EMCAL_NAME_TAG="${nametag}"
 export CB_EMCAL_GEN_FILE="${CB_EMCAL_NAME_TAG}.hepmc"
 
-export CB_EMCAL_SIM_FILE="sim_${CB_EMCAL_NAME_TAG}.root"
+export CB_EMCAL_SIM_FILE="sim_${CB_EMCAL_NAME_TAG}.edm4hep.root"
 export CB_EMCAL_REC_FILE="rec_${CB_EMCAL_NAME_TAG}.root"
 
 echo "CB_EMCAL_NUMEV = ${CB_EMCAL_NUMEV}"
@@ -81,9 +81,10 @@ fi
 ls -lh ${CB_EMCAL_GEN_FILE}
 
 # Run geant4 simulations
-npsim --runType batch \
+ddsim --runType batch \
       -v WARNING \
       --part.minimalKineticEnergy "0.5*MeV" \
+      --filter.tracker edep0 \
       --numberOfEvents ${CB_EMCAL_NUMEV} \
       --compactFile ${CB_EMCAL_COMPACT_PATH} \
       --inputFiles ${CB_EMCAL_GEN_FILE} \
diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster.py b/benchmarks/imaging_ecal/scripts/draw_cluster.py
index 02cc76832e35ee34230edb253f2720a4571ef311..e3b81460d2fb013abe3f7bffe40b4caec06be5d9 100755
--- a/benchmarks/imaging_ecal/scripts/draw_cluster.py
+++ b/benchmarks/imaging_ecal/scripts/draw_cluster.py
@@ -148,10 +148,12 @@ if __name__ == '__main__':
     df['eta'] = -np.log(np.tan(df['theta'].values/1000./2.))
 
     # Read all mc particles
-    dfallmcp = get_all_mcp(args.file, args.iev, 'mcparticles')
+    dfallmcp = get_all_mcp(args.file, args.iev, 'MCParticles')
     pdgbase = ROOT.TDatabasePDG()
     # Select decaying particles
-    dftemp = dfallmcp[dfallmcp['g4Parent'] == 1.0]
+    # FIXME g4Parent not in edm4hep
+    #dftemp = dfallmcp[dfallmcp['g4Parent'] == 1.0]
+    dftemp = dfallmcp
     if len(dftemp) > 0:
         dfdecaymcp = dftemp.copy()
         for iptl in [0, len(dfdecaymcp) - 1]:
@@ -169,7 +171,7 @@ if __name__ == '__main__':
         dfdecaymcp['eta'] = -np.log(np.tan(dfdecaymcp['theta'].values/1000./2.))
 
     # truth
-    dfmcp = get_mcp_simple(args.file, args.iev, 'mcparticles').iloc[0]
+    dfmcp = get_mcp_simple(args.file, args.iev, 'MCParticles').iloc[0]
     #pdgbase = ROOT.TDatabasePDG()
     inpart = pdgbase.GetParticle(int(dfmcp['pid']))
     if inpart:
diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py b/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py
index 42b80d4305f3bdc8927346196274dc3c4b2dc985..0a37a5222b2332112bd356f16afb814e03598fe6 100644
--- a/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py
+++ b/benchmarks/imaging_ecal/scripts/draw_cluster_layers.py
@@ -98,7 +98,7 @@ if __name__ == '__main__':
     df['eta'] = -np.log(np.tan(df['theta'].values/1000./2.))
 
     # truth
-    dfmcp = get_mcp_simple(args.file, args.iev, 'mcparticles').iloc[0]
+    dfmcp = get_mcp_simple(args.file, args.iev, 'MCParticles').iloc[0]
     pdgbase = ROOT.TDatabasePDG()
     inpart = pdgbase.GetParticle(int(dfmcp['pid']))
     print("Incoming particle = {}, pdgcode = {}, charge = {}, mass = {}"\
diff --git a/benchmarks/imaging_ecal/scripts/utils.py b/benchmarks/imaging_ecal/scripts/utils.py
index 0f4e40a239a6aced7430162e11f932a8bffab05a..2872bf36b712b5c7ddd52e8da1c0ef1edee6aec0 100644
--- a/benchmarks/imaging_ecal/scripts/utils.py
+++ b/benchmarks/imaging_ecal/scripts/utils.py
@@ -37,7 +37,7 @@ def load_root_macros(arg_macros):
 
 
 # read mc particles from root file
-def get_mcp_data(path, evnums=None, branch='mcparticles'):
+def get_mcp_data(path, evnums=None, branch='MCParticles'):
     f = ROOT.TFile(path)
     events = f.events
     if evnums is None:
@@ -55,13 +55,13 @@ def get_mcp_data(path, evnums=None, branch='mcparticles'):
         events.GetEntry(iev)
         # extract full mc particle data
         for part in getattr(events, branch):
-            dbuf[idb] = (iev, part.ps.x, part.ps.y, part.ps.z, part.pdgID, part.status)
+            dbuf[idb] = (iev, part.momentum.x, part.momentum.y, part.momentum.z, part.PDG, part.simulatorStatus)
             idb += 1
     return pd.DataFrame(data=dbuf[:idb], columns=['event', 'px', 'py', 'pz', 'pid', 'status'])
 
 
 # read mc particles from root file
-def get_mcp_simple(path, evnums=None, branch='mcparticles'):
+def get_mcp_simple(path, evnums=None, branch='MCParticles'):
     f = ROOT.TFile(path)
     events = f.events
     if evnums is None:
@@ -79,14 +79,14 @@ def get_mcp_simple(path, evnums=None, branch='mcparticles'):
         events.GetEntry(iev)
         # extract full mc particle data
         part = getattr(events, branch)[2]
-        dbuf[idb] = (iev, part.ps.x, part.ps.y, part.ps.z, part.pdgID, part.status)
+        dbuf[idb] = (iev, part.momentum.x, part.momentum.y, part.momentum.z, part.PDG, part.simulatorStatus)
         idb += 1
     return pd.DataFrame(data=dbuf[:idb], columns=['event', 'px', 'py', 'pz', 'pid', 'status'])
 
 #######################################
 # read all mc particles from root file
 #######################################
-def get_all_mcp(path, evnums=None, branch='mcparticles'):
+def get_all_mcp(path, evnums=None, branch='MCParticles'):
     f = ROOT.TFile(path)
     events = f.events
     if evnums is None:
@@ -94,7 +94,7 @@ def get_all_mcp(path, evnums=None, branch='mcparticles'):
     elif isinstance(evnums, int):
         evnums = [evnums]
 
-    dbuf = np.zeros(shape=(2000*len(evnums), 10))
+    dbuf = np.zeros(shape=(2000*len(evnums), 9))
     idb = 0
     for iev in evnums:
         if iev >= events.GetEntries():
@@ -104,10 +104,10 @@ def get_all_mcp(path, evnums=None, branch='mcparticles'):
         events.GetEntry(iev)
         # extract mc particle data
         for ptl in getattr(events, branch):
-            dbuf[idb] = (iev, ptl.ps.x, ptl.ps.y, ptl.ps.z, ptl.pdgID, ptl.status, ptl.g4Parent, ptl.ve.x, ptl.ve.y, ptl.ve.z)
+            dbuf[idb] = (iev, ptl.momentum.x, ptl.momentum.y, ptl.momentum.z, ptl.PDG, ptl.simulatorStatus, ptl.endpoint.x, ptl.endpoint.y, ptl.endpoint.z)
             idb += 1
     
-    return pd.DataFrame(data=dbuf[:idb], columns=['event', 'px', 'py', 'pz', 'pid', 'status', 'g4Parent', 'vex', 'vey', 'vez'])
+    return pd.DataFrame(data=dbuf[:idb], columns=['event', 'px', 'py', 'pz', 'pid', 'status', 'vex', 'vey', 'vez'])
 
 # read hits data from root file
 def get_hits_data(path, evnums=None, branch='RecoEcalBarreImaginglHits'):
diff --git a/benchmarks/imaging_shower_ML/options/imaging_ml_data.py b/benchmarks/imaging_shower_ML/options/imaging_ml_data.py
index 6bf9d7c1022972075c272db361a29df7ded24c41..12bc3ef49903d059c9ab20d12d696c15b98efbc8 100644
--- a/benchmarks/imaging_shower_ML/options/imaging_ml_data.py
+++ b/benchmarks/imaging_shower_ML/options/imaging_ml_data.py
@@ -34,7 +34,7 @@ sf = float(os.environ.get('JUGGLER_SAMP_FRAC', '1.0'))
 geo_service  = GeoSvc('GeoSvc', detectors=[f.strip() for f in kwargs['compact'].split(',')])
 podev = EICDataSvc('EventDataSvc', inputs=[f.strip() for f in kwargs['input'].split(',')])
 
-podin = PodioInput('PodioReader', collections=['mcparticles', 'EcalBarrelHits', 'EcalBarrelScFiHits'])
+podin = PodioInput('PodioReader', collections=['MCParticles', 'EcalBarrelHits', 'EcalBarrelScFiHits'])
 podout = PodioOutput('out', filename=kwargs['output'])
 
 # Central Barrel Ecal (Imaging Cal.)
@@ -127,7 +127,7 @@ podout.outputCommands = [
 #        'keep *',
         'drop *',
         'keep EcalBarrel*Reco',
-        'keep mcparticles*',
+        'keep MCParticles*',
         'keep EcalBarrel*ML',
         'keep *Corr',
 ]
diff --git a/benchmarks/imaging_shower_ML/scripts/check_edep_dists.py b/benchmarks/imaging_shower_ML/scripts/check_edep_dists.py
index 404855b1c6fac23883ad7b3537d67b4013edc123..76dca74439d1ee1086bf3dc6dd4a94a8f27a9b1d 100644
--- a/benchmarks/imaging_shower_ML/scripts/check_edep_dists.py
+++ b/benchmarks/imaging_shower_ML/scripts/check_edep_dists.py
@@ -63,7 +63,7 @@ if __name__ == '__main__':
     parser.add_argument('-m', '--macros', type=str, default='rootlogon.C', dest='macros',
                          help='root macros to load (accept multiple paths separated by \",\")')
     parser.add_argument('--branch', type=str, default='EcalBarrelImagingHitsReco', help='name of data branch (eic::CalorimeterHitCollection)')
-    parser.add_argument('--truth-branch', type=str, default='mcparticles', help='name of truth mc branch')
+    parser.add_argument('--truth-branch', type=str, default='MCParticles', help='name of truth mc branch')
     parser.add_argument('--edep-max', type=float, default=0., help='maximum edep (GeV) to plot')
     parser.add_argument('--edep-nbins', type=int, default=200, help='number of bins')
     parser.add_argument('--name-tag', type=str, default='test', help='name tag to save the file')
@@ -73,17 +73,17 @@ if __name__ == '__main__':
     os.makedirs(args.outdir, exist_ok=True)
     load_root_macros(args.macros)
 
-    # read data and mcparticles
+    # read data and MCParticles
     rdf = ROOT.RDataFrame("events", args.file)
 
     mc_branch = args.truth_branch
-    dfm = flatten_collection(rdf, mc_branch, ['genStatus', 'pdgID', 'ps.x', 'ps.y', 'ps.z', 'mass'])
+    dfm = flatten_collection(rdf, mc_branch, ['generatorStatus', 'PDG', 'momentum.x', 'momentum.y', 'momentum.z', 'mass'])
     dfm.rename(columns={c: c.replace(mc_branch + '.', '') for c in dfm.columns}, inplace=True)
     # selete incident particles
-    dfm = dfm[dfm['genStatus'].isin([0, 1])]
+    dfm = dfm[dfm['generatorStatus'].isin([0, 1])]
     # NOTE: assumed single particles
     dfm = dfm.groupby('event').first()
-    # p, theta, phi, pT, eta = cartesian_to_polar(*dfm[['ps.x', 'ps.y', 'ps.z']].values.T)
+    # p, theta, phi, pT, eta = cartesian_to_polar(*dfm[['momentum.x', 'momentum.y', 'momentum.z']].values.T)
 
     if args.sim:
         df = flatten_collection(rdf, args.branch, ['energyDeposit'])
@@ -111,12 +111,12 @@ if __name__ == '__main__':
 
     hist_vals, hist_cols = [], []
     pdgbase = ROOT.TDatabasePDG()
-    for pdgid in dfm['pdgID'].unique():
+    for pdgid in dfm['PDG'].unique():
         particle = pdgbase.GetParticle(int(pdgid))
         if not particle:
             print("Unknown pdgcode {}, they are ignored".format(int(pdgid)))
             continue
-        events_indices = dfm[dfm.loc[:, 'pdgID'] == pdgid].index.unique()
+        events_indices = dfm[dfm.loc[:, 'PDG'] == pdgid].index.unique()
         print("{} entries of particle {}".format(len(events_indices), particle.GetName()))
 
         dfe_part = dfe.loc[dfe['event'].isin(events_indices)]
diff --git a/benchmarks/imaging_shower_ML/scripts/prepare_tf_dataset.py b/benchmarks/imaging_shower_ML/scripts/prepare_tf_dataset.py
index d94d1cb6ad9f177625367a8aa134e8da5d917018..75bfb05cd091521e515fdf56841fa4bbb1e49589 100644
--- a/benchmarks/imaging_shower_ML/scripts/prepare_tf_dataset.py
+++ b/benchmarks/imaging_shower_ML/scripts/prepare_tf_dataset.py
@@ -73,7 +73,7 @@ if __name__ == '__main__':
     os.makedirs(args.outdir, exist_ok=True)
     load_root_macros(args.macros)
 
-    # read data and mcparticles
+    # read data and MCParticles
     rdf = ROOT.RDataFrame("events", args.file)
     df = flatten_collection(rdf, args.branch, ['layer', 'energy', 'position.x', 'position.y', 'position.z'])
     df.rename(columns={c: c.replace(args.branch + '.', '') for c in df.columns}, inplace=True)
@@ -85,13 +85,13 @@ if __name__ == '__main__':
     df.loc[:, 'rc'] = rc
     df.loc[:, 'eta'] = eta
 
-    dfm = flatten_collection(rdf, 'mcparticles', ['genStatus', 'pdgID', 'ps.x', 'ps.y', 'ps.z', 'mass'])
-    dfm.rename(columns={c: c.replace('mcparticles.', '') for c in dfm.columns}, inplace=True)
+    dfm = flatten_collection(rdf, 'MCParticles', ['generatorStatus', 'PDG', 'momentum.x', 'momentum.y', 'momentum.z', 'mass'])
+    dfm.rename(columns={c: c.replace('MCParticles.', '') for c in dfm.columns}, inplace=True)
     # selete incident particles
-    dfm = dfm[dfm['genStatus'].isin([0, 1])]
+    dfm = dfm[dfm['generatorStatus'].isin([0, 1])]
     # NOTE: assumed single particles
     dfm = dfm.groupby('event').first()
-    p, theta, phi, pT, eta = cartesian_to_polar(*dfm[['ps.x', 'ps.y', 'ps.z']].values.T)
+    p, theta, phi, pT, eta = cartesian_to_polar(*dfm[['momentum.x', 'momentum.y', 'momentum.z']].values.T)
     dfm.loc[:, 'p'] = p
     dfm.loc[:, 'theta'] = theta
     dfm.loc[:, 'phi'] = phi
diff --git a/benchmarks/imaging_shower_ML/sim_rec_tag.py b/benchmarks/imaging_shower_ML/sim_rec_tag.py
index 0cd6ce6bb08e330b7b004885c64d1beea0c641d1..e7425a61f8fd10ac7f7b4f3b9dd7602eed878b45 100755
--- a/benchmarks/imaging_shower_ML/sim_rec_tag.py
+++ b/benchmarks/imaging_shower_ML/sim_rec_tag.py
@@ -34,7 +34,7 @@ for mdir in ['gen_data', 'sim_data', 'rec_data', 'tag_data']:
     os.makedirs(mdir, exist_ok=True)
 
 gen_file = os.path.join('gen_data', '{nametag}_{pmin}_{pmax}.hepmc'.format(**kwargs))
-sim_file = os.path.join('sim_data', '{nametag}_{pmin}_{pmax}.root'.format(**kwargs))
+sim_file = os.path.join('sim_data', '{nametag}_{pmin}_{pmax}.edm4hep.root'.format(**kwargs))
 rec_file = os.path.join('rec_data', '{nametag}_{pmin}_{pmax}.root'.format(**kwargs))
 tag_dir = os.path.join('tag_data', '{nametag}_{pmin}_{pmax}'.format(**kwargs))
 
@@ -52,8 +52,9 @@ if 'sim' in procs:
             '--particles', args.particles]
     subprocess.run(gen_cmd)
     # simulation
-    sim_cmd = ['npsim',
+    sim_cmd = ['ddsim',
             '--part.minimalKineticEnergy', '1*TeV',
+            '--filter.tracker', 'edep0',
             '--numberOfEvents', '{}'.format(args.nev),
             '--runType', 'batch',
             # '--physics.list', args.physics_list,
diff --git a/benchmarks/rich/forward_hadrons.sh b/benchmarks/rich/forward_hadrons.sh
index d71bcece3dd0e97a10149aaf81e472961de06b57..08de28e9b2fede48ade6196228c21c5f3a4d6338 100644
--- a/benchmarks/rich/forward_hadrons.sh
+++ b/benchmarks/rich/forward_hadrons.sh
@@ -20,7 +20,7 @@ fi
 export JUGGLER_FILE_NAME_TAG="rich_forward_hadrons"
 export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
 
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
 export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
 
 echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
@@ -58,8 +58,9 @@ fi
 exit 0 
 
 
-npsim --runType batch \
+ddsim --runType batch \
       --part.minimalKineticEnergy 1000*GeV  \
+      --filter.tracker edep0 \
       -v WARNING \
       --numberOfEvents ${JUGGLER_N_EVENTS} \
       --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
diff --git a/benchmarks/rich/options/rich_reco.py b/benchmarks/rich/options/rich_reco.py
index 43a3d4bc39b101fd26b99776d5529df2739d7b86..d851004c79a071ce92ccbc25b441cff62d600eaf 100644
--- a/benchmarks/rich/options/rich_reco.py
+++ b/benchmarks/rich/options/rich_reco.py
@@ -23,12 +23,12 @@ from Configurables import Jug__Reco__PhotoMultiplierReco as PhotoMultiplierReco
 from Configurables import Jug__Reco__PhotoRingClusters as PhotoRingClusters
 
 qe_data = [(1.0, 0.25), (7.5, 0.25),]
-podioinput = PodioInput("PodioReader", collections=["mcparticles", "ForwardRICHHits"], OutputLevel=DEBUG)
+podioinput = PodioInput("PodioReader", collections=["MCParticles", "ForwardRICHHits"], OutputLevel=DEBUG)
 
 pmtdigi = PhotoMultiplierDigi(inputHitCollection="ForwardRICHHits", outputHitCollection="DigiForwardRICHHits",
                               quantumEfficiency=[(a*units.eV, b) for a, b in qe_data])
 pmtreco = PhotoMultiplierReco(inputHitCollection="DigiForwardRICHHits", outputHitCollection="RecoForwardRICHHits")
-richcluster = PhotoRingClusters(inputHitCollection="RecoForwardRICHHits", #inputTrackCollection="mcparticles",
+richcluster = PhotoRingClusters(inputHitCollection="RecoForwardRICHHits", #inputTrackCollection="MCParticles",
                                 outputClusterCollection="RICHClusters")
 
 out = PodioOutput("out", filename=output_rec_file)
diff --git a/benchmarks/track_finding/config.yml b/benchmarks/track_finding/config.yml
index a5c029010a2e7f8cf64b8070c4a5db07ad463e80..29be4a8be34fc7f530e81ab089fa0a3bf58c0a1e 100644
--- a/benchmarks/track_finding/config.yml
+++ b/benchmarks/track_finding/config.yml
@@ -1,5 +1,13 @@
+track_finding:compile:
+  extends: .compile_benchmark
+  stage: compile
+  script:
+    - compile_analyses.py --dir scripts track_finding
+
 track_finding:multiple_tracks:
   extends: .rec_benchmark
+  needs:
+    - ["track_finding:compile"]
   stage: run
   timeout: 24 hours
   script:
diff --git a/benchmarks/track_finding/multiple_tracks.sh b/benchmarks/track_finding/multiple_tracks.sh
index b1a2ffcfea9d08944ad9f82c197fec2e16009a20..914df3b10520be18edca90eed4de499409bd854f 100644
--- a/benchmarks/track_finding/multiple_tracks.sh
+++ b/benchmarks/track_finding/multiple_tracks.sh
@@ -48,14 +48,14 @@ fi
 export JUGGLER_FILE_NAME_TAG="multiple_tracks"
 export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
 
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
 export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
 
 echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
 echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
 
 ## generate the input events
-root -b -q "benchmarks/track_finding/scripts/gen_multiple_tracks.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+root -b -q "benchmarks/track_finding/scripts/gen_multiple_tracks.cxx+(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running script"
   exit 1
@@ -63,8 +63,9 @@ fi
 
 echo "Running geant4 simulation"
 ## run geant4 simulations
-npsim --runType batch \
+ddsim --runType batch \
   --part.minimalKineticEnergy 1000*GeV  \
+  --filter.tracker edep0 \
   -v WARNING \
   --numberOfEvents ${JUGGLER_N_EVENTS} \
   --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
@@ -91,7 +92,7 @@ rootls -t ${JUGGLER_REC_FILE}
 
 mkdir -p results/track_finding
 
-root -b -q "benchmarks/track_finding/scripts/rec_multiple_tracks.cxx(\"${JUGGLER_REC_FILE}\")"
+root -b -q "benchmarks/track_finding/scripts/rec_multiple_tracks.cxx+(\"${JUGGLER_REC_FILE}\")"
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running root script"
   exit 1
diff --git a/benchmarks/track_finding/options/track_reconstruction.py b/benchmarks/track_finding/options/track_reconstruction.py
index 884ccaf0811e6264fa9cd68db5dbfc6d89b1b908..88a34c61b65b40584dac61d008c87a86c7dc6634 100644
--- a/benchmarks/track_finding/options/track_reconstruction.py
+++ b/benchmarks/track_finding/options/track_reconstruction.py
@@ -22,11 +22,24 @@ input_sim_file  = str(os.environ["JUGGLER_SIM_FILE"])
 output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
 n_events = str(os.environ["JUGGLER_N_EVENTS"])
 
-geo_service  = GeoSvc("GeoSvc", detectors=["{}/{}.xml".format(detector_path,detector_name)], OutputLevel=WARNING)
+## note: old version of material map is called material-maps.XXX, new version is materials-map.XXX
+##       these names are somewhat inconsistent, and should probably all be renamed to 'material-map.XXX'
+##       FIXME
+if detector_version == 'acadia':
+    geo_service = GeoSvc("GeoSvc",
+                         detectors=["{}/{}.xml".format(detector_path,detector_name)],
+                         materials="config/material-maps.json",
+                         OutputLevel=WARNING)
+else:
+    geo_service = GeoSvc("GeoSvc",
+                         detectors=["{}/{}.xml".format(detector_path,detector_name)],
+                         materials="calibrations/materials-map.cbor",
+                         OutputLevel=WARNING)
 podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=WARNING)
 
 from Configurables import PodioInput
 
+from Configurables import Jug__Digi__SimTrackerHitsCollector as SimTrackerHitsCollector
 from Configurables import Jug__Digi__SiliconTrackerDigi as TrackerDigi
 
 from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
@@ -46,50 +59,104 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 
 algorithms = [ ]
 
-input_collections = ['mcparticles','TrackerEndcapHits','TrackerBarrelHits','VertexBarrelHits','GEMTrackerEndcapHits']
+tracker_endcap_collections = [
+    'TrackerEndcapHits1',
+    'TrackerEndcapHits2',
+    'TrackerEndcapHits3',
+    'TrackerEndcapHits4',
+    'TrackerEndcapHits5',
+    'TrackerEndcapHits6'
+]
+tracker_barrel_collections = [
+    'TrackerBarrelHits'
+]
+vertex_barrel_collections = [
+    'VertexBarrelHits'
+]
+gem_endcap_collections = [
+    'GEMTrackerEndcapHits1',
+    'GEMTrackerEndcapHits2',
+    'GEMTrackerEndcapHits3'
+]
+vertex_endcap_collections = [
+    'VertexEndcapHits'
+]
+mpgd_barrel_collections = [
+    'MPGDTrackerBarrelHits1',
+    'MPGDTrackerBarrelHits2'
+]
+
+input_collections = ['MCParticles'] + tracker_endcap_collections + tracker_barrel_collections + vertex_barrel_collections + gem_endcap_collections
+
 if 'acadia' in detector_version:
-    input_collections.append('VertexEndcapHits')
+    input_collections += vertex_endcap_collections
 else:
-    input_collections.append('MPGDTrackerBarrelHits')
-podioinput = PodioInput("PodioReader", 
-                        collections=input_collections)#, OutputLevel=DEBUG)
+    input_collections += mpgd_barrel_collections
+
+podioinput = PodioInput("PodioReader",
+                        collections=input_collections)
 algorithms.append( podioinput )
 
-trk_b_digi = TrackerDigi("trk_b_digi", 
-        inputHitCollection="TrackerBarrelHits",
-        outputHitCollection="TrackerBarrelRawHits",
+trk_b_coll = SimTrackerHitsCollector("trk_b_coll",
+        inputSimTrackerHits = tracker_barrel_collections,
+        outputSimTrackerHits = "TrackerBarrelAllHits")
+algorithms.append( trk_b_coll )
+trk_b_digi = TrackerDigi("trk_b_digi",
+        inputHitCollection = trk_b_coll.outputSimTrackerHits,
+        outputHitCollection = "TrackerBarrelRawHits",
         timeResolution=8)
 algorithms.append( trk_b_digi )
-trk_ec_digi = TrackerDigi("trk_ec_digi", 
-        inputHitCollection="TrackerEndcapHits",
-        outputHitCollection="TrackerEndcapRawHits",
+
+trk_ec_coll = SimTrackerHitsCollector("trk_ec_coll",
+        inputSimTrackerHits = tracker_endcap_collections,
+        outputSimTrackerHits = "TrackerEndcapAllHits")
+algorithms.append( trk_ec_coll )
+trk_ec_digi = TrackerDigi("trk_ec_digi",
+        inputHitCollection = trk_ec_coll.outputSimTrackerHits,
+        outputHitCollection = "TrackerEndcapRawHits",
         timeResolution=8)
 algorithms.append( trk_ec_digi )
 
-vtx_b_digi = TrackerDigi("vtx_b_digi", 
-        inputHitCollection="VertexBarrelHits",
-        outputHitCollection="VertexBarrelRawHits",
+vtx_b_coll = SimTrackerHitsCollector("vtx_b_coll",
+        inputSimTrackerHits = vertex_barrel_collections,
+        outputSimTrackerHits = "VertexBarrelAllHits")
+algorithms.append( vtx_b_coll )
+vtx_b_digi = TrackerDigi("vtx_b_digi",
+        inputHitCollection = vtx_b_coll.outputSimTrackerHits,
+        outputHitCollection = "VertexBarrelRawHits",
         timeResolution=8)
 algorithms.append( vtx_b_digi )
 
 if 'acadia' in detector_version:
-    vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
-            inputHitCollection="VertexEndcapHits",
-            outputHitCollection="VertexEndcapRawHits",
+    vtx_ec_coll = SimTrackerHitsCollector("vtx_ec_coll",
+            inputSimTrackerHits = vertex_endcap_collections,
+            outputSimTrackerHits = "VertexEndcapAllHits")
+    algorithms.append( vtx_ec_coll )
+    vtx_ec_digi = TrackerDigi("vtx_ec_digi",
+            inputHitCollection = vtx_ec_coll.outputSimTrackerHits,
+            outputHitCollection = "VertexEndcapRawHits",
             timeResolution=8)
     algorithms.append( vtx_ec_digi )
 else:
-    mm_b_digi = TrackerDigi("mm_b_digi", 
-            inputHitCollection="MPGDTrackerBarrelHits",
-            outputHitCollection="MPGDTrackerBarrelRawHits",
+    mm_b_coll = SimTrackerHitsCollector("mm_b_coll",
+            inputSimTrackerHits = mpgd_barrel_collections,
+            outputSimTrackerHits = "MPGDTrackerBarrelAllHits")
+    algorithms.append( mm_b_coll )
+    mm_b_digi = TrackerDigi("mm_b_digi",
+            inputHitCollection = mm_b_coll.outputSimTrackerHits,
+            outputHitCollection = "MPGDTrackerBarrelRawHits",
             timeResolution=8)
     algorithms.append( mm_b_digi )
 
+gem_ec_coll = SimTrackerHitsCollector("gem_ec_coll",
+        inputSimTrackerHits = gem_endcap_collections,
+        outputSimTrackerHits = "GEMTrackerEndcapAllHits")
+algorithms.append( gem_ec_coll )
 gem_ec_digi = TrackerDigi("gem_ec_digi",
-        inputHitCollection="GEMTrackerEndcapHits",
-        outputHitCollection="GEMTrackerEndcapRawHits",
+        inputHitCollection = gem_ec_coll.outputSimTrackerHits,
+        outputHitCollection = "GEMTrackerEndcapRawHits",
         timeResolution=10)
-algorithms.append(gem_ec_digi)
+algorithms.append( gem_ec_digi )
 
 # Tracker and vertex reconstruction
 trk_b_reco = TrackerHitReconstruction("trk_b_reco",
@@ -158,7 +225,7 @@ algorithms.append( sourcelinker )
 
 ## Track param init
 truth_trk_init = TrackParamTruthInit("truth_trk_init",
-        inputMCParticles="mcparticles",
+        inputMCParticles="MCParticles",
         outputInitialTrackParameters="InitTrackParams")
         #OutputLevel=DEBUG)
 algorithms.append( truth_trk_init )
@@ -204,7 +271,7 @@ out.outputCommands = ["keep *",
         "drop trajectories",
         "drop outputSourceLinks",
         "drop outputInitialTrackParameters",
-        "keep mcparticles"
+        "keep MCParticles"
         ]
 algorithms.append(out)
 
diff --git a/benchmarks/track_finding/scripts/gen_multiple_tracks.cxx b/benchmarks/track_finding/scripts/gen_multiple_tracks.cxx
index 8596e929d190d7b023f37f3fbbf7da07b88508cb..185bb741a23528d9d3bed93054128de1c7726b58 100644
--- a/benchmarks/track_finding/scripts/gen_multiple_tracks.cxx
+++ b/benchmarks/track_finding/scripts/gen_multiple_tracks.cxx
@@ -8,6 +8,7 @@
 #include <cmath>
 
 #include "TMath.h"
+#include "TRandom.h"
 
 #include "common_bench/particles.h"
 
diff --git a/benchmarks/track_finding/scripts/rec_multiple_tracks.cxx b/benchmarks/track_finding/scripts/rec_multiple_tracks.cxx
index 9b3e551205fed66a47e04e54e3ac2193202d8659..2f1fa6ad80525675a456603d846a6ce4020f0246 100644
--- a/benchmarks/track_finding/scripts/rec_multiple_tracks.cxx
+++ b/benchmarks/track_finding/scripts/rec_multiple_tracks.cxx
@@ -2,13 +2,15 @@
 #include "TCanvas.h"
 #include "TLegend.h"
 #include "TH1D.h"
+#include "THStack.h"
 #include "TProfile.h"
+#include "Math/Vector4D.h"
 
 #include <iostream>
 R__LOAD_LIBRARY(libJugBase.so)
 R__LOAD_LIBRARY(libeicd.so)
 R__LOAD_LIBRARY(libDD4pod.so)
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/TrackParametersCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
@@ -26,10 +28,10 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 };
 
 
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+std::vector<float> pt (std::vector<edm4hep::MCParticleData> 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));
+    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
   }
   return result;
 }
@@ -48,11 +50,11 @@ auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
   }
   return result;
 };
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto fourvec = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv);
   }
   return result;
@@ -87,11 +89,11 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root")
   ROOT::EnableImplicitMT();
   ROOT::RDataFrame df("events", fname);
 
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
+  auto df0 = df.Define("isThrown", "MCParticles.generatorStatus == 1")
+                 .Define("thrownParticles", "MCParticles[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
                  .Define("nThrown", "thrownParticles.size()")
-                 .Define("nProto", "outputProtoTracks.size()")
+                 .Define("nProto", "nProtoTracks")
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
                  .Define("theta0", "theta_thrown[0]")
diff --git a/benchmarks/track_fitting/config.yml b/benchmarks/track_fitting/config.yml
index 2363c6a73ffe82ebee142bc8b105509914fb8bdb..30a6b75c8a23640c5cabbd66e7283d05f13d6de2 100644
--- a/benchmarks/track_fitting/config.yml
+++ b/benchmarks/track_fitting/config.yml
@@ -1,6 +1,14 @@
+track_fitting:compile:
+  extends: .compile_benchmark
+  stage: compile
+  script:
+    - compile_analyses.py --dir scripts track_fitting
+
 track_fitting:single_tracks:
   extends: .rec_benchmark
   stage: run
+  needs:
+    - ["track_fitting:compile"]
   script:
     - bash benchmarks/track_fitting/single_tracks.sh
       
diff --git a/benchmarks/track_fitting/options/track_reconstruction.py b/benchmarks/track_fitting/options/track_reconstruction.py
index 9622ba8981e33c7e5045a42ff254bc20e3d6ac47..885aa2741a7781ad3bc79a1e574073226f005d38 100644
--- a/benchmarks/track_fitting/options/track_reconstruction.py
+++ b/benchmarks/track_fitting/options/track_reconstruction.py
@@ -22,18 +22,30 @@ input_sim_file  = str(os.environ["JUGGLER_SIM_FILE"])
 output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
 n_events = str(os.environ["JUGGLER_N_EVENTS"])
 
-geo_service  = GeoSvc("GeoSvc", detectors=["{}/{}.xml".format(detector_path,detector_name)], OutputLevel=WARNING)
+## note: old version of material map is called material-maps.XXX, new version is materials-map.XXX
+##       these names are somewhat inconsistent, and should probably all be renamed to 'material-map.XXX'
+##       FIXME
+if detector_version == 'acadia':
+    geo_service = GeoSvc("GeoSvc",
+                         detectors=["{}/{}.xml".format(detector_path,detector_name)],
+                         materials="config/material-maps.json",
+                         OutputLevel=WARNING)
+else:
+    geo_service = GeoSvc("GeoSvc",
+                         detectors=["{}/{}.xml".format(detector_path,detector_name)],
+                         materials="calibrations/materials-map.cbor",
+                         OutputLevel=WARNING)
 podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=WARNING)
 
 from Configurables import PodioInput
 
+from Configurables import Jug__Digi__SimTrackerHitsCollector as SimTrackerHitsCollector
 from Configurables import Jug__Digi__SiliconTrackerDigi as TrackerDigi
 
 from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
 from Configurables import Jug__Reco__TrackingHitsCollector2 as TrackingHitsCollector
 
 from Configurables import Jug__Reco__TrackerSourceLinker as TrackerSourceLinker
-#from Configurables import Jug__Reco__TrackingHitsSourceLinker as TrackingHitsSourceLinker
 from Configurables import Jug__Reco__TrackParamTruthInit as TrackParamTruthInit
 from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit
 from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit
@@ -46,50 +58,104 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 
 algorithms = [ ]
 
-input_collections = ['mcparticles','TrackerEndcapHits','TrackerBarrelHits','VertexBarrelHits','GEMTrackerEndcapHits']
+tracker_endcap_collections = [
+    'TrackerEndcapHits1',
+    'TrackerEndcapHits2',
+    'TrackerEndcapHits3',
+    'TrackerEndcapHits4',
+    'TrackerEndcapHits5',
+    'TrackerEndcapHits6'
+]
+tracker_barrel_collections = [
+    'TrackerBarrelHits'
+]
+vertex_barrel_collections = [
+    'VertexBarrelHits'
+]
+gem_endcap_collections = [
+    'GEMTrackerEndcapHits1',
+    'GEMTrackerEndcapHits2',
+    'GEMTrackerEndcapHits3'
+]
+vertex_endcap_collections = [
+    'VertexEndcapHits'
+]
+mpgd_barrel_collections = [
+    'MPGDTrackerBarrelHits1',
+    'MPGDTrackerBarrelHits2'
+]
+
+input_collections = ['MCParticles'] + tracker_endcap_collections + tracker_barrel_collections + vertex_barrel_collections + gem_endcap_collections
+
 if 'acadia' in detector_version:
-    input_collections.append('VertexEndcapHits')
+    input_collections += vertex_endcap_collections
 else:
-    input_collections.append('MPGDTrackerBarrelHits')
-podioinput = PodioInput("PodioReader", 
-                        collections=input_collections)#, OutputLevel=DEBUG)
+    input_collections += mpgd_barrel_collections
+
+podioinput = PodioInput("PodioReader",
+                        collections=input_collections)
 algorithms.append( podioinput )
 
-trk_b_digi = TrackerDigi("trk_b_digi", 
-        inputHitCollection="TrackerBarrelHits",
-        outputHitCollection="TrackerBarrelRawHits",
+trk_b_coll = SimTrackerHitsCollector("trk_b_coll",
+        inputSimTrackerHits = tracker_barrel_collections,
+        outputSimTrackerHits = "TrackerBarrelAllHits")
+algorithms.append( trk_b_coll )
+trk_b_digi = TrackerDigi("trk_b_digi",
+        inputHitCollection = trk_b_coll.outputSimTrackerHits,
+        outputHitCollection = "TrackerBarrelRawHits",
         timeResolution=8)
 algorithms.append( trk_b_digi )
-trk_ec_digi = TrackerDigi("trk_ec_digi", 
-        inputHitCollection="TrackerEndcapHits",
-        outputHitCollection="TrackerEndcapRawHits",
+
+trk_ec_coll = SimTrackerHitsCollector("trk_ec_coll",
+        inputSimTrackerHits = tracker_endcap_collections,
+        outputSimTrackerHits = "TrackerEndcapAllHits")
+algorithms.append( trk_ec_coll )
+trk_ec_digi = TrackerDigi("trk_ec_digi",
+        inputHitCollection = trk_ec_coll.outputSimTrackerHits,
+        outputHitCollection = "TrackerEndcapRawHits",
         timeResolution=8)
 algorithms.append( trk_ec_digi )
 
-vtx_b_digi = TrackerDigi("vtx_b_digi", 
-        inputHitCollection="VertexBarrelHits",
-        outputHitCollection="VertexBarrelRawHits",
+vtx_b_coll = SimTrackerHitsCollector("vtx_b_coll",
+        inputSimTrackerHits = vertex_barrel_collections,
+        outputSimTrackerHits = "VertexBarrelAllHits")
+algorithms.append( vtx_b_coll )
+vtx_b_digi = TrackerDigi("vtx_b_digi",
+        inputHitCollection = vtx_b_coll.outputSimTrackerHits,
+        outputHitCollection = "VertexBarrelRawHits",
         timeResolution=8)
 algorithms.append( vtx_b_digi )
 
 if 'acadia' in detector_version:
-    vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
-            inputHitCollection="VertexEndcapHits",
-            outputHitCollection="VertexEndcapRawHits",
+    vtx_ec_coll = SimTrackerHitsCollector("vtx_ec_coll",
+            inputSimTrackerHits = vertex_endcap_collections,
+            outputSimTrackerHits = "VertexEndcapAllHits")
+    algorithms.append( vtx_ec_coll )
+    vtx_ec_digi = TrackerDigi("vtx_ec_digi",
+            inputHitCollection = vtx_ec_coll.outputSimTrackerHits,
+            outputHitCollection = "VertexEndcapRawHits",
             timeResolution=8)
     algorithms.append( vtx_ec_digi )
 else:
-    mm_b_digi = TrackerDigi("mm_b_digi", 
-            inputHitCollection="MPGDTrackerBarrelHits",
-            outputHitCollection="MPGDTrackerBarrelRawHits",
+    mm_b_coll = SimTrackerHitsCollector("mm_b_coll",
+            inputSimTrackerHits = mpgd_barrel_collections,
+            outputSimTrackerHits = "MPGDTrackerBarrelAllHits")
+    algorithms.append( mm_b_coll )
+    mm_b_digi = TrackerDigi("mm_b_digi",
+            inputHitCollection = mm_b_coll.outputSimTrackerHits,
+            outputHitCollection = "MPGDTrackerBarrelRawHits",
             timeResolution=8)
     algorithms.append( mm_b_digi )
 
+gem_ec_coll = SimTrackerHitsCollector("gem_ec_coll",
+        inputSimTrackerHits = gem_endcap_collections,
+        outputSimTrackerHits = "GEMTrackerEndcapAllHits")
+algorithms.append( gem_ec_coll )
 gem_ec_digi = TrackerDigi("gem_ec_digi",
-        inputHitCollection="GEMTrackerEndcapHits",
-        outputHitCollection="GEMTrackerEndcapRawHits",
+        inputHitCollection = gem_ec_coll.outputSimTrackerHits,
+        outputHitCollection = "GEMTrackerEndcapRawHits",
         timeResolution=10)
-algorithms.append(gem_ec_digi)
+algorithms.append( gem_ec_digi )
 
 # Tracker and vertex reconstruction
 trk_b_reco = TrackerHitReconstruction("trk_b_reco",
@@ -149,7 +215,7 @@ algorithms.append( sourcelinker )
 
 ## Track param init
 truth_trk_init = TrackParamTruthInit("truth_trk_init",
-        inputMCParticles="mcparticles",
+        inputMCParticles="MCParticles",
         outputInitialTrackParameters="InitTrackParams",
         OutputLevel=DEBUG)
         #OutputLevel=DEBUG)
@@ -196,7 +262,7 @@ out.outputCommands = ["keep *",
         "drop trajectories",
         "drop outputSourceLinks",
         "drop outputInitialTrackParameters",
-        "keep mcparticles"
+        "keep MCParticles"
         ]
 algorithms.append(out)
 
diff --git a/benchmarks/track_fitting/scripts/gen_multiple_tracks.cxx b/benchmarks/track_fitting/scripts/gen_multiple_tracks.cxx
index f26840068ae6697d286e5f4c63f18e7d808798f7..a8251f09e11585047f42549aaf50402a7f320cd0 100644
--- a/benchmarks/track_fitting/scripts/gen_multiple_tracks.cxx
+++ b/benchmarks/track_fitting/scripts/gen_multiple_tracks.cxx
@@ -3,9 +3,11 @@
 #include "HepMC3/WriterAscii.h"
 #include "HepMC3/Print.h"
 
+#include "TRandom.h"
+
 #include <iostream>
-#include<random>
-#include<cmath>
+#include <random>
+#include <cmath>
 #include <math.h>
 #include <TMath.h>
 
diff --git a/benchmarks/track_fitting/scripts/gen_single_tracks.cxx b/benchmarks/track_fitting/scripts/gen_single_tracks.cxx
index 963810510115d37df8cff7e0f9c3e994a69858fc..de7c4c192c64e998bd185d6d99eacc07132d9174 100644
--- a/benchmarks/track_fitting/scripts/gen_single_tracks.cxx
+++ b/benchmarks/track_fitting/scripts/gen_single_tracks.cxx
@@ -3,6 +3,8 @@
 #include "HepMC3/WriterAscii.h"
 #include "HepMC3/Print.h"
 
+#include "TRandom.h"
+
 #include <iostream>
 #include <random>
 #include <cmath>
diff --git a/benchmarks/track_fitting/scripts/rec_multiple_tracks.cxx b/benchmarks/track_fitting/scripts/rec_multiple_tracks.cxx
index ada976a32d3f7af3c1f8762cb14bd729d8424c3e..99f0ded5554aa553b26ba7dc3eb68e6e16cdfc82 100644
--- a/benchmarks/track_fitting/scripts/rec_multiple_tracks.cxx
+++ b/benchmarks/track_fitting/scripts/rec_multiple_tracks.cxx
@@ -2,13 +2,15 @@
 #include "TCanvas.h"
 #include "TLegend.h"
 #include "TH1D.h"
+#include "THStack.h"
 #include "TProfile.h"
+#include "Math/Vector4D.h"
 
 #include <iostream>
 
 R__LOAD_LIBRARY(libeicd.so)
 R__LOAD_LIBRARY(libDD4pod.so)
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/TrackParametersCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
@@ -26,10 +28,10 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 };
 
 
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+std::vector<float> pt (std::vector<edm4hep::MCParticleData> 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));
+    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
   }
   return result;
 }
@@ -48,11 +50,11 @@ auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
   }
   return result;
 };
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto fourvec = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv);
   }
   return result;
@@ -84,8 +86,8 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root")
   ROOT::EnableImplicitMT();
   ROOT::RDataFrame df("events", fname);
 
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
+  auto df0 = df.Define("isThrown", "MCParticles.generatorStatus == 1")
+                 .Define("thrownParticles", "MCParticles[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
diff --git a/benchmarks/track_fitting/scripts/rec_single_tracks.cxx b/benchmarks/track_fitting/scripts/rec_single_tracks.cxx
index 9aad6734df65ecd7fc0dec19ab4484addfefe63a..4569358bd76c25fbf3280b149c7a25dd36039043 100644
--- a/benchmarks/track_fitting/scripts/rec_single_tracks.cxx
+++ b/benchmarks/track_fitting/scripts/rec_single_tracks.cxx
@@ -2,13 +2,15 @@
 #include "TCanvas.h"
 #include "TLegend.h"
 #include "TH1D.h"
+#include "THStack.h"
 #include "TProfile.h"
+#include "Math/Vector4D.h"
 
 #include <iostream>
 
 R__LOAD_LIBRARY(libeicd.so)
 R__LOAD_LIBRARY(libDD4pod.so)
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/TrackParametersCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
@@ -29,10 +31,10 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 };
 
 
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+std::vector<float> pt (std::vector<edm4hep::MCParticleData> 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));
+    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
   }
   return result;
 }
@@ -51,11 +53,11 @@ auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
   }
   return result;
 };
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto fourvec = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv);
   }
   return result;
@@ -87,8 +89,8 @@ int rec_single_tracks(const char* fname = "topside/rec_single_tracks.root")
   ROOT::EnableImplicitMT();
   ROOT::RDataFrame df("events", fname);
 
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
+  auto df0 = df.Define("isThrown", "MCParticles.generatorStatus == 1")
+                 .Define("thrownParticles", "MCParticles[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
diff --git a/benchmarks/track_fitting/single_tracks.sh b/benchmarks/track_fitting/single_tracks.sh
index 85a05d9740b48dec0c773e1b34c4c8de9cefecb0..19d166fa1661b4fec48031036d1e8a98e9d19791 100644
--- a/benchmarks/track_fitting/single_tracks.sh
+++ b/benchmarks/track_fitting/single_tracks.sh
@@ -51,14 +51,14 @@ export JUGGLER_N_EVENTS=$(expr ${JUGGLER_N_EVENTS} \* 1)
 export JUGGLER_FILE_NAME_TAG="single_tracks"
 export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
 
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
 export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
 
 echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
 echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
 
 ## generate the input events
-root -b -q "benchmarks/track_fitting/scripts/gen_single_tracks.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+root -b -q "benchmarks/track_fitting/scripts/gen_single_tracks.cxx+(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running script"
   exit 1
@@ -66,8 +66,9 @@ fi
 
 echo "Running geant4 simulation"
 ## run geant4 simulations
-npsim --runType batch \
+ddsim --runType batch \
   --part.minimalKineticEnergy 1000*GeV  \
+  --filter.tracker edep0 \
   -v WARNING \
   --numberOfEvents ${JUGGLER_N_EVENTS} \
   --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
@@ -95,7 +96,7 @@ rootls -t  ${JUGGLER_REC_FILE}
 
 mkdir -p results/track_fitting
 
-root -b -q "benchmarks/track_fitting/scripts/rec_single_tracks.cxx(\"${JUGGLER_REC_FILE}\")"
+root -b -q "benchmarks/track_fitting/scripts/rec_single_tracks.cxx+(\"${JUGGLER_REC_FILE}\")"
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running root script"
   exit 1
diff --git a/benchmarks/tracking/central_electrons.sh b/benchmarks/tracking/central_electrons.sh
index 1c0fe1bab7fec1a4613b41cd9e443f147cfd6b64..85065831e3fd30c4b549ccb05ac1fe2be2433ecf 100644
--- a/benchmarks/tracking/central_electrons.sh
+++ b/benchmarks/tracking/central_electrons.sh
@@ -60,7 +60,7 @@ export JUGGLER_N_EVENTS=$(expr ${JUGGLER_N_EVENTS} \* 1)
 export JUGGLER_FILE_NAME_TAG="central_electrons"
 export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
 
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
 export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
 
 echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
@@ -71,7 +71,7 @@ if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
 then
 
   ## generate the input events
-  root -b -q "benchmarks/tracking/scripts/gen_central_electrons.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+  root -b -q "benchmarks/tracking/scripts/gen_central_electrons.cxx+(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
   if [[ "$?" -ne "0" ]] ; then
     echo "ERROR running script"
     exit 1
@@ -79,8 +79,9 @@ then
 
   echo "Running geant4 simulation"
   ## run geant4 simulations
-  npsim --runType batch \
+  ddsim --runType batch \
     --part.minimalKineticEnergy 1000*GeV  \
+    --filter.tracker edep0 \
     -v WARNING \
     --numberOfEvents ${JUGGLER_N_EVENTS} \
     --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
@@ -106,17 +107,17 @@ fi
 
 mkdir -p results/tracking
 
-root -b -q "benchmarks/tracking/scripts/rec_central_electrons.cxx(\"${JUGGLER_REC_FILE}\")"
+root -b -q "benchmarks/tracking/scripts/rec_central_electrons.cxx+(\"${JUGGLER_REC_FILE}\")"
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running root script"
   exit 1
 fi
 
-#root -b -q "benchmarks/tracking/scripts/hits_central_electrons.cxx(\"${JUGGLER_SIM_FILE}\")"
-#if [[ "$?" -ne "0" ]] ; then
-#  echo "ERROR running root script"
-#  exit 1
-#fi
+root -b -q "benchmarks/tracking/scripts/hits_central_electrons.cxx+(\"${JUGGLER_REC_FILE}\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running root script"
+  exit 1
+fi
 
 root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
 if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then 
diff --git a/benchmarks/tracking/central_pions.sh b/benchmarks/tracking/central_pions.sh
index e006db0f7ae04d292658d0bddee699a1f1336619..1081f4fc9b87378eefb8c606755c90486dccca4f 100644
--- a/benchmarks/tracking/central_pions.sh
+++ b/benchmarks/tracking/central_pions.sh
@@ -59,7 +59,7 @@ fi
 export JUGGLER_FILE_NAME_TAG="central_pions"
 export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
 
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
 export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
 
 echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
@@ -70,7 +70,7 @@ if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
 then
 
   ## generate the input events
-  root -b -q "benchmarks/tracking/scripts/gen_central_pions.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+  root -b -q "benchmarks/tracking/scripts/gen_central_pions.cxx+(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
   if [[ "$?" -ne "0" ]] ; then
     echo "ERROR running script"
     exit 1
@@ -78,8 +78,9 @@ then
 
   echo "Running geant4 simulation"
   ## run geant4 simulations
-  npsim --runType batch \
+  ddsim --runType batch \
     --part.minimalKineticEnergy 1000*GeV  \
+    --filter.tracker edep0 \
     -v WARNING \
     --numberOfEvents ${JUGGLER_N_EVENTS} \
     --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
@@ -105,13 +106,13 @@ fi
 
 mkdir -p results/tracking
 
-root -b -q "benchmarks/tracking/scripts/rec_central_pions.cxx(\"${JUGGLER_REC_FILE}\")"
+root -b -q "benchmarks/tracking/scripts/rec_central_pions.cxx+(\"${JUGGLER_REC_FILE}\")"
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running root script"
   exit 1
 fi
 
-root -b -q "benchmarks/tracking/scripts/hits_central_pions.cxx(\"${JUGGLER_SIM_FILE}\")"
+root -b -q "benchmarks/tracking/scripts/hits_central_pions.cxx+(\"${JUGGLER_REC_FILE}\")"
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running root script"
   exit 1
diff --git a/benchmarks/tracking/config.yml b/benchmarks/tracking/config.yml
index 580f443937e0c947770178f455c043a69aeeb0c8..796a6f7e031d7779c7c6587bc61c94ae3d293cf7 100644
--- a/benchmarks/tracking/config.yml
+++ b/benchmarks/tracking/config.yml
@@ -1,36 +1,45 @@
+tracking_compile:
+  extends: .compile_benchmark
+  stage: compile
+  script:
+    - compile_analyses.py --dir scripts tracking
+
 tracking_central_electrons:
   extends: .rec_benchmark
   stage: run
-  timeout: 24 hours
+  needs:
+    - ["tracking_compile"]
   script:
     - bash benchmarks/tracking/central_electrons.sh
-      #allow_failure: true
       
 multiple_tracks:
   extends: .rec_benchmark
   stage: run
-  timeout: 24 hours
+  needs:
+    - ["tracking_compile"]
   script:
     - bash benchmarks/tracking/multiple_tracks.sh
 
 tracking_central_pions:
   extends: .rec_benchmark
   stage: run
-  timeout: 24 hours
+  needs:
+    - ["tracking_compile"]
   script:
     - bash benchmarks/tracking/central_pions.sh
-      #allow_failure: true
 
 tracking_truth_init_electrons:
   extends: .rec_benchmark
   stage: run
-  timeout: 24 hours
+  needs:
+    - ["tracking_compile"]
   script:
     - python benchmarks/tracking/run_tracking_benchmarks.py --nametag=truth_electron --particle=electron --etamin=-4 --etamax=4 -n 150
 
 tracking_truth_init_pions:
   extends: .rec_benchmark
   stage: run
-  timeout: 24 hours
+  needs:
+    - ["tracking_compile"]
   script:
     - python benchmarks/tracking/run_tracking_benchmarks.py --nametag=truth_pion --particle=pion+ --etamin=-4 --etamax=4 -n 150
diff --git a/benchmarks/tracking/multiple_tracks.sh b/benchmarks/tracking/multiple_tracks.sh
index 55cfa470160019bfe7ec466a0d269aeede81b433..b9aea549ccbd1e855236a0ee7ea8606d64d0b7d0 100644
--- a/benchmarks/tracking/multiple_tracks.sh
+++ b/benchmarks/tracking/multiple_tracks.sh
@@ -60,7 +60,7 @@ export JUGGLER_N_EVENTS=$(expr ${JUGGLER_N_EVENTS} \* 1)
 export JUGGLER_FILE_NAME_TAG="multiple_tracks"
 export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
 
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
 export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
 
 echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
@@ -71,7 +71,7 @@ if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
 then
 
   ## generate the input events
-  root -b -q "benchmarks/tracking/scripts/gen_multiple_tracks.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+  root -b -q "benchmarks/tracking/scripts/gen_multiple_tracks.cxx+(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
   if [[ "$?" -ne "0" ]] ; then
     echo "ERROR running script"
     exit 1
@@ -79,8 +79,9 @@ then
 
   echo "Running geant4 simulation"
   ## run geant4 simulations
-  npsim --runType batch \
+  ddsim --runType batch \
     --part.minimalKineticEnergy 1000*GeV  \
+    --filter.tracker edep0 \
     -v WARNING \
     --numberOfEvents ${JUGGLER_N_EVENTS} \
     --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
@@ -106,13 +107,13 @@ fi
 
 mkdir -p results/tracking
 
-root -b -q "benchmarks/tracking/scripts/rec_multiple_tracks.cxx(\"${JUGGLER_REC_FILE}\")"
+root -b -q "benchmarks/tracking/scripts/rec_multiple_tracks.cxx+(\"${JUGGLER_REC_FILE}\")"
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running root script"
   exit 1
 fi
 
-#root -b -q "benchmarks/tracking/scripts/hits_central_electrons.cxx(\"${JUGGLER_SIM_FILE}\")"
+#root -b -q "benchmarks/tracking/scripts/hits_central_electrons.cxx+(\"${JUGGLER_SIM_FILE}\")"
 #if [[ "$?" -ne "0" ]] ; then
 #  echo "ERROR running root script"
 #  exit 1
diff --git a/benchmarks/tracking/options/track_reconstruction.py b/benchmarks/tracking/options/track_reconstruction.py
index a79ccdf70e94363f0d7c179da450d97b3c0b62b8..729021d2bf8086042d980efd1ff1709a0f372d44 100644
--- a/benchmarks/tracking/options/track_reconstruction.py
+++ b/benchmarks/tracking/options/track_reconstruction.py
@@ -39,6 +39,7 @@ podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=W
 
 from Configurables import PodioInput
 
+from Configurables import Jug__Digi__SimTrackerHitsCollector as SimTrackerHitsCollector
 from Configurables import Jug__Digi__SiliconTrackerDigi as TrackerDigi
 
 from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
@@ -58,50 +59,104 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 
 algorithms = [ ]
 
-input_collections = ['mcparticles','TrackerEndcapHits','TrackerBarrelHits','VertexBarrelHits','GEMTrackerEndcapHits']
+tracker_endcap_collections = [
+    'TrackerEndcapHits1',
+    'TrackerEndcapHits2',
+    'TrackerEndcapHits3',
+    'TrackerEndcapHits4',
+    'TrackerEndcapHits5',
+    'TrackerEndcapHits6'
+]
+tracker_barrel_collections = [
+    'TrackerBarrelHits'
+]
+vertex_barrel_collections = [
+    'VertexBarrelHits'
+]
+gem_endcap_collections = [
+    'GEMTrackerEndcapHits1',
+    'GEMTrackerEndcapHits2',
+    'GEMTrackerEndcapHits3'
+]
+vertex_endcap_collections = [
+    'VertexEndcapHits'
+]
+mpgd_barrel_collections = [
+    'MPGDTrackerBarrelHits1',
+    'MPGDTrackerBarrelHits2'
+]
+
+input_collections = ['MCParticles'] + tracker_endcap_collections + tracker_barrel_collections + vertex_barrel_collections + gem_endcap_collections
+
 if 'acadia' in detector_version:
-    input_collections.append('VertexEndcapHits')
+    input_collections += vertex_endcap_collections
 else:
-    input_collections.append('MPGDTrackerBarrelHits')
-podioinput = PodioInput("PodioReader", 
-                        collections=input_collections)#, OutputLevel=DEBUG)
+    input_collections += mpgd_barrel_collections
+
+podioinput = PodioInput("PodioReader",
+                        collections=input_collections)
 algorithms.append( podioinput )
 
-trk_b_digi = TrackerDigi("trk_b_digi", 
-        inputHitCollection="TrackerBarrelHits",
-        outputHitCollection="TrackerBarrelRawHits",
+trk_b_coll = SimTrackerHitsCollector("trk_b_coll",
+        inputSimTrackerHits = tracker_barrel_collections,
+        outputSimTrackerHits = "TrackerBarrelAllHits")
+algorithms.append( trk_b_coll )
+trk_b_digi = TrackerDigi("trk_b_digi",
+        inputHitCollection = trk_b_coll.outputSimTrackerHits,
+        outputHitCollection = "TrackerBarrelRawHits",
         timeResolution=8)
 algorithms.append( trk_b_digi )
-trk_ec_digi = TrackerDigi("trk_ec_digi", 
-        inputHitCollection="TrackerEndcapHits",
-        outputHitCollection="TrackerEndcapRawHits",
+
+trk_ec_coll = SimTrackerHitsCollector("trk_ec_coll",
+        inputSimTrackerHits = tracker_endcap_collections,
+        outputSimTrackerHits = "TrackerEndcapAllHits")
+algorithms.append( trk_ec_coll )
+trk_ec_digi = TrackerDigi("trk_ec_digi",
+        inputHitCollection = trk_ec_coll.outputSimTrackerHits,
+        outputHitCollection = "TrackerEndcapRawHits",
         timeResolution=8)
 algorithms.append( trk_ec_digi )
 
-vtx_b_digi = TrackerDigi("vtx_b_digi", 
-        inputHitCollection="VertexBarrelHits",
-        outputHitCollection="VertexBarrelRawHits",
+vtx_b_coll = SimTrackerHitsCollector("vtx_b_coll",
+        inputSimTrackerHits = vertex_barrel_collections,
+        outputSimTrackerHits = "VertexBarrelAllHits")
+algorithms.append( vtx_b_coll )
+vtx_b_digi = TrackerDigi("vtx_b_digi",
+        inputHitCollection = vtx_b_coll.outputSimTrackerHits,
+        outputHitCollection = "VertexBarrelRawHits",
         timeResolution=8)
 algorithms.append( vtx_b_digi )
 
 if 'acadia' in detector_version:
-    vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
-            inputHitCollection="VertexEndcapHits",
-            outputHitCollection="VertexEndcapRawHits",
+    vtx_ec_coll = SimTrackerHitsCollector("vtx_ec_coll",
+            inputSimTrackerHits = vertex_endcap_collections,
+            outputSimTrackerHits = "VertexEndcapAllHits")
+    algorithms.append( vtx_ec_coll )
+    vtx_ec_digi = TrackerDigi("vtx_ec_digi",
+            inputHitCollection = vtx_ec_coll.outputSimTrackerHits,
+            outputHitCollection = "VertexEndcapRawHits",
             timeResolution=8)
     algorithms.append( vtx_ec_digi )
 else:
-    mm_b_digi = TrackerDigi("mm_b_digi", 
-            inputHitCollection="MPGDTrackerBarrelHits",
-            outputHitCollection="MPGDTrackerBarrelRawHits",
+    mm_b_coll = SimTrackerHitsCollector("mm_b_coll",
+            inputSimTrackerHits = mpgd_barrel_collections,
+            outputSimTrackerHits = "MPGDTrackerBarrelAllHits")
+    algorithms.append( mm_b_coll )
+    mm_b_digi = TrackerDigi("mm_b_digi",
+            inputHitCollection = mm_b_coll.outputSimTrackerHits,
+            outputHitCollection = "MPGDTrackerBarrelRawHits",
             timeResolution=8)
     algorithms.append( mm_b_digi )
 
+gem_ec_coll = SimTrackerHitsCollector("gem_ec_coll",
+        inputSimTrackerHits = gem_endcap_collections,
+        outputSimTrackerHits = "GEMTrackerEndcapAllHits")
+algorithms.append( gem_ec_coll )
 gem_ec_digi = TrackerDigi("gem_ec_digi",
-        inputHitCollection="GEMTrackerEndcapHits",
-        outputHitCollection="GEMTrackerEndcapRawHits",
+        inputHitCollection = gem_ec_coll.outputSimTrackerHits,
+        outputHitCollection = "GEMTrackerEndcapRawHits",
         timeResolution=10)
-algorithms.append(gem_ec_digi)
+algorithms.append( gem_ec_digi )
 
 # Tracker and vertex reconstruction
 trk_b_reco = TrackerHitReconstruction("trk_b_reco",
@@ -154,13 +209,12 @@ algorithms.append( trk_hit_col )
 sourcelinker = TrackerSourceLinker("sourcelinker",
         inputHitCollection=trk_hit_col.trackingHits,
         outputSourceLinks="TrackSourceLinks",
-        outputMeasurements="TrackMeasurements",
-        OutputLevel=VERBOSE)
+        outputMeasurements="TrackMeasurements")
 algorithms.append( sourcelinker )
 
 ## Track param init
 truth_trk_init = TrackParamTruthInit("truth_trk_init",
-        inputMCParticles="mcparticles",
+        inputMCParticles="MCParticles",
         outputInitialTrackParameters="InitTrackParams")
 algorithms.append( truth_trk_init )
 
@@ -176,13 +230,11 @@ parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
         inputTrajectories="trajectories",
         outputParticles="ReconstructedParticles",
         outputTrackParameters="outputTrackParameters")
-        #OutputLevel=DEBUG)
 algorithms.append( parts_from_fit )
 
 #trajs_from_fit = TrajectoryFromTrackFit("trajs_from_fit",
 #inputTrajectories = trk_find_alg.outputTrajectories,
 #outputTrajectoryParameters = "outputTrajectoryParameters")
-        #OutputLevel=DEBUG)
 #algorithms.append(trajs_from_fit)
 
 out = PodioOutput("out", filename=output_rec_file)
@@ -192,7 +244,7 @@ out.outputCommands = ["keep *",
         "drop trajectories",
         "drop outputSourceLinks",
         "drop outputInitialTrackParameters",
-        "keep mcparticles"
+        "keep MCParticles"
         ]
 algorithms.append(out)
 
diff --git a/benchmarks/tracking/run_tracking_benchmarks.py b/benchmarks/tracking/run_tracking_benchmarks.py
index 32a8a149f1b8ba7e3dabfb4741d0bcec1aa6059c..d4660c5e205650280133b3a9818042e81e184e49 100755
--- a/benchmarks/tracking/run_tracking_benchmarks.py
+++ b/benchmarks/tracking/run_tracking_benchmarks.py
@@ -40,7 +40,7 @@ option_script = os.path.join(script_dir, 'options', args.options)
 analysis_script = os.path.join(script_dir, 'scripts', 'tracking_performance.py')
 
 gen_file = 'gen_{nametag}_{pmin}_{pmax}.hepmc'.format(**kwargs)
-sim_file = 'sim_{nametag}_{pmin}_{pmax}.root'.format(**kwargs)
+sim_file = 'sim_{nametag}_{pmin}_{pmax}.edm4hep.root'.format(**kwargs)
 rec_file = 'rec_{nametag}_{pmin}_{pmax}.root'.format(**kwargs)
 
 procs = [p.strip() for p in args.process.split(',')]
@@ -56,8 +56,9 @@ if 'sim' in procs:
             '--particles', args.particles]
     subprocess.run(gen_cmd)
     # simulation
-    sim_cmd = ['npsim',
+    sim_cmd = ['ddsim',
             '--part.minimalKineticEnergy', '1*TeV',
+            '--filter.tracker', 'edep0',
             '--numberOfEvents', '{}'.format(args.nev),
             '--runType', 'batch',
             '--inputFiles', gen_file,
@@ -91,7 +92,7 @@ if 'rec' in procs:
 if 'ana' in procs:
     os.makedirs('results', exist_ok=True)
     ana_cmd = ['python', analysis_script, rec_file,
-               '--mc-collection', 'mcparticles',
+               '--mc-collection', 'MCParticles',
                '--tracking-collection', 'outputTrackParameters',
                '-o', 'results', '-t', args.nametag]
     return_code = subprocess.run(ana_cmd).returncode
diff --git a/benchmarks/tracking/scripts/gen_central_electrons.cxx b/benchmarks/tracking/scripts/gen_central_electrons.cxx
index 3e9d8c01034fe360ba166d707d7d209f9a40fec2..455c799dd598cc50ab70ba4be2d3ec5a54ae47a5 100644
--- a/benchmarks/tracking/scripts/gen_central_electrons.cxx
+++ b/benchmarks/tracking/scripts/gen_central_electrons.cxx
@@ -3,9 +3,11 @@
 #include "HepMC3/WriterAscii.h"
 #include "HepMC3/Print.h"
 
+#include "TRandom.h"
+
 #include <iostream>
-#include<random>
-#include<cmath>
+#include <random>
+#include <cmath>
 #include <math.h>
 #include <TMath.h>
 
diff --git a/benchmarks/tracking/scripts/gen_central_pions.cxx b/benchmarks/tracking/scripts/gen_central_pions.cxx
index abf691b5a984ad23758d94d900b9fa74e438b874..0e314c539844a700cb2a8041f7250a403a5c8a33 100644
--- a/benchmarks/tracking/scripts/gen_central_pions.cxx
+++ b/benchmarks/tracking/scripts/gen_central_pions.cxx
@@ -3,9 +3,11 @@
 #include "HepMC3/WriterAscii.h"
 #include "HepMC3/Print.h"
 
+#include "TRandom.h"
+
 #include <iostream>
-#include<random>
-#include<cmath>
+#include <random>
+#include <cmath>
 #include <math.h>
 #include <TMath.h>
 
diff --git a/benchmarks/tracking/scripts/gen_multiple_tracks.cxx b/benchmarks/tracking/scripts/gen_multiple_tracks.cxx
index f26840068ae6697d286e5f4c63f18e7d808798f7..a8251f09e11585047f42549aaf50402a7f320cd0 100644
--- a/benchmarks/tracking/scripts/gen_multiple_tracks.cxx
+++ b/benchmarks/tracking/scripts/gen_multiple_tracks.cxx
@@ -3,9 +3,11 @@
 #include "HepMC3/WriterAscii.h"
 #include "HepMC3/Print.h"
 
+#include "TRandom.h"
+
 #include <iostream>
-#include<random>
-#include<cmath>
+#include <random>
+#include <cmath>
 #include <math.h>
 #include <TMath.h>
 
diff --git a/benchmarks/tracking/scripts/hits_central_electrons.cxx b/benchmarks/tracking/scripts/hits_central_electrons.cxx
index 7e38cc88aea38933b9b0567bfc915198e72e779e..b36c350ccbc0c669f997c7ee045b3a3ddecba00f 100644
--- a/benchmarks/tracking/scripts/hits_central_electrons.cxx
+++ b/benchmarks/tracking/scripts/hits_central_electrons.cxx
@@ -2,15 +2,17 @@
 #include "TCanvas.h"
 #include "TLegend.h"
 #include "TH1D.h"
+#include "THStack.h"
 #include "TProfile.h"
+#include "Math/Vector4D.h"
 
 #include <iostream>
 
 R__LOAD_LIBRARY(libeicd.so)
 R__LOAD_LIBRARY(libDD4pod.so)
-#include "dd4pod/TrackerHitCollection.h"
-#include "dd4pod/TrackerHitData.h"
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/SimTrackerHitData.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/TrackParametersCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
@@ -28,10 +30,10 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 };
 
 
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+std::vector<float> pt (std::vector<edm4hep::MCParticleData> 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));
+    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
   }
   return result;
 }
@@ -50,11 +52,11 @@ auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
   }
   return result;
 };
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto fourvec = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv);
   }
   return result;
@@ -80,14 +82,14 @@ auto delta_p_over_p = [](const std::vector<double>& tracks, const std::vector<do
   return res;
 };
 
-int hits_central_electrons(const char* fname = "sim_central_electrons.root")
+int hits_central_electrons(const char* fname = "sim_central_electrons.edm4hep.root")
 {
 
   ROOT::EnableImplicitMT();
   ROOT::RDataFrame df("events", fname);
 
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
+  auto df0 = df.Define("isThrown", "MCParticles.generatorStatus == 1")
+                 .Define("thrownParticles", "MCParticles[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
@@ -103,21 +105,20 @@ int hits_central_electrons(const char* fname = "sim_central_electrons.root")
                  //.Define("delta_p_over_p1",delta_p_over_p, {"p_track1", "p_thrown"})
                  //.Define("delta_p_over_p2",delta_p_over_p, {"p_track2", "p_thrown"})
                  //.Define("N_VtxBarrelHits",[](std::vector<eic::TrackerHitData> hits) { return hits.size();},{"VertexBarrelRecHits"})
-                 //.Define("N_Hits",       [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"trackingHits"})
-                 .Define("N_BarrelHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelHits"})
-                 .Define("N_EndcapHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapHits"})
+                 .Define("N_BarrelHits", [](std::vector<edm4hep::SimTrackerHitData> hits) { return hits.size();}, {"TrackerBarrelAllHits"})
+                 .Define("N_EndcapHits", [](std::vector<edm4hep::SimTrackerHitData> hits) { return hits.size();}, {"TrackerEndcapAllHits"})
                  ;
 
-  auto hBarrel_x_vs_y = df0.Histo2D({"hBarrel_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "TrackerBarrelHits.position.x", "TrackerBarrelHits.position.y");
-  auto hEndcap_x_vs_y = df0.Histo2D({"hEndcap_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "TrackerEndcaplHits.position.x", "TrackerEndcapHits.position.y");
-  //auto hvtxBarrel_x_vs_y = df0.Histo2D({"hvtxBarrel_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "VertexBarrelHits.position.x", "VertexBarrelHits.position.y");
-  //auto hvtxEndcap_x_vs_y = df0.Histo2D({"hvtxEndcap_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "VertexEndcaplHits.position.x", "VertexEndcapHits.position.y");
+  auto hBarrel_x_vs_y = df0.Histo2D({"hBarrel_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "TrackerBarrelAllHits.position.x", "TrackerBarrelAllHits.position.y");
+  auto hEndcap_x_vs_y = df0.Histo2D({"hEndcap_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "TrackerEndcapAllHits.position.x", "TrackerEndcapAllHits.position.y");
+  //auto hvtxBarrel_x_vs_y = df0.Histo2D({"hvtxBarrel_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "VertexBarrelAllHits.position.x", "VertexBarrelAllHits.position.y");
+  //auto hvtxEndcap_x_vs_y = df0.Histo2D({"hvtxEndcap_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "VertexEndcaplAllHits.position.x", "VertexEndcapAllHits.position.y");
   //auto hAllHits_x_vs_y = df0.Histo2D({"hAllHitsx_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "allHits.position.x", "allHits.position.y");
 
-  auto hBarrel_x_vs_z = df0.Histo2D({"hBarrel_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "TrackerBarrelHits.position.z" , "TrackerBarrelHits.position.y");
-  auto hEndcap_x_vs_z = df0.Histo2D({"hEndcap_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "TrackerEndcaplHits.position.z", "TrackerEndcapHits.position.y");
-  //auto hvtxBarrel_x_vs_z = df0.Histo2D({"hvtxBarrel_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "VertexBarrelHits.position.z"  , "VertexBarrelHits.position.y" );
-  //auto hvtxEndcap_x_vs_z = df0.Histo2D({"hvtxEndcap_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "VertexEndcaplHits.position.z" , "VertexEndcapHits.position.y" );
+  auto hBarrel_x_vs_z = df0.Histo2D({"hBarrel_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "TrackerBarrelAllHits.position.z" , "TrackerBarrelAllHits.position.y");
+  auto hEndcap_x_vs_z = df0.Histo2D({"hEndcap_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "TrackerEndcapAllHits.position.z", "TrackerEndcapAllHits.position.y");
+  //auto hvtxBarrel_x_vs_z = df0.Histo2D({"hvtxBarrel_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "VertexBarrelAllHits.position.z"  , "VertexBarrelAllHits.position.y" );
+  //auto hvtxEndcap_x_vs_z = df0.Histo2D({"hvtxEndcap_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "VertexEndcaplAllHits.position.z" , "VertexEndcapAllHits.position.y" );
   //auto hAllHits_x_vs_z = df0.Histo2D({"hAllHitsx_vs_z","; z ; x ", 100, -900, 900,100, -900, 900 }, "allHits.position.z"           , "allHits.position.y"          );
 
   auto hBarrel_N_vs_theta = df0.Histo1D({"hBarrel_N_vs_theta", "; #theta [deg.]",   20, 0, 180 }, "theta0", "N_BarrelHits");
diff --git a/benchmarks/tracking/scripts/hits_central_pions.cxx b/benchmarks/tracking/scripts/hits_central_pions.cxx
index 0a70c5e487fa07404f524261ae9a4fe80f4d8b29..6263e1c167b15c8e446bec30a7e74c93e95caf55 100644
--- a/benchmarks/tracking/scripts/hits_central_pions.cxx
+++ b/benchmarks/tracking/scripts/hits_central_pions.cxx
@@ -2,15 +2,17 @@
 #include "TCanvas.h"
 #include "TLegend.h"
 #include "TH1D.h"
+#include "THStack.h"
 #include "TProfile.h"
+#include "Math/Vector4D.h"
 
 #include <iostream>
 
 R__LOAD_LIBRARY(libeicd.so)
 R__LOAD_LIBRARY(libDD4pod.so)
-#include "dd4pod/TrackerHitCollection.h"
-#include "dd4pod/TrackerHitData.h"
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/SimTrackerHitData.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/TrackParametersCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
@@ -28,10 +30,10 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 };
 
 
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+std::vector<float> pt (std::vector<edm4hep::MCParticleData> 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));
+    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
   }
   return result;
 }
@@ -50,11 +52,11 @@ auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
   }
   return result;
 };
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto fourvec = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv);
   }
   return result;
@@ -80,14 +82,14 @@ auto delta_p_over_p = [](const std::vector<double>& tracks, const std::vector<do
   return res;
 };
 
-int hits_central_pions(const char* fname = "sim_central_pions.root")
+int hits_central_pions(const char* fname = "sim_central_pions.edm4hep.root")
 {
 
   ROOT::EnableImplicitMT();
   ROOT::RDataFrame df("events", fname);
 
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
+  auto df0 = df.Define("isThrown", "MCParticles.generatorStatus == 1")
+                 .Define("thrownParticles", "MCParticles[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
@@ -103,11 +105,11 @@ int hits_central_pions(const char* fname = "sim_central_pions.root")
                  //.Define("delta_p_over_p1",delta_p_over_p, {"p_track1", "p_thrown"})
                  //.Define("delta_p_over_p2",delta_p_over_p, {"p_track2", "p_thrown"})
                  //.Define("N_VtxBarrelHits",[](std::vector<eic::TrackerHitData> hits) { return hits.size();},{"VertexBarrelRecHits"})
-                 .Define("N_BarrelHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelHits"})
-                 .Define("N_EndcapHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapHits"})
+                 .Define("N_BarrelHits", [](std::vector<edm4hep::SimTrackerHitData> hits) { return hits.size();}, {"TrackerBarrelAllHits"})
+                 .Define("N_EndcapHits", [](std::vector<edm4hep::SimTrackerHitData> hits) { return hits.size();}, {"TrackerEndcapAllHits"})
                  ;
 
-  auto hBarrel_x_vs_y = df0.Histo2D({"hBarrel_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "TrackerBarrelHits.position.x", "TrackerBarrelHits.position.y");
+  auto hBarrel_x_vs_y = df0.Histo2D({"hBarrel_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "TrackerBarrelAllHits.position.x", "TrackerBarrelAllHits.position.y");
 
   auto hBarrel_N_vs_theta = df0.Histo1D({"hBarrel_N_vs_theta", "; #theta [deg.]",   20, 0, 180 }, "theta0", "N_BarrelHits");
   auto hEndcap_N_vs_theta = df0.Histo1D({"hEndcap_N_vs_theta", "; #theta [deg.]",   20, 0, 180 }, "theta0", "N_EndcapHits");
diff --git a/benchmarks/tracking/scripts/rec_central_electrons.cxx b/benchmarks/tracking/scripts/rec_central_electrons.cxx
index 8c3708f1c52d9d08b65dd53e52188176f94512a1..ad241ebbb72ab670b0902ca717261056493db908 100644
--- a/benchmarks/tracking/scripts/rec_central_electrons.cxx
+++ b/benchmarks/tracking/scripts/rec_central_electrons.cxx
@@ -2,13 +2,15 @@
 #include "TCanvas.h"
 #include "TLegend.h"
 #include "TH1D.h"
+#include "THStack.h"
 #include "TProfile.h"
+#include "Math/Vector4D.h"
 
 #include <iostream>
 
 R__LOAD_LIBRARY(libeicd.so)
 R__LOAD_LIBRARY(libDD4pod.so)
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/TrackParametersCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
@@ -26,10 +28,10 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 };
 
 
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+std::vector<float> pt (std::vector<edm4hep::MCParticleData> 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));
+    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
   }
   return result;
 }
@@ -48,11 +50,11 @@ auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
   }
   return result;
 };
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto fourvec = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv);
   }
   return result;
@@ -84,8 +86,8 @@ int rec_central_electrons(const char* fname = "topside/rec_central_electrons.roo
   ROOT::EnableImplicitMT();
   ROOT::RDataFrame df("events", fname);
 
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
+  auto df0 = df.Define("isThrown", "MCParticles.generatorStatus == 1")
+                 .Define("thrownParticles", "MCParticles[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
diff --git a/benchmarks/tracking/scripts/rec_central_pions.cxx b/benchmarks/tracking/scripts/rec_central_pions.cxx
index 516ecb2b315252dd495878472f0d8b58f9fab00a..588a5a02f2b9b529727c759da7bffe3a0b1f03d1 100644
--- a/benchmarks/tracking/scripts/rec_central_pions.cxx
+++ b/benchmarks/tracking/scripts/rec_central_pions.cxx
@@ -2,13 +2,15 @@
 #include "TCanvas.h"
 #include "TLegend.h"
 #include "TH1D.h"
+#include "THStack.h"
 #include "TProfile.h"
+#include "Math/Vector4D.h"
 
 #include <iostream>
 
 R__LOAD_LIBRARY(libeicd.so)
 R__LOAD_LIBRARY(libDD4pod.so)
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/TrackParametersCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
@@ -26,10 +28,10 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 };
 
 
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+std::vector<float> pt (std::vector<edm4hep::MCParticleData> 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));
+    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
   }
   return result;
 }
@@ -48,11 +50,11 @@ auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
   }
   return result;
 };
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto fourvec = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv);
   }
   return result;
@@ -84,8 +86,8 @@ int rec_central_pions(const char* fname = "topside/rec_central_pions.root")
   ROOT::EnableImplicitMT();
   ROOT::RDataFrame df("events", fname);
 
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
+  auto df0 = df.Define("isThrown", "MCParticles.generatorStatus == 1")
+                 .Define("thrownParticles", "MCParticles[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
diff --git a/benchmarks/tracking/scripts/rec_multiple_tracks.cxx b/benchmarks/tracking/scripts/rec_multiple_tracks.cxx
index ada976a32d3f7af3c1f8762cb14bd729d8424c3e..99f0ded5554aa553b26ba7dc3eb68e6e16cdfc82 100644
--- a/benchmarks/tracking/scripts/rec_multiple_tracks.cxx
+++ b/benchmarks/tracking/scripts/rec_multiple_tracks.cxx
@@ -2,13 +2,15 @@
 #include "TCanvas.h"
 #include "TLegend.h"
 #include "TH1D.h"
+#include "THStack.h"
 #include "TProfile.h"
+#include "Math/Vector4D.h"
 
 #include <iostream>
 
 R__LOAD_LIBRARY(libeicd.so)
 R__LOAD_LIBRARY(libDD4pod.so)
-#include "dd4pod/Geant4ParticleCollection.h"
+#include "edm4hep/MCParticleCollection.h"
 #include "eicd/TrackParametersCollection.h"
 #include "eicd/ClusterCollection.h"
 #include "eicd/ClusterData.h"
@@ -26,10 +28,10 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 };
 
 
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+std::vector<float> pt (std::vector<edm4hep::MCParticleData> 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));
+    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
   }
   return result;
 }
@@ -48,11 +50,11 @@ auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
   }
   return result;
 };
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+auto fourvec = [](ROOT::VecOps::RVec<edm4hep::MCParticleData> 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);
+    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
     result.push_back(lv);
   }
   return result;
@@ -84,8 +86,8 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root")
   ROOT::EnableImplicitMT();
   ROOT::RDataFrame df("events", fname);
 
-  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
-                 .Define("thrownParticles", "mcparticles[isThrown]")
+  auto df0 = df.Define("isThrown", "MCParticles.generatorStatus == 1")
+                 .Define("thrownParticles", "MCParticles[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
diff --git a/benchmarks/tracking/scripts/tracking_performance.py b/benchmarks/tracking/scripts/tracking_performance.py
index 813f4d93b17909ca0c66925565a9c0e567f4c029..c39c97171e57ba22f4b32a81ee1a4095c987e0c3 100644
--- a/benchmarks/tracking/scripts/tracking_performance.py
+++ b/benchmarks/tracking/scripts/tracking_performance.py
@@ -47,12 +47,12 @@ def flatten_collection(rdf, collection, cols=None):
     dfp.loc[:, 'event'] = evns
     return dfp
 
-def thrown_particles_figure(rdf, save, mcbranch="mcparticles"):
+def thrown_particles_figure(rdf, save, mcbranch="MCParticles"):
     # define truth particle info
-    dft = flatten_collection(rdf, mcbranch, ['genStatus', 'pdgID', 'ps.x', 'ps.y', 'ps.z', 'mass'])
+    dft = flatten_collection(rdf, mcbranch, ['generatorStatus', 'PDG', 'momentum.x', 'momentum.y', 'momentum.z', 'mass'])
     dft.rename(columns={c: c.replace(mcbranch + '.', '') for c in dft.columns}, inplace=True)
     # select thrown particles
-    dft = dft[dft['genStatus'] == 1]
+    dft = dft[dft['generatorStatus'] == 1]
 
     # figure
     fig, axs = plt.subplots(2, 2, figsize=(16, 12), dpi=120)
@@ -62,8 +62,8 @@ def thrown_particles_figure(rdf, save, mcbranch="mcparticles"):
     get_pcharge = np.vectorize(lambda pid: pdgbase.GetParticle(int(pid)).Charge()/3.)
 
     # enumerate particle names
-    dft.loc[:, 'pname'] = get_pname(dft['pdgID'].values)
-    dft.loc[:, 'charge'] = get_pcharge(dft['pdgID'].values)
+    dft.loc[:, 'pname'] = get_pname(dft['PDG'].values)
+    dft.loc[:, 'charge'] = get_pcharge(dft['PDG'].values)
     penum = {pname: i for i, pname in enumerate(dft['pname'].unique())}
     dft.loc[:, 'pname_id'] = dft['pname'].map(penum)
 
@@ -78,7 +78,7 @@ def thrown_particles_figure(rdf, save, mcbranch="mcparticles"):
 
     # calculate kinematics
     get_4vecs = np.vectorize(lambda x, y, z, m: ROOT.Math.PxPyPzMVector(x, y, z, m))
-    fourvecs = get_4vecs(*dft[['ps.x', 'ps.y', 'ps.z', 'mass']].values.T)
+    fourvecs = get_4vecs(*dft[['momentum.x', 'momentum.y', 'momentum.z', 'mass']].values.T)
 
     dft.loc[:, 'p'] = [v.P() for v in fourvecs]
     dft.loc[:, 'theta'] = [v.Theta() for v in fourvecs]
@@ -112,7 +112,7 @@ if __name__ == '__main__':
     parser.add_argument('-t', '--nametag', type=str, default='tracking', help='Name tag for output files.')
     parser.add_argument('-m', '--macros', dest='macros', default='rootlogon.C',
             help='Macro files to be loaded by root, separated by \",\".')
-    parser.add_argument('--mc-collection', dest='mc', default='mcparticles',
+    parser.add_argument('--mc-collection', dest='mc', default='MCParticles',
             help='Collection name of MC particles truth info.')
     parser.add_argument('--tracking-collection', dest='coll', default='outputTrackParameters',
             help='Collection name of clusters to plot')