src/components/TrackerHitReconstruction.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Sylvester Joosten, Wouter Deconinck
0003
0004 #include <algorithm>
0005
0006 // Gaudi
0007 #include "Gaudi/Algorithm.h"
0008 //#include "GaudiKernel/ToolHandle.h"
0009 #include "Gaudi/Property.h"
0010 #include "GaudiKernel/RndmGenerators.h"
0011
0012 #include "DD4hep/DD4hepUnits.h"
0013 #include "DDRec/CellIDPositionConverter.h"
0014 #include "DDRec/Surface.h"
0015 #include "DDRec/SurfaceManager.h"
0016
0017 #include <k4FWCore/DataHandle.h>
0018 #include <k4Interface/IGeoSvc.h>
0019
0020 // Event Model related classes
0021 //#include "GaudiExamples/MyTrack.h"
0022 #include "edm4eic/RawTrackerHitCollection.h"
0023 #include "edm4eic/TrackerHitCollection.h"
0024
0025 #include "DD4hep/DD4hepUnits.h"
0026
0027 namespace {
0028 inline double get_resolution(const double pixel_size) {
0029   constexpr const double sqrt_12 = 3.4641016151;
0030   return pixel_size / sqrt_12;
0031 }
0032 inline double get_variance(const double pixel_size) {
0033   const double res = get_resolution(pixel_size);
0034   return res * res;
0035 }
0036 } // namespace
0037
0038 namespace Jug::Reco {
0039
0040   /** Tracker hit reconstruction.
0041    *
0042    * \ingroup reco
0043    */
0044   class TrackerHitReconstruction : public Gaudi::Algorithm {
0045   private:
0046     Gaudi::Property<float> m_timeResolution{this, "timeResolution", 10}; // in ns
0047     mutable DataHandle<edm4eic::RawTrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
0048                                                                    this};
0049     mutable DataHandle<edm4eic::TrackerHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
0050                                                                  this};
0051
0052     /// Pointer to the geometry service
0053     SmartIF<IGeoSvc> m_geoSvc;
0054     std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_converter;
0055
0056   public:
0057     //  ill-formed: using Gaudi::Algorithm::GaudiAlgorithm;
0058     TrackerHitReconstruction(const std::string& name, ISvcLocator* svcLoc) : Gaudi::Algorithm(name, svcLoc) {
0059       declareProperty("inputHitCollection", m_inputHitCollection, "");
0060       declareProperty("outputHitCollection", m_outputHitCollection, "");
0061     }
0062
0063     StatusCode initialize() override {
0064       if (Gaudi::Algorithm::initialize().isFailure()) {
0065         return StatusCode::FAILURE;
0066       }
0067       m_geoSvc = service("GeoSvc");
0068       if (!m_geoSvc) {
0069         error() << "Unable to locate Geometry Service. "
0070                 << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0071         return StatusCode::FAILURE;
0072       }
0073       m_converter = std::make_shared<const dd4hep::rec::CellIDPositionConverter>(*(m_geoSvc->getDetector()));
0074       return StatusCode::SUCCESS;
0075     }
0076
0077     StatusCode execute(const EventContext&) const override {
0078       constexpr auto mm = dd4hep::mm;
0079       // input collection
0080       const auto* const rawhits = m_inputHitCollection.get();
0081       // Create output collections
0082       auto* rec_hits = m_outputHitCollection.createAndPut();
0083
0084       debug() << " raw hits size : " << std::size(*rawhits) << endmsg;
0085       for (const auto& ahit : *rawhits) {
0086         // debug() << "cell ID : " << ahit.cellID() << endmsg;
0087         auto pos = m_converter->position(ahit.getCellID());
0088         auto dim = m_converter->cellDimensions(ahit.getCellID());
0089
0090         if (msgLevel(MSG::VERBOSE)) {
0091           size_t i = 0;
0092           for (const auto& p : {pos.x(), pos.y(), pos.z()}) {
0093             verbose() << "position " << i++ << " [mm]: " << p / mm << endmsg;
0094           }
0095           verbose() << "dimension size: " << std::size(dim) << endmsg;
0096           for (size_t j = 0; j < std::size(dim); ++j) {
0097             verbose() << " - dimension " << j << " size: " << dim[j] << endmsg;
0098           }
0099         }
0100         // Note about variance:
0101         //    The variance is used to obtain a diagonal covariance matrix.
0102         //    Note that the covariance matrix is written in DD4hep surface coordinates,
0103         //    *NOT* global position coordinates. This implies that:
0104         //      - XY segmentation: xx -> sigma_x, yy-> sigma_y, zz -> 0, tt -> 0
0105         //      - XZ segmentation: xx -> sigma_x, yy-> sigma_z, zz -> 0, tt -> 0
0106         //      - XYZ segmentation: xx -> sigma_x, yy-> sigma_y, zz -> sigma_z, tt -> 0
0107         //    This is properly in line with how we get the local coordinates for the hit
0108         //    in the TrackerSourceLinker.
0109         rec_hits->create(
0110           ahit.getCellID(), // Raw DD4hep cell ID
0111           edm4hep::Vector3f{static_cast<float>(pos.x() / mm), static_cast<float>(pos.y() / mm), static_cast<float>(pos.z() / mm)},                    // mm
0112           edm4eic::CovDiag3f{get_variance(dim[0] / mm), get_variance(dim[1] / mm), // variance (see note above)
0113           std::size(dim) > 2 ? get_variance(dim[2] / mm) : 0.},
0114           static_cast<float>(ahit.getTimeStamp() / 1000), // ns
0115           m_timeResolution,                            // in ns
0116           static_cast<float>(ahit.getCharge() / 1.0e6),   // Collected energy (GeV)
0117           0.0F
0118         );                                       // Error on the energy
0119       }
0120       return StatusCode::SUCCESS;
0121     }
0122   };
0123   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0124   DECLARE_COMPONENT(TrackerHitReconstruction)
0125
0126 } // namespace Jug::Reco