Move digitization code to the main branch.

src/algorithms/tracking/BTOFHitReconstruction.cc

@@ -0,0 +1,152 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2022 Whitney Armstrong, Sylvester Joosten, Wouter Deconinck, Dmitry Romanov
+
+#include "BTOFHitReconstruction.h"
+
+#include <Evaluator/DD4hepUnits.h>
+#include <Math/GenVector/Cartesian3D.h>
+#include <Math/GenVector/DisplacementVector3D.h>
+#include <edm4eic/CovDiag3f.h>
+#include <edm4eic/Cov3f.h>
+#include <edm4hep/Vector3f.h>
+#include <fmt/core.h>
+#include <spdlog/common.h>
+#include <stddef.h>
+#include <iterator>
+#include <utility>
+#include <vector>
+
+#include "TMatrixT.h"
+
+namespace eicrecon {
+
+void BTOFHitReconstruction::init(const dd4hep::rec::CellIDPositionConverter* converter, 
+		                 const dd4hep::Detector* detector,
+				 std::shared_ptr<spdlog::logger>& logger) {
+
+    m_log = logger;
+
+    m_converter = converter;
+    m_detector = detector;
+}
+
+dd4hep::rec::CellID BTOFHitReconstruction::getDetInfos(const dd4hep::rec::CellID& id) {
+    // retrieve segmentation class if that hasn't been done
+    if(! m_decoder) {
+        const auto det = m_converter -> findContext(id) -> element;
+        auto readout = m_converter -> findReadout(det);
+        auto seg = readout.segmentation();
+        m_decoder = seg.decoder();
+    }
+
+    // CellID for BarrelTOF is composed of 6 parts
+    // system, layer, module, sensor, x, y
+    // If we fix x and y to zero, what remains will be the detector information only
+    auto id_return = id;
+    m_decoder -> set(id_return, "x", 0);
+    m_decoder -> set(id_return, "y", 0);
+    return id_return;
+}
+
+std::unique_ptr<edm4eic::TrackerHitCollection> BTOFHitReconstruction::process(const edm4eic::RawTrackerHitCollection& TDCADC_hits) {
+    using dd4hep::mm;
+
+    auto rec_hits { std::make_unique<edm4eic::TrackerHitCollection>() };
+
+    // collection of ADC values from all sensors
+    std::unordered_map<dd4hep::rec::CellID, 
+	               std::vector<HitInfo>> hitsBySensors;
+
+    for (const auto& TDCADC_hit : TDCADC_hits) {
+
+        auto id = TDCADC_hit.getCellID();
+
+        // Get position and dimension
+        auto pos = m_converter->position(id);
+	// Get sensors info
+	auto detID = this -> getDetInfos(id);
+	hitsBySensors[detID].emplace_back(pos.x(), pos.y(), pos.z(), int(TDCADC_hit.getCharge()), int(TDCADC_hit.getTimeStamp()), id);
+
+    }
+
+    auto geoManager = m_detector -> world().volume() -> GetGeoManager();
+
+    // loop through each sensors for Hit information
+    TMatrixT<double> varLocal(3, 3);
+    for(int i = 0; i < 3; ++i)
+	    for(int j = 0; j < 3; ++j)
+		    varLocal[i][j] = 0;
+
+    for (const auto& sensor : hitsBySensors) {
+	// INSERT clustering algorithm for each sensors here
+	// Right now I just perform a simple average over all hits in a sensors
+	// Will be problematic near the edges, but it's just an illustration
+        double ave_x = 0, ave_y = 0, ave_z = 0;
+	double tot_charge = 0;
+	const auto& hits = sensor.second;
+	// find cellID for the cell with maximum ADC value within a sensor
+	// I don't know why you need cellID for reconstructed hits, but we'll do it anyway
+	auto id = hits[0].id;
+	auto curr_adc = hits[0].adc;
+	auto first_tdc = hits[0].tdc;
+	for(const auto& hit : hits) {
+            // weigh all hits by ADC value
+            ave_x += hit.adc*hit.x;
+	    ave_y += hit.adc*hit.y;
+	    ave_z += hit.adc*hit.z;
+
+	    tot_charge += hit.adc;
+	    if(hit.adc > curr_adc) {
+	        curr_adc = hit.adc;
+		id = hit.id;
+	    }
+	    first_tdc = std::min(first_tdc, hit.tdc);
+	}
+
+	ave_x /= tot_charge;
+	ave_y /= tot_charge;
+	ave_z /= tot_charge;
+
+	auto cellSize = m_converter -> cellDimensions(id);
+
+	// get rotation matrix
+	auto node = geoManager -> FindNode(hits[0].x, hits[0].y, hits[0].z);
+	auto currMatrix = geoManager -> GetCurrentMatrix();
+	auto rotMatrixElements = currMatrix -> GetRotationMatrix();
+
+	// rotMatrix transforms local coordinates to global coordinates
+	// see line 342 of https://root.cern.ch/doc/master/TGeoMatrix_8cxx_source.html#l00342
+	TMatrixT<double> rot(3, 3, rotMatrixElements);
+	TMatrixT<double> rotT(3, 3);
+	rotT.Transpose(rot);
+
+	varLocal[0][0] = cellSize[0]*cellSize[0] / mm / mm / 12.; // final division by 12 because I assumed uniform distribution
+	varLocal[1][1] = cellSize[1]*cellSize[1] / mm / mm / 12.; // Std. dev of uniformation = width/sqrt(12)
+	varLocal[2][2] = 0;
+
+	// transform variance. see https://robotics.stackexchange.com/questions/2556/how-to-rotate-covariance
+	auto varGlobal = rot*varLocal*rotT;
+
+
+	// adc to charge
+	float charge = tot_charge * m_cfg.c_slope + m_cfg.c_intercept;
+	// TDC to time
+	float time = first_tdc * m_cfg.t_slope + m_cfg.t_intercept;
+        // >oO trace
+        rec_hits->create(
+	    id,
+            edm4hep::Vector3f{static_cast<float>(ave_x / mm), 
+	                      static_cast<float>(ave_y / mm), 
+			      static_cast<float>(ave_z / mm)}, // mm
+            edm4eic::CovDiag3f{varGlobal[0][0], varGlobal[1][1], varGlobal[2][2]}, // should be the covariance of position
+            time, // ns
+            0.0F,                            // covariance of time
+            charge,   // total ADC sum
+            0.0F);                                       // Error on the energy
+
+    }
+
+    return std::move(rec_hits);
+}
+
+} // namespace eicrecon

src/algorithms/tracking/BTOFHitReconstruction.h

@@ -0,0 +1,44 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2022 Whitney Armstrong, Sylvester Joosten, Wouter Deconinck, Dmitry Romanov
+
+#pragma once
+
+#include <DDRec/CellIDPositionConverter.h>
+#include <DD4hep/Detector.h>
+#include <edm4eic/RawTrackerHitCollection.h>
+#include <edm4eic/TrackerHitCollection.h>
+#include <spdlog/logger.h>
+#include <memory>
+
+#include "BTOFHitReconstructionConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+    /**
+     * Produces edm4eic::TrackerHit with geometric info from edm4eic::RawTrackerHit
+     */
+    class BTOFHitReconstruction : public WithPodConfig<BTOFHitReconstructionConfig> {
+
+    public:
+        /// Once in a lifetime initialization
+        void init(const dd4hep::rec::CellIDPositionConverter* converter, const dd4hep::Detector* detector, std::shared_ptr<spdlog::logger>& logger);
+
+        /// Processes RawTrackerHit and produces a TrackerHit
+        std::unique_ptr<edm4eic::TrackerHitCollection> process(const edm4eic::RawTrackerHitCollection& TDCADC_hits);
+
+    private:
+	struct HitInfo{ double x, y, z; int adc, tdc; dd4hep::rec::CellID id;};
+	dd4hep::rec::CellID getDetInfos(const dd4hep::rec::CellID& id);
+        /** algorithm logger */
+        std::shared_ptr<spdlog::logger> m_log;
+
+        /// Cell ID position converter
+        const dd4hep::rec::CellIDPositionConverter* m_converter = nullptr; 
+
+	/// fetch sensor information from cellID
+	const dd4hep::DDSegmentation::BitFieldCoder* m_decoder = nullptr;
+
+	const dd4hep::Detector* m_detector = nullptr;
+    };
+}

src/algorithms/tracking/BTOFHitReconstructionConfig.h

@@ -0,0 +1,13 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2022 Whitney Armstrong, Sylvester Joosten, Wouter Deconinck, Dmitry Romanov
+
+#pragma once
+
+namespace eicrecon {
+    struct BTOFHitReconstructionConfig {
+        // parameters that convert ADC to EDep
+	double c_slope = 3.86976e-7, c_intercept=2.42716e-5;
+	// parameters that convert TDC to hit time (ns)
+	double t_slope = 0.0197305, t_intercept = 0.208047;
+    };
+}

src/detectors/BTOF/BTOF.cc

@@ -17,7 +17,7 @@
 #include "algorithms/pid_lut/PIDLookupConfig.h"
 #include "extensions/jana/JOmniFactoryGeneratorT.h"
 #include "factories/digi/SiliconTrackerDigi_factory.h"
-#include "factories/tracking/TrackerHitReconstruction_factory.h"
+#include "factories/tracking/BTOFHitReconstruction_factory.h"
 #include "global/pid_lut/PIDLookup_factory.h"
 #include "services/geometry/dd4hep/DD4hep_service.h"
 
@@ -45,12 +45,11 @@ void InitPlugin(JApplication *app) {
     ));
 
     // Convert raw digitized hits into hits with geometry info (ready for tracking)
-    app->Add(new JOmniFactoryGeneratorT<TrackerHitReconstruction_factory>(
+    app->Add(new JOmniFactoryGeneratorT<BTOFHitReconstruction_factory>(
         "TOFBarrelRecHit",
-        {"TOFBarrelRawHit"},    // Input data collection tags
+        {"TOFBarrelADCTDC"},    // Input data collection tags
         {"TOFBarrelRecHit"},     // Output data tag
         {
-            .timeResolution = 10,
         },
         app
     ));         // Hit reco default config for factories

src/detectors/BTOF/BTOFHitDigi.cc

@@ -161,8 +161,9 @@ std::unique_ptr<edm4eic::RawTrackerHitCollection> BTOFHitDigi::execute(const edm
                     glandau.GetPoint(j, x1, y1);
                     glandau.GetPoint(j + 1, x2, y2);
 
-                    if ((y1 >= (-thres[1]*adc_range/Vm) && y2 <= (-thres[1]*adc_range/Vm))) {
-                        intersectionX = x1 + (x2 - x1) * (thres[1] - y1) / (y2 - y1);
+		    double norm_threshold = -thres[1]*adc_range/Vm;
+                    if (y1 >= norm_threshold && y2 <= norm_threshold) {
+                        intersectionX = x1 + (x2 - x1) * (norm_threshold - y1) / (y2 - y1);
 
                         tdc = /*BTOFHitDigi::ToDigitalCode(*/ceil(intersectionX/0.02);//, tdc_bit);
 			for (; j < nBins - 1; j++) {
@@ -182,8 +183,11 @@ std::unique_ptr<edm4eic::RawTrackerHitCollection> BTOFHitDigi::execute(const edm
 
 
 
-                adc = round(-V);
-                rawhits -> create(neighbour, adc, tdc);
+		// limit the range of adc values
+                adc = std::min(static_cast<double>(adc_range), round(-V));
+		// only store valid hits
+		if(tdc < std::numeric_limits<int>::max() && adc > 0)
+                    rawhits -> create(neighbour, adc, tdc);
 
             }
 //-----------------------------------------------------------

src/detectors/BTOF/BarrelTOFNeighborFinder.cc

@@ -47,7 +47,9 @@ void BarrelTOFNeighborFinder::_findAllNeighborsInSensor(
     }
 
     try {
-      _converter -> position(testCell);
+      auto pos = _converter -> position(testCell);
+      if(testCell != _converter -> cellID(pos))
+        continue;
     } catch (const std::invalid_argument& err) {
       // Ignore CellID that is invalid
       continue;
@@ -106,12 +108,15 @@ BarrelTOFNeighborFinder::findAllNeighborInSensor(const dd4hep::rec::CellID& hitC
   return neighbors;
 }
 
-int BarrelTOFNeighborFinder::_getSensorID(const dd4hep::rec::CellID& hitCell) {
+const dd4hep::rec::CellID BarrelTOFNeighborFinder::_getSensorID(const dd4hep::rec::CellID& hitCell) {
   _initWithCell(hitCell);
-  // when x or y-index goes out of bound, sometimes the position will corresponds to a new sensor
-  // will fetch the new cellID here
-  auto newID = _converter->cellID(_converter->position(hitCell)); 
-  return _decoder->get(newID, "sensor"); 
+  // fix x-y, what you left with are ids that corresponds to sensor info
+  // cellID may change when position changes. 
+  auto sensorID = hitCell;//_converter -> cellID(_converter -> position(hitCell));
+  _decoder -> set(sensorID, "x", 0);
+  _decoder -> set(sensorID, "y", 0);
+
+  return sensorID; 
 }
 
 /****************************************

src/detectors/BTOF/BarrelTOFNeighborFinder.h

@@ -38,7 +38,7 @@ class BarrelTOFNeighborFinder {
   void _findAllNeighborsInSensor(dd4hep::rec::CellID hitCell,
                                  std::shared_ptr<std::vector<dd4hep::rec::CellID>>& ans,
                                  std::unordered_set<dd4hep::rec::CellID>& dp);
-  int _getSensorID(const dd4hep::rec::CellID& hitCell);
+  const dd4hep::rec::CellID _getSensorID(const dd4hep::rec::CellID& hitCell);
   // need to initialize the class with a cell from Barrel TOF
   void _initWithCell(const dd4hep::rec::CellID& hitCell);
 

src/factories/tracking/BTOFHitReconstruction_factory.h

@@ -0,0 +1,36 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2023 Wouter Deconinck
+
+#pragma once
+
+#include "algorithms/tracking/BTOFHitReconstruction.h"
+#include "services/geometry/dd4hep/DD4hep_service.h"
+#include "extensions/jana/JOmniFactory.h"
+
+namespace eicrecon {
+
+class BTOFHitReconstruction_factory :
+public JOmniFactory<BTOFHitReconstruction_factory, BTOFHitReconstructionConfig> {
+ private:
+    eicrecon::BTOFHitReconstruction m_algo;
+
+    PodioInput<edm4eic::RawTrackerHit> m_raw_hits_input {this};
+    PodioOutput<edm4eic::TrackerHit> m_rec_hits_output {this};
+
+    Service<DD4hep_service> m_geoSvc {this};
+
+public:
+    void Configure() {
+        m_algo.applyConfig(config());
+        m_algo.init(m_geoSvc().converter(), m_geoSvc().detector(), logger());
+    }
+
+    void ChangeRun(int64_t run_number) {
+    }
+
+    void Process(int64_t run_number, uint64_t event_number) {
+        m_rec_hits_output() = m_algo.process(*m_raw_hits_input());
+    }
+};
+
+} // eicrecon