src/todo/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 "JugBase/DataHandle.h"
0020 #include "JugBase/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
0057   public:
0058     //  ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
0059     TrackerHitReconstruction(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
0060       declareProperty("inputHitCollection", m_inputHitCollection, "");
0061       declareProperty("outputHitCollection", m_outputHitCollection, "");
0062     }
0063
0064     StatusCode initialize() override {
0065       if (GaudiAlgorithm::initialize().isFailure()) {
0066         return StatusCode::FAILURE;
0067       }
0068       m_geoSvc = service("GeoSvc");
0069       if (!m_geoSvc) {
0070         error() << "Unable to locate Geometry Service. "
0071                 << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0072         return StatusCode::FAILURE;
0073       }
0074       return StatusCode::SUCCESS;
0075     }
0076
0077     StatusCode execute() 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_geoSvc->cellIDPositionConverter()->position(ahit.getCellID());
0088         auto dim = m_geoSvc->cellIDPositionConverter()->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         edm4eic::TrackerHit hit{ahit.getCellID(), // Raw DD4hep cell ID
0110                              {static_cast<float>(pos.x() / mm), static_cast<float>(pos.y() / mm),
0111                               static_cast<float>(pos.z() / mm)},                    // mm
0112                              {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};                                       // Error on the energy
0118         rec_hits->push_back(hit);
0119       }
0120       return StatusCode::SUCCESS;
0121     }
0122   };
0123   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0124   DECLARE_COMPONENT(TrackerHitReconstruction)
0125
0126 } // namespace Jug::Reco