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