Merge branch 'LDMX-Software:trunk' into trunk

Tracking/include/Tracking/Reco/TruthSeedProcessor.h

@@ -249,6 +249,9 @@ class TruthSeedProcessor : public TrackingGeometryUser {
   // skip the recoil tracker
   bool skip_recoil_{false};
 
+  // Maximum track id for hit to be selected from target scoring plane
+  int max_track_id_{5};
+
   std::shared_ptr<LinPropagator> linpropagator_;
 
   // Track Extrapolator Tool :: TODO Use the real extrapolator!

Tracking/python/tracking.py

@@ -280,7 +280,9 @@ class TruthSeedProcessor(Producer):
     skip_tagger : bool 
         Ignore the tagger tracker(makes empty collections). 
     skip_recoil : bool 
-        Ignore the recoil tracker(makes empty collections)
+        Ignore the recoil tracker(makes empty collections).
+    max_track_id : double
+        Maximum track ID for a hit to be selected in the target scoring plane.
     """
     def __init__(self, instance_name = "TruthSeedProcessor"):
         super().__init__(instance_name, 'tracking::reco::TruthSeedProcessor',
@@ -299,3 +301,4 @@ def __init__(self, instance_name = "TruthSeedProcessor"):
         self.p_cut_ecal = -1. #MeV
         self.skip_tagger = False
         self.skip_recoil = False
+        self.max_track_id = 5

Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx

@@ -40,6 +40,7 @@ void TruthSeedProcessor::configure(framework::config::Parameters& parameters) {
   recoil_sp_ = parameters.getParameter<double>("recoil_sp", true);
   target_sp_ = parameters.getParameter<double>("tagger_sp", true);
   seedSmearing_ = parameters.getParameter<bool>("seedSmearing", false);
+  max_track_id_ = parameters.getParameter<int>("max_track_id", 5);
 
   ldmx_log(info) << "Seed Smearing is set to " << seedSmearing_;
 
@@ -556,8 +557,8 @@ void TruthSeedProcessor::produce(framework::Event& event) {
   std::vector<int> recoil_sh_idxs;
   std::unordered_map<int, std::vector<int>> recoil_sh_count_map;
 
-  // We are only interested in the beam electron
-  int idx_taggerhit = -1;
+  std::vector<int> tagger_sh_idxs;
+  std::unordered_map<int, std::vector<int>> tagger_sh_count_map;
 
   // Target scoring hits for Tagger will have Z<0, Recoil scoring hits will have
   // Z>0
@@ -573,16 +574,13 @@ void TruthSeedProcessor::produce(framework::Event& event) {
     if (zhit < 0.) {
       // Tagger selection cuts
       // Negative scoring plane hit, with momentum > p_cut
-      if (p_vec(2) < 0. || p_vec.norm() < p_cut_ || hit.getPdgID() != 11 ||
-          hit.getTrackID() != 1)
+      if (p_vec(2) < 0. || p_vec.norm() < p_cut_) continue;
 
-      {
-        continue;
-      }
+      // Check that the hit was left by a charged particle
+      if (abs(particleMap[hit.getTrackID()].getCharge()) < 1e-8) continue;
 
       if (p_vec.norm() > tagger_p_max) {
-        idx_taggerhit = i_sh;
-        tagger_p_max = p_vec.norm();
+        tagger_sh_count_map[hit.getTrackID()].push_back(i_sh);
       }
     }  // Tagger loop
 
@@ -616,6 +614,23 @@ void TruthSeedProcessor::produce(framework::Event& event) {
         });
   }
 
+  // Sort tagger hits.
+  for (auto& [_track_id, hit_indices] : tagger_sh_count_map) {
+    std::sort(
+        hit_indices.begin(), hit_indices.end(),
+        [&](const int idx1, int idx2) -> bool {
+          const ldmx::SimTrackerHit& hit1 = scoring_hits.at(idx1);
+          const ldmx::SimTrackerHit& hit2 = scoring_hits.at(idx2);
+
+          Acts::Vector3 phit1{hit1.getMomentum()[0], hit1.getMomentum()[1],
+                              hit1.getMomentum()[2]};
+          Acts::Vector3 phit2{hit2.getMomentum()[0], hit2.getMomentum()[1],
+                              hit2.getMomentum()[2]};
+
+          return phit1.norm() > phit2.norm();
+        });
+  }
+
   // Building of the event truth information and the truth seeds
   // TODO remove the truthtracks in the future as the truth seeds are enough
 
@@ -639,36 +654,28 @@ void TruthSeedProcessor::produce(framework::Event& event) {
   auto beamOriginSurface{Acts::Surface::makeShared<Acts::PerigeeSurface>(
       Acts::Vector3(beamOrigin_[0], beamOrigin_[1], beamOrigin_[2]))};
 
-  // I found a tagger scoring plane hit
-  if (idx_taggerhit != -1 && !skip_tagger_) {
-    const ldmx::SimTrackerHit& hit = scoring_hits.at(idx_taggerhit);
-    const ldmx::SimParticle& phit = particleMap[hit.getTrackID()];
-
-    if (hit_count_map_tagger[hit.getTrackID()].size() > n_min_hits_tagger_) {
-      ldmx::Track truth_tagger_track;
-      createTruthTrack(phit, hit, truth_tagger_track, targetSurface);
-      truth_tagger_track.setNhits(
-          hit_count_map_tagger[hit.getTrackID()].size());
-      // get track state at the generation point
-      // ldmx::TruthTrack::TrackState ts_beamOrigin(phit,"beam_origin");
-      // propagate track to target
-      // trackExtrapolator(truth_tagger_track, perigee_surface);
-      // truth_tagger_track.addTrackState(ts_beamOrigin);
-      tagger_truth_tracks.push_back(truth_tagger_track);
+  if (!skip_tagger_) {
+    for (const auto& [track_id, hit_indices] : tagger_sh_count_map) {
+      const ldmx::SimTrackerHit& hit = scoring_hits.at(hit_indices.at(0));
+      const ldmx::SimParticle& phit = particleMap[hit.getTrackID()];
+
+      if (hit_count_map_tagger[hit.getTrackID()].size() > n_min_hits_tagger_) {
+        ldmx::Track truth_tagger_track;
+        createTruthTrack(phit, hit, truth_tagger_track, targetSurface);
+        truth_tagger_track.setNhits(
+            hit_count_map_tagger[hit.getTrackID()].size());
+        tagger_truth_tracks.push_back(truth_tagger_track);
+
+        if (hit.getPdgID() == 11 && hit.getTrackID() < max_track_id_) {
+          ldmx::Track beamETruthSeed = TaggerFullSeed(
+              particleMap[hit.getTrackID()], hit.getTrackID(), hit,
+              hit_count_map_tagger, beamOriginSurface, targetUnboundSurface);
+          beam_electrons.push_back(beamETruthSeed);
+        }
+      }
     }
   }
 
-  // Form the tagger full seed.
-  // TODO This won't work for multiple electrons sample. Fix.
-
-  if (idx_taggerhit != -1 && !skip_tagger_) {
-    ldmx::Track beamETruthSeed = TaggerFullSeed(
-        particleMap[1], 1, scoring_hits.at(idx_taggerhit), hit_count_map_tagger,
-        beamOriginSurface, targetUnboundSurface);
-
-    beam_electrons.push_back(beamETruthSeed);
-  }
-
   // Recover the EcalScoring hits
   std::vector<ldmx::SimTrackerHit> ecal_spHits =
       event.getCollection<ldmx::SimTrackerHit>("EcalScoringPlaneHits");