Adding some helpers on the main tracking script

	modified:   benchmarks/tracking/central_electrons.sh
	new file:   benchmarks/tracking/scripts/hits_central_electrons.cxx
	modified:   options/tracker_reconstruction.py

benchmarks/tracking/central_electrons.sh

 #!/bin/bash
+
+function print_the_help {
+  echo "USAGE: ${0} [--rec-only]  "
+  echo "  OPTIONS: "
+  echo "            --rec-only   only run gaudi reconstruction part"
+  exit 
+}
+
+REC_ONLY=
+ANALYSIS_ONLY=
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    --rec-only)
+      REC_ONLY=1
+      shift # past value
+      ;;
+    --ana-only)
+      ANALYSIS_ONLY=1
+      shift # past value
+      ;;
+    *)    # unknown option
+      #POSITIONAL+=("$1") # save it in an array for later
+      echo "unknown option $1"
+      print_the_help
+      shift # past argument
+      ;;
+  esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+
 ./util/print_env.sh
 
 ## To run the reconstruction, we need the following global variables:
@@ -27,34 +67,40 @@ echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
 echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
 
 
-## generate the input events
-root -b -q "benchmarks/tracking/scripts/gen_central_electrons.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
-if [[ "$?" -ne "0" ]] ; then
-  echo "ERROR running script"
-  exit 1
-fi
+if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
+then
 
+  ## generate the input events
+  root -b -q "benchmarks/tracking/scripts/gen_central_electrons.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running script"
+    exit 1
+  fi
 
-## 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
+  ## run geant4 simulations
+  npsim --runType batch \
+    --part.minimalKineticEnergy 1000*GeV  \
+    -v WARNING \
+    --numberOfEvents ${JUGGLER_N_EVENTS} \
+    --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+    --inputFiles  ${JUGGLER_FILE_NAME_TAG}.hepmc \
+    --outputFile  ${JUGGLER_SIM_FILE}
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running script"
+    exit 1
+  fi
 fi
 
-# 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
-if [[ "$?" -ne "0" ]] ; then
-  echo "ERROR running juggler"
-  exit 1
-fi
 
+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 options/tracker_reconstruction.py
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running juggler"
+    exit 1
+  fi
+fi
 
 mkdir -p results/tracking
 
@@ -64,6 +110,12 @@ 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_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
 if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then 
   # file must be less than 10 MB to upload

benchmarks/tracking/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].psx * in[i].psx + in[i].psy * in[i].psy));
+  }
+  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].psx, in[i].psy, in[i].psz, 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_SiBarrelHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"SiTrackerBarrelHits"})
+                 .Define("N_SiEndcapHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"SiTrackerEndcapHits"})
+                 ;
+
+  auto hSiBarrel_x_vs_y = df0.Histo2D({"hSiBarrel_x_vs_y", "; x ; y ",   100, -900, 900,100, -900, 900 }, "SiTrackerBarrelHits.position.x", "SiTrackerBarrelHits.position.y");
+
+  auto hSiBarrel_N_vs_theta = df0.Histo1D({"hSiBarrel_N_vs_theta", "; #theta [deg.]",   20, 0, 180 }, "theta0", "N_SiBarrelHits");
+  auto hSiEndcap_N_vs_theta = df0.Histo1D({"hSiEndcap_N_vs_theta", "; #theta [deg.]",   20, 0, 180 }, "theta0", "N_SiEndcapHits");
+  //auto hVtxBarrel_N_vs_theta = df0.Histo1D({"hVtxBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_VtxBarrelHits");
+
+  auto hSiBarrel_Nhits  = df0.Histo1D({"hSiBarrel_Nhits", "; #theta [deg.]",   20, 0, 20 }, "N_SiBarrelHits");
+  auto hSiEndcap_Nhits  = df0.Histo1D({"hSiEndcap_Nhits", "; #theta [deg.]",   20, 0, 20 }, "N_SiEndcapHits");
+  //auto hVtxBarrel_Nhits = df0.Histo1D({"hVtxBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 },  "N_VtxBarrelHits");
+
+  auto hSiBarrel_Ntheta = df0.Histo1D({"hSiBarrel_Ntheta", "; #theta [deg.]",   20, 0, 180 }, "theta0");
+  auto hSiEndcap_Ntheta = df0.Histo1D({"hSiEndcap_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*) hSiBarrel_N_vs_theta->Clone();
+  auto h2 = (TH1D*) hSiBarrel_Ntheta->Clone();
+  h1->Divide(h2);
+  hs->Add(h1);
+  h1 = (TH1D*) hSiEndcap_N_vs_theta->Clone();
+  h2 = (TH1D*) hSiEndcap_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*) hSiBarrel_N_vs_theta->Clone();
+  h2 = (TH1D*) hSiBarrel_Ntheta->Clone();
+  //h1->Divide(h2);
+  hs->Add(h2);
+  h1 = (TH1D*) hSiEndcap_N_vs_theta->Clone();
+  h2 = (TH1D*) hSiEndcap_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*) hSiBarrel_Nhits->Clone();
+  hs->Add(h1);
+  h1 = (TH1D*) hSiEndcap_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();
+  hSiBarrel_x_vs_y->DrawCopy("colz");
+  c->SaveAs("results/tracking/hits_central_electrons_xy.png");
+  c->SaveAs("results/tracking/hits_central_electrons_xy.pdf");
+  return 0;
+}

options/tracker_reconstruction.py

@@ -73,11 +73,13 @@ ecal_digi = EMCalorimeterDigi("ecal_digi",
 ufsd_digi = UFSDTrackerDigi("ufsd_digi", 
         inputHitCollection="SiTrackerBarrelHits",
         outputHitCollection="SiTrackerBarrelRawHits",
-        timeResolution=8)
+        timeResolution=8,
+        OutputLevel=DEBUG)
 ufsd_digi2 = UFSDTrackerDigi("ufsd_digi2", 
         inputHitCollection="SiTrackerEndcapHits",
         outputHitCollection="SiTrackerEndcapRawHits",
-        timeResolution=8)
+        timeResolution=8,
+        OutputLevel=DEBUG)
 
 #vtx_digi = UFSDTrackerDigi("vtx_digi", 
 #        inputHitCollection="SiVertexBarrelHits",
@@ -88,15 +90,15 @@ ufsd_digi2 = UFSDTrackerDigi("ufsd_digi2",
 ecal_reco = EMCalReconstruction("ecal_reco", 
         inputHitCollection="RawEcalBarrelHits", 
         outputHitCollection="RecEcalBarrelHits",
-        minModuleEdep=0.0*units.MeV,
-        OutputLevel=DEBUG)
+        minModuleEdep=0.0*units.MeV)
+        #OutputLevel=DEBUG)
 
 simple_cluster = SimpleClustering("simple_cluster", 
         inputHitCollection="RecEcalBarrelHits", 
         outputClusters="SimpleClusters",
         minModuleEdep=1.0*units.MeV,
-        maxDistance=50.0*units.cm,
-        OutputLevel=DEBUG)
+        maxDistance=50.0*units.cm)
+        #OutputLevel=DEBUG)
 
 trk_barrel_reco = TrackerHitReconstruction("trk_barrel_reco",
         inputHitCollection="SiTrackerBarrelRawHits",
@@ -113,27 +115,26 @@ trk_endcap_reco = TrackerHitReconstruction("trk_endcap_reco",
 # Source linker 
 sourcelinker = TrackerSourceLinker("trk_srclinker",
         inputHitCollection="TrackerBarrelRecHits",
-        outputSourceLinks="BarrelTrackSourceLinks",
-        OutputLevel=DEBUG)
+        outputSourceLinks="BarrelTrackSourceLinks")
 
 trk_hits_srclnkr = Tracker2SourceLinker("trk_hits_srclnkr",
         TrackerBarrelHits="TrackerBarrelRecHits",
         TrackerEndcapHits="TrackerEndcapRecHits",
         outputMeasurements="lnker2Measurements",
         outputSourceLinks="lnker2Links",
-        allTrackerHits="linker2AllHits",
-        OutputLevel=DEBUG)
+        allTrackerHits="linker2AllHits")
+        #OutputLevel=DEBUG)
 
 ## Track param init
 truth_trk_init = TrackParamTruthInit("truth_trk_init",
         inputMCParticles="mcparticles",
-        outputInitialTrackParameters="InitTrackParams",
-        OutputLevel=DEBUG)
+        outputInitialTrackParameters="InitTrackParams")
+        #OutputLevel=DEBUG)
 
 clust_trk_init = TrackParamClusterInit("clust_trk_init",
         inputClusters="SimpleClusters",
-        outputInitialTrackParameters="InitTrackParamsFromClusters",
-        OutputLevel=DEBUG)
+        outputInitialTrackParameters="InitTrackParamsFromClusters")
+        #OutputLevel=DEBUG)
 
 #vtxcluster_trk_init = TrackParamVertexClusterInit("vtxcluster_trk_init",
 #        inputVertexHits="VertexBarrelRecHits",
@@ -147,38 +148,38 @@ trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
         inputSourceLinks = sourcelinker.outputSourceLinks,
         inputMeasurements = sourcelinker.outputMeasurements,
         inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
-        outputTrajectories="trajectories",
-        OutputLevel=DEBUG)
+        outputTrajectories="trajectories")
+        #OutputLevel=DEBUG)
 parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
         inputTrajectories="trajectories",
         outputParticles="ReconstructedParticles",
-        outputTrackParameters="outputTrackParameters",
-        OutputLevel=DEBUG)
+        outputTrackParameters="outputTrackParameters")
+        #OutputLevel=DEBUG)
 
 trk_find_alg1 = TrackFindingAlgorithm("trk_find_alg1",
         inputSourceLinks = trk_hits_srclnkr.outputSourceLinks,
         inputMeasurements = trk_hits_srclnkr.outputMeasurements,
         inputInitialTrackParameters= "InitTrackParamsFromClusters", 
-        outputTrajectories="trajectories1",
-        OutputLevel=DEBUG)
+        outputTrajectories="trajectories1")
+        #OutputLevel=DEBUG)
 parts_from_fit1 = ParticlesFromTrackFit("parts_from_fit1",
         inputTrajectories="trajectories1",
         outputParticles="ReconstructedParticles1",
-        outputTrackParameters="outputTrackParameters1",
-        OutputLevel=DEBUG)
+        outputTrackParameters="outputTrackParameters1")
+        #OutputLevel=DEBUG)
 
 trk_find_alg2 = TrackFindingAlgorithm("trk_find_alg2",
         inputSourceLinks = trk_hits_srclnkr.outputSourceLinks,
         inputMeasurements = trk_hits_srclnkr.outputMeasurements,
         inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
         #inputInitialTrackParameters= "InitTrackParamsFromVtxClusters", 
-        outputTrajectories="trajectories2",
-        OutputLevel=DEBUG)
+        outputTrajectories="trajectories2")
+        #OutputLevel=DEBUG)
 parts_from_fit2 = ParticlesFromTrackFit("parts_from_fit2",
         inputTrajectories="trajectories2",
         outputParticles="ReconstructedParticles2",
-        outputTrackParameters="outputTrackParameters2",
-        OutputLevel=DEBUG)
+        outputTrackParameters="outputTrackParameters2")
+        #OutputLevel=DEBUG)
 
 
 #types = []