Added average number of hits vs eta

benchmarks/track_fitting/options/track_reconstruction.py

@@ -3,7 +3,7 @@ from Gaudi.Configuration import *
 from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
 from GaudiKernel import SystemOfUnits as units
 
-detector_name = "topside"
+detector_name = "athena"
 if "JUGGLER_DETECTOR" in os.environ :
   detector_name = str(os.environ["JUGGLER_DETECTOR"])
 
@@ -124,16 +124,14 @@ trk_hit_col = TrackingHitsCollector("trk_hit_col",
             str(vtx_b_reco.outputHitCollection),
             str(vtx_ec_reco.outputHitCollection),
             str(gem_ec_reco.outputHitCollection) ],
-        trackingHits="trackingHits",
-        OutputLevel=DEBUG)
+        trackingHits="trackingHits")
 algorithms.append( trk_hit_col )
 
 # Hit Source linker 
 sourcelinker = TrackerSourceLinker("sourcelinker",
         inputHitCollection=trk_hit_col.trackingHits,
         outputSourceLinks="TrackSourceLinks",
-        outputMeasurements="TrackMeasurements",
-        OutputLevel=DEBUG)
+        outputMeasurements="TrackMeasurements")
 algorithms.append( sourcelinker )
 
 ## Track param init

benchmarks/track_fitting/scripts/gen_single_tracks.cxx

@@ -19,8 +19,13 @@ using namespace HepMC3;
 void gen_single_tracks(int n_events = 100, 
                      const char* out_fname = "single_tracks.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));
+  // Throw flat in cos(theta)
+  //double cos_theta_min = std::cos( 2.0*(M_PI/180.0));
+  //double cos_theta_max = std::cos(178.0*(M_PI/180.0));
+
+  // Throw flat in eta
+  double eta_min = -3.5;
+  double eta_max =  3.5;
 
   WriterAscii hepmc_output(out_fname);
   int events_parsed = 0;
@@ -41,13 +46,15 @@ void gen_single_tracks(int n_events = 100,
         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);
+    double p     = r1->Uniform(1.0, 10.0);
+    double phi   = r1->Uniform(0.0, 2.0 * M_PI);
+    double eta   = r1->Uniform(eta_min, eta_max);
+    double th    = 2.0*std::atan(std::exp(-eta));
+    //double costh = r1->Uniform(cos_theta_min, cos_theta_max);
+    //double th    = std::acos(costh);
+    double px    = p * std::cos(phi) * std::sin(th);
+    double py    = p * std::sin(phi) * std::sin(th);
+    double 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);

benchmarks/track_fitting/scripts/rec_single_tracks.cxx

@@ -14,6 +14,9 @@ R__LOAD_LIBRARY(libDD4pod.so)
 #include "eicd/ClusterData.h"
 #include "eicd/TrackerHitCollection.h"
 
+#include "common_bench/util.h"
+namespace cb = common_bench;
+
 using ROOT::RDataFrame;
 using namespace ROOT::VecOps;
 
@@ -89,24 +92,20 @@ int rec_single_tracks(const char* fname = "topside/rec_single_tracks.root")
                  .Define("thrownP", fourvec, {"thrownParticles"})
                  .Define("p_thrown", momentum, {"thrownP"})
                  .Define("theta_thrown", theta, {"thrownP"})
+                 .Define("eta_thrown",cb::eta , {"thrownP"})
                  .Define("theta0", "theta_thrown[0]")
+                 .Define("eta0", "eta_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<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_vs_eta = df0.Histo2D({"h_nTracks_vs_eta", "; #eta; N tracks ", 50,-4,4,10, 0, 10}, "eta0","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");
@@ -120,6 +119,10 @@ int rec_single_tracks(const char* fname = "topside/rec_single_tracks.root")
   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 hNhits_vs_eta    = df0.Histo1D({"hNhits_vs_eta", "; #eta ",    50, -4, 4 }, "eta0", "N_Hits");
+  auto hBarrel_N_vs_eta = df0.Histo1D({"hBarrel_N_vs_eta", "; #eta ", 50, -4, 4 }, "eta0", "N_BarrelHits");
+  auto hEndcap_N_vs_eta = df0.Histo1D({"hEndcap_N_vs_eta", "; #eta ", 50, -4, 4 }, "eta0", "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");
@@ -130,6 +133,9 @@ int rec_single_tracks(const char* fname = "topside/rec_single_tracks.root")
   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 hHits_Neta   = df0.Histo1D({"hHits_Neta", "; #eta [deg.]",   50, -4, 4 }, "eta0");
+  auto hBarrel_Neta = df0.Histo1D({"hBarrel_Neta", "; #eta [deg.]", 50, -4, 4 }, "eta0");
+  auto hEndcap_Neta = df0.Histo1D({"hEndcap_Neta", "; #eta [deg.]", 50, -4, 4 }, "eta0");
   // -----------------------------------------------
   auto c = new TCanvas();
 
@@ -197,6 +203,33 @@ int rec_single_tracks(const char* fname = "topside/rec_single_tracks.root")
   c->SaveAs("results/track_fitting/rec_single_tracks_n_hits_vs_theta.png");
   c->SaveAs("results/track_fitting/rec_single_tracks_n_hits_vs_theta.pdf");
 
+  // -----------------------------------------------
+  c  = new TCanvas();
+  hs = new THStack("n_hits","; #eta  ");
+
+  h1 = (TH1D*) hBarrel_N_vs_eta->Clone();
+  h2 = (TH1D*) hBarrel_Neta->Clone();
+  h1->SetLineColor(4);
+  h1->Divide(h2);
+  hs->Add(h1);
+
+  h1 = (TH1D*) hEndcap_N_vs_eta->Clone();
+  h2 = (TH1D*) hEndcap_Neta->Clone();
+  h1->Divide(h2);
+  h1->SetLineColor(2);
+  hs->Add(h1);
+
+  h1 = (TH1D*) hNhits_vs_eta->Clone();
+  h2 = (TH1D*) hHits_Neta->Clone();
+  h1->Divide(h2);
+  h1->SetLineColor(1);
+  hs->Add(h1);
+
+  hs->Draw("nostack, hist");
+  c->BuildLegend();
+  c->SaveAs("results/track_fitting/rec_single_tracks_n_hits_vs_eta.png");
+  c->SaveAs("results/track_fitting/rec_single_tracks_n_hits_vs_eta.pdf");
+
   // -----------------------------------------------
   c  = new TCanvas();
   hs = new THStack("theta","; #theta  ");