From 2d110c855b4897553101358292fc8a98cdd53fb5 Mon Sep 17 00:00:00 2001
From: Zhenyu Ye <yezhenyu@10-8-239-141.vlan877.client.wireless.uic.edu>
Date: Thu, 4 Nov 2021 18:12:17 -0500
Subject: [PATCH] add BTOF benchmark

---
 benchmarks/tof/central_electrons.sh           | 130 +++++++++
 benchmarks/tof/config.yml                     |  14 +
 benchmarks/tof/multiple_tracks.sh             | 107 ++++++++
 .../central_electron_reconstruction.py        | 239 +++++++++++++++++
 .../central_electron_reconstruction.py.org    | 193 ++++++++++++++
 .../tof/options/track_reconstruction.py       | 239 +++++++++++++++++
 .../tof/scripts/gen_central_electrons.cxx     |  83 ++++++
 .../tof/scripts/gen_multiple_tracks.cxx       |  95 +++++++
 .../tof/scripts/rec_central_electrons.cxx     | 249 ++++++++++++++++++
 .../tof/scripts/rec_multiple_tracks.cxx       | 248 +++++++++++++++++
 10 files changed, 1597 insertions(+)
 create mode 100644 benchmarks/tof/central_electrons.sh
 create mode 100644 benchmarks/tof/config.yml
 create mode 100644 benchmarks/tof/multiple_tracks.sh
 create mode 100644 benchmarks/tof/options/central_electron_reconstruction.py
 create mode 100644 benchmarks/tof/options/central_electron_reconstruction.py.org
 create mode 100644 benchmarks/tof/options/track_reconstruction.py
 create mode 100644 benchmarks/tof/scripts/gen_central_electrons.cxx
 create mode 100644 benchmarks/tof/scripts/gen_multiple_tracks.cxx
 create mode 100644 benchmarks/tof/scripts/rec_central_electrons.cxx
 create mode 100644 benchmarks/tof/scripts/rec_multiple_tracks.cxx

diff --git a/benchmarks/tof/central_electrons.sh b/benchmarks/tof/central_electrons.sh
new file mode 100644
index 00000000..5839c187
--- /dev/null
+++ b/benchmarks/tof/central_electrons.sh
@@ -0,0 +1,130 @@
+#!/bin/bash
+
+function print_the_help {
+  echo "USAGE: ${0} [--rec-only]  "
+  echo "  OPTIONS: "
+  echo "            --rec-only   only run gaudi reconstruction part"
+  exit 
+}
+
+REC_ONLY=
+ANALYSIS_ONLY=
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    --rec-only)
+      REC_ONLY=1
+      shift # past value
+      ;;
+    --ana-only)
+      ANALYSIS_ONLY=1
+      shift # past value
+      ;;
+    *)    # unknown option
+      #POSITIONAL+=("$1") # save it in an array for later
+      echo "unknown option $1"
+      print_the_help
+      shift # past argument
+      ;;
+  esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+
+print_env.sh
+
+## To run the reconstruction, we need the following global variables:
+## - JUGGLER_INSTALL_PREFIX:   Install prefix for Juggler (simu/recon)
+## - JUGGLER_DETECTOR:         the detector package we want to use for this benchmark
+## - JUGGLER_DETECTOR_VERSION: the detector package we want to use for this benchmark
+## - DETECTOR_PATH:            full path to the detector definitions
+##
+## You can ready options/env.sh for more in-depth explanations of the variables
+## and how they can be controlled.
+export DETECTOR_PATH=${DETECTOR_PATH}
+
+if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
+  export JUGGLER_N_EVENTS=100
+fi
+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_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+
+if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
+then
+
+  ## generate the input events
+  root -b -q "benchmarks/tof/scripts/gen_central_electrons.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running script"
+    exit 1
+  fi
+
+  echo "Running geant4 simulation"
+  ## run geant4 simulations
+  npsim --runType batch \
+    --part.minimalKineticEnergy 1000*GeV  \
+    -v WARNING \
+    --numberOfEvents ${JUGGLER_N_EVENTS} \
+    --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+    --inputFiles  ${JUGGLER_FILE_NAME_TAG}.hepmc \
+    --outputFile  ${JUGGLER_SIM_FILE}
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running script"
+    exit 1
+  fi
+fi
+
+rootls -t ${JUGGLER_SIM_FILE}
+
+if [[ -z "${ANALYSIS_ONLY}" ]] ;
+then
+  # Need to figure out how to pass file name to juggler from the commandline
+  xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv gaudirun.py benchmarks/tof/options/central_electron_reconstruction.py
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running juggler"
+    exit 1
+  fi
+fi
+
+mkdir -p results/tof
+
+root -b -q "benchmarks/tof/scripts/rec_central_electrons.cxx(\"${JUGGLER_REC_FILE}\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running root script"
+  exit 1
+fi
+
+#root -b -q "benchmarks/tof/scripts/hits_central_electrons.cxx(\"${JUGGLER_SIM_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 
+  # file must be less than 10 MB to upload
+  if [[ "${root_filesize}" -lt "10000000" ]] ; then 
+    cp ${JUGGLER_REC_FILE} results/.
+  fi
+fi
+
+cp tracking_geometry.obj results/.
+
diff --git a/benchmarks/tof/config.yml b/benchmarks/tof/config.yml
new file mode 100644
index 00000000..569c6414
--- /dev/null
+++ b/benchmarks/tof/config.yml
@@ -0,0 +1,14 @@
+tof:multiple_tracks:
+  extends: .rec_benchmark
+  stage: run
+  timeout: 24 hours
+  script:
+    - bash benchmarks/tof/multiple_tracks.sh
+
+tof:collect:
+  stage: collect
+  needs:
+    - ["tof:multiple_tracks"]
+  script:
+    - echo "Done collecting artifacts."
+
diff --git a/benchmarks/tof/multiple_tracks.sh b/benchmarks/tof/multiple_tracks.sh
new file mode 100644
index 00000000..d5790802
--- /dev/null
+++ b/benchmarks/tof/multiple_tracks.sh
@@ -0,0 +1,107 @@
+#!/bin/bash
+
+function print_the_help {
+  echo "USAGE: ${0}   "
+  echo "  OPTIONS: "
+  exit 
+}
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    #--ana-only)
+    #  ANALYSIS_ONLY=1
+    #  shift # past value
+    #  ;;
+    *)    # unknown option
+      #POSITIONAL+=("$1") # save it in an array for later
+      echo "unknown option $1"
+      print_the_help
+      shift # past argument
+      ;;
+  esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+
+print_env.sh
+
+## To run the reconstruction, we need the following global variables:
+## - JUGGLER_INSTALL_PREFIX:   Install prefix for Juggler (simu/recon)
+## - JUGGLER_DETECTOR:         the detector package we want to use for this benchmark
+## - JUGGLER_DETECTOR_VERSION: the detector package we want to use for this benchmark
+## - DETECTOR_PATH:            full path to the detector definitions
+
+if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
+  export JUGGLER_N_EVENTS=100
+fi
+#export JUGGLER_N_EVENTS=$(expr ${JUGGLER_N_EVENTS} \* 1)
+
+export 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_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/tof/scripts/gen_multiple_tracks.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script"
+  exit 1
+fi
+
+echo "Running geant4 simulation"
+## run geant4 simulations
+npsim --runType batch \
+  --part.minimalKineticEnergy 1000*GeV  \
+  -v WARNING \
+  --numberOfEvents ${JUGGLER_N_EVENTS} \
+  --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+  --inputFiles  ${JUGGLER_FILE_NAME_TAG}.hepmc \
+  --outputFile  ${JUGGLER_SIM_FILE}
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script"
+  exit 1
+fi
+
+rootls -t ${JUGGLER_SIM_FILE}
+
+if [[ -z "${ANALYSIS_ONLY}" ]] ;
+then
+  # Need to figure out how to pass file name to juggler from the commandline
+  xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv gaudirun.py benchmarks/tof/options/track_reconstruction.py
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running juggler"
+    exit 1
+  fi
+fi
+
+rootls -t ${JUGGLER_REC_FILE}
+
+mkdir -p results/tof
+
+root -b -q "benchmarks/tof/scripts/rec_multiple_tracks.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 
+  # file must be less than 10 MB to upload
+  if [[ "${root_filesize}" -lt "10000000" ]] ; then 
+    cp ${JUGGLER_REC_FILE} results/.
+  fi
+fi
+
diff --git a/benchmarks/tof/options/central_electron_reconstruction.py b/benchmarks/tof/options/central_electron_reconstruction.py
new file mode 100644
index 00000000..d4fbac2d
--- /dev/null
+++ b/benchmarks/tof/options/central_electron_reconstruction.py
@@ -0,0 +1,239 @@
+from Gaudi.Configuration import *
+
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
+
+detector_name = "topside"
+if "JUGGLER_DETECTOR" in os.environ :
+  detector_name = str(os.environ["JUGGLER_DETECTOR"])
+
+detector_path = ""
+if "DETECTOR_PATH" in os.environ :
+  detector_path = str(os.environ["DETECTOR_PATH"])
+
+# todo add checks
+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)
+podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=WARNING)
+
+from Configurables import PodioInput
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__TrackerHitCollection_dd4pod__TrackerHitCollection_ as TrkCopier
+
+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__TrackParamTruthInit as TrackParamTruthInit
+from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit
+from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit
+
+#from Configurables import Jug__Reco__ConformalXYPeakProtoTracks as ConformalXYPeakProtoTracks
+from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
+from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
+
+from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
+
+
+algorithms = [ ]
+
+podioinput = PodioInput("PodioReader",
+                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","BarrelTOFHits","ForwardTOFHits","BackwardTOFHits"])
+                        #collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits","BarrelTOFHits","ForwardTOFHits","BackwardTOFHits"])
+algorithms.append( podioinput )
+
+## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
+copier = MCCopier("MCCopier", 
+        inputCollection="mcparticles", 
+        outputCollection="mcparticles2") 
+algorithms.append( copier )
+
+trkcopier = TrkCopier("TrkCopier", 
+        inputCollection="TrackerBarrelHits", 
+        outputCollection="TrackerBarrelHits2") 
+algorithms.append( trkcopier )
+
+trk_b_digi = TrackerDigi("trk_b_digi", 
+        inputHitCollection="TrackerBarrelHits",
+        outputHitCollection="TrackerBarrelRawHits",
+        timeResolution=8)
+algorithms.append( trk_b_digi )
+trk_ec_digi = TrackerDigi("trk_ec_digi", 
+        inputHitCollection="TrackerEndcapHits",
+        outputHitCollection="TrackerEndcapRawHits",
+        timeResolution=8)
+algorithms.append( trk_ec_digi )
+
+vtx_b_digi = TrackerDigi("vtx_b_digi", 
+        inputHitCollection="VertexBarrelHits",
+        outputHitCollection="VertexBarrelRawHits",
+        timeResolution=8)
+algorithms.append( vtx_b_digi )
+
+vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
+        inputHitCollection="VertexEndcapHits",
+        outputHitCollection="VertexEndcapRawHits",
+        timeResolution=8)
+algorithms.append( vtx_ec_digi )
+
+#gem_ec_digi = TrackerDigi("gem_ec_digi",
+#        inputHitCollection="GEMTrackerEndcapHits",
+#        outputHitCollection="GEMTrackerEndcapRawHits",
+#        timeResolution=10)
+#algorithms.append(gem_ec_digi)
+
+tof_b_digi = TrackerDigi("tof_b_digi",
+        inputHitCollection="BarrelTOFHits",
+        outputHitCollection="BarrelTOFRawHits",
+        timeResolution=20)
+algorithms.append(tof_b_digi)
+
+tof_forward_digi = TrackerDigi("tof_forward_digi",
+        inputHitCollection="ForwardTOFHits",
+        outputHitCollection="ForwardTOFRawHits",
+        timeResolution=20)
+algorithms.append(tof_forward_digi)
+
+tof_backward_digi = TrackerDigi("tof_backward_digi",
+        inputHitCollection="BackwardTOFHits",
+        outputHitCollection="BackwardTOFRawHits",
+        timeResolution=20)
+algorithms.append(tof_backward_digi)
+
+# Tracker and vertex reconstruction
+trk_b_reco = TrackerHitReconstruction("trk_b_reco",
+        inputHitCollection = trk_b_digi.outputHitCollection,
+        outputHitCollection="TrackerBarrelRecHits")
+algorithms.append( trk_b_reco )
+
+trk_ec_reco = TrackerHitReconstruction("trk_ec_reco",
+        inputHitCollection = trk_ec_digi.outputHitCollection,
+        outputHitCollection="TrackerEndcapRecHits")
+algorithms.append( trk_ec_reco )
+
+vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
+        inputHitCollection = vtx_b_digi.outputHitCollection,
+        outputHitCollection="VertexBarrelRecHits")
+algorithms.append( vtx_b_reco )
+
+vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
+        inputHitCollection = vtx_ec_digi.outputHitCollection,
+        outputHitCollection="VertexEndcapRecHits")
+algorithms.append( vtx_ec_reco )
+
+#gem_ec_reco = TrackerHitReconstruction("gem_ec_reco",
+#        inputHitCollection=gem_ec_digi.outputHitCollection,
+#        outputHitCollection="GEMTrackerEndcapRecHits")
+#algorithms.append(gem_ec_reco)
+
+tof_b_reco = TrackerHitReconstruction("tof_b_reco",
+        inputHitCollection=tof_b_digi.outputHitCollection,
+        outputHitCollection="BarrelTOFRecHits")
+algorithms.append(tof_b_reco)
+
+tof_forward_reco = TrackerHitReconstruction("tof_forward_reco",
+        inputHitCollection=tof_forward_digi.outputHitCollection,
+        outputHitCollection="ForwardTOFRecHits")
+algorithms.append(tof_forward_reco)
+
+tof_backward_reco = TrackerHitReconstruction("tof_backward_reco",
+        inputHitCollection=tof_backward_digi.outputHitCollection,
+        outputHitCollection="BackwardTOFRecHits")
+algorithms.append(tof_backward_reco)
+
+trk_hit_col = TrackingHitsCollector("trk_hit_col",
+        inputTrackingHits=[
+            str(trk_b_reco.outputHitCollection),
+            str(trk_ec_reco.outputHitCollection),
+            str(vtx_b_reco.outputHitCollection),
+            str(vtx_ec_reco.outputHitCollection),
+#            str(gem_ec_reco.outputHitCollection),
+            str(tof_b_reco.outputHitCollection),
+            str(tof_forward_reco.outputHitCollection),
+            str(tof_backward_reco.outputHitCollection) ],
+        trackingHits="trackingHits",
+        OutputLevel=DEBUG)
+algorithms.append( trk_hit_col )
+
+# track finding
+#conformal_find = ConformalXYPeakProtoTracks("conformal_find",
+#        inputTrackerHits=trk_hit_col.trackingHits,
+#        outputProtoTracks="outputProtoTracks",
+#        nProtoTracks="nProtoTracks",
+#        nPhiBins=90,
+#        OutputLevel=DEBUG)
+#algorithms.append( conformal_find )
+
+# Hit Source linker 
+sourcelinker = TrackerSourceLinker("sourcelinker",
+        inputHitCollection=trk_hit_col.trackingHits,
+        outputSourceLinks="TrackSourceLinks",
+        outputMeasurements="TrackMeasurements",        OutputLevel=DEBUG)
+algorithms.append( sourcelinker )
+
+## Track param init
+truth_trk_init = TrackParamTruthInit("truth_trk_init",
+        inputMCParticles="mcparticles",
+        outputInitialTrackParameters="InitTrackParams")
+        #OutputLevel=DEBUG)
+algorithms.append( truth_trk_init )
+
+# Tracking algorithms
+trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
+        inputSourceLinks = sourcelinker.outputSourceLinks,
+        inputMeasurements = sourcelinker.outputMeasurements,
+        inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
+        outputTrajectories="trajectories")
+        #OutputLevel=DEBUG)
+algorithms.append( trk_find_alg )
+
+parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
+        inputTrajectories="trajectories",
+        outputParticles="ReconstructedParticles",
+        outputTrackParameters="outputTrackParameters")
+        #OutputLevel=DEBUG)
+algorithms.append( parts_from_fit )
+
+#types = []
+## this printout is useful to check that the type information is passed to python correctly
+#print("---------------------------------------\n")
+#print("---\n# List of input and output types by class")
+#for configurable in sorted([ PodioInput, EICDataSvc, PodioOutput,
+#                             TrackerHitReconstruction,ExampleCaloDigi, 
+#                             UFSDTrackerDigi, TrackerSourceLinker,
+#                             PodioOutput],
+#                           key=lambda c: c.getType()):
+#    print("\"{}\":".format(configurable.getType()))
+#    props = configurable.getDefaultProperties()
+#    for propname, prop in sorted(props.items()):
+#        print(" prop name: {}".format(propname))
+#        if isinstance(prop, DataHandleBase):
+#            types.append(prop.type())
+#            print("  {}: \"{}\"".format(propname, prop.type()))
+#print("---")
+
+out = PodioOutput("out", filename=output_rec_file)
+out.outputCommands = ["keep *", 
+        "drop BarrelTrackSourceLinks", 
+        "drop InitTrackParams",
+        "drop trajectories",
+        "drop outputSourceLinks",
+        "drop outputInitialTrackParameters",
+        "drop mcparticles"
+        ]
+algorithms.append(out)
+
+ApplicationMgr(
+    TopAlg = algorithms,
+    EvtSel = 'NONE',
+    EvtMax   = n_events,
+    ExtSvc = [podioevent,geo_service],
+    OutputLevel=WARNING
+ )
+
diff --git a/benchmarks/tof/options/central_electron_reconstruction.py.org b/benchmarks/tof/options/central_electron_reconstruction.py.org
new file mode 100644
index 00000000..7086d8a9
--- /dev/null
+++ b/benchmarks/tof/options/central_electron_reconstruction.py.org
@@ -0,0 +1,193 @@
+from Gaudi.Configuration import *
+
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
+
+detector_name = "topside"
+if "JUGGLER_DETECTOR" in os.environ :
+  detector_name = str(os.environ["JUGGLER_DETECTOR"])
+
+detector_path = ""
+if "DETECTOR_PATH" in os.environ :
+  detector_path = str(os.environ["DETECTOR_PATH"])
+
+# todo add checks
+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)
+podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=WARNING)
+
+from Configurables import PodioInput
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__TrackerHitCollection_dd4pod__TrackerHitCollection_ as TrkCopier
+
+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__TrackParamTruthInit as TrackParamTruthInit
+from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit
+from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit
+
+from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
+from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
+
+from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
+
+
+algorithms = [ ]
+
+podioinput = PodioInput("PodioReader", 
+                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits"])#, OutputLevel=DEBUG)
+algorithms.append( podioinput )
+
+## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
+copier = MCCopier("MCCopier", 
+        inputCollection="mcparticles", 
+        outputCollection="mcparticles2") 
+algorithms.append( copier )
+
+trkcopier = TrkCopier("TrkCopier", 
+        inputCollection="TrackerBarrelHits", 
+        outputCollection="TrackerBarrelHits2") 
+algorithms.append( trkcopier )
+
+trk_b_digi = TrackerDigi("trk_b_digi", 
+        inputHitCollection="TrackerBarrelHits",
+        outputHitCollection="TrackerBarrelRawHits",
+        timeResolution=8)
+algorithms.append( trk_b_digi )
+trk_ec_digi = TrackerDigi("trk_ec_digi", 
+        inputHitCollection="TrackerEndcapHits",
+        outputHitCollection="TrackerEndcapRawHits",
+        timeResolution=8)
+algorithms.append( trk_ec_digi )
+
+vtx_b_digi = TrackerDigi("vtx_b_digi", 
+        inputHitCollection="VertexBarrelHits",
+        outputHitCollection="VertexBarrelRawHits",
+        timeResolution=8)
+algorithms.append( vtx_b_digi )
+
+vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
+        inputHitCollection="VertexEndcapHits",
+        outputHitCollection="VertexEndcapRawHits",
+        timeResolution=8)
+algorithms.append( vtx_ec_digi )
+
+gem_ec_digi = TrackerDigi("gem_ec_digi",
+        inputHitCollection="GEMTrackerEndcapHits",
+        outputHitCollection="GEMTrackerEndcapRawHits",
+        timeResolution=10)
+algorithms.append(gem_ec_digi)
+
+# Tracker and vertex reconstruction
+trk_b_reco = TrackerHitReconstruction("trk_b_reco",
+        inputHitCollection = trk_b_digi.outputHitCollection,
+        outputHitCollection="TrackerBarrelRecHits")
+algorithms.append( trk_b_reco )
+
+trk_ec_reco = TrackerHitReconstruction("trk_ec_reco",
+        inputHitCollection = trk_ec_digi.outputHitCollection,
+        outputHitCollection="TrackerEndcapRecHits")
+algorithms.append( trk_ec_reco )
+
+vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
+        inputHitCollection = vtx_b_digi.outputHitCollection,
+        outputHitCollection="VertexBarrelRecHits")
+algorithms.append( vtx_b_reco )
+
+vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
+        inputHitCollection = vtx_ec_digi.outputHitCollection,
+        outputHitCollection="VertexEndcapRecHits")
+algorithms.append( vtx_ec_reco )
+
+gem_ec_reco = TrackerHitReconstruction("gem_ec_reco",
+        inputHitCollection=gem_ec_digi.outputHitCollection,
+        outputHitCollection="GEMTrackerEndcapRecHits")
+algorithms.append(gem_ec_reco)
+
+trk_hit_col = TrackingHitsCollector("trk_hit_col",
+        inputTrackingHits=[
+            str(trk_b_reco.outputHitCollection),
+            str(trk_ec_reco.outputHitCollection),
+            str(vtx_b_reco.outputHitCollection),
+            str(vtx_ec_reco.outputHitCollection),
+            str(gem_ec_reco.outputHitCollection) ],
+        trackingHits="trackingHits",
+        OutputLevel=DEBUG)
+algorithms.append( trk_hit_col )
+
+# Hit Source linker 
+sourcelinker = TrackerSourceLinker("sourcelinker",
+        inputHitCollection=trk_hit_col.trackingHits,
+        outputSourceLinks="TrackSourceLinks",
+        outputMeasurements="TrackMeasurements",
+        OutputLevel=DEBUG)
+algorithms.append( sourcelinker )
+
+## Track param init
+truth_trk_init = TrackParamTruthInit("truth_trk_init",
+        inputMCParticles="mcparticles",
+        outputInitialTrackParameters="InitTrackParams")
+        #OutputLevel=DEBUG)
+algorithms.append( truth_trk_init )
+
+# Tracking algorithms
+trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
+        inputSourceLinks = sourcelinker.outputSourceLinks,
+        inputMeasurements = sourcelinker.outputMeasurements,
+        inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
+        outputTrajectories="trajectories")
+        #OutputLevel=DEBUG)
+algorithms.append( trk_find_alg )
+
+parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
+        inputTrajectories="trajectories",
+        outputParticles="ReconstructedParticles",
+        outputTrackParameters="outputTrackParameters")
+        #OutputLevel=DEBUG)
+algorithms.append( parts_from_fit )
+
+#types = []
+## this printout is useful to check that the type information is passed to python correctly
+#print("---------------------------------------\n")
+#print("---\n# List of input and output types by class")
+#for configurable in sorted([ PodioInput, EICDataSvc, PodioOutput,
+#                             TrackerHitReconstruction,ExampleCaloDigi, 
+#                             UFSDTrackerDigi, TrackerSourceLinker,
+#                             PodioOutput],
+#                           key=lambda c: c.getType()):
+#    print("\"{}\":".format(configurable.getType()))
+#    props = configurable.getDefaultProperties()
+#    for propname, prop in sorted(props.items()):
+#        print(" prop name: {}".format(propname))
+#        if isinstance(prop, DataHandleBase):
+#            types.append(prop.type())
+#            print("  {}: \"{}\"".format(propname, prop.type()))
+#print("---")
+
+out = PodioOutput("out", filename=output_rec_file)
+out.outputCommands = ["keep *", 
+        "drop BarrelTrackSourceLinks", 
+        "drop InitTrackParams",
+        "drop trajectories",
+        "drop outputSourceLinks",
+        "drop outputInitialTrackParameters",
+        "drop mcparticles"
+        ]
+algorithms.append(out)
+
+ApplicationMgr(
+    TopAlg = algorithms,
+    EvtSel = 'NONE',
+    EvtMax   = n_events,
+    ExtSvc = [podioevent,geo_service],
+    OutputLevel=WARNING
+ )
+
diff --git a/benchmarks/tof/options/track_reconstruction.py b/benchmarks/tof/options/track_reconstruction.py
new file mode 100644
index 00000000..0fa1fb6e
--- /dev/null
+++ b/benchmarks/tof/options/track_reconstruction.py
@@ -0,0 +1,239 @@
+from Gaudi.Configuration import *
+
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
+
+detector_name = "topside"
+if "JUGGLER_DETECTOR" in os.environ :
+  detector_name = str(os.environ["JUGGLER_DETECTOR"])
+
+detector_path = ""
+if "DETECTOR_PATH" in os.environ :
+  detector_path = str(os.environ["DETECTOR_PATH"])
+
+# todo add checks
+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)
+podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=WARNING)
+
+from Configurables import PodioInput
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__TrackerHitCollection_dd4pod__TrackerHitCollection_ as TrkCopier
+
+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__TrackParamTruthInit as TrackParamTruthInit
+from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit
+from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit
+
+from Configurables import Jug__Reco__ConformalXYPeakProtoTracks as ConformalXYPeakProtoTracks
+from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
+from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
+
+from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
+
+
+algorithms = [ ]
+
+podioinput = PodioInput("PodioReader",
+                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","BarrelTOFHits","ForwardTOFHits","BackwardTOFHits"])
+                        #collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits","BarrelTOFHits","ForwardTOFHits","BackwardTOFHits"])
+algorithms.append( podioinput )
+
+## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
+copier = MCCopier("MCCopier", 
+        inputCollection="mcparticles", 
+        outputCollection="mcparticles2") 
+algorithms.append( copier )
+
+trkcopier = TrkCopier("TrkCopier", 
+        inputCollection="TrackerBarrelHits", 
+        outputCollection="TrackerBarrelHits2") 
+algorithms.append( trkcopier )
+
+trk_b_digi = TrackerDigi("trk_b_digi", 
+        inputHitCollection="TrackerBarrelHits",
+        outputHitCollection="TrackerBarrelRawHits",
+        timeResolution=8)
+algorithms.append( trk_b_digi )
+trk_ec_digi = TrackerDigi("trk_ec_digi", 
+        inputHitCollection="TrackerEndcapHits",
+        outputHitCollection="TrackerEndcapRawHits",
+        timeResolution=8)
+algorithms.append( trk_ec_digi )
+
+vtx_b_digi = TrackerDigi("vtx_b_digi", 
+        inputHitCollection="VertexBarrelHits",
+        outputHitCollection="VertexBarrelRawHits",
+        timeResolution=8)
+algorithms.append( vtx_b_digi )
+
+vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
+        inputHitCollection="VertexEndcapHits",
+        outputHitCollection="VertexEndcapRawHits",
+        timeResolution=8)
+algorithms.append( vtx_ec_digi )
+
+#gem_ec_digi = TrackerDigi("gem_ec_digi",
+#        inputHitCollection="GEMTrackerEndcapHits",
+#        outputHitCollection="GEMTrackerEndcapRawHits",
+#        timeResolution=10)
+#algorithms.append(gem_ec_digi)
+
+tof_b_digi = TrackerDigi("tof_b_digi",
+        inputHitCollection="BarrelTOFHits",
+        outputHitCollection="BarrelTOFRawHits",
+        timeResolution=20)
+algorithms.append(tof_b_digi)
+
+tof_forward_digi = TrackerDigi("tof_forward_digi",
+        inputHitCollection="ForwardTOFHits",
+        outputHitCollection="ForwardTOFRawHits",
+        timeResolution=20)
+algorithms.append(tof_forward_digi)
+
+tof_backward_digi = TrackerDigi("tof_backward_digi",
+        inputHitCollection="BackwardTOFHits",
+        outputHitCollection="BackwardTOFRawHits",
+        timeResolution=20)
+algorithms.append(tof_backward_digi)
+
+# Tracker and vertex reconstruction
+trk_b_reco = TrackerHitReconstruction("trk_b_reco",
+        inputHitCollection = trk_b_digi.outputHitCollection,
+        outputHitCollection="TrackerBarrelRecHits")
+algorithms.append( trk_b_reco )
+
+trk_ec_reco = TrackerHitReconstruction("trk_ec_reco",
+        inputHitCollection = trk_ec_digi.outputHitCollection,
+        outputHitCollection="TrackerEndcapRecHits")
+algorithms.append( trk_ec_reco )
+
+vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
+        inputHitCollection = vtx_b_digi.outputHitCollection,
+        outputHitCollection="VertexBarrelRecHits")
+algorithms.append( vtx_b_reco )
+
+vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
+        inputHitCollection = vtx_ec_digi.outputHitCollection,
+        outputHitCollection="VertexEndcapRecHits")
+algorithms.append( vtx_ec_reco )
+
+#gem_ec_reco = TrackerHitReconstruction("gem_ec_reco",
+#        inputHitCollection=gem_ec_digi.outputHitCollection,
+#        outputHitCollection="GEMTrackerEndcapRecHits")
+#algorithms.append(gem_ec_reco)
+
+tof_b_reco = TrackerHitReconstruction("tof_b_reco",
+        inputHitCollection=tof_b_digi.outputHitCollection,
+        outputHitCollection="BarrelTOFRecHits")
+algorithms.append(tof_b_reco)
+
+tof_forward_reco = TrackerHitReconstruction("tof_forward_reco",
+        inputHitCollection=tof_forward_digi.outputHitCollection,
+        outputHitCollection="ForwardTOFRecHits")
+algorithms.append(tof_forward_reco)
+
+tof_backward_reco = TrackerHitReconstruction("tof_backward_reco",
+        inputHitCollection=tof_backward_digi.outputHitCollection,
+        outputHitCollection="BackwardTOFRecHits")
+algorithms.append(tof_backward_reco)
+
+trk_hit_col = TrackingHitsCollector("trk_hit_col",
+        inputTrackingHits=[
+            str(trk_b_reco.outputHitCollection),
+            str(trk_ec_reco.outputHitCollection),
+            str(vtx_b_reco.outputHitCollection),
+            str(vtx_ec_reco.outputHitCollection),
+#            str(gem_ec_reco.outputHitCollection),
+            str(tof_b_reco.outputHitCollection),
+            str(tof_forward_reco.outputHitCollection),
+            str(tof_backward_reco.outputHitCollection) ],
+        trackingHits="trackingHits",
+        OutputLevel=DEBUG)
+algorithms.append( trk_hit_col )
+
+# track finding
+conformal_find = ConformalXYPeakProtoTracks("conformal_find",
+        inputTrackerHits=trk_hit_col.trackingHits,
+        outputProtoTracks="outputProtoTracks",
+        nProtoTracks="nProtoTracks",
+        nPhiBins=90,
+        OutputLevel=DEBUG)
+algorithms.append( conformal_find )
+
+# Hit Source linker 
+sourcelinker = TrackerSourceLinker("sourcelinker",
+        inputHitCollection=trk_hit_col.trackingHits,
+        outputSourceLinks="TrackSourceLinks",
+        outputMeasurements="TrackMeasurements")
+algorithms.append( sourcelinker )
+
+## Track param init
+truth_trk_init = TrackParamTruthInit("truth_trk_init",
+        inputMCParticles="mcparticles",
+        outputInitialTrackParameters="InitTrackParams")
+        #OutputLevel=DEBUG)
+algorithms.append( truth_trk_init )
+
+# Tracking algorithms
+trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
+        inputSourceLinks = sourcelinker.outputSourceLinks,
+        inputMeasurements = sourcelinker.outputMeasurements,
+        inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
+        outputTrajectories="trajectories")
+        #OutputLevel=DEBUG)
+algorithms.append( trk_find_alg )
+
+parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
+        inputTrajectories="trajectories",
+        outputParticles="ReconstructedParticles",
+        outputTrackParameters="outputTrackParameters")
+        #OutputLevel=DEBUG)
+algorithms.append( parts_from_fit )
+
+#types = []
+## this printout is useful to check that the type information is passed to python correctly
+#print("---------------------------------------\n")
+#print("---\n# List of input and output types by class")
+#for configurable in sorted([ PodioInput, EICDataSvc, PodioOutput,
+#                             TrackerHitReconstruction,ExampleCaloDigi, 
+#                             UFSDTrackerDigi, TrackerSourceLinker,
+#                             PodioOutput],
+#                           key=lambda c: c.getType()):
+#    print("\"{}\":".format(configurable.getType()))
+#    props = configurable.getDefaultProperties()
+#    for propname, prop in sorted(props.items()):
+#        print(" prop name: {}".format(propname))
+#        if isinstance(prop, DataHandleBase):
+#            types.append(prop.type())
+#            print("  {}: \"{}\"".format(propname, prop.type()))
+#print("---")
+
+out = PodioOutput("out", filename=output_rec_file)
+out.outputCommands = ["keep *", 
+        "drop BarrelTrackSourceLinks", 
+        "drop InitTrackParams",
+        "drop trajectories",
+        "drop outputSourceLinks",
+        "drop outputInitialTrackParameters",
+        "drop mcparticles"
+        ]
+algorithms.append(out)
+
+ApplicationMgr(
+    TopAlg = algorithms,
+    EvtSel = 'NONE',
+    EvtMax   = n_events,
+    ExtSvc = [podioevent,geo_service],
+    OutputLevel=WARNING
+ )
+
diff --git a/benchmarks/tof/scripts/gen_central_electrons.cxx b/benchmarks/tof/scripts/gen_central_electrons.cxx
new file mode 100644
index 00000000..3e9d8c01
--- /dev/null
+++ b/benchmarks/tof/scripts/gen_central_electrons.cxx
@@ -0,0 +1,83 @@
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include <iostream>
+#include<random>
+#include<cmath>
+#include <math.h>
+#include <TMath.h>
+
+using namespace HepMC3;
+
+/** Generate electrons in the central region.
+ *  This is for testing detectors in the "barrel" region.
+ */
+void gen_central_electrons(int n_events = 100, 
+                     const char* out_fname = "central_electrons.hepmc")
+{
+  double cos_theta_min = std::cos( 140.0*(M_PI/180.0));
+  double cos_theta_max = std::cos(170.0*(M_PI/180.0));
+
+  WriterAscii hepmc_output(out_fname);
+  int events_parsed = 0;
+  GenEvent evt(Units::GEV, Units::MM);
+
+  // Random number generator
+  TRandom *r1 = new TRandom();
+
+  for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
+    // FourVector(px,py,pz,e,pdgid,status)
+    // type 4 is beam
+    // pdgid 11 - electron
+    // pdgid 111 - pi0
+    // pdgid 2212 - proton
+    GenParticlePtr p1 =
+        std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+    GenParticlePtr p2 = std::make_shared<GenParticle>(
+        FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
+
+    // Define momentum
+    Double_t p     = r1->Uniform(1.0, 10.0);
+    Double_t phi   = r1->Uniform(0.0, 2.0 * M_PI);
+    Double_t costh = r1->Uniform(cos_theta_min, cos_theta_max);
+    Double_t th    = std::acos(costh);
+    Double_t px    = p * std::cos(phi) * std::sin(th);
+    Double_t py    = p * std::sin(phi) * std::sin(th);
+    Double_t pz    = p * std::cos(th);
+    // Generates random vectors, uniformly distributed over the surface of a
+    // sphere of given radius, in this case momentum.
+    // r1->Sphere(px, py, pz, p);
+
+    //std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
+
+    // type 1 is final state
+    // pdgid 11 - electron 0.510 MeV/c^2
+    GenParticlePtr p3 = std::make_shared<GenParticle>(
+        FourVector(
+            px, py, pz,
+            sqrt(p*p + (0.000511 * 0.000511))),
+        11, 1);
+
+    GenVertexPtr v1 = std::make_shared<GenVertex>();
+    v1->add_particle_in(p1);
+    v1->add_particle_in(p2);
+
+    v1->add_particle_out(p3);
+    evt.add_vertex(v1);
+
+    if (events_parsed == 0) {
+      std::cout << "First event: " << std::endl;
+      Print::listing(evt);
+    }
+
+    hepmc_output.write_event(evt);
+    if (events_parsed % 10000 == 0) {
+      std::cout << "Event: " << events_parsed << std::endl;
+    }
+    evt.clear();
+  }
+  hepmc_output.close();
+  std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}
diff --git a/benchmarks/tof/scripts/gen_multiple_tracks.cxx b/benchmarks/tof/scripts/gen_multiple_tracks.cxx
new file mode 100644
index 00000000..58f74acc
--- /dev/null
+++ b/benchmarks/tof/scripts/gen_multiple_tracks.cxx
@@ -0,0 +1,95 @@
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include <iostream>
+#include <random>
+#include <cmath>
+
+#include "TMath.h"
+
+#include "common_bench/particles.h"
+
+
+using namespace HepMC3;
+
+/** Generate multiple electrons/positron tracks in the central region.
+ *  This is for testing detectors in the "barrel" region.
+ */
+void gen_multiple_tracks(int n_events = 100, 
+                         const char* out_fname = "multiple_tracks.hepmc",
+                         int n_parts = 1)
+{
+  double cos_theta_min = std::cos( 10.0*(M_PI/180.0));
+  double cos_theta_max = std::cos(170.0*(M_PI/180.0));
+
+  WriterAscii hepmc_output(out_fname);
+  int events_parsed = 0;
+  GenEvent evt(Units::GEV, Units::MM);
+
+  // Random number generator
+  TRandom *r1 = new TRandom();
+
+  for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
+    // FourVector(px,py,pz,e,pdgid,status)
+    // type 4 is beam
+    // pdgid 11 - electron
+    // pdgid 111 - pi0
+    // pdgid 2212 - proton
+    for (int ip = 0; ip < n_parts; ip++) {
+      GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+      GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
+
+      // Define momentum
+      Double_t p     = 0.5; //r1->Uniform(1.0, 10.0);
+      Double_t phi   = r1->Uniform(0.0, 2.0 * M_PI);
+      Double_t costh = r1->Uniform(-0.5, 0.5); //r1->Uniform(cos_theta_min, cos_theta_max);
+      Double_t th    = std::acos(costh);
+      Double_t px    = p * std::cos(phi) * std::sin(th);
+      Double_t py    = p * std::sin(phi) * std::sin(th);
+      Double_t pz    = p * std::cos(th);
+      // Generates random vectors, uniformly distributed over the surface of a
+      // sphere of given radius, in this case momentum.
+      // r1->Sphere(px, py, pz, p);
+
+      // std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
+
+      // type 1 is final state
+      // pdgid 11 - electron 0.510 MeV/c^2
+      GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.000511 * 0.000511))),
+                                                        ((ip % 2 == 0) ? 11 : -11), 1);
+      GenParticlePtr p4 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.938272 * 0.938272))),
+                                                            ((ip % 2 == 0) ? 2212 : -2212), 1);
+      GenParticlePtr p5 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.493677 * 0.493677))),
+                                                            ((ip % 2 == 0) ? 321 : -321), 1);
+      GenParticlePtr p6 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.139570 * 0.139570))),
+                                                            ((ip % 2 == 0) ? 211 : -211), 1);
+
+
+      GenVertexPtr v1 = std::make_shared<GenVertex>();
+      v1->add_particle_in(p1);
+      v1->add_particle_in(p2);
+
+      v1->add_particle_out(p3);
+      v1->add_particle_out(p4);
+      v1->add_particle_out(p5);
+      v1->add_particle_out(p6);
+
+      evt.add_vertex(v1);
+    }
+
+    if (events_parsed == 0) {
+      std::cout << "First event: " << std::endl;
+      Print::listing(evt);
+    }
+
+    hepmc_output.write_event(evt);
+    if (events_parsed % 10000 == 0) {
+      std::cout << "Event: " << events_parsed << std::endl;
+    }
+    evt.clear();
+  }
+  hepmc_output.close();
+  std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}
diff --git a/benchmarks/tof/scripts/rec_central_electrons.cxx b/benchmarks/tof/scripts/rec_central_electrons.cxx
new file mode 100644
index 00000000..c8c00737
--- /dev/null
+++ b/benchmarks/tof/scripts/rec_central_electrons.cxx
@@ -0,0 +1,249 @@
+#include "ROOT/RDataFrame.hxx"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TH1D.h"
+#include "TProfile.h"
+
+#include <iostream>
+
+R__LOAD_LIBRARY(libeicd.so)
+R__LOAD_LIBRARY(libDD4pod.so)
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/TrackParametersCollection.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.h"
+#include "eicd/TrackerHitCollection.h"
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
+  std::vector<double> result;
+  for (size_t i = 0; i < in.size(); ++i) {
+    result.push_back(std::abs(1.0/(in[i].qOverP)));
+  }
+  return result;
+};
+
+
+std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+  std::vector<float> result;
+  for (size_t i = 0; i < in.size(); ++i) {
+    result.push_back(std::sqrt(in[i].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y));
+  }
+  return result;
+}
+
+auto momentum = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
+  std::vector<double> result;
+  for (size_t i = 0; i < in.size(); ++i) {
+   result.push_back(in[i].P());
+  }
+  return result;
+};
+auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
+  std::vector<double> result;
+  for (size_t i = 0; i < in.size(); ++i) {
+   result.push_back(in[i].Theta()*180/M_PI);
+  }
+  return result;
+};
+auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+  std::vector<ROOT::Math::PxPyPzMVector> result;
+  ROOT::Math::PxPyPzMVector lv;
+  for (size_t i = 0; i < in.size(); ++i) {
+    lv.SetCoordinates(in[i].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
+    result.push_back(lv);
+  }
+  return result;
+};
+
+auto delta_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
+  std::vector<double> res;
+  for (const auto& p1 : thrown) {
+    for (const auto& p2 : tracks) {
+      res.push_back(p1 - p2);
+    }
+  }
+  return res;
+};
+
+auto delta_p_over_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
+  std::vector<double> res;
+  for (const auto& p1 : thrown) {
+    for (const auto& p2 : tracks) {
+      res.push_back((p1 - p2)/p1);
+    }
+  }
+  return res;
+};
+
+int rec_central_electrons(const char* fname = "topside/rec_central_electrons.root")
+{
+
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame df("events", fname);
+
+  auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
+                 .Define("thrownParticles", "mcparticles2[isThrown]")
+                 .Define("thrownP", fourvec, {"thrownParticles"})
+                 .Define("p_thrown", momentum, {"thrownP"})
+                 .Define("theta_thrown", theta, {"thrownP"})
+                 .Define("theta0", "theta_thrown[0]")
+                 .Define("nTracks", "outputTrackParameters.size()")
+                 .Define("p_track", p_track, {"outputTrackParameters"})
+                 //.Define("p_track1", p_track, {"outputTrackParameters1"})
+                 //.Define("p_track2", p_track, {"outputTrackParameters2"})
+                 .Define("delta_p0",delta_p, {"p_track", "p_thrown"})
+                 //.Define("delta_p1",delta_p, {"p_track1", "p_thrown"})
+                 //.Define("delta_p2",delta_p, {"p_track2", "p_thrown"})
+                 .Define("delta_p_over_p0",delta_p_over_p, {"p_track", "p_thrown"})
+                 //.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<eic::TrackerHitData> hits) { return hits.size();}, {"trackingHits"})
+                 .Define("N_BarrelHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelRecHits"})
+                 .Define("N_EndcapHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapRecHits"})
+                 ;
+
+  auto h_nTracks_vs_theta = df0.Histo2D({"h_nTracks_vs_theta", "; #theta; N tracks ", 40,0,180,10, 0, 10}, "theta0","nTracks");
+
+  auto h_nTracks = df0.Histo1D({"h_nTracks", "; N tracks ", 10, 0, 10}, "nTracks");
+  auto h_pTracks = df0.Histo1D({"h_pTracks", "; GeV/c ", 100, 0, 10}, "p_track");
+
+  auto h_delta_p0  = df0.Histo1D({"h_delta_p0", "Truth Track Init; GeV/c ",  100, -10, 10}, "delta_p0");
+
+  auto h_delta_p0_over_p = df0.Histo1D({"h_delta_p0_over_p",  "Truth Track Init; delta p/p ",  100, -0.1, 0.1}, "delta_p_over_p0");
+
+  auto hNhits_vs_theta = df0.Histo1D({"hNhits_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_Hits");
+  auto hBarrel_N_vs_theta = df0.Histo1D({"hBarrel_N_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_BarrelHits");
+  auto hEndcap_N_vs_theta = df0.Histo1D({"hEndcap_N_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_EndcapHits");
+  //auto hVtxBarrel_N_vs_theta = df0.Histo1D({"hVtxBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_VtxBarrelHits");
+
+  auto hHits_Nhits  = df0.Histo1D({"hHits_Nhits", "; #theta [deg.]",       20, 0, 20 }, "N_Hits");
+  auto hBarrel_Nhits  = df0.Histo1D({"hBarrel_Nhits", "; #theta [deg.]",   20, 0, 20 }, "N_BarrelHits");
+  auto hEndcap_Nhits  = df0.Histo1D({"hEndcap_Nhits", "; #theta [deg.]",   20, 0, 20 }, "N_EndcapHits");
+  //auto hVtxBarrel_Nhits = df0.Histo1D({"hVtxBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 },  "N_VtxBarrelHits");
+
+  auto hHits_Ntheta = df0.Histo1D({"hHits_Ntheta", "; #theta [deg.]",       40, 0, 180 }, "theta0");
+  auto hBarrel_Ntheta = df0.Histo1D({"hBarrel_Ntheta", "; #theta [deg.]",   40, 0, 180 }, "theta0");
+  auto hEndcap_Ntheta = df0.Histo1D({"hEndcap_Ntheta", "; #theta [deg.]",   40, 0, 180 }, "theta0");
+  //auto hVtxBarrel_Ntheta = df0.Histo1D({"hVtxBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
+
+  // -----------------------------------------------
+  auto c = new TCanvas();
+
+  h_nTracks->DrawCopy();
+  c->SaveAs("results/tof/rec_central_electrons_nTracks.png");
+  c->SaveAs("results/tof/rec_central_electrons_nTracks.pdf");
+
+  // -----------------------------------------------
+  h_pTracks->DrawCopy();
+  c->SaveAs("results/tof/rec_central_electrons_pTracks.png");
+  c->SaveAs("results/tof/rec_central_electrons_pTracks.pdf");
+
+  // -----------------------------------------------
+  c = new TCanvas();
+  THStack * hs = new THStack("hs_delta_p","; GeV/c ");
+  TH1D* h1 = (TH1D*) h_delta_p0->Clone();
+  hs->Add(h1);
+  //h1 = (TH1D*) h_delta_p1->Clone();
+  //h1->SetLineColor(2);
+  //hs->Add(h1);
+  //h1 = (TH1D*) h_delta_p2->Clone();
+  //h1->SetLineColor(4);
+  //h1->SetFillStyle(3001);
+  //h1->SetFillColor(4);
+  //hs->Add(h1);
+  hs->Draw("nostack");
+  c->BuildLegend();
+  c->SaveAs("results/tof/rec_central_electrons_delta_p.png");
+  c->SaveAs("results/tof/rec_central_electrons_delta_p.pdf");
+
+  // -----------------------------------------------
+  c  = new TCanvas();
+  hs = new THStack("hs_delta_p_over_p","; delta p/p ");
+  h1 = (TH1D*) h_delta_p0_over_p->Clone();
+  hs->Add(h1);
+  hs->Draw("nostack");
+  c->BuildLegend();
+  c->SaveAs("results/tof/rec_central_electrons_delta_p_over_p.png");
+  c->SaveAs("results/tof/rec_central_electrons_delta_p_over_p.pdf");
+
+  // -----------------------------------------------
+  c  = new TCanvas();
+  hs = new THStack("n_hits","; #theta  ");
+
+  h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
+  auto h2 = (TH1D*) hBarrel_Ntheta->Clone();
+  h1->SetLineColor(4);
+  h1->Divide(h2);
+  hs->Add(h1);
+
+  h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
+  h2 = (TH1D*) hEndcap_Ntheta->Clone();
+  h1->Divide(h2);
+  h1->SetLineColor(2);
+  hs->Add(h1);
+
+  h1 = (TH1D*) hNhits_vs_theta->Clone();
+  h2 = (TH1D*) hHits_Ntheta->Clone();
+  h1->Divide(h2);
+  h1->SetLineColor(1);
+  hs->Add(h1);
+
+  hs->Draw("nostack, hist");
+  c->BuildLegend();
+  c->SaveAs("results/tof/rec_central_electrons_n_hits_vs_theta.png");
+  c->SaveAs("results/tof/rec_central_electrons_n_hits_vs_theta.pdf");
+
+  // -----------------------------------------------
+  c  = new TCanvas();
+  hs = new THStack("theta","; #theta  ");
+
+  h2 = (TH1D*) hBarrel_Ntheta->Clone();
+  h2->SetLineColor(4);
+  hs->Add(h2);
+
+  h2 = (TH1D*) hEndcap_Ntheta->Clone();
+  h2->SetLineColor(2);
+  hs->Add(h2);
+
+  h2 = (TH1D*) hHits_Ntheta->Clone();
+  h2->SetLineColor(1);
+  hs->Add(h2);
+
+  hs->Draw("nostack hist");
+  c->BuildLegend();
+  c->SaveAs("results/tof/rec_central_electrons_theta.png");
+  c->SaveAs("results/tof/rec_central_electrons_theta.pdf");
+
+  // -----------------------------------------------
+  c  = new TCanvas();
+  hs = new THStack("hits","; hits  ");
+  h1 = (TH1D*) hBarrel_Nhits->Clone();
+  h1->SetLineColor(4);
+  hs->Add(h1);
+  h1 = (TH1D*) hEndcap_Nhits->Clone();
+  h1->SetLineColor(2);
+  hs->Add(h1);
+  h1 = (TH1D*) hHits_Nhits->Clone();
+  h1->SetLineColor(2);
+  hs->Add(h1);
+  //h1 = (TH1D*) hVtxBarrel_Nhits->Clone();
+  //h1->SetLineColor(4);
+  //h1->SetFillStyle(3001);
+  //h1->SetFillColor(4);
+  //hs->Add(h1);
+  //hs->Draw("nostack hist");
+  c->BuildLegend();
+  c->SaveAs("results/tof/rec_central_electrons_nhits.png");
+  c->SaveAs("results/tof/rec_central_electrons_nhits.pdf");
+
+  // -----------------------------------------------
+  c  = new TCanvas();
+  h_nTracks_vs_theta->DrawCopy("colz");
+  c->SaveAs("results/tof/rec_central_electrons_nTracks_vs_theta.png");
+  c->SaveAs("results/tof/rec_central_electrons_nTracks_vs_theta.pdf");
+  return 0;
+}
diff --git a/benchmarks/tof/scripts/rec_multiple_tracks.cxx b/benchmarks/tof/scripts/rec_multiple_tracks.cxx
new file mode 100644
index 00000000..d4d3756f
--- /dev/null
+++ b/benchmarks/tof/scripts/rec_multiple_tracks.cxx
@@ -0,0 +1,248 @@
+#include "ROOT/RDataFrame.hxx"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TH1D.h"
+#include "TProfile.h"
+
+#include <iostream>
+R__LOAD_LIBRARY(libJugBase.so)
+R__LOAD_LIBRARY(libeicd.so)
+R__LOAD_LIBRARY(libDD4pod.so)
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/TrackParametersCollection.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.h"
+#include "eicd/TrackerHitCollection.h"
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
+  std::vector<double> result;
+  for (size_t i = 0; i < in.size(); ++i) {
+    result.push_back(std::abs(1.0/(in[i].qOverP)));
+  }
+  return result;
+};
+
+
+std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
+  std::vector<float> result;
+  for (size_t i = 0; i < in.size(); ++i) {
+    result.push_back(std::sqrt(in[i].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y));
+  }
+  return result;
+}
+
+auto momentum = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
+  std::vector<double> result;
+  for (size_t i = 0; i < in.size(); ++i) {
+   result.push_back(in[i].P());
+  }
+  return result;
+};
+auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
+  std::vector<double> result;
+  for (size_t i = 0; i < in.size(); ++i) {
+   result.push_back(in[i].Theta()*180/M_PI);
+  }
+  return result;
+};
+auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+  std::vector<ROOT::Math::PxPyPzMVector> result;
+  ROOT::Math::PxPyPzMVector lv;
+  for (size_t i = 0; i < in.size(); ++i) {
+    lv.SetCoordinates(in[i].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
+    result.push_back(lv);
+  }
+  return result;
+};
+
+auto delta_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
+  std::vector<double> res;
+  for (const auto& p1 : thrown) {
+    for (const auto& p2 : tracks) {
+      res.push_back(p1 - p2);
+    }
+  }
+  return res;
+};
+
+auto delta_p_over_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
+  std::vector<double> res;
+  for (const auto& p1 : thrown) {
+    for (const auto& p2 : tracks) {
+      res.push_back((p1 - p2)/p1);
+    }
+  }
+  return res;
+};
+//gInterpreter->GenerateDictionary("vector<vector<float> >", "vector")
+//gInterpreter->GenerateDictionary("vector<unsigned long>", "vector")
+//gInterpreter->GenerateDictionary("vector<vector<unsigned long> >", "vector")
+
+int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root")
+{
+
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame df("events", fname);
+
+  auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
+                 .Define("thrownParticles", "mcparticles2[isThrown]")
+                 .Define("thrownP", fourvec, {"thrownParticles"})
+                 .Define("nThrown", "thrownParticles.size()")
+                 .Define("nProto", "outputProtoTracks.size()")
+                 .Define("p_thrown", momentum, {"thrownP"})
+                 .Define("theta_thrown", theta, {"thrownP"})
+                 .Define("theta0", "theta_thrown[0]")
+                 .Define("nTracks", "outputTrackParameters.size()")
+                 .Define("p_track", p_track, {"outputTrackParameters"})
+                 .Define("delta_p0",delta_p, {"p_track", "p_thrown"})
+                 .Define("delta_p_over_p0",delta_p_over_p, {"p_track", "p_thrown"})
+                 .Define("N_Hits",       [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"trackingHits"})
+                 .Define("N_BarrelHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelRecHits"})
+                 .Define("N_EndcapHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapRecHits"})
+                 .Define("N_BarrelTOFHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"BarrelTOFRecHits"})
+                 .Define("N_ForwardTOFHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"ForwardTOFRecHits"})
+                 .Define("N_BackwardTOFHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"BackwardTOFRecHits"})
+                 ;
+
+  auto h_nTracks_vs_theta = df0.Histo2D({"h_nTracks_vs_theta", "; #theta; N tracks ", 40,0,180,10, 0, 10}, "theta0","nTracks");
+
+  auto h_nTracks = df0.Histo1D({"h_nTracks", "; N tracks ", 10, 0, 10}, "nTracks");
+  auto h_pTracks = df0.Histo1D({"h_pTracks", "; GeV/c ", 100, 0, 10}, "p_track");
+
+  auto h_delta_p0  = df0.Histo1D({"h_delta_p0", "Truth Track Init; GeV/c ",  100, -10, 10}, "delta_p0");
+
+  auto h_nProtoTracks = df0.Histo1D({"h_nProtoTracks", "; n ", 10, 0, 10}, "nProto");
+  auto h_nProtoTracks2 = df0.Histo1D({"h_nProtoTracks2", "; n ", 10, 0, 10}, "nProtoTracks");
+  auto h_nThrown      = df0.Histo1D({"h_nThrown",      "; n ", 10, 0, 10}, "nThrown");
+
+  auto h_delta_p0_over_p = df0.Histo1D({"h_delta_p0_over_p",  "Truth Track Init; delta p/p ",  100, -0.1, 0.1}, "delta_p_over_p0");
+
+  auto hNhits_vs_theta = df0.Histo1D({"hNhits_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_Hits");
+  auto hBarrel_N_vs_theta = df0.Histo1D({"hBarrel_N_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_BarrelHits");
+  auto hEndcap_N_vs_theta = df0.Histo1D({"hEndcap_N_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_EndcapHits");
+    auto hBarrelTOF_N_vs_theta = df0.Histo1D({"hBarrelTOF_N_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_BarrelTOFHits");
+    auto hForwardTOF_N_vs_theta = df0.Histo1D({"hForwardTOF_N_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_ForwardTOFHits");
+    auto hBackwardTOF_N_vs_theta = df0.Histo1D({"hBackwardTOF_N_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_BackwardTOFHits");
+
+  auto hHits_Nhits  = df0.Histo1D({"hHits_Nhits", "; #theta [deg.]",       20, 0, 20 }, "N_Hits");
+  auto hBarrel_Nhits  = df0.Histo1D({"hBarrel_Nhits", "; #theta [deg.]",   20, 0, 20 }, "N_BarrelHits");
+  auto hEndcap_Nhits  = df0.Histo1D({"hEndcap_Nhits", "; #theta [deg.]",   20, 0, 20 }, "N_EndcapHits");
+
+  auto hHits_Ntheta = df0.Histo1D({"hHits_Ntheta", "; #theta [deg.]",       40, 0, 180 }, "theta0");
+  auto hBarrel_Ntheta = df0.Histo1D({"hBarrel_Ntheta", "; #theta [deg.]",   40, 0, 180 }, "theta0");
+  auto hEndcap_Ntheta = df0.Histo1D({"hEndcap_Ntheta", "; #theta [deg.]",   40, 0, 180 }, "theta0");
+
+  // -----------------------------------------------
+  auto c = new TCanvas();
+
+  h_nThrown->SetLineColor(2);
+  h_nThrown->DrawCopy();
+  h_nProtoTracks->DrawCopy("same");
+  h_nProtoTracks2->SetLineColor(4);
+  h_nProtoTracks2->DrawCopy("same");
+  c->SaveAs("results/tof/rec_multiple_tracks_nProtoTracks.png");
+  c->SaveAs("results/tof/rec_multiple_tracks_nProtoTracks.pdf");
+
+  h_nTracks->DrawCopy();
+  c->SaveAs("results/tof/rec_multiple_tracks_nTracks.png");
+  c->SaveAs("results/tof/rec_multiple_tracks_nTracks.pdf");
+
+  // -----------------------------------------------
+  h_pTracks->DrawCopy();
+  c->SaveAs("results/tof/rec_multiple_tracks_pTracks.png");
+  c->SaveAs("results/tof/rec_multiple_tracks_pTracks.pdf");
+
+  // -----------------------------------------------
+  c = new TCanvas();
+  THStack * hs = new THStack("hs_delta_p","; GeV/c ");
+  TH1D* h1 = (TH1D*) h_delta_p0->Clone();
+  hs->Add(h1);
+  hs->Draw("nostack");
+  c->BuildLegend();
+  c->SaveAs("results/tof/rec_multiple_tracks_delta_p.png");
+  c->SaveAs("results/tof/rec_multiple_tracks_delta_p.pdf");
+
+  // -----------------------------------------------
+  c  = new TCanvas();
+  hs = new THStack("hs_delta_p_over_p","; delta p/p ");
+  h1 = (TH1D*) h_delta_p0_over_p->Clone();
+  hs->Add(h1);
+  hs->Draw("nostack");
+  c->BuildLegend();
+  c->SaveAs("results/tof/rec_multiple_tracks_delta_p_over_p.png");
+  c->SaveAs("results/tof/rec_multiple_tracks_delta_p_over_p.pdf");
+
+  // -----------------------------------------------
+  c  = new TCanvas();
+  hs = new THStack("n_hits","; #theta  ");
+
+  h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
+  auto h2 = (TH1D*) hBarrel_Ntheta->Clone();
+  h1->SetLineColor(4);
+  h1->Divide(h2);
+  hs->Add(h1);
+
+  h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
+  h2 = (TH1D*) hEndcap_Ntheta->Clone();
+  h1->Divide(h2);
+  h1->SetLineColor(2);
+  hs->Add(h1);
+
+  h1 = (TH1D*) hNhits_vs_theta->Clone();
+  h2 = (TH1D*) hHits_Ntheta->Clone();
+  h1->Divide(h2);
+  h1->SetLineColor(1);
+  hs->Add(h1);
+
+  hs->Draw("nostack, hist");
+  c->BuildLegend();
+  c->SaveAs("results/tof/rec_multiple_tracks_n_hits_vs_theta.png");
+  c->SaveAs("results/tof/rec_multiple_tracks_n_hits_vs_theta.pdf");
+
+  // -----------------------------------------------
+  c  = new TCanvas();
+  hs = new THStack("theta","; #theta  ");
+
+  h2 = (TH1D*) hBarrel_Ntheta->Clone();
+  h2->SetLineColor(4);
+  hs->Add(h2);
+
+  h2 = (TH1D*) hEndcap_Ntheta->Clone();
+  h2->SetLineColor(2);
+  hs->Add(h2);
+
+  h2 = (TH1D*) hHits_Ntheta->Clone();
+  h2->SetLineColor(1);
+  hs->Add(h2);
+
+  hs->Draw("nostack hist");
+  c->BuildLegend();
+  c->SaveAs("results/tof/rec_multiple_tracks_theta.png");
+  c->SaveAs("results/tof/rec_multiple_tracks_theta.pdf");
+
+  // -----------------------------------------------
+  c  = new TCanvas();
+  hs = new THStack("hits","; hits  ");
+  h1 = (TH1D*) hBarrel_Nhits->Clone();
+  h1->SetLineColor(4);
+  hs->Add(h1);
+  h1 = (TH1D*) hEndcap_Nhits->Clone();
+  h1->SetLineColor(2);
+  hs->Add(h1);
+  h1 = (TH1D*) hHits_Nhits->Clone();
+  h1->SetLineColor(2);
+  hs->Add(h1);
+  c->BuildLegend();
+  c->SaveAs("results/tof/rec_multiple_tracks_nhits.png");
+  c->SaveAs("results/tof/rec_multiple_tracks_nhits.pdf");
+
+  // -----------------------------------------------
+  c  = new TCanvas();
+  h_nTracks_vs_theta->DrawCopy("colz");
+  c->SaveAs("results/tof/rec_multiple_tracks_nTracks_vs_theta.png");
+  c->SaveAs("results/tof/rec_multiple_tracks_nTracks_vs_theta.pdf");
+  return 0;
+}
-- 
GitLab