benchmarks/imaging_ecal/scripts/epi_separation.py

@@ -42,7 +42,7 @@ if __name__ == '__main__':
     parser.add_argument('pifile', type=str, help='path to root file (pions)')
     parser.add_argument('--prof', type=str, default='profile.csv', help='path to electron profile')
     parser.add_argument('--plot-dir', type=str, default='./plots', dest='outdir', help='output directory')
-    parser.add_argument('-b', '--branch-name', type=str, default='EcalBarrelClustersLayers', dest='branch',
+    parser.add_argument('-b', '--branch-name', type=str, default='EcalBarrelImagingClustersLayers', dest='branch',
                         help='branch name in the root file (outputLayerCollection from ImagingClusterReco)')
     parser.add_argument('-m', '--macros', type=str, default='rootlogon.C', dest='macros',
                         help='root macros to load (accept multiple paths separated by \",\")')

benchmarks/imaging_ecal/scripts/utils.py

@@ -3,6 +3,10 @@
 
     Author: Chao Peng (ANL)
     Date: 04/30/2021
+
+    Added all mc particle info to obtain decaying particles
+    Author: Jihee Kim (ANL)
+    Date: 08/06/2021
 '''
 
 import os
@@ -51,7 +55,7 @@ def get_mcp_data(path, evnums=None, branch='mcparticles2'):
         events.GetEntry(iev)
         # extract full mc particle data
         for part in getattr(events, branch):
-            dbuf[idb] = (iev, part.psx, part.psy, part.psz, part.pdgID, part.status)
+            dbuf[idb] = (iev, part.ps.x, part.ps.y, part.ps.z, part.pdgID, part.status)
             idb += 1
     return pd.DataFrame(data=dbuf[:idb], columns=['event', 'px', 'py', 'pz', 'pid', 'status'])
 
@@ -75,13 +79,38 @@ def get_mcp_simple(path, evnums=None, branch='mcparticles2'):
         events.GetEntry(iev)
         # extract full mc particle data
         part = getattr(events, branch)[2]
-        dbuf[idb] = (iev, part.psx, part.psy, part.psz, part.pdgID, part.status)
+        dbuf[idb] = (iev, part.ps.x, part.ps.y, part.ps.z, part.pdgID, part.status)
         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='mcparticles2'):
+    f = ROOT.TFile(path)
+    events = f.events
+    if evnums is None:
+        evnums = np.arange(events.GetEntries())
+    elif isinstance(evnums, int):
+        evnums = [evnums]
+
+    dbuf = np.zeros(shape=(2000*len(evnums), 10))
+    idb = 0
+    for iev in evnums:
+        if iev >= events.GetEntries():
+            print('Error: event {:d} is out of range (0 - {:d})'.format(iev, events.GetEntries() - 1))
+            continue
+
+        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)
+            idb += 1
+    
+    return pd.DataFrame(data=dbuf[:idb], columns=['event', 'px', 'py', 'pz', 'pid', 'status', 'g4Parent', 'vex', 'vey', 'vez'])
 
 # read hits data from root file
-def get_hits_data(path, evnums=None, branch='RecoEcalBarrelHits'):
+def get_hits_data(path, evnums=None, branch='RecoEcalBarreImaginglHits'):
     f = ROOT.TFile(path)
     events = f.events
     if evnums is None:
@@ -107,7 +136,7 @@ def get_hits_data(path, evnums=None, branch='RecoEcalBarrelHits'):
 
 
 # read layers data from root file
-def get_layers_data(path, evnums=None, branch="EcalBarrelClustersLayers"):
+def get_layers_data(path, evnums=None, branch="EcalBarrelImagingClustersLayers"):
     f = ROOT.TFile(path)
     events = f.events
     if evnums is None:
@@ -134,7 +163,7 @@ def get_layers_data(path, evnums=None, branch="EcalBarrelClustersLayers"):
 
 
 # read clusters data from root file
-def get_clusters_data(path, evnums=None, branch='EcalBarrelClustersReco'):
+def get_clusters_data(path, evnums=None, branch='EcalBarrelImagingClustersReco'):
     f = ROOT.TFile(path)
     events = f.events
     if evnums is None:

benchmarks/imaging_shower_ML/scripts/prepare_tf_dataset.py

@@ -49,7 +49,7 @@ def get_mcp_simple(path, evnums=None, branch='mcparticles2'):
         events.GetEntry(iev)
         # extract full mc particle data
         part = getattr(events, branch)[2]
-        dbuf[idb] = (iev, part.psx, part.psy, part.psz, part.pdgID, part.status)
+        dbuf[idb] = (iev, part.ps.x, part.ps.y, part.ps.z, part.pdgID, part.status)
         idb += 1
     return pd.DataFrame(data=dbuf[:idb], columns=['event', 'px', 'py', 'pz', 'pid', 'status'])
 

benchmarks/track_finding/config.yml

+track_finding:multiple_tracks:
+  extends: .rec_benchmark
+  stage: run
+  timeout: 24 hours
+  script:
+    - bash benchmarks/track_finding/multiple_tracks.sh

benchmarks/track_finding/multiple_tracks.sh

+#!/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_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/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
+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
+  gaudirun.py benchmarks/tracking/options/track_reconstruction.py
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running juggler"
+    exit 1
+  fi
+fi
+
+mkdir -p results/track_finding
+
+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
+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
+

benchmarks/tracking/options/tracker_reconstruction.py → benchmarks/track_finding/options/track_reconstruction.py

@@ -29,7 +29,7 @@ from Configurables import Jug__Digi__UFSDTrackerDigi as TrackerDigi
 from Configurables import Jug__Digi__EMCalorimeterDigi as EMCalorimeterDigi
 
 from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerHitReconstruction
-from Configurables import Jug__Reco__TrackingHitsCollector as TrackingHitsCollector
+from Configurables import Jug__Reco__TrackingHitsCollector2 as TrackingHitsCollector
 
 from Configurables import Jug__Reco__TrackerSourceLinker as TrackerSourceLinker
 from Configurables import Jug__Reco__TrackerSourcesLinker as TrackerSourcesLinker
@@ -48,7 +48,7 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 algorithms = [ ]
 
 podioinput = PodioInput("PodioReader", 
-                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits"])#, OutputLevel=DEBUG)
+                        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.
@@ -85,6 +85,12 @@ vtx_ec_digi = TrackerDigi("vtx_ec_digi",
         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,
@@ -106,17 +112,25 @@ vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
         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",
-        trackerBarrelHits=trk_b_reco.outputHitCollection,
-        trackerEndcapHits=trk_ec_reco.outputHitCollection,
-        vertexBarrelHits=vtx_b_reco.outputHitCollection,
-        vertexEndcapHits=vtx_ec_reco.outputHitCollection,
-        outputHitCollection="trackingHits")
+        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.outputHitCollection,
+        inputHitCollection=trk_hit_col.trackingHits,
         outputSourceLinks="TrackSourceLinks",
         outputMeasurements="TrackMeasurements",
         OutputLevel=DEBUG)

benchmarks/track_finding/scripts/gen_central_particles.cxx

+#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 particles in the central region.
+ */
+void gen_central_particles(int n_events = 100, 
+                     const char* out_fname = "central_electrons.hepmc")
+{
+  double cos_theta_min = std::cos( 10.0*(M_PI/180.0));
+  double cos_theta_max = std::cos(170.0*(M_PI/180.0));
+
+  const auto& partmap  = common_bench::particleMap;
+
+  WriterAscii hepmc_output(out_fname);
+  GenEvent evt(Units::GEV, Units::MM);
+
+  int events_parsed = 0;
+
+  // 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, partmap.at(2212).mass), 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 + std::pow(partmap.at(11).mass,2))),
+        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;
+}

benchmarks/track_finding/scripts/gen_multiple_tracks.cxx

+#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 = 2)
+{
+  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     = 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))),
+                                                        ((ip % 2 == 0) ? 11 : -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;
+}

benchmarks/track_finding/scripts/hits_central_electrons.cxx

+#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/TrackerHitCollection.h"
+#include "dd4pod/TrackerHitData.h"
+#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].E());
+  }
+  return result;
+};
+auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
+  std::vector<double> result;
+  for (size_t i = 0; i < in.size(); ++i) {
+   result.push_back(in[i].Theta()*180/M_PI);
+  }
+  return result;
+};
+auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+  std::vector<ROOT::Math::PxPyPzMVector> result;
+  ROOT::Math::PxPyPzMVector lv;
+  for (size_t i = 0; i < in.size(); ++i) {
+    lv.SetCoordinates(in[i].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
+    result.push_back(lv);
+  }
+  return result;
+};
+
+auto 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 hits_central_electrons(const char* fname = "sim_central_electrons.root")
+{
+
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame df("events", fname);
+
+  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
+                 .Define("thrownParticles", "mcparticles[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<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"})
+                 ;
+
+  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 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 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");
+  auto hEndcap_N_vs_theta = df0.Histo1D({"hEndcap_N_vs_theta", "; #theta [deg.]",   20, 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 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 hBarrel_Ntheta = df0.Histo1D({"hBarrel_Ntheta", "; #theta [deg.]",   20, 0, 180 }, "theta0");
+  auto hEndcap_Ntheta = df0.Histo1D({"hEndcap_Ntheta", "; #theta [deg.]",   20, 0, 180 }, "theta0");
+  //auto hVtxBarrel_Ntheta = df0.Histo1D({"hVtxBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
+
+  auto c = new TCanvas();
+  auto hs = new THStack("n_hits","; #theta  ");
+  auto h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
+  auto h2 = (TH1D*) hBarrel_Ntheta->Clone();
+  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*) hVtxBarrel_vs_theta->Clone();
+  //h1->SetLineColor(4);
+  //h1->SetFillStyle(3001);
+  //h1->SetFillColor(4);
+  //hs->Add(h1);
+  hs->Draw("nostack, hist");
+  c->BuildLegend();
+  c->SaveAs("results/tracking/hits_central_electrons_n_hits_vs_theta.png");
+  c->SaveAs("results/tracking/hits_central_electrons_n_hits_vs_theta.pdf");
+
+  c  = new TCanvas();
+  hs = new THStack("theta","; #theta  ");
+  h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
+  h2 = (TH1D*) hBarrel_Ntheta->Clone();
+  //h1->Divide(h2);
+  hs->Add(h2);
+  h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
+  h2 = (TH1D*) hEndcap_Ntheta->Clone();
+  //h1->Divide(h2);
+  h1->SetLineColor(2);
+  h2->SetLineColor(2);
+  hs->Add(h2);
+  //h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
+  //h1->SetLineColor(4);
+  //h1->SetFillStyle(3001);
+  //h1->SetFillColor(4);
+  //hs->Add(h1);
+  hs->Draw("nostack hist");
+  c->BuildLegend();
+  c->SaveAs("results/tracking/hits_central_electrons_theta.png");
+  c->SaveAs("results/tracking/hits_central_electrons_theta.pdf");
+
+  c  = new TCanvas();
+  hs = new THStack("hits","; hits  ");
+  h1 = (TH1D*) hBarrel_Nhits->Clone();
+  hs->Add(h1);
+  h1 = (TH1D*) hEndcap_Nhits->Clone();
+  h1->SetLineColor(2);
+  h2->SetLineColor(2);
+  hs->Add(h2);
+  //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/tracking/hits_central_electrons_nhits.png");
+  c->SaveAs("results/tracking/hits_central_electrons_nhits.pdf");
+
+  c = new TCanvas();
+  hBarrel_x_vs_y->DrawCopy("colz");
+  c->SaveAs("results/tracking/hits_central_electrons_trkBarrel_xy.png");
+  c->SaveAs("results/tracking/hits_central_electrons_trkBarrel_xy.pdf");
+
+  c = new TCanvas();
+  hBarrel_x_vs_y->DrawCopy("colz");
+  hEndcap_x_vs_y->DrawCopy("colz same");
+  //hvtxBarrel_x_vs_y->DrawCopy("colz same");
+  //hvtxEndcap_x_vs_y->DrawCopy("colz same");
+
+  c->SaveAs("results/tracking/hits_central_electrons_Hits_xy.png");
+  c->SaveAs("results/tracking/hits_central_electrons_Hits_xy.pdf");
+
+  //hAllHits_x_vs_z->DrawCopy("colz");
+  hBarrel_x_vs_z->DrawCopy("colz");
+  hEndcap_x_vs_z->DrawCopy("colz same");
+  c->SaveAs("results/tracking/hits_central_electrons_Hits_xz.png");
+  c->SaveAs("results/tracking/hits_central_electrons_Hits_xz.pdf");
+  return 0;
+}

benchmarks/track_finding/scripts/hits_central_pions.cxx

+#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/TrackerHitCollection.h"
+#include "dd4pod/TrackerHitData.h"
+#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].E());
+  }
+  return result;
+};
+auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
+  std::vector<double> result;
+  for (size_t i = 0; i < in.size(); ++i) {
+   result.push_back(in[i].Theta()*180/M_PI);
+  }
+  return result;
+};
+auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+  std::vector<ROOT::Math::PxPyPzMVector> result;
+  ROOT::Math::PxPyPzMVector lv;
+  for (size_t i = 0; i < in.size(); ++i) {
+    lv.SetCoordinates(in[i].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
+    result.push_back(lv);
+  }
+  return result;
+};
+
+auto 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 hits_central_pions(const char* fname = "sim_central_pions.root")
+{
+
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame df("events", fname);
+
+  auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
+                 .Define("thrownParticles", "mcparticles[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_BarrelHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelHits"})
+                 .Define("N_EndcapHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapHits"})
+                 ;
+
+  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_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");
+  //auto hVtxBarrel_N_vs_theta = df0.Histo1D({"hVtxBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_VtxBarrelHits");
+
+  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 hBarrel_Ntheta = df0.Histo1D({"hBarrel_Ntheta", "; #theta [deg.]",   20, 0, 180 }, "theta0");
+  auto hEndcap_Ntheta = df0.Histo1D({"hEndcap_Ntheta", "; #theta [deg.]",   20, 0, 180 }, "theta0");
+  //auto hVtxBarrel_Ntheta = df0.Histo1D({"hVtxBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
+
+  auto c = new TCanvas();
+  auto hs = new THStack("n_hits","; #theta  ");
+  auto h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
+  auto h2 = (TH1D*) hBarrel_Ntheta->Clone();
+  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*) hVtxBarrel_vs_theta->Clone();
+  //h1->SetLineColor(4);
+  //h1->SetFillStyle(3001);
+  //h1->SetFillColor(4);
+  //hs->Add(h1);
+  hs->Draw("nostack, hist");
+  c->BuildLegend();
+  c->SaveAs("results/tracking/hits_central_pions_n_hits_vs_theta.png");
+  c->SaveAs("results/tracking/hits_central_pions_n_hits_vs_theta.pdf");
+
+  c  = new TCanvas();
+  hs = new THStack("theta","; #theta  ");
+  h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
+  h2 = (TH1D*) hBarrel_Ntheta->Clone();
+  //h1->Divide(h2);
+  hs->Add(h2);
+  h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
+  h2 = (TH1D*) hEndcap_Ntheta->Clone();
+  //h1->Divide(h2);
+  h1->SetLineColor(2);
+  h2->SetLineColor(2);
+  hs->Add(h2);
+  //h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
+  //h1->SetLineColor(4);
+  //h1->SetFillStyle(3001);
+  //h1->SetFillColor(4);
+  //hs->Add(h1);
+  hs->Draw("nostack hist");
+  c->BuildLegend();
+  c->SaveAs("results/tracking/hits_central_pions_theta.png");
+  c->SaveAs("results/tracking/hits_central_pions_theta.pdf");
+
+  c  = new TCanvas();
+  hs = new THStack("hits","; hits  ");
+  h1 = (TH1D*) hBarrel_Nhits->Clone();
+  hs->Add(h1);
+  h1 = (TH1D*) hEndcap_Nhits->Clone();
+  h1->SetLineColor(2);
+  h2->SetLineColor(2);
+  hs->Add(h2);
+  //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/tracking/hits_central_pions_nhits.png");
+  c->SaveAs("results/tracking/hits_central_pions_nhits.pdf");
+
+  c = new TCanvas();
+  hBarrel_x_vs_y->DrawCopy("colz");
+  c->SaveAs("results/tracking/hits_central_pions_xy.png");
+  c->SaveAs("results/tracking/hits_central_pions_xy.pdf");
+  return 0;
+}

benchmarks/track_finding/scripts/rec_multiple_tracks.cxx

+#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_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("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"})
+                 ;
+
+  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 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_nTracks->DrawCopy();
+  c->SaveAs("results/track_finding/rec_multiple_tracks_nTracks.png");
+  c->SaveAs("results/track_finding/rec_multiple_tracks_nTracks.pdf");
+
+  // -----------------------------------------------
+  h_pTracks->DrawCopy();
+  c->SaveAs("results/track_finding/rec_multiple_tracks_pTracks.png");
+  c->SaveAs("results/track_finding/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/track_finding/rec_multiple_tracks_delta_p.png");
+  c->SaveAs("results/track_finding/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/track_finding/rec_multiple_tracks_delta_p_over_p.png");
+  c->SaveAs("results/track_finding/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/track_finding/rec_multiple_tracks_n_hits_vs_theta.png");
+  c->SaveAs("results/track_finding/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/track_finding/rec_multiple_tracks_theta.png");
+  c->SaveAs("results/track_finding/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/track_finding/rec_multiple_tracks_nhits.png");
+  c->SaveAs("results/track_finding/rec_multiple_tracks_nhits.pdf");
+
+  // -----------------------------------------------
+  c  = new TCanvas();
+  h_nTracks_vs_theta->DrawCopy("colz");
+  c->SaveAs("results/track_finding/rec_multiple_tracks_nTracks_vs_theta.png");
+  c->SaveAs("results/track_finding/rec_multiple_tracks_nTracks_vs_theta.pdf");
+  return 0;
+}

benchmarks/tracking/central_electrons.sh

@@ -97,7 +97,7 @@ rootls -t ${JUGGLER_SIM_FILE}
 if [[ -z "${ANALYSIS_ONLY}" ]] ;
 then
   # Need to figure out how to pass file name to juggler from the commandline
-  gaudirun.py benchmarks/tracking/options/tracker_reconstruction.py
+  xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv gaudirun.py benchmarks/tracking/options/track_reconstruction.py
   if [[ "$?" -ne "0" ]] ; then
     echo "ERROR running juggler"
     exit 1

benchmarks/tracking/central_pions.sh

@@ -96,7 +96,7 @@ 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/tracking/options/tracker_reconstruction.py
+  xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv gaudirun.py benchmarks/tracking/options/track_reconstruction.py
   if [[ "$?" -ne "0" ]] ; then
     echo "ERROR running juggler"
     exit 1

benchmarks/tracking/config.yml

@@ -26,25 +26,11 @@ tracking_truth_init_electrons:
   stage: run
   timeout: 24 hours
   script:
-    - python benchmarks/tracking/run_tracking_benchmarks.py --options truth_seeded_tracking.py --nametag=truth_electron --particle=electron --etamin=-3 --etamax=3 -n 100
+    - python benchmarks/tracking/run_tracking_benchmarks.py --nametag=truth_electron --particle=electron --etamin=-3 --etamax=3 -n 100
 
 tracking_truth_init_pions:
   extends: .rec_benchmark
   stage: run
   timeout: 24 hours
   script:
-    - python benchmarks/tracking/run_tracking_benchmarks.py --options truth_seeded_tracking.py --nametag=truth_pion --particle=pion+ --etamin=-3 --etamax=3 -n 100
-
-      #tracking_imclust_init_electrons:
-      #  extends: .rec_benchmark
-      #  stage: run
-      #  timeout: 24 hours
-      #  script:
-      #    - python benchmarks/tracking/run_tracking_benchmarks.py --options imagingcluster_seeded_tracking.py --nametag=imclust_electron --particle=electron --etamin=-3 --etamax=3 -n 100
-      #
-      #tracking_imclust_init_pions:
-      #  extends: .rec_benchmark
-      #  stage: run
-      #  timeout: 24 hours
-      #  script:
-      #    - python benchmarks/tracking/run_tracking_benchmarks.py --options imagingcluster_seeded_tracking.py --nametag=imclust_pion --particle=pion+ --etamin=-3 --etamax=3 -n 100
+    - python benchmarks/tracking/run_tracking_benchmarks.py --nametag=truth_pion --particle=pion+ --etamin=-3 --etamax=3 -n 100

benchmarks/tracking/multiple_tracks.sh

@@ -97,7 +97,7 @@ rootls -t ${JUGGLER_SIM_FILE}
 if [[ -z "${ANALYSIS_ONLY}" ]] ;
 then
   # Need to figure out how to pass file name to juggler from the commandline
-  gaudirun.py benchmarks/tracking/options/tracker_reconstruction.py
+  gaudirun.py benchmarks/tracking/options/track_reconstruction.py
   if [[ "$?" -ne "0" ]] ; then
     echo "ERROR running juggler"
     exit 1

benchmarks/tracking/options/imagingcluster_seeded_tracking.py

-from Gaudi.Configuration import *
-
-from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
-from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
-
-detector_name = str(os.environ.get("JUGGLER_DETECTOR", "athena"))
-detector_path = str(os.environ.get("JUGGLER_DETECTOR_PATH", "."))
-compact_path = str(os.environ.get("JUGGLER_COMPACT_PATH", "{}.xml".format(os.path.join(detector_path, detector_name))))
-
-# CAL reconstruction
-# get sampling fractions from system environment variable, 1.0 by default
-cb_ecal_sf = float(os.environ.get("CB_ECAL_SAMP_FRAC", 0.01324))
-
-# todo add checks
-input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
-output_rec = str(os.environ["JUGGLER_REC_FILE"])
-n_events = int(os.environ["JUGGLER_N_EVENTS"])
-
-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__Reco__TrackerHitReconstruction as TrackerReco
-
-from Configurables import Jug__Reco__TrackerSourceLinker as TrackerSourceLinker
-from Configurables import Jug__Reco__TrackerSourcesLinker as TrackerSourcesLinker
-#from Configurables import Jug__Reco__TrackingHitsSourceLinker as TrackingHitsSourceLinker
-from Configurables import Jug__Reco__TrackParamImagingClusterInit as TrackParamImagingClusterInit
-
-from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
-from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
-from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
-
-from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
-from Configurables import Jug__Reco__ImagingPixelReco as ImCalPixelReco
-from Configurables import Jug__Reco__ImagingTopoCluster as ImagingCluster
-from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
-
-
-sim_colls = [
-    "mcparticles",
-    "TrackerEndcapHits",
-    "TrackerBarrelHits",
-    "VertexBarrelHits",
-    "VertexEndcapHits",
-    "EcalBarrelHits"
-]
-podin = PodioInput("PodioReader", collections=sim_colls)
-podout = PodioOutput("out", filename=output_rec)
-
-## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
-mccopier = MCCopier("MCCopier",
-        inputCollection="mcparticles",
-        outputCollection="mcparticles2")
-
-trk_b_digi = TrackerDigi("trk_b_digi",
-        inputHitCollection="TrackerBarrelHits",
-        outputHitCollection="TrackerBarrelRawHits",
-        timeResolution=8)
-trk_ec_digi = TrackerDigi("trk_ec_digi",
-        inputHitCollection="TrackerEndcapHits",
-        outputHitCollection="TrackerEndcapRawHits",
-        timeResolution=8)
-
-vtx_b_digi = TrackerDigi("vtx_b_digi",
-        inputHitCollection="VertexBarrelHits",
-        outputHitCollection="VertexBarrelRawHits",
-        timeResolution=8)
-
-vtx_ec_digi = TrackerDigi("vtx_ec_digi",
-        inputHitCollection="VertexEndcapHits",
-        outputHitCollection="VertexEndcapRawHits",
-        timeResolution=8)
-
-# Tracker and vertex reconstruction
-trk_b_reco = TrackerReco("trk_b_reco",
-        inputHitCollection = trk_b_digi.outputHitCollection,
-        outputHitCollection="TrackerBarrelRecHits")
-
-trk_ec_reco = TrackerReco("trk_ec_reco",
-        inputHitCollection = trk_ec_digi.outputHitCollection,
-        outputHitCollection="TrackerEndcapRecHits")
-
-vtx_b_reco = TrackerReco("vtx_b_reco",
-        inputHitCollection = vtx_b_digi.outputHitCollection,
-        outputHitCollection="VertexBarrelRecHits")
-
-vtx_ec_reco = TrackerReco("vtx_ec_reco",
-        inputHitCollection = vtx_ec_digi.outputHitCollection,
-        outputHitCollection="VertexEndcapRecHits")
-
-sourcelinker = TrackerSourcesLinker("trk_srcslnkr",
-        inputHitCollections=["VertexBarrelRecHits", "TrackerBarrelRecHits"],
-        outputSourceLinks="TrackerSourceLinks",
-        outputMeasurements="TrackerMeasurements",
-        OutputLevel=DEBUG)
-
-# Central Barrel Ecal (Imaging Cal.)
-cb_ecal_daq = dict(
-        dynamicRangeADC=3*MeV,
-        capacityADC=8192,
-        pedestalMean=400,
-        pedestalSigma=20)   # about 6 keV
-
-cb_ecal_digi = CalHitDigi("cb_ecal_digi",
-        inputHitCollection="EcalBarrelHits",
-        outputHitCollection="EcalBarrelHitsDigi",
-        energyResolutions=[0., 0.02, 0.],   # 2% flat resolution
-        **cb_ecal_daq)
-
-cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
-        inputHitCollection=cb_ecal_digi.outputHitCollection,
-        outputHitCollection="EcalBarrelHitsReco",
-        thresholdFactor=3,  # about 20 keV
-        readoutClass="EcalBarrelHits",  # readout class
-        layerField="layer",             # field to get layer id
-        sectorField="module",           # field to get sector id
-        **cb_ecal_daq)
-
-cb_ecal_cl = ImagingCluster("cb_ecal_cl",
-        inputHitCollection=cb_ecal_reco.outputHitCollection,
-        outputProtoClusterCollection="EcalBarrelProtoClusters",
-        localDistXY=[2.*mm, 2*mm],  # same layer
-        layerDistEtaPhi=[10*mrad, 10*mrad],     # adjacent layer
-        neighbourLayersRange=2,                 # id diff for adjacent layer
-        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="EcalBarrelClusters",
-        outputInfoCollection="EcalBarrelClustersInfo",
-        outputLayerCollection="EcalBarrelLayers")
-
-## Track param init
-imclust_trk_init = TrackParamImagingClusterInit("imclust_trk_init",
-        inputClusters="EcalBarrelClusters",
-        outputInitialTrackParameters="InitTrackParams",
-        OutputLevel=DEBUG)
-
-# Tracking algorithms
-trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
-        inputSourceLinks = sourcelinker.outputSourceLinks,
-        inputMeasurements = sourcelinker.outputMeasurements,
-        inputInitialTrackParameters= imclust_trk_init.outputInitialTrackParameters,
-        outputTrajectories="trajectories",
-        OutputLevel=DEBUG)
-
-parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
-        inputTrajectories=trk_find_alg.outputTrajectories,
-        outputParticles="ReconstructedParticles",
-        outputTrackParameters="outputTrackParameters",
-        OutputLevel=DEBUG)
-
-podout.outputCommands = [
-    "keep *",
-    "drop InitTrackParams",
-    "drop trajectories",
-    "drop outputSourceLinks",
-    "drop outputInitialTrackParameters",
-    "drop mcparticles",
-]
-
-ApplicationMgr(
-    TopAlg = [podin, mccopier,
-              trk_b_digi, trk_ec_digi, vtx_b_digi, vtx_ec_digi,
-              trk_b_reco, trk_ec_reco, vtx_b_reco, vtx_ec_reco,
-              cb_ecal_digi, cb_ecal_reco, cb_ecal_cl, cb_ecal_clreco,
-              sourcelinker,
-              imclust_trk_init,
-              trk_find_alg, parts_from_fit,
-              podout
-              ],
-    EvtSel = 'NONE',
-    EvtMax   = n_events,
-    ExtSvc = [podioevent, geo_service],
-    OutputLevel=WARNING
- )
-

benchmarks/tracking/options/truth_seeded_tracking.py → benchmarks/tracking/options/track_reconstruction.py

 from Gaudi.Configuration import *
 
 from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
-from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
+from GaudiKernel import SystemOfUnits as units
 
-detector_name = str(os.environ.get("JUGGLER_DETECTOR", "athena"))
-detector_path = str(os.environ.get("JUGGLER_DETECTOR_PATH", "."))
-compact_path = str(os.environ.get("JUGGLER_COMPACT_PATH", "{}.xml".format(os.path.join(detector_path, detector_name))))
+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_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
-output_rec = str(os.environ["JUGGLER_REC_FILE"])
-n_events = int(os.environ["JUGGLER_N_EVENTS"])
+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=[compact_path], OutputLevel=WARNING)
-podioevent = EICDataSvc("EventDataSvc", inputs=input_sims, OutputLevel=WARNING)
+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__ExampleCaloDigi as ExampleCaloDigi
 from Configurables import Jug__Digi__UFSDTrackerDigi as TrackerDigi
-from Configurables import Jug__Reco__TrackerHitReconstruction as TrackerReco
+from Configurables import Jug__Digi__EMCalorimeterDigi as EMCalorimeterDigi
+
+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__TrackerSourcesLinker as TrackerSourcesLinker
 #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
 
 from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
 from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
-from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction
+
+from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
+
 
+algorithms = [ ]
 
-sim_colls = [
-    "mcparticles",
-    "TrackerEndcapHits",
-    "TrackerBarrelHits",
-    "VertexBarrelHits",
-    "VertexEndcapHits",
-    "EcalBarrelHits"
-]
-podin = PodioInput("PodioReader", collections=sim_colls)
-podout = PodioOutput("out", filename=output_rec)
+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.
-mccopier = MCCopier("MCCopier",
-        inputCollection="mcparticles",
-        outputCollection="mcparticles2")
+copier = MCCopier("MCCopier", 
+        inputCollection="mcparticles", 
+        outputCollection="mcparticles2") 
+algorithms.append( copier )
 
-trk_b_digi = TrackerDigi("trk_b_digi",
+trkcopier = TrkCopier("TrkCopier", 
+        inputCollection="TrackerBarrelHits", 
+        outputCollection="TrackerBarrelHits2") 
+algorithms.append( trkcopier )
+
+trk_b_digi = TrackerDigi("trk_b_digi", 
         inputHitCollection="TrackerBarrelHits",
         outputHitCollection="TrackerBarrelRawHits",
         timeResolution=8)
-trk_ec_digi = TrackerDigi("trk_ec_digi",
+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",
+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",
+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 = TrackerReco("trk_b_reco",
+trk_b_reco = TrackerHitReconstruction("trk_b_reco",
         inputHitCollection = trk_b_digi.outputHitCollection,
         outputHitCollection="TrackerBarrelRecHits")
+algorithms.append( trk_b_reco )
 
-trk_ec_reco = TrackerReco("trk_ec_reco",
+trk_ec_reco = TrackerHitReconstruction("trk_ec_reco",
         inputHitCollection = trk_ec_digi.outputHitCollection,
         outputHitCollection="TrackerEndcapRecHits")
+algorithms.append( trk_ec_reco )
 
-vtx_b_reco = TrackerReco("vtx_b_reco",
+vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
         inputHitCollection = vtx_b_digi.outputHitCollection,
         outputHitCollection="VertexBarrelRecHits")
+algorithms.append( vtx_b_reco )
 
-vtx_ec_reco = TrackerReco("vtx_ec_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 )
 
-sourcelinker = TrackerSourcesLinker("trk_srcslnkr",
-        inputHitCollections=["VertexBarrelRecHits", "TrackerBarrelRecHits"],
-        outputSourceLinks="TrackerSourceLinks",
-        outputMeasurements="TrackerMeasurements",
+# 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)
+        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= truth_trk_init.outputInitialTrackParameters,
-        outputTrajectories="trajectories",
-        OutputLevel=DEBUG)
+        inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
+        outputTrajectories="trajectories")
+        #OutputLevel=DEBUG)
+algorithms.append( trk_find_alg )
 
 parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
-        inputTrajectories=trk_find_alg.outputTrajectories,
+        inputTrajectories="trajectories",
         outputParticles="ReconstructedParticles",
-        outputTrackParameters="outputTrackParameters",
-        OutputLevel=DEBUG)
-
-podout.outputCommands = [
-    "keep *",
-    "drop InitTrackParams",
-    "drop trajectories",
-    "drop outputSourceLinks",
-    "drop outputInitialTrackParameters",
-    "drop mcparticles",
-]
+        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 = [podin, mccopier,
-              trk_b_digi, trk_ec_digi, vtx_b_digi, vtx_ec_digi,
-              trk_b_reco, trk_ec_reco, vtx_b_reco, vtx_ec_reco,
-              sourcelinker,
-              truth_trk_init,
-              trk_find_alg, parts_from_fit,
-              podout
-              ],
+    TopAlg = algorithms,
     EvtSel = 'NONE',
     EvtMax   = n_events,
-    ExtSvc = [podioevent, geo_service],
+    ExtSvc = [podioevent,geo_service],
     OutputLevel=WARNING
  )
 

benchmarks/tracking/run_tracking_benchmarks.py

@@ -28,7 +28,7 @@ parser.add_argument('--etamax', type=float, default=4, help='Maximum pseudorapid
 parser.add_argument('--pmin', type=float, default=5.0, help='Minimum momentum of particles (GeV).')
 parser.add_argument('--pmax', type=float, default=5.0, help='Maximum momentum of particles (GeV).')
 parser.add_argument('--compact', type=str, default=default_compact, help='Path to detector compact file.')
-parser.add_argument('--options', type=str, default='truth_seeded_tracking.py', help='Path to reconstruction options script.')
+parser.add_argument('--options', type=str, default='track_reconstruction.py', help='Path to reconstruction options script.')
 parser.add_argument('--uploadSizeLimit', type=float, default=10, help='Upper limit of file size (Mb) to be uploaded.')
 
 args = parser.parse_args()

benchmarks/tracking/scripts/hits_central_electrons.cxx

@@ -31,7 +31,7 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
   std::vector<float> result;
   for (size_t i = 0; i < in.size(); ++i) {
-    result.push_back(std::sqrt(in[i].psx * in[i].psx + in[i].psy * in[i].psy));
+    result.push_back(std::sqrt(in[i].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y));
   }
   return result;
 }
@@ -54,7 +54,7 @@ auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
   std::vector<ROOT::Math::PxPyPzMVector> result;
   ROOT::Math::PxPyPzMVector lv;
   for (size_t i = 0; i < in.size(); ++i) {
-    lv.SetCoordinates(in[i].psx, in[i].psy, in[i].psz, in[i].mass);
+    lv.SetCoordinates(in[i].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
     result.push_back(lv);
   }
   return result;

benchmarks/tracking/scripts/hits_central_pions.cxx

@@ -31,7 +31,7 @@ auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
 std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
   std::vector<float> result;
   for (size_t i = 0; i < in.size(); ++i) {
-    result.push_back(std::sqrt(in[i].psx * in[i].psx + in[i].psy * in[i].psy));
+    result.push_back(std::sqrt(in[i].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y));
   }
   return result;
 }
@@ -54,7 +54,7 @@ auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
   std::vector<ROOT::Math::PxPyPzMVector> result;
   ROOT::Math::PxPyPzMVector lv;
   for (size_t i = 0; i < in.size(); ++i) {
-    lv.SetCoordinates(in[i].psx, in[i].psy, in[i].psz, in[i].mass);
+    lv.SetCoordinates(in[i].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
     result.push_back(lv);
   }
   return result;