Merge pull request #27 from ornlneutronimaging/add_time_packet_func

Add time packet func

sophiread/FastSophiread/benchmarks/benchmark_raw2events.cpp

@@ -19,6 +19,8 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <https://www.gnu.org/licenses/>.
  */
+#include <spdlog/spdlog.h>
+
 #include <chrono>
 #include <fstream>
 #include <iostream>
@@ -41,14 +43,14 @@ std::vector<char> readTimepix3RawFile(const std::string& filepath) {
 
   // Check if file is open successfully
   if (!file.is_open()) {
-    std::cerr << "Failed to open file: " << filepath << std::endl;
+    spdlog::error("Failed to open file: {}", filepath);
     exit(EXIT_FAILURE);
   }
 
   // Get the size of the file
   std::streamsize fileSize = file.tellg();
   file.seekg(0, std::ios::beg);
-  std::cout << "File size (bytes): " << fileSize << std::endl;
+  spdlog::info("File size (bytes): {}", fileSize);
 
   // Create a vector to store the data
   std::vector<char> fileData(fileSize);
@@ -65,7 +67,7 @@ std::vector<char> readTimepix3RawFile(const std::string& filepath) {
 int main(int argc, char* argv[]) {
   // sanity check
   if (argc < 2) {
-    std::cerr << "Usage: " << argv[0] << " <input file>" << std::endl;
+    spdlog::error("Usage: {} <input file>", argv[0]);
     return 1;
   }
 
@@ -75,13 +77,23 @@ int main(int argc, char* argv[]) {
   auto raw_data = readTimepix3RawFile(in_tpx3);
   auto end = std::chrono::high_resolution_clock::now();
   auto elapsed = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
-  std::cout << "Read raw data: " << elapsed / 1e6 << " s" << std::endl;
+  spdlog::info("Read raw data: {} s", elapsed / 1e6);
 
   // single thread processing
   // -- run
-  std::cout << "\nSingle thread processing..." << std::endl;
+  spdlog::info("Single thread processing...");
   start = std::chrono::high_resolution_clock::now();
+  // locate all the TPX3H (chip dataset) in the raw data
   auto batches = findTPX3H(raw_data);
+
+  // extract all tdc and gdc timestamps
+  unsigned long tdc_timestamp = 0;
+  unsigned long long int gdc_timestamp = 0;
+  for (auto& tpx3 : batches) {
+    extractTGDC(tpx3, raw_data, tdc_timestamp, gdc_timestamp);
+  }
+
+  // process all batches to get neutron events
   auto abs_alg = std::make_unique<ABS>(5.0, 1, 75);
   int total_events = 0;
   for (auto& tpx3 : batches) {
@@ -102,18 +114,25 @@ int main(int argc, char* argv[]) {
     auto hits = tpx3.hits;
     n_hits += hits.size();
   }
-  std::cout << "Number of hits: " << n_hits << std::endl;
-  std::cout << "Number of events: " << total_events << std::endl;
+  spdlog::info("Number of hits: {}", n_hits);
+  spdlog::info("Number of events: {}", total_events);
   elapsed = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
-  std::cout << "Single thread processing: " << elapsed / 1e6 << " s" << std::endl;
+  spdlog::info("Single thread processing: {} s", elapsed / 1e6);
   auto speed = n_hits / (elapsed / 1e6);
-  std::cout << "Single thread processing speed: " << speed << " hits/s" << std::endl;
+  spdlog::info("Single thread processing speed: {} hits/s", speed);
 
   // multi-thread processing
   // -- run
-  std::cout << "\nMulti-thread processing..." << std::endl;
+  spdlog::info("Multi-thread processing...");
   start = std::chrono::high_resolution_clock::now();
+  // locate all the TPX3H (chip dataset) in the raw data, single thread
   auto batches_mt = findTPX3H(raw_data);
+  // extract all tdc and gdc timestamps, single thread
+  tdc_timestamp = 0;
+  gdc_timestamp = 0;
+  for (auto& tpx3 : batches_mt) {
+    extractTGDC(tpx3, raw_data, tdc_timestamp, gdc_timestamp);
+  }
   // use tbb parallel_for to process batches
   tbb::parallel_for(tbb::blocked_range<size_t>(0, batches_mt.size()), [&](const tbb::blocked_range<size_t>& r) {
     auto abs_alg_mt = std::make_unique<ABS>(5.0, 1, 75);
@@ -135,9 +154,25 @@ int main(int argc, char* argv[]) {
     auto hits = tpx3.hits;
     n_hits += hits.size();
   }
-  std::cout << "Number of hits: " << n_hits << std::endl;
+  spdlog::info("Number of hits: {}", n_hits);
   elapsed = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
-  std::cout << "Multi-thread processing: " << elapsed / 1e6 << " s" << std::endl;
+  spdlog::info("Multi-thread processing: {} s", elapsed / 1e6);
   speed = n_hits / (elapsed / 1e6);
-  std::cout << "Multi-thread processing speed: " << speed << " hits/s" << std::endl;
+  spdlog::info("Multi-thread processing speed: {} hits/s", speed);
+
+  // sanity check: hit.getTOF() should be smaller than 666,667 clock, which is
+  //               equivalent to 16.67 ms
+  int n_bad_hits = 0;
+  for (const auto& tpx3 : batches_mt) {
+    for (const auto& hit : tpx3.hits) {
+      auto tof_ms = hit.getTOF_ns() * 1e-6;
+      if (tof_ms > 16.67) {
+        spdlog::error("TOF: {} ms", tof_ms);
+        n_bad_hits++;
+      } else {
+        spdlog::debug("TOF: {} ms", tof_ms);
+      }
+    }
+  }
+  spdlog::info("bad/total hits: {}/{}", n_bad_hits, n_hits);
 }

sophiread/FastSophiread/include/tpx3_fast.h

@@ -36,10 +36,14 @@ struct TPX3 {
   const int num_packets;       // number of packets in the dataset batch (time packet and data packet)
   const int chip_layout_type;  // data source (sub-chip ID)
   std::vector<Hit> hits;       // hits extracted from the dataset batch
+  std::vector<unsigned long> tdcs;           // tdc extracted from the dataset batch
+  std::vector<unsigned long long int> gdcs;  // gdc extracted from the dataset batch
 
   TPX3(std::size_t index, int num_packets, int chip_layout_type)
       : index(index), num_packets(num_packets), chip_layout_type(chip_layout_type) {
     hits.reserve(num_packets);
+    tdcs.reserve(num_packets);
+    gdcs.reserve(num_packets);
   };
 
   void emplace_back(const char* packet, const unsigned long long tdc, const unsigned long long gdc) {
@@ -52,6 +56,14 @@ std::vector<TPX3> findTPX3H(ForwardIt first, ForwardIt last);
 std::vector<TPX3> findTPX3H(const std::vector<char>& raw_bytes);
 std::vector<TPX3> findTPX3H(char* raw_bytes, std::size_t size);
 
+template <typename ForwardIter>
+void extractTGDC(TPX3& tpx3h, ForwardIter bytes_begin, ForwardIter bytes_end, unsigned long& tdc_timestamp,
+                 unsigned long long int& gdc_timestamp);
+void extractTGDC(TPX3& tpx3h, const std::vector<char>& raw_bytes, unsigned long& tdc_timestamp,
+                 unsigned long long int& gdc_timestamp);
+void extractTGDC(TPX3& tpx3h, char* raw_bytes, std::size_t size, unsigned long& tdc_timestamp,
+                 unsigned long long int& gdc_timestamp);
+
 template <typename ForwardIter>
 void extractHits(TPX3& tpx3h, ForwardIter bytes_begin, ForwardIter bytes_end);
 void extractHits(TPX3& tpx3h, const std::vector<char>& raw_bytes);

sophiread/FastSophiread/src/hit.cpp

@@ -30,13 +30,16 @@
  * @param gdc
  * @param chip_layout_type
  */
-Hit::Hit(const char *packet, const unsigned long long tdc, const unsigned long long gdc, const int chip_layout_type) {
+Hit::Hit(const char *packet, const unsigned long long TDC_timestamp, const unsigned long long GDC_timestamp,
+         const int chip_layout_type) {
+  // local variables
   unsigned short pixaddr, dcol, spix, pix;
   unsigned short *spider_time;
   unsigned short *nTOT;    // bytes 2,3, raw time over threshold
   unsigned int *nTOA;      // bytes 3,4,5,6, raw time of arrival
   unsigned int *npixaddr;  // bytes 4,5,6,7
   unsigned int spidertime = 0;
+
   // timing information
   spider_time = (unsigned short *)(&packet[0]);  // Spider time  (16 bits)
   nTOT = (unsigned short *)(&packet[2]);         // ToT          (10 bits)
@@ -46,14 +49,9 @@ Hit::Hit(const char *packet, const unsigned long long tdc, const unsigned long l
   m_toa = (*nTOA >> 6) & 0x3FFF;
   spidertime = 16384 * (*spider_time) + m_toa;
 
-  // rename variables for clarity
-  unsigned long long int GDC_timestamp = gdc;
-  unsigned long long TDC_timestamp = tdc;
-
   // convert spidertime to global timestamp
   unsigned long SPDR_LSB30 = 0;
   unsigned long SPDR_MSB18 = 0;
-  //   unsigned long long SPDR_timestamp = 0;
 
   SPDR_LSB30 = GDC_timestamp & 0x3FFFFFFF;
   SPDR_MSB18 = (GDC_timestamp >> 30) & 0x3FFFF;
@@ -66,7 +64,7 @@ Hit::Hit(const char *packet, const unsigned long long tdc, const unsigned long l
   // tof calculation
   // TDC packets not always arrive before corresponding data packets
   if (m_spidertime < TDC_timestamp) {
-    m_tof = m_spidertime - TDC_timestamp + 666667;  // 1E9 / 60.0 is approximately 16666667
+    m_tof = m_spidertime - TDC_timestamp + 666667;
   } else {
     m_tof = m_spidertime - TDC_timestamp;
   }

sophiread/FastSophiread/src/tpx3_fast.cpp

@@ -77,21 +77,30 @@ std::vector<TPX3> findTPX3H(const std::vector<char> &raw_bytes) {
  */
 std::vector<TPX3> findTPX3H(char *raw_bytes, std::size_t size) { return findTPX3H(raw_bytes, raw_bytes + size); }
 
+/**
+ * @brief Extract the TDC and GDC timestamps from a TPX3H (chip dataset).
+ *
+ * @tparam ForwardIter
+ * @param tpx3h
+ * @param bytes_begin
+ * @param bytes_end
+ * @param tdc_timestamp: global-tracked TDC timestamp
+ * @param gdc_timestamp: global-tracked GDC timestamp
+ */
 template <typename ForwardIter>
-void extractHits(TPX3 &tpx3h, ForwardIter bytes_begin, ForwardIter bytes_end) {
+void extractTGDC(TPX3 &tpx3h, ForwardIter bytes_begin, ForwardIter bytes_end, unsigned long &tdc_timestamp,
+                 unsigned long long int &gdc_timestamp) {
   // Define the local variables
   // -- TDC
   unsigned long *tdclast;
   unsigned long long mytdc = 0;
   unsigned long TDC_MSB16 = 0;
   unsigned long TDC_LSB32 = 0;
-  unsigned long TDC_timestamp = 0;
   // -- GDC
   unsigned long *gdclast;
   unsigned long long mygdc = 0;
   unsigned long Timer_LSB32 = 0;
   unsigned long Timer_MSB16 = 0;
-  unsigned long long int GDC_timestamp = 0;  // 48-bit
 
   // Move to the first packet
   auto bytes_iter = bytes_begin;
@@ -110,28 +119,101 @@ void extractHits(TPX3 &tpx3h, ForwardIter bytes_begin, ForwardIter bytes_end) {
     if (char_array[7] == 0x6F) {  // TDC data packets
       tdclast = (unsigned long *)(&char_array[0]);
       mytdc = (((*tdclast) >> 12) & 0xFFFFFFFF);  // rick: 0x3fffffff, get 32-bit tdc
-      TDC_LSB32 = GDC_timestamp & 0xFFFFFFFF;
-      TDC_MSB16 = (GDC_timestamp >> 32) & 0xFFFF;
+      TDC_LSB32 = gdc_timestamp & 0xFFFFFFFF;
+      TDC_MSB16 = (gdc_timestamp >> 32) & 0xFFFF;
       if (mytdc < TDC_LSB32) {
         TDC_MSB16++;
       }
-      TDC_timestamp = (TDC_MSB16 << 32) & 0xFFFF00000000;
-      TDC_timestamp = TDC_timestamp | mytdc;
+      tdc_timestamp = (TDC_MSB16 << 32) & 0xFFFF00000000;
+      tdc_timestamp = tdc_timestamp | mytdc;
     } else if ((char_array[7] & 0xF0) == 0x40) {  // GDC data packet
+
       gdclast = (unsigned long *)(&char_array[0]);
       mygdc = (((*gdclast) >> 16) & 0xFFFFFFFFFFF);
+
       if (((mygdc >> 40) & 0xF) == 0x4) {
         Timer_LSB32 = mygdc & 0xFFFFFFFF;  // 32-bit
       } else if (((mygdc >> 40) & 0xF) == 0x5) {
+        // Serval sometimes report 0 GDC during experiment, so we need to check
+        // if the GDC is 0, if so, we use the previous GDC
+        auto gdc_tmp = gdc_timestamp;
+
         Timer_MSB16 = mygdc & 0xFFFF;  // 16-bit
-        GDC_timestamp = Timer_MSB16;
-        GDC_timestamp = (GDC_timestamp << 32) & 0xFFFF00000000;
-        GDC_timestamp = GDC_timestamp | Timer_LSB32;
+        gdc_tmp = Timer_MSB16;
+        gdc_tmp = (gdc_tmp << 32) & 0xFFFF00000000;
+        gdc_tmp = gdc_tmp | Timer_LSB32;
+
+        if (gdc_tmp != 0) {
+          gdc_timestamp = gdc_tmp;
+        }
       }
+
     } else if ((char_array[7] & 0xF0) == 0xb0) {  // data packet
+      tpx3h.tdcs.emplace_back(tdc_timestamp);
+      tpx3h.gdcs.emplace_back(gdc_timestamp);
+    }
+  }
+}
+
+/**
+ * @brief Extract the TDC and GDC timestamps from a TPX3H (chip dataset).
+ *
+ * @param tpx3h
+ * @param raw_bytes
+ * @param tdc_timestamp: global-tracked TDC timestamp
+ * @param gdc_timestamp: global-tracked GDC timestamp
+ */
+void extractTGDC(TPX3 &tpx3h, const std::vector<char> &raw_bytes, unsigned long &tdc_timestamp,
+                 unsigned long long int &gdc_timestamp) {
+  extractTGDC(tpx3h, raw_bytes.cbegin(), raw_bytes.cend(), tdc_timestamp, gdc_timestamp);
+}
+
+/**
+ * @brief Extract the TDC and GDC timestamps from a TPX3H (chip dataset).
+ *
+ * @param tpx3h
+ * @param raw_bytes
+ * @param size
+ * @param tdc_timestamp: global-tracked TDC timestamp
+ * @param gdc_timestamp: global-tracked GDC timestamp
+ */
+void extractTGDC(TPX3 &tpx3h, char *raw_bytes, std::size_t size, unsigned long &tdc_timestamp,
+                 unsigned long long int &gdc_timestamp) {
+  extractTGDC(tpx3h, raw_bytes, raw_bytes + size, tdc_timestamp, gdc_timestamp);
+}
+
+/**
+ * @brief Extract all hits from a TPX3H (chip dataset).
+ *
+ * @tparam ForwardIter
+ * @param tpx3h
+ * @param bytes_begin
+ * @param bytes_end
+ */
+template <typename ForwardIter>
+void extractHits(TPX3 &tpx3h, ForwardIter bytes_begin, ForwardIter bytes_end) {
+  // Move to the first packet
+  auto bytes_iter = bytes_begin;
+  std::advance(bytes_iter, tpx3h.index);
+
+  // Loop over all packets
+  auto hit_idx = 0;
+  for (auto j = 0; j < tpx3h.num_packets; ++j) {
+    if (std::next(bytes_iter, 8) >= bytes_end) {
+      continue;
+    }
+
+    bytes_iter = std::next(bytes_iter, 8);
+    const char *char_array = &(*bytes_iter);
+
+    // extract the data from the data packet
+    if ((char_array[7] & 0xF0) == 0xb0) {  // data packet
       // NOTE: we are implicitly calling the Hit constructor directly within the
       //       vector for speed.
-      tpx3h.emplace_back(char_array, TDC_timestamp, GDC_timestamp);
+      if ( tpx3h.tdcs[hit_idx] != 0 && tpx3h.gdcs[hit_idx] != 0){
+         tpx3h.emplace_back(char_array, tpx3h.tdcs[hit_idx], tpx3h.gdcs[hit_idx]);
+      }
+      ++hit_idx;
     }
   }
 }

sophiread/FastSophiread/tests/test_hit.cpp

@@ -44,7 +44,7 @@ class HitTest : public ::testing::Test {
  protected:
  protected:
   char packet[8] = {0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08};
-  unsigned long long tdc = 1000;
+  unsigned long long tdc = 8411155;
   unsigned long long gdc = 2000;
   int chip_layout_type = 0;
   Hit hit;
@@ -63,12 +63,12 @@ TEST_F(HitTest, CheckSpidertimens) {
 }
 
 TEST_F(HitTest, CheckTOF) {
-  unsigned long long expectedTOF = 8410156;  // Expected TOF
+  unsigned long long expectedTOF = 1;  // Expected TOF
   ASSERT_EQ(expectedTOF, hit.getTOF());
 }
 
 TEST_F(HitTest, CheckTOFns) {
-  double expectedTOFns = 210253900;  // Expected TOFns
+  double expectedTOFns = 25;  // Expected TOFns
   ASSERT_DOUBLE_EQ(expectedTOFns, hit.getTOF_ns());
 }