add noise hits to reconstructed hits by randomly generating cell id

src/algorithms/tracking/TrackerHitReconstruction.cc

@@ -15,6 +15,8 @@
 #include <iterator>
 #include <utility>
 #include <vector>
+#include <random>
+#include <bitset> // For displaying bits (optional)
 
 namespace eicrecon {
 
@@ -41,6 +43,8 @@ std::unique_ptr<edm4eic::TrackerHitCollection> TrackerHitReconstruction::process
 
     auto rec_hits { std::make_unique<edm4eic::TrackerHitCollection>() };
 
+
+    // add hits from DD4hep
     for (const auto& raw_hit : raw_hits) {
 
         auto id = raw_hit.getCellID();
@@ -85,6 +89,74 @@ std::unique_ptr<edm4eic::TrackerHitCollection> TrackerHitReconstruction::process
 
     }
 
+    //------------------------------------------------
+    // Example code of adding random noise hits
+    // Shujie Li, 08.2024
+    // use the first raw hits to obtain the volume ID of the detector system
+    // this also make sure the hit container has at least one hit
+    int num_hits=10; // total number of noise hits in this detector volume. Should move this to the config file
+    std::vector<uint64_t> noise_ids;
+
+    // generate valid random cell id
+    for (const auto& hit0 : raw_hits) {
+        uint64_t id0,vol_id;
+        id0 = hit0.getCellID();
+        // vol ID is the last 8 bits in Si tracker. Need to make it more flexible
+        // may want to predefine the layer/module ID as well to speed up the radom ID generation
+        vol_id = id0 & 0xFF;
+        // std::cout<<"::: "<<raw_hits.size()<<"/"<<id0<<"=="<<std::bitset<8>(id0)<<"  ::"<<vol_id<<std::endl;
+
+        // Setup random number generator
+        std::random_device rd; // Obtain a random number from hardware
+        std::mt19937_64 eng(rd()); // Seed the generator
+        std::uniform_int_distribution<uint64_t> distr; // Define the range for 64-bit integers
+        dd4hep::Position pos;
+        int nn=0;
+        int i=0;
+        while (i < num_hits) {
+            uint64_t randomNum = distr(eng); // Generate a random 64-bit number
+            randomNum = (randomNum & ~uint64_t(0xFF)) | vol_id; // Clear the last 8 bits and set them to 'vol ID'
+            try {
+                pos = m_converter->position(randomNum);
+                // std::cout<<" ::: converter position "<<pos.x()<<std::endl;
+                noise_ids.push_back(randomNum);
+                i++;
+            } catch(std::exception &e) {
+                // std::cout<<"::: cell ID error caught"<<std::endl;
+                nn++;
+            }
+            std::cout <<std::bitset<8>(vol_id) <<":"<< i<<"/"<<nn<<":::"<<randomNum<<"   ";
+        }
+        break;
+    }
+
+    // generate noise hits from ids.
+       for (auto id : noise_ids) {
+        // Get position and dimension
+        auto pos = m_converter->position(id);
+        auto dim = m_converter->cellDimensions(id);
+        // >oO trace
+        if(m_log->level() == spdlog::level::trace) {
+            m_log->trace("Noise hits inserted: position x={:.2f} y={:.2f} z={:.2f} [mm]: ", pos.x()/ mm, pos.y()/ mm, pos.z()/ mm);
+            m_log->trace("dimension size: {}", dim.size());
+            for (size_t j = 0; j < std::size(dim); ++j) {
+                m_log->trace(" - dimension {:<5} size: {:.2}",  j, dim[j]);
+            }
+        }
+ #if EDM4EIC_VERSION_MAJOR >= 7
+       auto rec_hit =
+ #endif
+       rec_hits->create(
+            id, // Raw DD4hep cell ID
+            edm4hep::Vector3f{static_cast<float>(pos.x() / mm), static_cast<float>(pos.y() / mm), static_cast<float>(pos.z() / mm)}, // mm
+            edm4eic::CovDiag3f{get_variance(dim[0] / mm), get_variance(dim[1] / mm), // variance (see note above)
+            std::size(dim) > 2 ? get_variance(dim[2] / mm) : 0.},
+                static_cast<float>(0), // ns
+            m_cfg.timeResolution,                            // in ns
+            static_cast<float>(0),   // Collected energy (GeV)
+            0.0F);                                       // Error on the energy
+
+    }
     return std::move(rec_hits);
 }