feat: Bucket silicon raw hits in time

.gitignore

@@ -49,6 +49,8 @@ install/*
 *.out
 *.m_app
 
+# Apple stuff
+.DS_Store
 
 ##### i g n o r i n g   I D E s #####
 

src/algorithms/digi/SiliconTrackerDigi.cc

@@ -9,7 +9,7 @@
 #include <edm4hep/Vector3d.h>
 #include <edm4hep/Vector3f.h>
 #include <fmt/core.h>
-#include <algorithm>
+#include <vector>
 #include <cmath>
 #include <cstdint>
 #include <gsl/pointers>
@@ -28,20 +28,42 @@ void SiliconTrackerDigi::init() {
     };
 }
 
+// Logic: timeResolution is the smearing of the time.
+// Smeared time marks the beginning of a bucket. If another hit falls within the same bucket,
+// the energy is summed up. If not, a new hit is created.
 
 void SiliconTrackerDigi::process(
-        const SiliconTrackerDigi::Input& input,
-        const SiliconTrackerDigi::Output& output) const {
+    const SiliconTrackerDigi::Input& input,
+    const SiliconTrackerDigi::Output& output) const {
 
     const auto [sim_hits] = input;
     auto [raw_hits,associations] = output;
 
     // A map of unique cellIDs with temporary structure RawHit
-    std::unordered_map<std::uint64_t, edm4eic::MutableRawTrackerHit> cell_hit_map;
-
+    std::unordered_map<std::uint64_t, std::vector<edm4eic::MutableRawTrackerHit>> cell_hit_map;
+
+    auto bucket_width = (std::int32_t) (m_cfg.integrationWindow * 1e3);
+
+    // Some detectors have a steady pulse of hits - generate empty time buckets for those
+    // Conundrum: We can't know the number of buckets without knowing the timeSlice length
+    // Assuming this is only used for MAPS, and that the time slice length is <= 1ms,
+    // we need at most 1ms / 2000ns = 500 buckets
+    // Let's use 1000, and pad with over and underflow to allow smearing out of range
+    if ( m_cfg.prepopulate ) {
+        for (auto cell : cell_hit_map) {
+            auto& hits = cell.second;
+            for (int i = -1; i <= 1001; i++) {
+                hits.push_back(
+                    edm4eic::MutableRawTrackerHit{
+                    cell.first,
+                    0,
+                    0
+                } );
+            }
+        }
+    }
 
     for (const auto& sim_hit : *sim_hits) {
-
         // time smearing
         double time_smearing = m_gauss();
         double result_time = sim_hit.getTime() + time_smearing;
@@ -58,52 +80,74 @@ void SiliconTrackerDigi::process(
         debug("   time smearing: {:.4f}, resulting time = {:.4f} [ns]", time_smearing, result_time);
         debug("   hit_time_stamp: {} [~ps]", hit_time_stamp);
 
-
         if (sim_hit.getEDep() < m_cfg.threshold) {
             debug("  edep is below threshold of {:.2f} [keV]", m_cfg.threshold / dd4hep::keV);
             continue;
         }
 
-        if (cell_hit_map.count(sim_hit.getCellID()) == 0) {
-            // This cell doesn't have hits
-            cell_hit_map[sim_hit.getCellID()] = {
+        bool bucket_found = false;
+        if (cell_hit_map.count(sim_hit.getCellID()) == 1) {
+            // Update an existing hit?
+            for (auto& hit : cell_hit_map[sim_hit.getCellID()]) {
+                auto existing_time = hit.getTimeStamp();
+                // TODO: edge cases?
+                if ( hit_time_stamp >= existing_time && hit_time_stamp <= existing_time + bucket_width ) {
+                    // There is already a hit within the same time window
+                    debug("  Hit already exists in cell ID={}, within the same time bucket. Time stamp: {}, bucket from {} to {}",
+                         sim_hit.getCellID(), hit.getTimeStamp(), existing_time, existing_time + bucket_width);
+                    // sum deposited energy
+                    auto charge = hit.getCharge();
+                    hit.setCharge(charge + (std::int32_t) std::llround(sim_hit.getEDep() * 1e6));
+                    bucket_found = true;
+                    break;
+                } // time bucket found
+            } // loop over existing hits
+        } // cellID found
+
+        if (!bucket_found) {
+            // There is no hit in the same time bucket
+            debug("  No pre-existing hit in cell ID={} in the same time bucket. Time stamp: {}",
+                        sim_hit.getCellID(), sim_hit.getTime());
+
+            // Create a new hit
+            // Note: time uncertainty in the TrackerHitReconstruction is set independently
+            // It would probably be better to move it into the RawTrackerHit class
+            // (same for spatial uncertainty, actually)
+            // Note 2: It's possible to not fall into a bucket but still be close enough to one or
+            // more that uncertainties overlap. Cannot be avoided in the current setup.
+            // It could lead to ambiguity which bucket is chosen for a third hit in this area.
+            // In reality, this is probably more like dead time; revisit later.
+            cell_hit_map[sim_hit.getCellID()].push_back(
+                edm4eic::MutableRawTrackerHit{
                 sim_hit.getCellID(),
                 (std::int32_t) std::llround(sim_hit.getEDep() * 1e6),
-                hit_time_stamp  // ns->ps
-            };
-        } else {
-            // There is previous values in the cell
-            auto& hit = cell_hit_map[sim_hit.getCellID()];
-            debug("  Hit already exists in cell ID={}, prev. hit time: {}", sim_hit.getCellID(), hit.getTimeStamp());
-
-            // keep earliest time for hit
-            auto time_stamp = hit.getTimeStamp();
-            hit.setTimeStamp(std::min(hit_time_stamp, hit.getTimeStamp()));
-
-            // sum deposited energy
-            auto charge = hit.getCharge();
-            hit.setCharge(charge + (std::int32_t) std::llround(sim_hit.getEDep() * 1e6));
-        }
-    }
-
-    for (auto item : cell_hit_map) {
-        raw_hits->push_back(item.second);
-
-        for (const auto& sim_hit : *sim_hits) {
-          if (item.first == sim_hit.getCellID()) {
-            // set association
-            auto hitassoc = associations->create();
-            hitassoc.setWeight(1.0);
-            hitassoc.setRawHit(item.second);
+                hit_time_stamp
+            } );
+        } // bucket found
+    } // loop over sim hits
+
+    // Make associations
+    for (auto cell : cell_hit_map) {
+        for (auto& hit : cell.second) {
+            raw_hits->push_back(hit);
+
+            for (const auto& sim_hit : *sim_hits) {
+                if (cell.first == sim_hit.getCellID()) {
+                    // set association
+                    auto hitassoc = associations->create();
+                    hitassoc.setWeight(1.0);
+                    hitassoc.setRawHit(hit);
 #if EDM4EIC_VERSION_MAJOR >= 6
-            hitassoc.setSimHit(sim_hit);
+                    hitassoc.setSimHit(sim_hit);
 #else
-            hitassoc.addToSimHits(sim_hit);
+                    hitassoc.addToSimHits(sim_hit);
 #endif
-          }
-        }
+                } // if cellID matches
+            } // sim_hits
+        } //hits
+    } // cell_hit_map
 
-    }
-}
+
+} // process
 
 } // namespace eicrecon

src/algorithms/digi/SiliconTrackerDigiConfig.h

@@ -4,14 +4,18 @@
 #pragma once
 
 #include <DD4hep/DD4hepUnits.h>
+#include <edm4eic/unit_system.h>
+
 
 namespace eicrecon {
 
   struct SiliconTrackerDigiConfig {
     // sub-systems should overwrite their own
     // NB: be aware of thresholds in npsim! E.g. https://github.com/eic/npsim/pull/9/files
     double threshold  = 0 * dd4hep::keV;
-    double timeResolution = 8;   /// TODO 8 of what units??? Same TODO in juggler. Probably [ns]
+    double timeResolution    = 2000 * edm4eic::unit::ns; // Source for 8? This is arbitrary --> however, it should be >0, maybe large, to not distort the fitter by default. Choosing 2us for now.
+    double integrationWindow = 2000 * edm4eic::unit::ns; // Source for 8? This is arbitrary --> however, it should be >0, maybe large, to not distort the fitter by default. Choosing 2us for now.
+    bool prepopulate = false; // prepopulate cellHit map with empty hits, for MAPS
   };
 
 } // eicrecon

src/algorithms/tracking/TrackerHitReconstructionConfig.h

@@ -2,9 +2,10 @@
 // Copyright (C) 2022 Whitney Armstrong, Sylvester Joosten, Wouter Deconinck, Dmitry Romanov
 
 #pragma once
+#include <edm4eic/unit_system.h>
 
 namespace eicrecon {
     struct TrackerHitReconstructionConfig {
-        float timeResolution = 10;
+    double timeResolution    = 2000 * edm4eic::unit::ns; // It should be >0, maybe large, to not distort the fitter by default. Choosing 2us for now.
     };
 }

src/detectors/B0TRK/B0TRK.cc

@@ -19,6 +19,8 @@ void InitPlugin(JApplication *app) {
     using namespace eicrecon;
 
     // Digitization
+    auto B0TrackerTimeResolution    =   0.02 * dd4hep::ns; // 15-20 ps
+    auto B0TrackerIntegrationWindow =  0.750 * dd4hep::ns; // 750 ps shaping time
     app->Add(new JOmniFactoryGeneratorT<SiliconTrackerDigi_factory>(
         "B0TrackerRawHits",
         {
@@ -30,7 +32,8 @@ void InitPlugin(JApplication *app) {
         },
         {
             .threshold = 10.0 * dd4hep::keV,
-            .timeResolution = 8,
+            .timeResolution = B0TrackerTimeResolution,
+            .integrationWindow = B0TrackerIntegrationWindow,
         },
         app
     ));
@@ -41,7 +44,7 @@ void InitPlugin(JApplication *app) {
         {"B0TrackerRawHits"},
         {"B0TrackerRecHits"},
         {
-            .timeResolution = 8,
+            .timeResolution = B0TrackerTimeResolution,
         },
         app
     ));

src/detectors/BTOF/BTOF.cc

@@ -30,6 +30,8 @@ void InitPlugin(JApplication *app) {
     using namespace eicrecon;
 
     // Digitization
+    auto BTOFBarrelTimeResolution    = 0.02 * dd4hep::ns; // 20 ps bin width
+    auto BTOFBarrelIntegrationWindow =  5.0 * dd4hep::ns; // shaping time
     app->Add(new JOmniFactoryGeneratorT<SiliconTrackerDigi_factory>(
         "TOFBarrelRawHit",
         {
@@ -41,7 +43,8 @@ void InitPlugin(JApplication *app) {
         },
         {
             .threshold = 6.0 * dd4hep::keV,
-            .timeResolution = 0.025,    // [ns]
+            .timeResolution = BTOFBarrelTimeResolution,
+            .integrationWindow = BTOFBarrelIntegrationWindow,
         },
         app
     ));
@@ -52,7 +55,7 @@ void InitPlugin(JApplication *app) {
         {"TOFBarrelRawHit"},    // Input data collection tags
         {"TOFBarrelRecHit"},     // Output data tag
         {
-            .timeResolution = 10,
+            .timeResolution = BTOFBarrelTimeResolution,
         },
         app
     ));         // Hit reco default config for factories

src/detectors/BTRK/BTRK.cc

@@ -1,7 +1,5 @@
 // Copyright 2022, Dmitry Romanov
 // Subject to the terms in the LICENSE file found in the top-level directory.
-//
-//
 
 #include <Evaluator/DD4hepUnits.h>
 #include <JANA/JApplication.h>
@@ -19,6 +17,8 @@ void InitPlugin(JApplication *app) {
     using namespace eicrecon;
 
     // Digitization
+    auto BTRKBarrelTimeResolution    = 2000 * dd4hep::ns;
+    auto BTRKBarrelIntegrationWindow = 2000 * dd4hep::ns; // shaping time
     app->Add(new JOmniFactoryGeneratorT<SiliconTrackerDigi_factory>(
         "SiBarrelRawHits",
         {
@@ -30,17 +30,21 @@ void InitPlugin(JApplication *app) {
         },
         {
             .threshold = 0.54 * dd4hep::keV,
+            .timeResolution = BTRKBarrelTimeResolution,
+            .integrationWindow = BTRKBarrelIntegrationWindow,
+            .prepopulate = true, // for MAPS, initialize digitization with a "pulse" of empty hits
         },
         app
     ));
 
-
     // Convert raw digitized hits into hits with geometry info (ready for tracking)
     app->Add(new JOmniFactoryGeneratorT<TrackerHitReconstruction_factory>(
         "SiBarrelTrackerRecHits",
         {"SiBarrelRawHits"},
         {"SiBarrelTrackerRecHits"},
-        {}, // default config
+        {
+            .timeResolution = BTRKBarrelTimeResolution,
+        },
         app
     ));
 

src/detectors/BVTX/BVTX.cc

@@ -19,6 +19,8 @@ void InitPlugin(JApplication *app) {
     using namespace eicrecon;
 
     // Digitization
+    auto SiBarrelVertexTimeResolution    = 2000 * dd4hep::ns;
+    auto SiBarrelVertexIntegrationWindow = 2000 * dd4hep::ns;
     app->Add(new JOmniFactoryGeneratorT<SiliconTrackerDigi_factory>(
         "SiBarrelVertexRawHits",
         {
@@ -30,6 +32,9 @@ void InitPlugin(JApplication *app) {
         },
         {
             .threshold = 0.54 * dd4hep::keV,
+            .timeResolution    = SiBarrelVertexTimeResolution,
+            .integrationWindow = SiBarrelVertexIntegrationWindow,
+            .prepopulate = true, // for MAPS, initialize digitization with a "pulse" of empty hits
         },
         app
     ));
@@ -39,7 +44,9 @@ void InitPlugin(JApplication *app) {
         "SiBarrelVertexRecHits",
         {"SiBarrelVertexRawHits"},
         {"SiBarrelVertexRecHits"},
-        {}, // default config
+        {
+            .timeResolution = SiBarrelVertexTimeResolution,
+        },
         app
     ));
 

src/detectors/ECTOF/ECTOF.cc

@@ -19,6 +19,8 @@ void InitPlugin(JApplication *app) {
     using namespace eicrecon;
 
     // Digitization
+    auto TOFEndcapTimeResolution    =  20.0 * dd4hep::ns; // 20 ps bin width
+    auto TOFEndcapIntegrationWindow =   5.0 * dd4hep::ns; // shaping time
     app->Add(new JOmniFactoryGeneratorT<SiliconTrackerDigi_factory>(
       "TOFEndcapRawHits",
       {
@@ -29,7 +31,9 @@ void InitPlugin(JApplication *app) {
       },
       {
         .threshold = 6.0 * dd4hep::keV,
-        .timeResolution = 0.025,
+        .timeResolution = TOFEndcapTimeResolution,
+        .integrationWindow = TOFEndcapIntegrationWindow,
+
       },
       app
     ));
@@ -40,7 +44,7 @@ void InitPlugin(JApplication *app) {
       {"TOFEndcapRawHits"},     // Input data collection tags
       {"TOFEndcapRecHits"},     // Output data tag
       {
-        .timeResolution = 0.025,
+        .timeResolution = TOFEndcapTimeResolution,
       },
       app
     ));

src/detectors/ECTRK/ECTRK.cc

@@ -19,6 +19,8 @@ void InitPlugin(JApplication *app) {
     using namespace eicrecon;
 
     // Digitization
+    auto SiEndcapTrackerTimeResolution    = 2000. * dd4hep::ns;
+    auto SiEndcapTrackerIntegrationWindow = 2000. * dd4hep::ns;
     app->Add(new JOmniFactoryGeneratorT<SiliconTrackerDigi_factory>(
         "SiEndcapTrackerRawHits",
         {
@@ -30,6 +32,9 @@ void InitPlugin(JApplication *app) {
         },
         {
             .threshold = 0.54 * dd4hep::keV,
+            .timeResolution = SiEndcapTrackerTimeResolution,
+            .integrationWindow = SiEndcapTrackerIntegrationWindow,
+            .prepopulate = true, // for MAPS, initialize digitization with a "pulse" of empty hits
         },
         app
     ));
@@ -39,7 +44,9 @@ void InitPlugin(JApplication *app) {
         "SiEndcapTrackerRecHits",
         {"SiEndcapTrackerRawHits"},
         {"SiEndcapTrackerRecHits"},
-        {}, // default config
+        {
+            .timeResolution = SiEndcapTrackerTimeResolution,
+        },
         app
     ));