diff --git a/benchmarks/track_finding/options/track_reconstruction.py b/benchmarks/track_finding/options/track_reconstruction.py
index b657a008f9e2cde7fea75385da4b4c287c862ac8..328c0e0ba311a5c469f8891588674595a922b803 100644
--- a/benchmarks/track_finding/options/track_reconstruction.py
+++ b/benchmarks/track_finding/options/track_reconstruction.py
@@ -38,6 +38,7 @@ 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__EMCalReconstruction as EMCalReconstruction
@@ -48,7 +49,7 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 algorithms = [ ]
 
 podioinput = PodioInput("PodioReader", 
-                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits"])#, OutputLevel=DEBUG)
+                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits"])
 algorithms.append( podioinput )
 
 ## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
@@ -128,6 +129,14 @@ trk_hit_col = TrackingHitsCollector("trk_hit_col",
         OutputLevel=DEBUG)
 algorithms.append( trk_hit_col )
 
+# Hit Source linker 
+conformal_find = ConformalXYPeakProtoTracks("conformal_find",
+        inputTrackerHits=trk_hit_col.trackingHits,
+        outputProtoTracks="outputProtoTracks",
+        nProtoTracks="nProtoTracks",
+        OutputLevel=DEBUG)
+algorithms.append( conformal_find )
+
 # Hit Source linker 
 sourcelinker = TrackerSourceLinker("sourcelinker",
         inputHitCollection=trk_hit_col.trackingHits,
diff --git a/benchmarks/track_finding/scripts/gen_central_particles.cxx b/benchmarks/track_finding/scripts/gen_central_particles.cxx
deleted file mode 100644
index fda5133aad4c3f05f18611eb5cbd335ebc5a9bcd..0000000000000000000000000000000000000000
--- a/benchmarks/track_finding/scripts/gen_central_particles.cxx
+++ /dev/null
@@ -1,88 +0,0 @@
-#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;
-}
diff --git a/benchmarks/track_finding/scripts/gen_multiple_tracks.cxx b/benchmarks/track_finding/scripts/gen_multiple_tracks.cxx
index a9c0999f4bd23a844b6fb8822b0db7e019d53142..8596e929d190d7b023f37f3fbbf7da07b88508cb 100644
--- a/benchmarks/track_finding/scripts/gen_multiple_tracks.cxx
+++ b/benchmarks/track_finding/scripts/gen_multiple_tracks.cxx
@@ -19,7 +19,7 @@ using namespace HepMC3;
  */
 void gen_multiple_tracks(int n_events = 100, 
                          const char* out_fname = "multiple_tracks.hepmc",
-                         int n_parts = 2)
+                         int n_parts = 5)
 {
   double cos_theta_min = std::cos( 10.0*(M_PI/180.0));
   double cos_theta_max = std::cos(170.0*(M_PI/180.0));
diff --git a/benchmarks/track_finding/scripts/hits_central_electrons.cxx b/benchmarks/track_finding/scripts/hits_central_electrons.cxx
deleted file mode 100644
index 7e38cc88aea38933b9b0567bfc915198e72e779e..0000000000000000000000000000000000000000
--- a/benchmarks/track_finding/scripts/hits_central_electrons.cxx
+++ /dev/null
@@ -1,216 +0,0 @@
-#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;
-}
diff --git a/benchmarks/track_finding/scripts/hits_central_pions.cxx b/benchmarks/track_finding/scripts/hits_central_pions.cxx
deleted file mode 100644
index 0a70c5e487fa07404f524261ae9a4fe80f4d8b29..0000000000000000000000000000000000000000
--- a/benchmarks/track_finding/scripts/hits_central_pions.cxx
+++ /dev/null
@@ -1,190 +0,0 @@
-#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;
-}
diff --git a/benchmarks/track_finding/scripts/rec_multiple_tracks.cxx b/benchmarks/track_finding/scripts/rec_multiple_tracks.cxx
index f66932cd0790434f0b2619e4604909c9cdfa9c43..fdd02666a18ccb94380fdb9924be1322c59f1159 100644
--- a/benchmarks/track_finding/scripts/rec_multiple_tracks.cxx
+++ b/benchmarks/track_finding/scripts/rec_multiple_tracks.cxx
@@ -87,6 +87,7 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root")
   auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
                  .Define("thrownParticles", "mcparticles2[isThrown]")
                  .Define("thrownP", fourvec, {"thrownParticles"})
+                 .Define("nThrown", "thrownParticles.size()")
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
                  .Define("theta0", "theta_thrown[0]")
@@ -106,6 +107,9 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root")
 
   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}, "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");
@@ -123,6 +127,12 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root")
   // -----------------------------------------------
   auto c = new TCanvas();
 
+  h_nThrown->SetLineColor(2);
+  h_nThrown->DrawCopy();
+  h_nProtoTracks->DrawCopy("same");
+  c->SaveAs("results/track_finding/rec_multiple_tracks_nProtoTracks.png");
+  c->SaveAs("results/track_finding/rec_multiple_tracks_nProtoTracks.pdf");
+
   h_nTracks->DrawCopy();
   c->SaveAs("results/track_finding/rec_multiple_tracks_nTracks.png");
   c->SaveAs("results/track_finding/rec_multiple_tracks_nTracks.pdf");