From c02ddfec83326a620fc341c68db42d8d4f9f08f1 Mon Sep 17 00:00:00 2001
From: Whitney Armstrong <warmstrong@anl.gov>
Date: Sat, 7 Aug 2021 01:10:43 +0000
Subject: [PATCH] Modified tracking

---
 benchmarks/tracking/central_electrons.sh      |  13 +-
 .../options/tracker_reconstruction.py         | 124 ++++++------------
 .../scripts/gen_central_electrons.cxx         |   4 +-
 .../scripts/hits_central_electrons.cxx        |  30 ++++-
 .../scripts/rec_central_electrons.cxx         |  79 ++++++-----
 options/tracker_reconstruction.py             |  10 +-
 6 files changed, 124 insertions(+), 136 deletions(-)

diff --git a/benchmarks/tracking/central_electrons.sh b/benchmarks/tracking/central_electrons.sh
index 31d18ca0..58fe7c13 100644
--- a/benchmarks/tracking/central_electrons.sh
+++ b/benchmarks/tracking/central_electrons.sh
@@ -54,6 +54,7 @@ export DETECTOR_PATH=${DETECTOR_PATH}
 if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
   export JUGGLER_N_EVENTS=100
 fi
+export JUGGLER_N_EVENTS=$(expr ${JUGGLER_N_EVENTS} \* 1)
 
 
 export JUGGLER_FILE_NAME_TAG="central_electrons"
@@ -96,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/tracker_reconstruction.py
   if [[ "$?" -ne "0" ]] ; then
     echo "ERROR running juggler"
     exit 1
@@ -111,11 +112,11 @@ if [[ "$?" -ne "0" ]] ; then
   exit 1
 fi
 
-root -b -q "benchmarks/tracking/scripts/hits_central_electrons.cxx(\"${JUGGLER_SIM_FILE}\")"
-if [[ "$?" -ne "0" ]] ; then
-  echo "ERROR running root script"
-  exit 1
-fi
+#root -b -q "benchmarks/tracking/scripts/hits_central_electrons.cxx(\"${JUGGLER_SIM_FILE}\")"
+#if [[ "$?" -ne "0" ]] ; then
+#  echo "ERROR running root script"
+#  exit 1
+#fi
 
 root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
 if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then 
diff --git a/benchmarks/tracking/options/tracker_reconstruction.py b/benchmarks/tracking/options/tracker_reconstruction.py
index 94ffd88e..bcb4755f 100644
--- a/benchmarks/tracking/options/tracker_reconstruction.py
+++ b/benchmarks/tracking/options/tracker_reconstruction.py
@@ -29,6 +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__TrackerSourceLinker as TrackerSourceLinker
 from Configurables import Jug__Reco__TrackerSourcesLinker as TrackerSourcesLinker
@@ -44,143 +45,105 @@ from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction
 from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 
 
+algorithms = [ ]
 
 podioinput = PodioInput("PodioReader", 
-                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","EcalBarrelHits"])#, OutputLevel=DEBUG)
+                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits"])#, OutputLevel=DEBUG)
+algorithms.append( podioinput )
 
 ## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
 copier = MCCopier("MCCopier", 
         inputCollection="mcparticles", 
         outputCollection="mcparticles2") 
+algorithms.append( copier )
+
 trkcopier = TrkCopier("TrkCopier", 
         inputCollection="TrackerBarrelHits", 
         outputCollection="TrackerBarrelHits2") 
-
-ecal_digi = EMCalorimeterDigi("ecal_digi", 
-        inputHitCollection="EcalBarrelHits", 
-        outputHitCollection="RawEcalBarrelHits")
+algorithms.append( trkcopier )
 
 trk_b_digi = TrackerDigi("trk_b_digi", 
         inputHitCollection="TrackerBarrelHits",
         outputHitCollection="TrackerBarrelRawHits",
         timeResolution=8)
+algorithms.append( trk_b_digi )
 trk_ec_digi = TrackerDigi("trk_ec_digi", 
         inputHitCollection="TrackerEndcapHits",
         outputHitCollection="TrackerEndcapRawHits",
         timeResolution=8)
+algorithms.append( trk_ec_digi )
 
 vtx_b_digi = TrackerDigi("vtx_b_digi", 
         inputHitCollection="VertexBarrelHits",
         outputHitCollection="VertexBarrelRawHits",
         timeResolution=8)
+algorithms.append( vtx_b_digi )
 
 vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
         inputHitCollection="VertexEndcapHits",
         outputHitCollection="VertexEndcapRawHits",
         timeResolution=8)
-
-ecal_reco = EMCalReconstruction("ecal_reco", 
-        inputHitCollection="RawEcalBarrelHits", 
-        outputHitCollection="RecEcalBarrelHits",
-        minModuleEdep=0.0*units.MeV,
-        OutputLevel=WARNING)
-
-#simple_cluster = SimpleClustering("simple_cluster", 
-#        inputHitCollection="RecEcalBarrelHits", 
-#        outputClusters="SimpleClusters",
-#        minModuleEdep=1.0*units.MeV,
-#        maxDistance=50.0*units.cm,
-#        OutputLevel=WARNING)
+algorithms.append( vtx_ec_digi )
 
 # Tracker and vertex reconstruction
 trk_b_reco = TrackerHitReconstruction("trk_b_reco",
         inputHitCollection = trk_b_digi.outputHitCollection,
         outputHitCollection="TrackerBarrelRecHits")
+algorithms.append( trk_b_reco )
 
 trk_ec_reco = TrackerHitReconstruction("trk_ec_reco",
         inputHitCollection = trk_ec_digi.outputHitCollection,
         outputHitCollection="TrackerEndcapRecHits")
+algorithms.append( trk_ec_reco )
 
 vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
         inputHitCollection = vtx_b_digi.outputHitCollection,
         outputHitCollection="VertexBarrelRecHits")
+algorithms.append( vtx_b_reco )
 
 vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
         inputHitCollection = vtx_ec_digi.outputHitCollection,
         outputHitCollection="VertexEndcapRecHits")
+algorithms.append( vtx_ec_reco )
+
+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")
+algorithms.append( trk_hit_col )
 
 # Hit Source linker 
 sourcelinker = TrackerSourceLinker("sourcelinker",
-        inputHitCollection=trk_b_reco.outputHitCollection,
-        outputSourceLinks="BarrelTrackSourceLinks",
-        outputMeasurements="BarrelTrackMeasurements",
+        inputHitCollection=trk_hit_col.outputHitCollection,
+        outputSourceLinks="TrackSourceLinks",
+        outputMeasurements="TrackMeasurements",
         OutputLevel=DEBUG)
-
-#trk_hits_srclnkr = TrackerSourcesLinker("trk_srcslnkr",
-#        ITrackerBarrelHits = vtx_b_reco.outputHitCollection,
-#        ITrackerEndcapHits = vtx_ec_reco.outputHitCollection,
-#        OTrackerBarrelHits = trk_b_reco.outputHitCollection,
-#        OTrackerEndcapHits = trk_ec_reco.outputHitCollection,
-#        outputSourceLinks="TrackerMeasurements",
-#        OutputLevel=DEBUG)
+algorithms.append( sourcelinker )
 
 ## Track param init
 truth_trk_init = TrackParamTruthInit("truth_trk_init",
         inputMCParticles="mcparticles",
-        outputInitialTrackParameters="InitTrackParams",
-        OutputLevel=DEBUG)
-
-#clust_trk_init = TrackParamClusterInit("clust_trk_init",
-#        inputClusters="SimpleClusters",
-#        outputInitialTrackParameters="InitTrackParamsFromClusters",
-#        OutputLevel=DEBUG)
-
-#vtxcluster_trk_init = TrackParamVertexClusterInit("vtxcluster_trk_init",
-#        inputVertexHits="VertexBarrelRecHits",
-#        inputClusters="SimpleClusters",
-#        outputInitialTrackParameters="InitTrackParamsFromVtxClusters",
-#        maxHitRadius=40.0*units.mm,
-#        OutputLevel=WARNING)
+        outputInitialTrackParameters="InitTrackParams")
+        #OutputLevel=DEBUG)
+algorithms.append( truth_trk_init )
 
 # Tracking algorithms
 trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
         inputSourceLinks = sourcelinker.outputSourceLinks,
         inputMeasurements = sourcelinker.outputMeasurements,
         inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
-        outputTrajectories="trajectories",
-        OutputLevel=DEBUG)
+        outputTrajectories="trajectories")
+        #OutputLevel=DEBUG)
+algorithms.append( trk_find_alg )
 
 parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
         inputTrajectories="trajectories",
         outputParticles="ReconstructedParticles",
-        outputTrackParameters="outputTrackParameters",
-        OutputLevel=DEBUG)
-
-#trk_find_alg1 = TrackFindingAlgorithm("trk_find_alg1",
-#        inputSourceLinks = sourcelinker.outputSourceLinks,
-#        inputMeasurements = sourcelinker.outputMeasurements,
-#        inputInitialTrackParameters= "InitTrackParamsFromClusters", 
-#        outputTrajectories="trajectories1",
-#        OutputLevel=DEBUG)
-#
-#parts_from_fit1 = ParticlesFromTrackFit("parts_from_fit1",
-#        inputTrajectories="trajectories1",
-#        outputParticles="ReconstructedParticles1",
-#        outputTrackParameters="outputTrackParameters1",
-#        OutputLevel=DEBUG)
-
-#trk_find_alg2 = TrackFindingAlgorithm("trk_find_alg2",
-#        inputSourceLinks = trk_hits_srclnkr.outputSourceLinks,
-#        inputMeasurements = trk_hits_srclnkr.outputMeasurements,
-#        inputInitialTrackParameters= "InitTrackParamsFromVtxClusters", 
-#        outputTrajectories="trajectories2",
-#        OutputLevel=WARNING)
-#parts_from_fit2 = ParticlesFromTrackFit("parts_from_fit2",
-#        inputTrajectories="trajectories2",
-#        outputParticles="ReconstructedParticles2",
-#        outputTrackParameters="outputTrackParameters2",
-#        OutputLevel=WARNING)
-
+        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
@@ -209,21 +172,10 @@ out.outputCommands = ["keep *",
         "drop outputInitialTrackParameters",
         "drop mcparticles"
         ]
+algorithms.append(out)
 
 ApplicationMgr(
-    TopAlg = [podioinput, 
-              copier, trkcopier,
-              ecal_digi, 
-              trk_b_digi, trk_ec_digi, vtx_b_digi, vtx_ec_digi, 
-              trk_b_reco, trk_ec_reco, vtx_b_reco, vtx_ec_reco, 
-              ecal_reco, 
-              sourcelinker, #trk_hits_srclnkr,
-              truth_trk_init, #clust_trk_init, vtxcluster_trk_init, 
-              trk_find_alg, parts_from_fit,
-              #trk_find_alg1, parts_from_fit1,
-              #trk_find_alg2, parts_from_fit2,
-              out
-              ],
+    TopAlg = algorithms,
     EvtSel = 'NONE',
     EvtMax   = n_events,
     ExtSvc = [podioevent,geo_service],
diff --git a/benchmarks/tracking/scripts/gen_central_electrons.cxx b/benchmarks/tracking/scripts/gen_central_electrons.cxx
index e07d320f..3e9d8c01 100644
--- a/benchmarks/tracking/scripts/gen_central_electrons.cxx
+++ b/benchmarks/tracking/scripts/gen_central_electrons.cxx
@@ -17,8 +17,8 @@ using namespace HepMC3;
 void gen_central_electrons(int n_events = 100, 
                      const char* out_fname = "central_electrons.hepmc")
 {
-  double cos_theta_min = std::cos( 80.0*(M_PI/180.0));
-  double cos_theta_max = std::cos(100.0*(M_PI/180.0));
+  double cos_theta_min = std::cos( 140.0*(M_PI/180.0));
+  double cos_theta_max = std::cos(170.0*(M_PI/180.0));
 
   WriterAscii hepmc_output(out_fname);
   int events_parsed = 0;
diff --git a/benchmarks/tracking/scripts/hits_central_electrons.cxx b/benchmarks/tracking/scripts/hits_central_electrons.cxx
index 6f3402b8..a0d76fd4 100644
--- a/benchmarks/tracking/scripts/hits_central_electrons.cxx
+++ b/benchmarks/tracking/scripts/hits_central_electrons.cxx
@@ -103,11 +103,22 @@ int hits_central_electrons(const char* fname = "sim_central_electrons.root")
                  //.Define("delta_p_over_p1",delta_p_over_p, {"p_track1", "p_thrown"})
                  //.Define("delta_p_over_p2",delta_p_over_p, {"p_track2", "p_thrown"})
                  //.Define("N_VtxBarrelHits",[](std::vector<eic::TrackerHitData> hits) { return hits.size();},{"VertexBarrelRecHits"})
+                 //.Define("N_Hits",       [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"trackingHits"})
                  .Define("N_BarrelHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelHits"})
                  .Define("N_EndcapHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapHits"})
                  ;
 
   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");
@@ -184,7 +195,22 @@ int hits_central_electrons(const char* fname = "sim_central_electrons.root")
 
   c = new TCanvas();
   hBarrel_x_vs_y->DrawCopy("colz");
-  c->SaveAs("results/tracking/hits_central_electrons_xy.png");
-  c->SaveAs("results/tracking/hits_central_electrons_xy.pdf");
+  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/tracking/scripts/rec_central_electrons.cxx b/benchmarks/tracking/scripts/rec_central_electrons.cxx
index 22a9e0ed..89a4521b 100644
--- a/benchmarks/tracking/scripts/rec_central_electrons.cxx
+++ b/benchmarks/tracking/scripts/rec_central_electrons.cxx
@@ -101,43 +101,48 @@ int rec_central_electrons(const char* fname = "topside/rec_central_electrons.roo
                  //.Define("delta_p_over_p1",delta_p_over_p, {"p_track1", "p_thrown"})
                  //.Define("delta_p_over_p2",delta_p_over_p, {"p_track2", "p_thrown"})
                  //.Define("N_VtxBarrelHits",[](std::vector<eic::TrackerHitData> hits) { return hits.size();},{"VertexBarrelRecHits"})
+                 .Define("N_Hits",       [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"trackingHits"})
                  .Define("N_BarrelHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelRecHits"})
                  .Define("N_EndcapHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapRecHits"})
                  ;
 
+  auto h_nTracks_vs_theta = df0.Histo2D({"h_nTracks_vs_theta", "; #theta; N tracks ", 40,0,180,10, 0, 10}, "theta0","nTracks");
+
   auto h_nTracks = df0.Histo1D({"h_nTracks", "; N tracks ", 10, 0, 10}, "nTracks");
   auto h_pTracks = df0.Histo1D({"h_pTracks", "; GeV/c ", 100, 0, 10}, "p_track");
 
   auto h_delta_p0  = df0.Histo1D({"h_delta_p0", "Truth Track Init; GeV/c ",  100, -10, 10}, "delta_p0");
-  //auto h_delta_p1 = df0.Histo1D({"h_delta_p1", "Ecal Cluster Init; GeV/c ", 100, -10, 10}, "delta_p1");
-  //auto h_delta_p2 = df0.Histo1D({"h_delta_p2", "Ecal Cluster, innner vtx hit Init; GeV/c ", 100, -10, 10}, "delta_p2");
 
   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 h_delta_p1_over_p = df0.Histo1D({"h_delta_p1_over_p", "Ecal Cluster Init; delta p/p ", 100, -0.1, 0.1}, "delta_p_over_p1");
-  //auto h_delta_p2_over_p = df0.Histo1D({"h_delta_p2_over_p", "Ecal Cluster, innner vtx hit Init; delta p/p ", 100, -0.1, 0.1}, "delta_p_over_p2");
 
-  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 hNhits_vs_theta = df0.Histo1D({"hNhits_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_Hits");
+  auto hBarrel_N_vs_theta = df0.Histo1D({"hBarrel_N_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_BarrelHits");
+  auto hEndcap_N_vs_theta = df0.Histo1D({"hEndcap_N_vs_theta", "; #theta [deg.]",   40, 0, 180 }, "theta0", "N_EndcapHits");
   //auto hVtxBarrel_N_vs_theta = df0.Histo1D({"hVtxBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_VtxBarrelHits");
 
+  auto hHits_Nhits  = df0.Histo1D({"hHits_Nhits", "; #theta [deg.]",       20, 0, 20 }, "N_Hits");
   auto hBarrel_Nhits  = df0.Histo1D({"hBarrel_Nhits", "; #theta [deg.]",   20, 0, 20 }, "N_BarrelHits");
   auto hEndcap_Nhits  = df0.Histo1D({"hEndcap_Nhits", "; #theta [deg.]",   20, 0, 20 }, "N_EndcapHits");
   //auto hVtxBarrel_Nhits = df0.Histo1D({"hVtxBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 },  "N_VtxBarrelHits");
 
-  auto 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 hHits_Ntheta = df0.Histo1D({"hHits_Ntheta", "; #theta [deg.]",       40, 0, 180 }, "theta0");
+  auto hBarrel_Ntheta = df0.Histo1D({"hBarrel_Ntheta", "; #theta [deg.]",   40, 0, 180 }, "theta0");
+  auto hEndcap_Ntheta = df0.Histo1D({"hEndcap_Ntheta", "; #theta [deg.]",   40, 0, 180 }, "theta0");
   //auto hVtxBarrel_Ntheta = df0.Histo1D({"hVtxBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
 
+  // -----------------------------------------------
   auto c = new TCanvas();
 
   h_nTracks->DrawCopy();
   c->SaveAs("results/tracking/rec_central_electrons_nTracks.png");
   c->SaveAs("results/tracking/rec_central_electrons_nTracks.pdf");
 
+  // -----------------------------------------------
   h_pTracks->DrawCopy();
   c->SaveAs("results/tracking/rec_central_electrons_pTracks.png");
   c->SaveAs("results/tracking/rec_central_electrons_pTracks.pdf");
 
+  // -----------------------------------------------
   c = new TCanvas();
   THStack * hs = new THStack("hs_delta_p","; GeV/c ");
   TH1D* h1 = (TH1D*) h_delta_p0->Clone();
@@ -155,74 +160,76 @@ int rec_central_electrons(const char* fname = "topside/rec_central_electrons.roo
   c->SaveAs("results/tracking/rec_central_electrons_delta_p.png");
   c->SaveAs("results/tracking/rec_central_electrons_delta_p.pdf");
 
+  // -----------------------------------------------
   c  = new TCanvas();
   hs = new THStack("hs_delta_p_over_p","; delta p/p ");
   h1 = (TH1D*) h_delta_p0_over_p->Clone();
   hs->Add(h1);
-  //h1 = (TH1D*) h_delta_p1_over_p->Clone();
-  //h1->SetLineColor(2);
-  //hs->Add(h1);
-  //h1 = (TH1D*) h_delta_p2_over_p->Clone();
-  //h1->SetLineColor(4);
-  //h1->SetFillStyle(3001);
-  //h1->SetFillColor(4);
-  //hs->Add(h1);
   hs->Draw("nostack");
   c->BuildLegend();
   c->SaveAs("results/tracking/rec_central_electrons_delta_p_over_p.png");
   c->SaveAs("results/tracking/rec_central_electrons_delta_p_over_p.pdf");
 
+  // -----------------------------------------------
   c  = new TCanvas();
   hs = new THStack("n_hits","; #theta  ");
+
   h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
   auto h2 = (TH1D*) hBarrel_Ntheta->Clone();
+  h1->SetLineColor(4);
   h1->Divide(h2);
   hs->Add(h1);
+
   h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
   h2 = (TH1D*) hEndcap_Ntheta->Clone();
   h1->Divide(h2);
   h1->SetLineColor(2);
   hs->Add(h1);
-  //h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
-  //h1->SetLineColor(4);
-  //h1->SetFillStyle(3001);
-  //h1->SetFillColor(4);
-  //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/tracking/rec_central_electrons_n_hits_vs_theta.png");
   c->SaveAs("results/tracking/rec_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);
+  h2->SetLineColor(4);
   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);
+
+  h2 = (TH1D*) hHits_Ntheta->Clone();
+  h2->SetLineColor(1);
+  hs->Add(h2);
+
   hs->Draw("nostack hist");
   c->BuildLegend();
   c->SaveAs("results/tracking/rec_central_electrons_theta.png");
   c->SaveAs("results/tracking/rec_central_electrons_theta.pdf");
 
+  // -----------------------------------------------
   c  = new TCanvas();
   hs = new THStack("hits","; hits  ");
   h1 = (TH1D*) hBarrel_Nhits->Clone();
+  h1->SetLineColor(4);
   hs->Add(h1);
   h1 = (TH1D*) hEndcap_Nhits->Clone();
   h1->SetLineColor(2);
-  h2->SetLineColor(2);
-  hs->Add(h2);
+  hs->Add(h1);
+  h1 = (TH1D*) hHits_Nhits->Clone();
+  h1->SetLineColor(2);
+  hs->Add(h1);
   //h1 = (TH1D*) hVtxBarrel_Nhits->Clone();
   //h1->SetLineColor(4);
   //h1->SetFillStyle(3001);
@@ -233,6 +240,10 @@ int rec_central_electrons(const char* fname = "topside/rec_central_electrons.roo
   c->SaveAs("results/tracking/rec_central_electrons_nhits.png");
   c->SaveAs("results/tracking/rec_central_electrons_nhits.pdf");
 
-
+  // -----------------------------------------------
+  c  = new TCanvas();
+  h_nTracks_vs_theta->DrawCopy("colz");
+  c->SaveAs("results/tracking/rec_central_electrons_nTracks_vs_theta.png");
+  c->SaveAs("results/tracking/rec_central_electrons_nTracks_vs_theta.pdf");
   return 0;
 }
diff --git a/options/tracker_reconstruction.py b/options/tracker_reconstruction.py
index 19ca4255..a9f0347d 100644
--- a/options/tracker_reconstruction.py
+++ b/options/tracker_reconstruction.py
@@ -19,7 +19,7 @@ if "JUGGLER_DETECTOR_PATH" in os.environ :
 
 geo_service  = GeoSvc("GeoSvc",
         detectors=["{}/{}.xml".format(detector_path, detector_name)])
-podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=DEBUG)
+podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file])#, OutputLevel=DEBUG)
 
 from Configurables import PodioInput
 from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
@@ -73,13 +73,11 @@ ecal_digi = EMCalorimeterDigi("ecal_digi",
 ufsd_digi = UFSDTrackerDigi("ufsd_digi", 
         inputHitCollection="SiTrackerBarrelHits",
         outputHitCollection="SiTrackerBarrelRawHits",
-        timeResolution=8,
-        OutputLevel=DEBUG)
+        timeResolution=8)
 ufsd_digi2 = UFSDTrackerDigi("ufsd_digi2", 
         inputHitCollection="SiTrackerEndcapHits",
         outputHitCollection="SiTrackerEndcapRawHits",
-        timeResolution=8,
-        OutputLevel=DEBUG)
+        timeResolution=8)
 
 #vtx_digi = UFSDTrackerDigi("vtx_digi", 
 #        inputHitCollection="SiVertexBarrelHits",
@@ -232,7 +230,7 @@ ApplicationMgr(
     EvtSel = 'NONE',
     EvtMax   = n_events,
     ExtSvc = [podioevent,geo_service],
-    OutputLevel=DEBUG
+    OutputLevel=INFO
  )
 
 
-- 
GitLab