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/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;
+}