algorithms/tracking/MPGDHitReconstruction.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2026 Whitney Armstrong, Sylvester Joosten, Wouter Deconinck, Dmitry Romanov, Yann Bedfer
0003
0004 #include "MPGDHitReconstruction.h"
0005
0006 #include <DD4hep/Detector.h>
0007 #include <DD4hep/IDDescriptor.h>
0008 #include <DD4hep/Objects.h>
0009 #include <DD4hep/Readout.h>
0010 #include <DDSegmentation/BitFieldCoder.h>
0011 #include <Evaluator/DD4hepUnits.h>
0012 #include <Math/GenVector/Cartesian3D.h>
0013 #include <Math/GenVector/DisplacementVector3D.h>
0014 // Access "algorithms:GeoSvc"
0015 #include <algorithms/geo.h>
0016 #include <algorithms/logger.h>
0017 #include <edm4eic/CovDiag3f.h>
0018 #include <edm4hep/Vector3f.h>
0019 #include <algorithm>
0020 #include <array>
0021 #include <cstddef>
0022 #include <cstdint>
0023 #include <iterator>
0024 #include <stdexcept>
0025 #include <tuple>
0026 #include <vector>
0027
0028 using namespace dd4hep;
0029
0030 namespace eicrecon {
0031
0032 void MPGDHitReconstruction::init() {
0033
0034   // ***** Parse IDDescriptor
0035   // Access id decoder
0036   const dd4hep::Detector* detector = m_geo.detector();
0037   if (m_cfg.readout.empty()) {
0038     throw std::runtime_error("Readout is empty");
0039   }
0040   try {
0041     m_id_dec = detector->readout(m_cfg.readout).idSpec().decoder();
0042   } catch (...) {
0043     critical(R"(Failed to load ID decoder for "{}" readout.)", m_cfg.readout);
0044     throw std::runtime_error("Failed to load ID decoder");
0045   }
0046   parseIDDescriptor();
0047 }
0048
0049 void MPGDHitReconstruction::process(const Input& input, const Output& output) const {
0050   using dd4hep::mm;
0051
0052   const auto [raw_hits] = input;
0053   auto [rec_hits]       = output;
0054
0055   // Reorder input raw_hits
0056   // Sort them in ascending order of channel# on a per SUBVOLUME * (p|n) basis.
0057   int nRawHits = raw_hits->size();
0058   std::vector<size_t> idcs(nRawHits);
0059   size_t idx(0);
0060   std::generate(idcs.begin(), idcs.end(), [&] { return idx++; });
0061   std::sort(idcs.begin(), idcs.end(), [&](int idxa, int idxb) {
0062     CellID cIDa = raw_hits->at(idxa).getCellID(), vIDa = cIDa & m_subVolBits;
0063     CellID cIDb = raw_hits->at(idxb).getCellID(), vIDb = cIDb & m_subVolBits;
0064     if (vIDa != vIDb)
0065       return vIDa < vIDb;
0066     int pna = (vIDa & m_pStripBit) ? 0 : 1, pnb = (vIDb & m_pStripBit) ? 0 : 1;
0067     if (pna != pnb)
0068       return pna < pnb;
0069     std::int16_t hIDa = pna == 0 ? cIDa >> m_coordOffsets[0] : cIDa >> m_coordOffsets[1];
0070     std::int16_t hIDb = pnb == 0 ? cIDb >> m_coordOffsets[0] : cIDb >> m_coordOffsets[1];
0071     return hIDa < hIDb;
0072   });
0073
0074   CellID prvID; // CellID of previous
0075   // Current cluster (cc): in the making or to be stored
0076   int currentPN;
0077   size_t currentNDims;
0078   Position clusPos; // cc: weighted position
0079   std::vector<double> clusDim(
0080       3); // cc: uncertainty = resolution along measurement axis, else weighted dimension
0081   double clusCharge(0); // cc: sum of charges
0082   double sW(0);         // cc: sum of Weights (= "clusCharge", as of 2025/12)
0083   double clusChMx(0);
0084   int clusChMxIdx(-1);  // cc: Max. charge and index of: determine timing
0085   int clusFirstIdx(-1); // cc: Index of 1st RawHit contributing
0086   // Loop on ordered raw_hits + a last iteration to store last cluster
0087   int jdx;
0088   for (jdx = 0, prvID = 0, currentPN = -1, sW = 0; jdx <= nRawHits; jdx++) {
0089     bool lastHit = jdx == nRawHits;
0090     bool newCluster;
0091     CellID cID;
0092     int idx;
0093     if (!lastHit) {
0094       idx                 = idcs[jdx];
0095       const auto& raw_hit = raw_hits->at(idx);
0096       cID                 = raw_hit.getCellID();
0097       if ((cID & m_subVolBits) != (prvID & m_subVolBits) ||
0098           (cID & m_pStripBit) != (prvID & m_pStripBit)) {
0099         newCluster = true;
0100       } else {
0101         // Cluster continuation? Require channel#+1
0102         std::int16_t prvHID, hID;
0103         if (cID & m_pStripBit) {
0104           prvHID = prvID >> m_coordOffsets[0];
0105           hID    = cID >> m_coordOffsets[0];
0106         } else {
0107           prvHID = prvID >> m_coordOffsets[1];
0108           hID    = cID >> m_coordOffsets[1];
0109         }
0110         newCluster = hID != prvHID + 1;
0111       }
0112       if (prvID && !newCluster) {
0113         // ***** ADD HIT TO CURRENT CLUSTER
0114         double weight = raw_hit.getCharge();
0115         clusPos += m_converter->position(cID) * weight;
0116         sW += weight;
0117         std::vector<double> dim = m_converter->cellDimensions(cID);
0118         size_t nDims            = std::size(dim);
0119         for (size_t i = 0; i < nDims; i++)
0120           clusDim[i] += dim[i] * weight;
0121         clusCharge += weight;
0122         if (weight > clusChMx) {
0123           clusChMxIdx = idx;
0124           clusChMx    = weight;
0125         }
0126         if (idx < clusFirstIdx) {
0127           clusFirstIdx = idx;
0128         }
0129         if (nDims != currentNDims) {
0130           critical(
0131               R"(RawHits from "{}" readout have different #dims: cellID 0x{:0>16x} = %u, cellID 0x{:0>16x} = %u.)",
0132               m_cfg.readout, prvID, currentNDims, cID, nDims);
0133           throw std::runtime_error("Inconsistent RawHits");
0134         }
0135         prvID = cID;
0136         continue; // ***** TRY AND ADD YET ANOTHER HIT
0137       }
0138     } else
0139       newCluster = false;
0140     if (prvID) {
0141       // ***** STORE (and log) CURRENT CLUSTER
0142       // Timing = timing of max. charge.
0143       const auto& clusHit = raw_hits->at(clusChMxIdx);
0144       double clusTime     = clusHit.getTimeStamp();
0145       // ID = ID of max. charge by convention; it's not very meaningful in fact.
0146       CellID clusID = clusHit.getCellID();
0147       // Averaging
0148       clusPos /= sW * mm;
0149       // Variance
0150       for (int i = 0; i < (int)currentNDims; i++) {
0151         if (i == currentPN) { // Measurement axis
0152           clusDim[i] = m_cfg.stripResolutions[currentPN];
0153         } else {
0154           constexpr const double sqrt_12 = 3.4641016151;
0155           clusDim[i] /= sW * sqrt_12;
0156         }
0157         clusDim[i] /= mm;
0158         clusDim[i] *= clusDim[i];
0159       }
0160       // >oO trace
0161       if (level() == algorithms::LogLevel::kDebug) {
0162         debug("cellID = 0x{:08x}, 0x{:08x}", clusID >> 32, (std::uint32_t)clusID);
0163         debug("position x={:.2f} y={:.2f} z={:.2f} [mm]: ", clusPos.x(), clusPos.y(), clusPos.z());
0164         debug("dimension size: {}", currentNDims);
0165         for (size_t j = 0; j < currentNDims; ++j) {
0166           debug(" - dimension {:<5} size: {:.2}", j, clusDim[j]);
0167         }
0168       }
0169
0170       // Note about variance:
0171       //    The variance is used to obtain a diagonal covariance matrix.
0172       //    Note that the covariance matrix is written in DD4hep surface coordinates,
0173       //    *NOT* global position coordinates. This implies that:
0174       //      - XY segmentation: xx -> sigma_x, yy-> sigma_y, zz -> 0, tt -> 0
0175       //      - XZ segmentation: xx -> sigma_x, yy-> sigma_z, zz -> 0, tt -> 0
0176       //      - XYZ segmentation: xx -> sigma_x, yy-> sigma_y, zz -> sigma_z, tt -> 0
0177       //    This is properly in line with how we get the local coordinates for the hit
0178       //    in the TrackerSourceLinker.
0179       // CartesianGridUV:
0180       // - Covariance matrix is along U and V. I(Y.B.) don't know whether this
0181       //  is what is meant by "DD4hep surface coordinates" above. And don't
0182       //  know whether it's in line w/ TrackerSourceLinker.
0183       // - But I do know that the matrix has to be rotated along X,Y (i.e.
0184       //  the reference frame local to the DD4hep volume) lest hit-to-track
0185       //  association fail completely. At present (2025/12), this is done in
0186       //  TrackerMeasurementFromHits. It cannot be done here, since it would
0187       //  result in a non-diagonal matrix (and a violently non-diagonal one:
0188       //  precision along one of U or V being very large for a strip hit).
0189       auto rec_hit = rec_hits->create(
0190           clusID, // Raw DD4hep cell ID of largest cluster's hit.
0191           edm4hep::Vector3f{static_cast<float>(clusPos.x()), static_cast<float>(clusPos.y()),
0192                             static_cast<float>(clusPos.z())}, // mm
0193           edm4eic::CovDiag3f{clusDim[0],
0194                              clusDim[1], // variance (see note above)
0195                              currentNDims > 2 ? clusDim[2] : 0.},
0196           static_cast<float>(clusTime / 1000.0),  // ns
0197           m_cfg.timeResolution,                   // in ns
0198           static_cast<float>(clusCharge / 1.0e6), // Collected energy (GeV)
0199           0.0F);                                  // Error on the energy
0200       // ********** REC <- RAW ASSOCIATION
0201       // - In EDM4eic, there's room for ONLY ONE ASSOCIATED RAW.
0202       // - Here we NEED MORE.
0203       // - Temporarily, simply put the earliest one. So that subsequent RawHits
0204       //  that should also be associated can easily be guessed.
0205       // - But imho(Y.B), "EDM4eic::TrackerHit" should be modified.
0206       const auto& firstHit = raw_hits->at(clusFirstIdx);
0207       rec_hit.setRawHit(firstHit);
0208     }
0209     if (newCluster) {
0210       const auto& raw_hit     = raw_hits->at(idx);
0211       double weight           = raw_hit.getCharge();
0212       clusPos                 = m_converter->position(cID) * weight;
0213       sW                      = weight;
0214       std::vector<double> dim = m_converter->cellDimensions(cID);
0215       size_t nDims            = std::size(dim);
0216       for (size_t i = 0; i < nDims; i++)
0217         clusDim[i] = dim[i] * weight;
0218       clusCharge   = weight;
0219       clusChMx     = weight;
0220       clusChMxIdx  = idx;
0221       currentPN    = (cID & m_pStripBit) ? 0 : 1;
0222       currentNDims = nDims;
0223       prvID        = cID;
0224       clusFirstIdx = idx;
0225     }
0226   }
0227 }
0228 void MPGDHitReconstruction::parseIDDescriptor() {
0229   // Parse IDDescriptor: Retrieve cellIDs of relevant fields.
0230   // (As an illustration, here is the IDDescriptor of CyMBaL (as of 2025/12):
0231   // <id>system:8,layer:4,module:12,sensor:2,strip:28:4,phi:-16,z:-16</id>.)
0232
0233   // - "m_(p|n)stripBit".
0234   debug(R"(Parsing IDDescriptor for "{}" readout)", m_cfg.readout);
0235   const char fieldName[] = "strip";
0236   CellID stripBits[2]    = {0, 0};
0237   try {
0238     m_id_dec->get(~((CellID)0x0), fieldName);
0239   } catch (const std::runtime_error& error) {
0240     critical(R"(No field "{}" in IDDescriptor of readout "{}".)", fieldName, m_cfg.readout);
0241     throw std::runtime_error("Invalid IDDescriptor");
0242   }
0243   const BitFieldElement& fieldElement = (*m_id_dec)[fieldName];
0244   // Coordinates are assigned specific bits by convention
0245   FieldID bits[2] = {0x1, 0x2};
0246   for (int coord = 0; coord < 2; coord++) {
0247     stripBits[coord] = bits[coord] << fieldElement.offset();
0248   }
0249   // CellIDs derived from above
0250   m_pStripBit = stripBits[0];
0251   m_nStripBit = stripBits[1];
0252   // Require coordinate fields to start @ bits 32 and 48: this ensures that
0253   // these fields fit into std::int16_t and is taken advantage of by debug
0254   // messages to cleanly separate coordinates from the rest of CellID.
0255   m_coordOffsets[0] = m_coordOffsets[1] = 64;
0256   for (int pn = 0; pn < 2; pn++) {
0257     std::string coordName;
0258     if (m_cfg.readout == "MPGDBarrelHits") {
0259       coordName = pn ? "z" : "phi";
0260     } else if (m_cfg.readout == "OuterMPGDBarrelHits") {
0261       coordName = pn ? "v" : "u";
0262     } else {
0263       coordName = pn ? "y" : "x";
0264     }
0265     try {
0266       m_id_dec->index(coordName);
0267     } catch (const std::runtime_error& error) {
0268       critical(R"(No "{}" field in IDDescriptor of readout "{}")", coordName, m_cfg.readout);
0269       throw std::runtime_error("Invalid IDDescriptor");
0270     }
0271     const BitFieldElement& fieldElement = (*m_id_dec)[coordName];
0272     int offset                          = fieldElement.offset();
0273     if (offset < m_coordOffsets[pn])
0274       m_coordOffsets[pn] = offset;
0275   }
0276   if (m_coordOffsets[0] != 32 || m_coordOffsets[1] != 48) {
0277     critical(R"(Coordinate fields in IDDescriptor of readout "{}" do not start @ bits 32 and 48.)",
0278              m_cfg.readout);
0279     throw std::runtime_error("Invalid IDDescriptor");
0280   }
0281   for (int i = 0; i < 32; i++)
0282     // Now that offsets have been checked, we can safely define "m_subVolBits"
0283     // to span [0,32[.
0284     m_subVolBits |= ((CellID)0x1) << i;
0285 }
0286
0287 } // namespace eicrecon