EIC code displayed by LXR master/​EICrecon/​src/​algorithms/​digi/​MPGDTrackerDigi.cc	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-07-03 07:52:35

 // SPDX-License-Identifier: LGPL-3.0-or-later
 // Copyright (C) 2022 - 2025 Whitney Armstrong, Wouter Deconinck, Sylvester Joosten, Dmitry Romanov, Yann Bedfer
 
 /*
   Digitization specific to MPGDs.
   - What's special in MPGDs is their 2D-strip readout: i.e. simultaneous
    registering along two sets of coordinate strips.
   - The fact imposes strong constraints on digitization, stemming from DD4hep's
    provision that one and only one readout surface be associated to a sensitive
    volume. The one-and-only-one rule is enforced through to ACTS's TrackState,
    forbidding simple solutions exploiting MultiSegmentation to create two strip
    hits on a same surface.
   - Present solution is based on multiple sensitive SUBVOLUMES. FIVE of them:
     + TWO persistent ones, acting as READOUT SURFACES (i.e. to which raw hits
      are assigned, in both RealData and MC): very thin and sitting very close
      to mid-plane,
     + THREE HELPER ones: TWO RADIATORS, thick, serving temporarily to collect
      energy deposit and ONE REFERENCE, thin and sitting exactly at mid-plane.
      (The REFERENCE is not essential, but found to help rationalize the source
      code.)
   - Each SUBVOLUME has a distinct ID, w/ a distinctive "STRIP" FIELD, but all
    share a common XML <readout> thanks to a MultiSegmentation based on the
    "STRIP" FIELD.
   - Several subHits are created (up to five, one per individual SUBVOLUME)
    corresponding to a SINGLE I[ONISATION]+A[MPLIFICATION] process, inducing in
    turn CORRELATED signals on two coordinates.
   - These subHits are EXTENDED so as to span the full sensitive volume and be
    on the sensitive surface of the REFERENCE SUBVOLUME (its mid-plane), and
    hence reconstitute the proper hit position of the TRAVERSING particle.
     This, for particles that are determined to be indeed traversing.
   - We then want to preserve the I+A correlation when accumulating hits on the
    readout channels. This could be obtained by accumulating directly subHits (
    as is done in "SiliconTrackerDigi", as of commit #151496af). For such a
    scheme to preserve CORRELATION, we would have to collect all related subHits
    (i.e. those deemed to originate from a given I+A) and apply to them a common
    smearing (in amplitude and time). Here, we go one step further and COALESCE
    the related subHits. Our guideline can be formulated as: reproduce what we
    would get would there be a single volume. This comes with a cost of extra
    complication in the source code, but a limited one. This should go on w/
    the common smearing simulating I+A, independent smearings simulating CHARGE
    SHARING and CHARGE SPREADING. Then accumulation channel by channel. Then,
    possibly, apply smearing simulating FE electronics.
   - In any case, we therefore have a DOUBLE LOOP ON INPUT SimTrackerHits, to
    collect those subHits arising from the same I+A. I.e. SAME P[ARTICLE],
    M[ODULE] (a particle can fire two overlapping modules, not to be mixed) and
    O[RIGIN] (direct or via secondary).
   - The double loop is optimized for speed (by keeping track of the start index
    of the second loop) in order to not waste time on high multiplicity events.
   - We assume subHits originating from the same ionization process to come in
    sequences unbroken by any intertwining.
   - A number of checks are conducted before undertaking subHit COALESCENCE:
     I) SubHits should be TRAVERSING, i.e. exiting and then entering through
     opposite walls, as opposed to exiting or entering through the edge or
     dying/being-born within their SUBVOLUME, see methods "(c|b)Traversing".
     II) SubHits have SAME PMO, see "samePMO".
     III) SubHits extrapolate to a common point, see "outInDistance".
     All these checks are done to be on the safe side. A subset of them, hinging
    on (III), should be sufficient. But there are so many cases to take into
    account that it's difficult to settle on a minimal subset.
   - Evaluation:
     + MIPs are found to be properly handled, whether a cutoff on deposited
      energy is applied or not: hits are COALESCED when needed. They are
      assigned to mid-plane, except for few cases, see "flagUnexpected".
     + For lower energy particles, there are at times ambiguous cases where
      seemingly related subHits turn out to not be COALESCED, because they fail
      to pass above-mentioned checks due to their erratic trajectories.
       Note that even if they are genuinely related, this is not dramatic:
      subHits will still be directly accumulated on their common (within pitch
      precision) cell, only that we miss the correlation.
     + Special cases of i) particle exiting and reEntering the same SUBVOLUME
      in the cylindrical setup, ii) SimHit positioned outside the SUBVOLUME (
      cylindrical setup once again) it is assigned to, are catered for. In case
      (i), subHits are COALESCED but not EXTENDED.
   - Imperfections:
     + In order to validate (III), one has to decide on a TOLERANCE. The ideal
      would be to do this based on multiscattering. Instead, what's used are
      guesstimates, and a recipe to RELAX TOLERANCE for low energy stuff.
   - SIMULATION
     In addition to DIGITIZATION proper, the method involves a SIMULATION of
    some of the inner workings of the MPGD, viz. AMPLIFICATION, SIGNAL SHARING
    and SPREADING.
     + This SIMULATION relies on parameterizations obtained from beam tests.
     + Present version is simplistic: 2-hit clusters (except when on the edge),
      with identical timing and amplitude, and possibly, BELOW-THRESHOLD charge.
   - DIGITIZATION:
     + It involves the production of 2-hit clusters, called here CLUSTERIZATION.
     + It otherwise follows the standard steps. Remains to agree on the handling
      of the discrimination threshold, see Issue #1722.
  */
 
 #include "MPGDTrackerDigi.h"
 
 #include <DD4hep/Alignments.h>
 #include <DD4hep/DetElement.h>
 #include <DD4hep/IDDescriptor.h>
 #include <DD4hep/Objects.h>
 #include <DD4hep/Readout.h>
 #include <DD4hep/Shapes.h>
 #include <DD4hep/VolumeManager.h>
 #include <DD4hep/detail/SegmentationsInterna.h>
 #include <DDSegmentation/BitFieldCoder.h>
 #include <DDSegmentation/CartesianGridUV.h>
 #include <DDSegmentation/CartesianGridXY.h>
 #include <DDSegmentation/CylindricalGridPhiZ.h>
 #include <DDSegmentation/MultiSegmentation.h>
 #include <DDSegmentation/Segmentation.h>
 #include <Evaluator/DD4hepUnits.h>
 #include <Math/GenVector/Cartesian3D.h>
 #include <Math/GenVector/DisplacementVector3D.h>
 #include <Parsers/Primitives.h>
 #include <TGeoMatrix.h>
 #include <TMath.h>
 // Access "algorithms:GeoSvc"
 #include <algorithms/geo.h>
 #include <algorithms/logger.h>
 #include <edm4eic/unit_system.h>
 #include <edm4hep/MCParticleCollection.h>
 #include <edm4hep/Vector3d.h>
 #include <edm4hep/Vector3f.h>
 #include <fmt/format.h>
 #include <podio/detail/Link.h>
 #include <podio/detail/LinkCollectionImpl.h>
 #include <algorithm>
 #include <array>
 #include <cmath>
 #include <cstdint>
 #include <gsl/pointers>
 #include <gsl/util>
 #include <initializer_list>
 #include <iterator>
 #include <map>
 #include <memory>
 #include <random>
 #include <stdexcept>
 #include <tuple>
 #include <unordered_map>
 #include <utility>
 #include <vector>
 
 #include "algorithms/digi/MPGDTrackerDigiConfig.h"
 
 using namespace dd4hep;
 
 namespace eicrecon {
 
 void MPGDTrackerDigi::init() {
   // Access id decoder
   m_detector = algorithms::GeoSvc::instance().detector();
 
   if (m_cfg.readout.empty()) {
     throw std::runtime_error("Readout is empty");
   }
   try {
     m_seg    = m_detector->readout(m_cfg.readout).segmentation();
     m_id_dec = m_detector->readout(m_cfg.readout).idSpec().decoder();
   } catch (const std::runtime_error&) {
     critical(R"(Failed to load ID decoder for "{}" readout.)", m_cfg.readout);
     throw std::runtime_error("Failed to load ID decoder");
   }
 
   // IDDescriptor and Segmentation
   parseIDDescriptor();
   parseSegmentation();
 
   // Ordering of SUBVOLUMES (based on "STRIP" FIELD)
   m_stripRank = [&](CellID vID) {
     CellID sID = vID & m_stripBits;
     for (int rank = 0; rank < 5; rank++) {
       if (sID == m_stripIDs[rank]) {
         return rank;
       }
     }
     return -1;
   };
   m_orientation = [&](CellID vID, CellID vJD) {
     int ranki = m_stripRank(vID);
     int rankj = m_stripRank(vJD);
     if (rankj > ranki) {
       return +1;
     }
     if (rankj < ranki) {
       return -1;
     }
     return 0;
   };
   m_isUpstream = [](int orientation, unsigned int status) {
     // Outgoing particle exits...
     bool isUpstream =
         (orientation < 0 && (status & 0x2)) ||           // ...lower wall
         (orientation > 0 && (status & 0x8)) ||           // ...upper wall
         (orientation == 0 && (status & 0x102) == 0x102); // ...lower wall and can reEnter
     return isUpstream;
   };
   m_isDownstream = [](int orientation, unsigned int status) {
     // Incoming particle enters...
     bool isDownstream =
         (orientation > 0 && (status & 0x1)) ||           // ...lower wall
         (orientation < 0 && (status & 0x4)) ||           // ...upper wall
         (orientation == 0 && (status & 0x101) == 0x101); // ...lower wall and can be reEntering
     return isDownstream;
   };
   // RELAXED TOLERANCE
   m_toleranceFactor = [](double P) {
     int factor = 0;
     if (P < 1 * dd4hep::MeV) {
       factor = 4;
     } else if (P < 10 * dd4hep::MeV) {
       factor = 2;
     } else {
       factor = 1;
     }
     return factor;
   };
 
   // CLUSTERIZATION
   m_binToPosition = [](FieldID bin, double cellSize, double offset) {
     return bin * cellSize + offset;
   };
 }
 
 // Interfaces
 void getLocalPosMom(const edm4hep::SimTrackerHit& sim_hit, const TGeoHMatrix& toModule,
                     double* lpos, double* lmom);
 bool cExtrapolate(const double* lpos, const double* lmom, // Input subHit
                   double rT,                              // Target radius
                   double* lext);                          // Extrapolated position @ <rT>
 double getRef2Cur(DetElement refVol, DetElement curVol);
 bool bExtrapolate(const double* lpos, const double* lmom, // Input subHit
                   double zT,                              // Target Z
                   double* lext);                          // Extrapolated position @ <zT>
 std::string inconsistency(const edm4hep::EventHeader& event, unsigned int status, CellID cID,
                           const double* lpos, const double* lmom);
 std::string oddity(const edm4hep::EventHeader& event, unsigned int status, double dist, CellID cID,
                    const double* lpos, const double* lmom, CellID cJD, const double* lpoj,
                    const double* lmoj);
 double outInDistance(int shape, int orientation, double lintos[][3], double louts[][3],
                      double* lmom, double* lmoj);
 void flagUnexpected(const edm4hep::EventHeader& event, int shape, double expected,
                     const edm4hep::SimTrackerHit& sim_hit, double* lpini, double* lpend,
                     double* lpos, double* lmom);
 
 void MPGDTrackerDigi::process(const MPGDTrackerDigi::Input& input,
                               const MPGDTrackerDigi::Output& output) const {
 
   const auto [headers, sim_hits] = input;
 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
   auto [raw_hits, links, associations] = output;
 #else
   auto [raw_hits, associations] = output;
 #endif
 
   // local random generator
   auto seed = m_uid.getUniqueID(*headers, name());
   std::default_random_engine generator(seed);
   std::normal_distribution<double> gaussian;
 
   // Maps of unique cellIDs with temporary structure RawHit
   std::unordered_map<std::uint64_t, edm4eic::MutableRawTrackerHit> cell_hit_maps[2];
   // A map of strip cellIDs with vector of contributing cellIDs
   std::map<std::uint64_t, std::vector<std::uint64_t>> stripID2cIDs;
   // Prepare for strip segmentation
   const VolumeManager& volman = m_detector->volumeManager();
 
   // Reference to event, to be used to document error messages
   // (N.B.: I don't know how to properly handle these "headers": may there
   // be more than one? none?...)
   const edm4hep::EventHeader& header = headers->at(0);
 
   // *************** LOOP ON sim_hits
   size_t sim_size = sim_hits->size();
   for (int idx = 0; idx < (int)sim_size; idx++) {
     const edm4hep::SimTrackerHit& sim_hit = sim_hits->at(idx);
 
     // ***** TIME SMEARING
     // - Simplistic treatment.
     // - A more realistic one would have to distinguish a smearing common to
     //  both coordinates of the 2D-strip readout (mainly due to the drifting of
     //  the leading primary electrons of the I+A process) from other smearing
     //  effects, specific to each coordinate.
     double time_smearing = gaussian(generator) * m_cfg.timeResolution;
 
     // ***** REFERENCE SUBVOLUME
     CellID refID      = sim_hit.getCellID() & m_moduleBits;
     DetElement refVol = volman.lookupDetElement(refID);
     // ***** COALESCE ALL MUTUALLY CONSISTENT SUBHITS
     //       EXTEND TRAVERSING SUBHITS
     // - Needed because we want to preserve the correlation between 'p' and
     //  'n' strip hits resulting from a given I+A process (which is lost when
     //  one accumulates hits independently based on cellID).
     double lpos[3], eDep, time;
     std::vector<std::uint64_t> cIDs;
     const auto& shape = refVol.solid();
     if (std::string_view{shape.type()} == "TGeoTubeSeg") {
       // ********** TUBE GEOMETRY
       if (!cCoalesceExtend(input, idx, cIDs, lpos, eDep, time))
         continue;
     } else if (std::string_view{shape.type()} == "TGeoBBox") {
       // ********** BOX GEOMETRY
       if (!bCoalesceExtend(input, idx, cIDs, lpos, eDep, time))
         continue;
     } else {
       critical(R"(Bad input data: CellID {:x} has invalid shape "{}")", refID, shape.type());
       throw std::runtime_error(R"(Inconsistency: Inappropriate SimHits fed to "MPGDTrackerDigi".)");
     }
 
     // ***** 2D-position on sensitive surface
     double surfPos[2];
     Position locPos;
     if (std::string_view{shape.type()} == "TGeoTubeSeg") {
       // Sensitive surface radius = REFERENCE VOLUME radius
       const Tube& tRef = refVol.solid();
       double R         = (tRef.rMin() + tRef.rMax()) / 2;
       double phi       = atan2(lpos[1], lpos[0]);
       surfPos[0]       = phi * R;
       surfPos[1]       = lpos[2];
       locPos           = Position(R * cos(phi), R * sin(phi), lpos[2]);
     } else {
       locPos = Position(lpos[0], lpos[1], lpos[2]);
       if (m_gridAngle != 0.0) { // Transform to strip frame
         dd4hep::DDSegmentation::Vector3D position;
         position.X = lpos[0];
         position.Y = lpos[1];
         position   = RotationZ(m_gridAngle) * position;
         surfPos[0] = position.X;
         surfPos[1] = position.Y;
       } else {
         surfPos[0] = lpos[0];
         surfPos[1] = lpos[1];
       }
     }
 
     // ********** LOOP ON p|n STRIPS
     for (int pn = 0; pn < 2; pn++) {
       // ***** CLUSTERIZATION
       // Cluster = (CellID, energy fraction)
       Cluster cluster;
       int status = get2HitCluster(refID, locPos, surfPos, pn, generator, cluster);
       if (status != 0) {
         CellID vIDs = (std::uint32_t)refID;
         CellID hIDs = refID >> 32;
         error(R"(SimHit (= 0x{:08x}, 0x{:08x}, {:.2f},{:.2f} cm) beyond limits of {}Strips.)", hIDs,
               vIDs, surfPos[0] / cm, surfPos[1] / cm, (pn != 0) ? 'n' : 'p');
       }
       // ***** DEBUGGING INFO
       if (level() >= algorithms::LogLevel::kDebug) {
         std::string sCellID = (pn != 0) ? "cellIDn" : "cellIDp";
         if (pn == 0) {
           debug("  =>=>");
         }
         for (auto clusterHit : cluster) {
           CellID cID = clusterHit.first;
           CellID hID = cID >> 32;
           CellID vID = (std::uint32_t)cID;
           debug("Hit {} = 0x{:08x}, 0x{:08x}", sCellID, hID, vID);
           debug("  edep = {:.0f} [eV]", clusterHit.second * eDep / eV);
         }
       }
       // ***** APPLY THRESHOLD / STORE (sim_hit -> stripIDs) / ACCUMULATION
       // (Note: Threshold is applied first. See issue #1722 in
       // "https://github.com/eic/EICrecon/issues/1722".)
       std::unordered_map<std::uint64_t, edm4eic::MutableRawTrackerHit>& cell_hit_map =
           gsl::at(cell_hit_maps, pn);
       double g = m_cfg.gain;
       for (auto clusterHit : cluster) {
         CellID cID = clusterHit.first;
         double f   = clusterHit.second;
         // Threshold
         // - Let's apply same threshold to the two components of the cluster.
         // - This, so that cluster position be preserved at Hit Reconstruction
         //  time (when clustering is performed).
         // => This has the obvious drawback of creating BELOW-THRESHOLD RawHits.
         if (eDep < m_cfg.threshold) {
           debug("  eDep {:.2f} is below threshold of {:.2f} [keV]", eDep, m_cfg.threshold / keV);
           continue;
         }
         stripID2cIDs[cID]   = cIDs;
         double result_time  = time + time_smearing;
         auto hit_time_stamp = (std::int32_t)(result_time * 1e3);
         if (!cell_hit_map.contains(cID)) {
           // This cell doesn't have hits
           cell_hit_map[cID] = {
               cID, (std::int32_t)std::llround(eDep * 1e6 * f * g),
               hit_time_stamp // ns->ps
           };
         } else {
           // There is previous values in the cell
           auto& hit = cell_hit_map[cID];
           debug("  Hit already exists in cell ID={}, prev. hit time: {}", cID, hit.getTimeStamp());
 
           // keep earliest time for hit
           hit.setTimeStamp(std::min(hit_time_stamp, hit.getTimeStamp()));
 
           // sum deposited energy
           auto charge = hit.getCharge();
           hit.setCharge(charge + (std::int32_t)std::llround(eDep * 1e6 * f * g));
         }
       }
     } // End loop on strip = p,n
   } // End loop on sim_hit's
 
   // ***** RawHit INSTANTIATION AND RawHit<-SimHits ASSOCIATION:
   for (auto& cell_hit_map : cell_hit_maps) {
     for (auto item : cell_hit_map) {
       raw_hits->push_back(item.second);
       CellID stripID = item.first;
       const auto is  = stripID2cIDs.find(stripID);
       if (is == stripID2cIDs.end()) {
         error(R"(Inconsistency: CellID {:x} not found in "stripID2cIDs" map)", stripID);
         throw std::runtime_error(R"(Inconsistency in the handling of "stripID2cIDs" map)");
       }
       std::vector<std::uint64_t> cIDs = is->second;
       for (CellID cID : cIDs) {
         for (const auto& sim_hit : *sim_hits) {
           if (sim_hit.getCellID() == cID) {
 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
             // create link
             auto link = links->create();
             link.setFrom(item.second);
             link.setTo(sim_hit);
             link.setWeight(1.0);
 #endif
             // set association
             auto hitassoc = associations->create();
             hitassoc.setWeight(1.0);
             hitassoc.setRawHit(item.second);
             hitassoc.setSimHit(sim_hit);
           }
         }
       }
     }
   }
 }
 
 //  The whole of MPGDTrackerDigi relies on a number of assumptions concerning
 // the structure of the MPGD detector and its subdivision into SUBVOLUMES.
 // These assumptions imply in turn a specific segmentation scheme. In
 // particular, a MultiSegmentation allowing to navigate among the SUBVOLUMES.
 // This MultiSegmentation, based on a "strip" discriminator, cf. IDDescriptor,
 // can be embedded in yet another MultiSegmentation (case of CyMBaL).
 //  On the other hand, the CLUSTERIZATION procedure requires accessing
 // segmentation parameters (pitch, offset, etc...), which in turn requires a
 // navigation through the multiple levels of MultiSegmentation. In order to
 // simplify this navigation, we decide to impose the following strict schemes:
 // - CyMBaL (identified by it <readout> name, viz.: "MPGDBarrelHits"):
 //   + MultiSegmentation discriminating on "sector",
 //   + MultiSegmentation discriminating on "strip",
 //   + CylindricalGridPhiZ, with 'p' = "phi" and 'n' = "z".
 // - OuterBarrel (""OuterMPGDBarrelHits"):
 //   + MultiSegmentation discriminating on "strip",
 //   + CartesianGridUV, with 'p' = "u" and 'n' = "v".
 // - Else:
 //   + MultiSegmentation discriminating on "strip",
 //   + CartesianGridXY, with 'p' = "x" and 'n' = "y".
 //  Additional restrictions:
 // I) Concerning parameters:
 //  Resolution sigma must be small enough compared to pitch that smearing does
 // not bring us beyond closest neighbor.
 // II) Concerning IDDescriptor:
 //  The two coordinate fields span [32,48[ and [48,64[, cf. "(in|de)crementID".
 // => Let's check.
 void MPGDTrackerDigi::parseIDDescriptor() {
   // Parse IDDescriptor: Retrieve CellIDs of relevant fields.
   // (As an illustration, here is the IDDescriptor of CyMBaL (as of 2025/11):
   // <id>system:8,layer:4,module:12,sensor:2,strip:28:4,phi:-16,z:-16</id>.)
 
   // - "m_volumeBits" (Volume = CellID excluding channel specification)
   // - "m_stripBits".
   // - "m_sensorStripBits"
   debug(R"(Parsing IDDescriptor for "{}" readout)", m_cfg.readout);
   CellID sensorBits = 0;
   for (int field = 0; field < 5; field++) {
     const char* fieldName = m_fieldNames[field];
     CellID fieldID        = 0;
     try {
       fieldID = m_id_dec->get(~((CellID)0x0), fieldName);
     } catch (const std::runtime_error& error) {
       critical(R"(No field "{}" in IDDescriptor of readout "{}".)", fieldName, m_cfg.readout);
       throw std::runtime_error("Invalid IDDescriptor");
     }
     const BitFieldElement& fieldElement = (*m_id_dec)[fieldName];
     CellID fieldBits                    = fieldID << fieldElement.offset();
     m_volumeBits |= fieldBits;
     // - "m_stripBits".
     if (std::string_view{fieldName} == "strip") {
       m_stripBits = fieldBits;
       // SUBVOLUMES are assigned specific bits by convention
       CellID bits[5] = {0x3, 0x1, 0, 0x2, 0x4};
       for (int subVolume = 0; subVolume < 5; subVolume++) {
         m_stripIDs[subVolume] = bits[subVolume] << fieldElement.offset();
       }
     }
     // - "m_sensorStripBits"
     else if (m_cfg.readout == "MPGDBarrelHits" && std::string_view{fieldName} == "sensor") {
       sensorBits     = fieldBits;
       m_sensorOffset = fieldElement.offset();
     }
   }
   // CellIDs derived from above
   m_moduleBits = m_volumeBits & ~m_stripBits;
   m_pStripBit  = m_stripIDs[1];
   m_nStripBit  = m_stripIDs[3];
   // Sensor|strip mask: Will serve to index the parameter map.
   if (m_cfg.readout == "MPGDBarrelHits") {
     m_sensorStripBits = sensorBits | m_stripBits;
   } else {
     m_sensorStripBits = m_stripBits;
   }
   // Get coordinate field and index.
   // Require coordinate fields to start @ bits 32 and 48. This is taken
   // advantage of by "(in|de)crementID" (and by debug messages, to cleanly
   // separate coordinates from the rest of CellID).
   int coordOffsets[2] = {64, 64};
   for (int pn = 0; pn < 2; pn++) {
     std::string coordName;
     if (m_cfg.readout == "MPGDBarrelHits") {
       coordName = pn ? "z" : "phi";
     } else if (m_cfg.readout == "OuterMPGDBarrelHits") {
       coordName = pn ? "v" : "u";
     } else {
       coordName = pn ? "y" : "x";
     }
     try {
       m_stripIndices[pn] = m_id_dec->index(coordName);
     } catch (const std::runtime_error& error) {
       critical(R"(No "{}" field in IDDescriptor of readout "{}".)", coordName, m_cfg.readout);
       throw std::runtime_error("Invalid IDDescriptor");
     }
     const BitFieldElement& fieldElement = (*m_id_dec)[coordName];
     int offset                          = fieldElement.offset();
     m_stripIncs[pn]                     = ((CellID)0x1) << offset;
     if (offset < coordOffsets[pn])
       coordOffsets[pn] = offset;
   }
   if (coordOffsets[0] != 32 || coordOffsets[1] != 48) {
     critical(R"(Coordinate fields in IDDescriptor of readout "{}" do not start @ bits 32 and 48.)",
              m_cfg.readout);
     throw std::runtime_error("Invalid IDDescriptor");
   }
 }
 void MPGDTrackerDigi::parseSegmentation() {
   debug(R"(Find valid "MultiSegmentation" for "{}" readout.)", m_cfg.readout);
 
   // Local function: Check restriction (I).
   std::function<void(int, double)> checkResolutionVsPitch = [&](int pn, double pitch) {
     double sigma = m_cfg.stripResolutions[pn];
     if (m_truncation * sigma > pitch) {
       critical(R"(stripResolutions[{}] (= {}) too large for pitch (={}) of "{}" readout.)", pn,
                sigma, pitch, m_cfg.readout);
       throw std::runtime_error("Space resolution parameter too large");
     }
     return;
   };
   using Segmentation               = dd4hep::DDSegmentation::Segmentation;
   const Segmentation* segmentation = m_seg->segmentation;
   // Retrieve <segmentation> parameters: pitch, offset, min, max, index.
   unsigned int required = 0, fulfilled = 0;
   if (segmentation->type() == "MultiSegmentation") {
     using MultiSegmentation = dd4hep::DDSegmentation::MultiSegmentation;
     const auto* multiSeg    = dynamic_cast<const MultiSegmentation*>(segmentation);
     StripParameters pars;
     if (m_cfg.readout == "MPGDBarrelHits") {
       // ********** SPECIFIC CASE: "MPGDBarrelHits"
       required              = 0x3 | m_pStripBit | m_nStripBit; // 0x3 sectors | stripBits
       unsigned int innerBit = 0x0, outerBit = 0x1;
       for (unsigned int sectorBit : {innerBit, outerBit}) {
         CellID sensorID               = ((CellID)sectorBit) << m_sensorOffset;
         const Segmentation& sectorSeg = multiSeg->subsegmentation(sensorID);
         if (multiSeg->type() != "MultiSegmentation")
           continue;
         const auto& stripMultiSeg = dynamic_cast<const MultiSegmentation&>(sectorSeg);
         for (CellID stripID : {m_pStripBit, m_nStripBit}) {
           const Segmentation& stripSeg = stripMultiSeg.subsegmentation(stripID);
           if (stripSeg.type() != "CylindricalGridPhiZ") {
             critical(
                 R"(Segmentation type for "{}" readout = "{}", whereas expected = "CylindricalGridPhiZ".)",
                 m_cfg.readout, stripSeg.type());
             continue;
           }
           const dd4hep::DDSegmentation::CylindricalGridPhiZ& gridPhiZ =
               dynamic_cast<const dd4hep::DDSegmentation::CylindricalGridPhiZ&>(stripSeg);
           fulfilled |= sectorBit == innerBit ? 0x1 : 0x2;
           fulfilled |= stripID;
           // Parameters -> StripParameters "pars"
           double radius = gridPhiZ.radius();
           int pn;
           std::string stripName;
           if (stripID == m_pStripBit) {
             pars.pitch  = gridPhiZ.gridSizePhi() * radius;
             pars.offset = gridPhiZ.offsetPhi() * radius;
             pn          = 0;
           } else {
             pars.pitch  = gridPhiZ.gridSizeZ();
             pars.offset = gridPhiZ.offsetZ();
             pn          = 1;
           }
           pars.index = m_stripIndices[pn];
           // Check restriction (I)
           checkResolutionVsPitch(pn, pars.pitch);
           int nStrips = m_cfg.stripNumbers[pn];
           // The center of the Cartesian grid is in the center of the detector, see, e.g.:
           //  "https://github.com/HEP-FCC/FCCDetectors/blob/main/doc/DD4hepInFCCSW.md".
           // I(Y.B.) recommend an offset = pitch/2, so that layout be symmetric.
           pars.max = (nStrips / 2 - .5) * pars.pitch + pars.offset;
           pars.min = (-nStrips / 2 - .5) * pars.pitch + pars.offset;
           // "pars" stored in "m_stripParameters" map.
           CellID sensorStripID             = sensorID | stripID;
           m_stripParameters[sensorStripID] = pars;
         }
       }
     } else if (m_cfg.readout == "OuterMPGDBarrelHits") {
       // ********** SPECIFIC CASE: "OuterMPGDBarrelHits"
       required = m_pStripBit | m_nStripBit;
       double gridAngles[2];
       for (unsigned int stripID : {m_pStripBit, m_nStripBit}) {
         const dd4hep::DDSegmentation::Segmentation& stripSeg = multiSeg->subsegmentation(stripID);
         if (stripSeg.type() != "CartesianGridUV") {
           critical(
               R"(Segmentation type for "{}" readout = "{}", whereas expected = "CartesianGridUV".)",
               m_cfg.readout, stripSeg.type());
           continue;
         }
         const dd4hep::DDSegmentation::CartesianGridUV& gridUV =
             dynamic_cast<const dd4hep::DDSegmentation::CartesianGridUV&>(stripSeg);
         fulfilled |= stripID;
         // Parameters -> StripParameters "pars"
         int pn;
         std::string stripName;
         if (stripID == m_pStripBit) {
           pars.pitch  = gridUV.gridSizeU();
           pars.offset = gridUV.offsetU();
           pn          = 0;
         } else {
           pars.pitch  = gridUV.gridSizeV();
           pars.offset = gridUV.offsetV();
           pn          = 1;
         }
         pars.index     = m_stripIndices[pn];
         gridAngles[pn] = gridUV.gridAngle();
         // Check restriction (I)
         checkResolutionVsPitch(pn, pars.pitch);
         int nStrips = m_cfg.stripNumbers[pn];
         pars.max    = (nStrips / 2 - .5) * pars.pitch + pars.offset;
         pars.min    = (-nStrips / 2 - .5) * pars.pitch + pars.offset;
         // "pars" stored in "m_stripParameters" map.
         m_stripParameters[stripID] = pars;
       }
       if (fabs(gridAngles[0] - gridAngles[1]) > 1e-6) {
         critical(
             R"(Inconsistent gridAngles of "CartesianGridUV" for readout "{}": 'p' = {:.6f}, 'n' = {:.6f}.)",
             m_cfg.readout, gridAngles[0], gridAngles[1]);
         fulfilled = 0;
       }
       m_gridAngle = gridAngles[0];
     } else {
       // ********** DEFAULT: "CartesianGridXY" Segmentation assumed
       required = m_pStripBit | m_nStripBit;
       for (unsigned int stripID : {m_pStripBit, m_nStripBit}) {
         const dd4hep::DDSegmentation::Segmentation& stripSeg = multiSeg->subsegmentation(stripID);
         if (stripSeg.type() != "CartesianGridXY") {
           critical(
               R"(Segmentation type for "{}" readout = "{}", whereas expected = "CartesianGridXY".)",
               m_cfg.readout, stripSeg.type());
           continue;
         }
         const dd4hep::DDSegmentation::CartesianGridXY& gridXY =
             dynamic_cast<const dd4hep::DDSegmentation::CartesianGridXY&>(stripSeg);
         fulfilled |= stripID;
         // Parameters -> StripParameters "pars"
         int pn;
         std::string stripName;
         if (stripID == m_pStripBit) {
           pars.pitch  = gridXY.gridSizeX();
           pars.offset = gridXY.offsetX();
           pn          = 0;
         } else {
           pars.pitch  = gridXY.gridSizeY();
           pars.offset = gridXY.offsetY();
           pn          = 1;
         }
         pars.index = m_stripIndices[pn];
         // Check restriction (I)
         checkResolutionVsPitch(pn, pars.pitch);
         int nStrips = m_cfg.stripNumbers[pn];
         pars.max    = (nStrips / 2 - .5) * pars.pitch + pars.offset;
         pars.min    = (-nStrips / 2 - .5) * pars.pitch + pars.offset;
         // "pars" stored in "m_stripParameters" map.
         m_stripParameters[stripID] = pars;
       }
     }
   } else {
     critical(
         R"(Segmentation type for "{}" readout = "{}", whereas expected = "MultiSegmentation".)",
         m_cfg.readout, segmentation->type());
   }
   if (fulfilled != required) {
     critical(R"(Error retrieving Segmentation parameters for "{}" readout.)", m_cfg.readout);
     throw std::runtime_error("Error retrieving Segmentation parameters");
   }
 }
 
 // ***** COALESCE subHits with same PMO
 //       EXTEND hits (subHits or coalesced hits) to full sensitive volume
 // - Input = Elementary subHit, specified as index into collection of SimHits.
 // - Output = Coalesced/extended hit, specified by:
 //  + list of cellIDs of elementary subHits contributing,
 //  + local position,
 //  + EDep,
 //  + time.
 //   Given that all segmentation classes foreseen for MPGDs ("CartesianGrid.."
 //  for Outer and EndCaps, "CylindricalGridPhiZ" for "CyMBaL") disregard the
 //  _global_ position argument to "dd4hep::Segmentation::cellID", we need
 //  the _local_ position and only that.
 // - Also returned: updated index.
 bool MPGDTrackerDigi::cCoalesceExtend(const Input& input, int& idx,
                                       std::vector<std::uint64_t>& cIDs, double* lpos, double& eDep,
                                       double& time) const {
   const auto [headers, sim_hits]        = input;
   const edm4hep::EventHeader& header    = headers->at(0);
   const edm4hep::SimTrackerHit& sim_hit = sim_hits->at(idx);
   CellID vID                            = sim_hit.getCellID() & m_volumeBits;
   CellID refID                          = vID & m_moduleBits; // => The REFERENCE SUBVOLUME
   const VolumeManager& volman           = m_detector->volumeManager();
   DetElement refVol                     = volman.lookupDetElement(refID);
   // TGeoHMatrix: In order to avoid a "dangling-reference" warning, let's take
   // a copy of the matrix instead of a reference to it.
   const TGeoHMatrix toRefVol = refVol.nominal().worldTransformation();
   double lmom[3];
   getLocalPosMom(sim_hit, toRefVol, lpos, lmom);
   const double edmm = edm4eic::unit::mm, ed2dd = dd4hep::mm / edmm;
   using dd4hep::mm;
   using edm4eic::unit::eV, edm4eic::unit::GeV;
   // Hit in progress
   eDep = sim_hit.getEDep();
   time = sim_hit.getTime();
   // Get VOLUME parameters
   const Tube& tRef = refVol.solid();
   double dZ        = tRef.dZ();
   // phi?
   // In "https://root.cern.ch/root/html534/guides/users-guide/Geometry.html"
   // TGeoTubeSeg: "phi1 is converted to [0,360] (but still expressed in
   // radian, as far as I can tell) and phi2 > phi1."
   // => Convert it to [-pi,+pi].
   double startPhi = tRef.startPhi() * radian;
   startPhi -= 2 * TMath::Pi();
   double endPhi = tRef.endPhi() * radian;
   endPhi -= 2 * TMath::Pi();
   // Get current SUBVOLUME
   DetElement curVol = volman.lookupDetElement(vID);
   const Tube& tCur  = curVol.solid();
   double rMin = tCur.rMin(), rMax = tCur.rMax();
   // Is TRAVERSING?
   double lintos[2][3], louts[2][3], lpini[3], lpend[3], lmend[3];
   std::copy(std::begin(lmom), std::end(lmom), std::begin(lmend));
   unsigned int status =
       cTraversing(lpos, lmom, sim_hit.getPathLength() * ed2dd, sim_hit.isProducedBySecondary(),
                   rMin, rMax, dZ, startPhi, endPhi, lintos, louts, lpini, lpend);
   if (status & m_inconsistency) { // Inconsistency => Drop current "sim_hit"
     error(inconsistency(header, status, sim_hit.getCellID(), lpos, lmom));
     return false;
   }
   cIDs.push_back(sim_hit.getCellID());
   std::vector<int> subHitList;
   if (level() >= algorithms::LogLevel::kDebug) {
     subHitList.push_back(idx);
   }
   // Continuations?
   bool isContinuation  = status & (m_intoLower | m_intoUpper);
   bool hasContinuation = status & (m_outLower | m_outUpper);
   bool canReEnter      = status & m_canReEnter;
   int rank             = m_stripRank(vID);
   if (!canReEnter) {
     if (rank == 0 && (status & m_intoLower))
       isContinuation = false;
     if (rank == 0 && (status & m_outLower))
       hasContinuation = false;
   }
   if (rank == 4 && (status & m_intoUpper))
     isContinuation = false;
   if (rank == 4 && (status & m_outUpper))
     hasContinuation = false;
   if (hasContinuation) {
     // ***** LOOP OVER HITS
     int jdx;
     CellID vIDPrv   = vID;
     size_t sim_size = sim_hits->size();
     for (jdx = idx + 1; jdx < (int)sim_size; jdx++) {
       const edm4hep::SimTrackerHit& sim_hjt = sim_hits->at(jdx);
       CellID vJD                            = sim_hjt.getCellID() & m_volumeBits;
       // Particle may start inward and re-enter, being then outward-going.
       // => Orientation has to be evaluated w.r.t. previous vID.
       int orientation = m_orientation(vIDPrv, vJD);
       bool isUpstream = m_isUpstream(orientation, status);
       bool pmoStatus  = samePMO(sim_hit, sim_hjt);
       if (!pmoStatus || !isUpstream) {
         if ((pmoStatus && !isUpstream) && !sim_hit.isProducedBySecondary()) {
           // Bizarre, except if it's a low energy stuff (when it then can be a
           // looping particle). If it's not let's flag the case, for debugging.
           double P = sqrt(lmom[0] * lmom[0] + lmom[1] * lmom[1] + lmom[2] * lmom[2]);
           if (P > 10 * dd4hep::MeV)
             debug(inconsistency(header, 0, sim_hit.getCellID(), lpos, lmom));
         }
         break;
       }
       // Get 'j' radii
       curVol           = volman.lookupDetElement(vJD);
       const Tube& tubj = curVol.solid();
       rMin             = tubj.rMin();
       rMax             = tubj.rMax();
       double lpoj[3], lmoj[3];
       getLocalPosMom(sim_hjt, toRefVol, lpoj, lmoj);
       // Is TRAVERSING through the (quasi-)common wall?
       double ljns[2][3], lovts[2][3], lpjni[3], lpfnd[3];
       status =
           cTraversing(lpoj, lmoj, sim_hjt.getPathLength() * ed2dd, sim_hit.isProducedBySecondary(),
                       rMin, rMax, dZ, startPhi, endPhi, ljns, lovts, lpjni, lpfnd);
       if (status & m_inconsistency) { // Inconsistency => Drop current "sim_hjt"
         error(inconsistency(header, status, sim_hjt.getCellID(), lpoj, lmoj));
         break;
       }
       // ij-Compatibility: status
       bool jsDownstream = m_isDownstream(orientation, status);
       if (!jsDownstream)
         break;
       // ij-Compatibility: close exit/entrance-distance
       double dist = outInDistance(0, orientation, ljns, louts, lmom, lmoj);
       // RELAXED TOLERANCE for low energy stuff
       double P          = sqrt(lmom[0] * lmom[0] + lmom[1] * lmom[1] + lmom[2] * lmom[2]);
       double tolerance  = m_toleranceFactor(P) * 25 * dd4hep::um;
       bool isCompatible = dist > 0 && dist < tolerance;
       if (!isCompatible) {
         if (!sim_hit.isProducedBySecondary())
           debug(oddity(header, status, dist, sim_hit.getCellID(), lpos, lmom,
                        /* */ sim_hjt.getCellID(), lpoj, lmoj));
         break;
       }
       // ***** UPDATE
       vIDPrv = vJD;
       eDep += sim_hjt.getEDep();
       for (int i = 0; i < 3; i++) { // Update end point position/momentum.
         lpend[i] = lpfnd[i];
         lmend[i] = lmoj[i];
       }
       // ***** BOOK-KEEPING
       cIDs.push_back(sim_hjt.getCellID());
       if (level() >= algorithms::LogLevel::kDebug) {
         subHitList.push_back(jdx);
       }
       // ***** CONTINUATION?
       hasContinuation = status & 0xa;
       canReEnter      = status & 0x100;
       if (!canReEnter && m_stripRank(vJD) == 4)
         hasContinuation = false;
       if (!hasContinuation) {
         jdx++;
         break;
       } else { // Update outgoing position/momentum for next iteration.
         for (int i = 0; i < 3; i++) {
           louts[0][i] = lovts[0][i];
           louts[1][i] = lovts[1][i];
         }
       }
     }
     idx = jdx - 1;
   }
   // ***** EXTENSION?...
   if (sim_hit.isProducedBySecondary() && cIDs.size() < 2)
     if (denyExtension(sim_hit, tCur.rMax() - tCur.rMin())) {
       isContinuation = hasContinuation = false;
     }
   for (int io = 0; io < 2; io++) { // ...into/out-of
     if ((io == 0 && !isContinuation) || (io == 1 && !hasContinuation))
       continue;
     int direction = io ? +1 : -1;
     extendHit(refID, cIDs, direction, lpini, lmom, lpend, lmend);
   }
   // ***** FLAG CASES W/ UNEXPECTED OUTCOME
   flagUnexpected(header, 0, (tRef.rMin() + tRef.rMax()) / 2, sim_hit, lpini, lpend, lpos, lmom);
   // ***** UPDATE (local position <lpos>, DoF)
   double DoF2 = 0, dir = 0;
   for (int i = 0; i < 3; i++) {
     double neu = (lpini[i] + lpend[i]) / 2, alt = lpos[i];
     lpos[i]  = neu;
     double d = neu - alt;
     dir += d * lmom[i];
     DoF2 += d * d;
   }
   // Update time by ToF from original subHit to extended/COALESCED.
   time += ((dir > 0) ? 1 : ((dir < 0) ? -1 : 0)) * sqrt(DoF2) / dd4hep::c_light;
   if (level() >= algorithms::LogLevel::kDebug) {
     debug("--------------------");
     printSubHitList(input, subHitList);
     debug("  =");
     // Print position, eDep and time of coalesced/extended hit
     Position locPos(lpos[0], lpos[1], lpos[2]); // Simplification: strip surface = REFERENCE surface
     Position globPos = refVol.nominal().localToWorld(locPos);
     debug("  position  = ({:7.2f},{:7.2f},{:7.2f}) [mm]", globPos.X() / mm, globPos.Y() / mm,
           globPos.Z() / mm);
     debug("  edep = {:.0f} [eV]", eDep / eV);
     debug("  time = {:.2f} [ns]", time);
   }
   return true;
 }
 bool MPGDTrackerDigi::bCoalesceExtend(const Input& input, int& idx,
                                       std::vector<std::uint64_t>& cIDs, double* lpos, double& eDep,
                                       double& time) const {
   const auto [headers, sim_hits]        = input;
   const edm4hep::EventHeader& header    = headers->at(0);
   const edm4hep::SimTrackerHit& sim_hit = sim_hits->at(idx);
   CellID vID                            = sim_hit.getCellID() & m_volumeBits;
   CellID refID                          = vID & m_moduleBits; // => The REFERENCE SUBVOLUME
   const VolumeManager& volman           = m_detector->volumeManager();
   DetElement refVol                     = volman.lookupDetElement(refID);
   // TGeoHMatrix: In order to avoid a "dangling-reference" warning, let's take
   // a copy of the matrix instead of a reference to it.
   const TGeoHMatrix toRefVol = refVol.nominal().worldTransformation();
   double lmom[3];
   getLocalPosMom(sim_hit, toRefVol, lpos, lmom);
   const double edmm = edm4eic::unit::mm, ed2dd = dd4hep::mm / edmm;
   using dd4hep::mm;
   using edm4eic::unit::eV, edm4eic::unit::GeV;
   // Hit in progress
   eDep = sim_hit.getEDep();
   time = sim_hit.getTime();
   // Get VOLUME parameters
   const Box& bRef = refVol.solid(); // REFERENCE SUBVOLUME
   double dX = bRef.x(), dY = bRef.y();
   // Get current SUBVOLUME
   DetElement curVol = volman.lookupDetElement(vID);
   const Box& bCur   = curVol.solid();
   double dZ         = bCur.z();
   double ref2Cur    = getRef2Cur(refVol, curVol);
   // Is TRAVERSING?
   double lintos[2][3], louts[2][3], lpini[3], lpend[3], lmend[3];
   std::copy(std::begin(lmom), std::end(lmom), std::begin(lmend));
   unsigned int status =
       bTraversing(lpos, lmom, ref2Cur, sim_hit.getPathLength() * ed2dd,
                   sim_hit.isProducedBySecondary(), dZ, dX, dY, lintos, louts, lpini, lpend);
   if (status & m_inconsistency) { // Inconsistency => Drop current "sim_hit"
     error(inconsistency(header, status, sim_hit.getCellID(), lpos, lmom));
     return false;
   }
   cIDs.push_back(sim_hit.getCellID());
   std::vector<int> subHitList;
   if (level() >= algorithms::LogLevel::kDebug) {
     subHitList.push_back(idx);
   }
   // Continuations?
   int rank             = m_stripRank(vID);
   bool isContinuation  = status & (m_intoLower | m_intoUpper);
   bool hasContinuation = status & (m_outLower | m_outUpper);
   if ((rank == 0 && (status & m_intoLower)) || (rank == 4 && (status & m_intoUpper)))
     isContinuation = false;
   if ((rank == 0 && (status & m_outLower)) || (rank == 4 && (status & m_outUpper)))
     hasContinuation = false;
   if (hasContinuation) {
     // ***** LOOP OVER SUBHITS
     int jdx;
     CellID vIDPrv   = vID;
     size_t sim_size = sim_hits->size();
     for (jdx = idx + 1; jdx < (int)sim_size; jdx++) {
       const edm4hep::SimTrackerHit& sim_hjt = sim_hits->at(jdx);
       CellID vJD                            = sim_hjt.getCellID() & m_volumeBits;
       int orientation                       = m_orientation(vIDPrv, vJD);
       bool isUpstream                       = m_isUpstream(orientation, status);
       bool pmoStatus                        = samePMO(sim_hit, sim_hjt);
       if (!pmoStatus || !isUpstream) {
         if ((pmoStatus && !isUpstream) && !sim_hit.isProducedBySecondary()) {
           // Bizarre: let's flag the case for debugging, if not low energy.
           double P = sqrt(lmom[0] * lmom[0] + lmom[1] * lmom[1] + lmom[2] * lmom[2]);
           if (P > 10 * dd4hep::MeV)
             debug(inconsistency(header, 0, sim_hit.getCellID(), lpos, lmom));
         }
         break;
       }
       // Get 'j' Z
       curVol          = volman.lookupDetElement(vJD); // 'j' SUBVOLUME
       const Box& boxj = curVol.solid();
       dZ              = boxj.z();
       double ref2j    = getRef2Cur(refVol, curVol);
       // Is TRAVERSING through the (quasi)-common border?
       double lpoj[3], lmoj[3];
       getLocalPosMom(sim_hjt, toRefVol, lpoj, lmoj);
       double ljns[2][3], lovts[2][3], lpjni[3], lpfnd[3];
       status = bTraversing(lpoj, lmoj, ref2j, sim_hjt.getPathLength() * ed2dd,
                            sim_hit.isProducedBySecondary(), dZ, dX, dY, ljns, lovts, lpjni, lpfnd);
       if (status & m_inconsistency) { // Inconsistency => Drop current "sim_hjt"
         error(inconsistency(header, status, sim_hjt.getCellID(), lpoj, lmoj));
         break;
       }
       // ij-Compatibility: status
       bool jsDownstream = m_isDownstream(orientation, status);
       if (!jsDownstream)
         break;
       // ij-Compatibility: close exit/entrance-distance
       double dist = outInDistance(1, orientation, ljns, louts, lmom, lmoj);
       // RELAXED TOLERANCE for low energy stuff
       double P          = sqrt(lmom[0] * lmom[0] + lmom[1] * lmom[1] + lmom[2] * lmom[2]);
       double tolerance  = m_toleranceFactor(P) * 25 * dd4hep::um;
       bool isCompatible = dist > 0 && dist < tolerance;
       if (!isCompatible) {
         if (!sim_hit.isProducedBySecondary())
           debug(oddity(header, status, dist, sim_hit.getCellID(), lpos, lmom,
                        /* */ sim_hjt.getCellID(), lpoj, lmoj));
         break;
       }
       // ***** UPDATE
       vIDPrv = vJD;
       eDep += sim_hjt.getEDep();
       for (int i = 0; i < 3; i++) { // Update end point position/momentum.
         lpend[i] = lpfnd[i];
         lmend[i] = lmoj[i];
       }
       // ***** BOOK-KEEPING
       cIDs.push_back(sim_hjt.getCellID());
       if (level() >= algorithms::LogLevel::kDebug) {
         subHitList.push_back(jdx);
       }
       // ***** CONTINUATION?
       hasContinuation = status & 0xa;
       if (!hasContinuation) {
         jdx++;
         break;
       } else { // Update outgoing position/momentum for next iteration.
         for (int i = 0; i < 3; i++) {
           louts[0][i] = lovts[0][i];
           louts[1][i] = lovts[1][i];
         }
       }
     }
     idx = jdx - 1;
   }
   // ***** EXTENSION?...
   if (sim_hit.isProducedBySecondary() && cIDs.size() < 2)
     if (denyExtension(sim_hit, bCur.z())) {
       isContinuation = hasContinuation = false;
     }
   for (int io = 0; io < 2; io++) { // ...into/out-of
     if ((io == 0 && !isContinuation) || (io == 1 && !hasContinuation))
       continue;
     int direction = io ? +1 : -1;
     extendHit(refID, cIDs, direction, lpini, lmom, lpend, lmend);
   }
   // ***** FLAG CASES W/ UNEXPECTED OUTCOME
   flagUnexpected(header, 1, 0, sim_hit, lpini, lpend, lpos, lmom);
   // ***** UPDATE (local position <lpos>, DoF)
   double DoF2 = 0, dir = 0;
   for (int i = 0; i < 3; i++) {
     double neu = (lpini[i] + lpend[i]) / 2, alt = lpos[i];
     lpos[i]  = neu;
     double d = neu - alt;
     dir += d * lmom[i];
     DoF2 += d * d;
   }
   // Update time by ToF from original subHit to extended/COALESCED.
   time += ((dir > 0) ? 1 : ((dir < 0) ? -1 : 0)) * sqrt(DoF2) / dd4hep::c_light;
   if (level() >= algorithms::LogLevel::kDebug) {
     debug("--------------------");
     printSubHitList(input, subHitList);
     debug("  =");
     // Print position, eDep and time of coalesced/extended hit
     Position locPos(lpos[0], lpos[1], lpos[2]); // Simplification: strip surface = REFERENCE surface
     Position globPos = refVol.nominal().localToWorld(locPos);
     debug("  position  = ({:7.2f},{:7.2f},{:7.2f}) [mm]", globPos.X() / mm, globPos.Y() / mm,
           globPos.Z() / mm);
     debug("  edep = {:.0f} [eV]", eDep / eV);
     debug("  time = {:.2f} [ns]", time);
   }
   return true;
 }
 void MPGDTrackerDigi::printSubHitList(const Input& input, std::vector<int>& subHitList) const {
   const auto [headers, sim_hits] = input;
   const double edmm              = edm4eic::unit::mm;
   using edm4eic::unit::eV, edm4eic::unit::GeV;
   int ldx = 0;
   for (int kdx : subHitList) {
     const edm4hep::SimTrackerHit& sim_hp = sim_hits->at(kdx);
     CellID cIDk                          = sim_hp.getCellID();
     CellID hIDk = cIDk >> 32, vIDk = cIDk & m_volumeBits;
     if (ldx == 0) {
       debug("Hit cellID{:d} = 0x{:08x}, 0x{:08x}", ldx++, hIDk, vIDk);
     } else {
       debug("  + cellID{:d} = 0x{:08x}, 0x{:08x}", ldx++, hIDk, vIDk);
     }
     debug("  position  = ({:7.2f},{:7.2f},{:7.2f}) [mm]", sim_hp.getPosition().x / edmm,
           sim_hp.getPosition().y / edmm, sim_hp.getPosition().z / edmm);
     debug("  xy_radius = {:.2f}",
           std::hypot(sim_hp.getPosition().x, sim_hp.getPosition().y) / edmm);
     debug("  momentum  = ({:.2f}, {:.2f}, {:.2f}) [GeV]", sim_hp.getMomentum().x / GeV,
           sim_hp.getMomentum().y / GeV, sim_hp.getMomentum().z / GeV);
     debug("  edep = {:.0f} [eV]", sim_hp.getEDep() / eV);
     debug("  time = {:.2f} [ns]", sim_hp.getTime());
   }
 }
 
 void getLocalPosMom(const edm4hep::SimTrackerHit& sim_hit, const TGeoHMatrix& toModule,
                     double* lpos, double* lmom) {
   const edm4hep::Vector3d& pos = sim_hit.getPosition();
   // Length: Inputs are in EDM4eic units. Let's move to DD4hep units.
   const double edmm = edm4eic::unit::mm, ed2dd = dd4hep::mm / edmm;
   const double gpos[3]         = {pos.x * ed2dd, pos.y * ed2dd, pos.z * ed2dd};
   const edm4hep::Vector3f& mom = sim_hit.getMomentum();
   const double gmom[3]         = {mom.x, mom.y, mom.z};
   toModule.MasterToLocal(gpos, lpos);
   toModule.MasterToLocalVect(gmom, lmom);
 }
 
 // ******************** TRAVERSING?
 // Particle can be born/dead (then its position is not (entrance+exit)/2).
 // Or it can exit through the edge.
 // - Returned Status code, see header.
 //     Also, for internal use:
 //     0x10: Enters through edge
 //     0x20: Exits  through edge
 // - <lintos>/<louts>: Positions @ lower/upper wall upon Enter-/Exit-ing (when endorsed by <status>)
 // - <lpini>/<lpend>: Positions of extrema
 // - TOLERANCE? For MIPs, a tolerance of 1 µM works fine. But for lower energy,
 //  looks like we need something somewhat larger. The ideal would be to base
 //  the value on Molière width. Here, I use a somewhat arbitrary built-in.
 //  - If particle found to reach wall, w/in tolerance, assign end points to
 //   walls (instead of <lpos>+/-path/2). It will make so that the eventual
 //   extrapolated position falls exactly at mid-plane, even in the case of a
 //   low energy particle, where path may be affected by multiscattering. This,
 //   provided that particle is not a secondary.
 unsigned int MPGDTrackerDigi::cTraversing(const double* lpos, const double* lmom, double path,
                                           bool isSecondary,         // Input subHit
                                           double rMin, double rMax, // Current instance of SUBVOLUME
                                           double dZ, double startPhi,
                                           double endPhi, // Module parameters
                                           double lintos[][3], double louts[][3], double* lpini,
                                           double* lpend) const {
   unsigned int status = 0;
   double Mx = lpos[0], My = lpos[1], Mz = lpos[2], M2 = Mx * Mx + My * My;
   double Px = lmom[0], Py = lmom[1], Pz = lmom[2];
   // Intersection w/ the edge in phi
   double tIn = 0, tOut = 0;
   for (double phi : {startPhi, endPhi}) {
     // M+t*P = 0 + t'*U. t = (My*Ux-Mx*Uy)/(Px*Uy-Py*Ux);
     double Ux = cos(phi), Uy = sin(phi);
     double D = Px * Uy - Py * Ux;
     if (D) { // If P not // to U
       double t  = (My * Ux - Mx * Uy) / D;
       double Ex = Mx + t * Px, Ey = My + t * Py, Ez = Mz + t * Pz;
       double rE = sqrt(Ex * Ex + Ey * Ey), phiE = atan2(Ey, Ex);
       // The above does not distinguish between phi and phi+pi.
       // => Have to explicitly discard the latter.
       if (rMin < rE && rE < rMax && fabs(Ez) < dZ && fabs(phiE - phi) < 1) {
         if (t < 0) {
           status |= 0x10;
           tIn = t;
         } else {
           status |= 0x20;
           tOut = t;
         }
       }
     }
   }
   // Intersection w/ the edge in Z
   double zLow = -dZ, zUp = +dZ;
   for (double Z : {zLow, zUp}) {
     // Mz+t*Pz = Z
     if (Pz) {
       double t  = (Z - Mz) / Pz;
       double Ex = Mx + t * Px, Ey = My + t * Py, rE = sqrt(Ex * Ex + Ey * Ey);
       double phi = atan2(Ey, Ex);
       if (rMin < rE && rE < rMax && startPhi < phi && phi < endPhi) {
         if (t < 0) {
           if (!(status & 0x10) || ((status & 0x10) && t > tIn)) {
             status |= 0x10;
             tIn = t;
           }
         } else if (t > 0) {
           if (!(status & 0x20) || ((status & 0x20) && t < tOut)) {
             status |= 0x20;
             tOut = t;
           }
         }
       }
     }
   }
   // Intersection w/ tube walls
   double ts[3 /* rMin/rMax/edge */][2 /* In/Out */] = {
       {0, 0}, {0, 0}, {tIn, tOut}}; // Up to two intersections
   double a = Px * Px + Py * Py, b = Px * Mx + Py * My;
   for (int lu = 0; lu < 2; lu++) { // rMin/rMax
     double R;
     unsigned int statGene;
     if (lu == 1) {
       R        = rMax;
       statGene = 0x4;
     } else {
       R        = rMin;
       statGene = 0x1;
     }
     double c = M2 - R * R;
     if (!a) { // P is // to Z. Yet no intersect w/ Z edge.
       if ((status & 0x30) != 0x30)
         status |= 0x1000;
       continue;      // Inconsistency
     } else if (!c) { // Hit is on wall: inconsistency.
       status |= 0x2000;
       continue;
     } else {
       double det = b * b - a * c;
       if (det < 0) {
         if (lu == 1) { // No intersection w/ outer wall: inconsistency.
           status |= 0x4000;
           continue;
         }
       } else {
         double sqdet = sqrt(det);
         for (int is = 0; is < 2; is++) {
           int s     = 1 - 2 * is;
           double t  = (-b + s * sqdet) / a;
           double Ix = Mx + t * Px, Iy = My + t * Py, Iz = Mz + t * Pz, phi = atan2(Iy, Ix);
           if (fabs(Iz) > dZ || phi < startPhi || endPhi < phi)
             continue;
           if (t < 0) {
             // Two rMin intersects in same back/forward direction may happen
             // (one and and only one of them may then be hidden by edge).
             // => Have to allow wall intersect to coexist w/ edge intersect.
             //   This only for rMin, but for simplicity's sake...
             if (status & statGene) { // Two <0 ts: can only happen when rMin
               double tPrv = ts[lu][0];
               if (t > tPrv) {
                 ts[lu][0] = t;
                 ts[lu][1] = tPrv; // Current is actually IN, previous is OUT despite being <0
               }
               status |= statGene << 1;
             } else {
               ts[lu][0] = t;
               status |= statGene;
             }
           } else {                        // (if t > 0)
             if (status & statGene << 1) { // Two >0 ts: can only happen when rMin
               double tPrv = ts[lu][1];
               if (t < tPrv) {
                 ts[lu][1] = t;
                 ts[lu][0] = tPrv; // Current is actually OUT, previous is IN despite being >0
               }
               status |= statGene;
             } else {
               ts[lu][1] = t;
               status |= statGene << 1;
             }
           }
         }
       }
     }
   }
   // Combine w/ edge in/out, based on "t".
   // - A priori, wall crossing (conditioned by w/in edges) and edge crossing (
   //  conditioned by rMin<rE<rMax) are mutually exclusive...
   // - ...This, except when particle re-enters through the lower wall...
   // =>
   //  - No reEntrance: Let's double-check that the above holds, canceling wall
   //   crossing and raising an inconsistency status bit when not.
   //  - Else:
   //    - If doesReEnter, reEntrance is a mere transient trip outside the
   //     SUBVOLUME. => Cancel it.
   //    - Else disregard edge, possibly raising an inconsistency status bit.
   // Hit lies outside?
   // - ...Or when hit position lies outside volume (it can happen, thanks to
   //  limited precision).
   //  => If this is the case, let's swap their exit vs. entrance status. This
   //  will prevent inconsistencies from showing up below.
   // Can reEnter: does reEnter?
   bool canReEnter = (status & 0x3) == 0x3;
   double norm     = sqrt(a + Pz * Pz);
   if (canReEnter) {
     bool doesReEnter = true;
     for (int i12 = 0; i12 < 2; i12++) {
       double t               = ts[0][i12];
       int s                  = t > 0 ? +1 : -1;
       const double tolerance = 20 * dd4hep::um;
       if (path / 2 - s * t * norm < tolerance)
         doesReEnter = false;
     }
     if (doesReEnter) {
       status &= ~0x3;
       canReEnter = false;
     }
   }
   double rHit = sqrt(M2);
   if (rHit < rMin && !canReEnter) {
     unsigned int statvs = status;
     if (statvs & 0x1) {
       status &= ~0x1;
       status |= 0x2;
       ts[0][1] = -ts[0][0];
     }
     if (statvs & 0x2) {
       status &= ~0x2;
       status |= 0x1;
       ts[0][0] = -ts[0][1];
     }
   } else if (rHit > rMax) {
     unsigned int statvs = status;
     if (statvs & 0x4) {
       status &= ~0x4;
       status |= 0x8;
       ts[1][1] = -ts[1][0];
     }
     if (statvs & 0x8) {
       status &= ~0x8;
       status |= 0x4;
       ts[1][0] = -ts[1][1];
     }
   }
   for (int lu = 0; lu < 2; lu++) { // rMin/rMax
     unsigned int statGene = lu ? 0x4 : 0x1;
     if (status & statGene) {
       double t = ts[lu][0];
       if (t < 0) {
         if (status & 0x10) {
           if (lu == 1 || !canReEnter) { // No reEntrance:
             status |= 0x10000;          //   Inconsistency
             status &= ~statGene;        //   Cancel wall crossing
           } else {                      // ReEntrance: disregard edge crossing
             if (t < tIn)
               status |= 0x10000; // Inconsistency
           }
         }
       } else { // if (t > 0)
         if (status & 0x20) {
           if (lu == 1 || !canReEnter) {
             status |= 0x20000;
             status &= ~statGene;
           } else {
             if (t > tOut)
               status |= 0x20000;
           }
         }
       }
     }
     if (status & statGene << 1) {
       double t = ts[lu][1];
       if (t < 0) {
         if (status & 0x10) {
           if (lu == 1 || !canReEnter) {
             status |= 0x40000;
             status &= ~(statGene << 1);
           } else {
             if (t < tIn)
               status |= 0x40000;
           }
         }
       } else { // if (t > 0)
         if (status & 0x20) {
           if (lu == 1 || !canReEnter) {
             status |= 0x80000;
             status &= ~(statGene << 1);
           } else {
             if (t > tOut)
               status |= 0x80000;
           }
         }
       }
     }
   }
   // Is particle born/dead prior to entering/exiting?
   // - sim_hit must have been assigned the mean position: (entrance+exit)/2
   // - Let's then check entrance/exit against sim_hit's position +/- path/2.
   //   When the latter is too short, it means that the particle firing the
   //  hit gets born or dies in the SUBVOLUME ("dying" taken here in the broad
   //  sense of undergoing a discrete physics process).
   // => We remove the corresponding bit in the <status> pattern.
   // - Note that we not only require that the path be long enough, but also
   //  that it matches exactly distances to entrance/exit.
   if (canReEnter)
     status |= 0x100; // Remember that particle can re-enter.
   double at           = path / 2 / norm;
   unsigned int statws = 0;
   for (int is = 0; is < 2; is++) {
     int s     = 1 - 2 * is;
     double Ix = s * at * Px, Iy = s * at * Py, Iz = s * at * Pz;
     for (int lu = 0; lu < 2; lu++) { // Lower/upper wall
       unsigned int statvs = lu ? 0x4 : 0x1;
       for (int io = 0; io < 2; io++) {
         statvs <<= io;
         if (status & statvs) {
           double t = ts[lu][io];
           if (t * s < 0)
             continue;
           double dIx = t * Px - Ix, dIy = t * Py - Iy, dIz = t * Pz - Iz;
           double dist = sqrt(dIx * dIx + dIy * dIy + dIz * dIz);
           // RELAXED TOLERANCE for low energy stuff
           double tolerance = m_toleranceFactor(norm) * 20 * dd4hep::um;
           if (dist < tolerance)
             statws |= statvs;
         }
       }
     }
   }
   if (!(statws & 0x5)) /* No entrance */
     status &= ~0x5;
   if (!(statws & 0xa)) /* No exit */
     status &= ~0xa;
   // ***** End points
   // Assign end points to walls, if not a secondary and provided it's not a
   // reEntrance case, which case is more difficult to handle and we leave aside.
   if (((status & 0x5) == 0x1 || (status & 0x5) == 0x4) && !isSecondary) {
     double tIn = (status & 0x1) ? ts[0][0] : ts[1][0];
     lpini[0]   = Mx + tIn * Px;
     lpini[1]   = My + tIn * Py;
     lpini[2]   = Mz + tIn * Pz;
   } else {
     lpini[0] = Mx - at * Px;
     lpini[1] = My - at * Py;
     lpini[2] = Mz - at * Pz;
   }
   if (((status & 0xa) == 0x2 || (status & 0xa) == 0x8) && !isSecondary) {
     double tOut = (status & 0x2) ? ts[0][1] : ts[1][1];
     lpend[0]    = Mx + tOut * Px;
     lpend[1]    = My + tOut * Py;
     lpend[2]    = Mz + tOut * Pz;
   } else {
     lpend[0] = Mx + at * Px;
     lpend[1] = My + at * Py;
     lpend[2] = Mz + at * Pz;
   }
   // End points when on the walls
   for (int lu = 0; lu < 2; lu++) {
     unsigned int statvs = lu ? 0x4 : 0x1;
     double tIn = ts[lu][0], tOut = ts[lu][1];
     if (status & statvs) {
       lintos[lu][0] = Mx + tIn * Px;
       lintos[lu][1] = My + tIn * Py;
       lintos[lu][2] = Mz + tIn * Pz;
     }
     statvs <<= 1;
     if (status & statvs) {
       louts[lu][0] = Mx + tOut * Px;
       louts[lu][1] = My + tOut * Py;
       louts[lu][2] = Mz + tOut * Pz;
     }
   }
   return status;
 }
 unsigned int MPGDTrackerDigi::bTraversing(const double* lpos, const double* lmom, double ref2Cur,
                                           double path,
                                           bool isSecondary,     // Input subHit
                                           double dZ,            // Current instance of SUBVOLUME
                                           double dX, double dY, // Module parameters
                                           double lintos[][3], double louts[][3], double* lpini,
                                           double* lpend) const {
   unsigned int status = 0;
   double Mx = lpos[0], My = lpos[1], Mxy[2] = {Mx, My};
   double Px = lmom[0], Py = lmom[1], Pxy[2] = {Px, Py};
   double Mz = lpos[2] + ref2Cur, Pz = lmom[2];
   // Intersection w/ the edge in X,Y
   double tIn = 0, tOut = 0;
   double xyLow[2] = {-dX, -dY}, xyUp[2] = {+dX, +dY};
   for (int xy = 0; xy < 2; xy++) {
     int yx       = 1 - xy;
     double a_Low = xyLow[xy], a_Up = xyUp[xy], Ma = Mxy[xy], Pa = Pxy[xy];
     double b_Low = xyLow[yx], b_Up = xyUp[yx], Mb = Mxy[yx], Pb = Pxy[yx];
     for (double A : {a_Low, a_Up}) {
       // Ma+t*Pa = A
       if (Pa) {
         double t  = (A - Ma) / Pa;
         double Eb = Mb + t * Pb, Ez = Mz + t * Pz;
         if (b_Low < Eb && Eb < b_Up && fabs(Ez) < dZ) {
           if (t < 0) {
             if (!(status & 0x10) || ((status & 0x10) && t > tIn)) {
               status |= 0x10;
               tIn = t;
             }
           } else if (t > 0) {
             if (!(status & 0x20) || ((status & 0x20) && t < tOut)) {
               status |= 0x20;
               tOut = t;
             }
           }
         }
       }
     }
   }
   // Intersection w/ box walls
   for (int lu = 0; lu < 2; lu++) {
     int s                 = 2 * lu - 1;
     double Z              = s * dZ;
     unsigned int statGene = lu ? 0x4 : 0x1;
     // Mz+t*Pz = Z
     if (Pz) {
       double t = (Z - Mz) / Pz;
       if (t < 0) {
         if (!(status & 0x10) || ((status & 0x10) && t > tIn)) {
           status |= statGene;
           tIn = t;
         }
       } else if (t > 0) {
         if (!(status & 0x20) || ((status & 0x20) && t < tOut)) {
           status |= statGene << 1;
           tOut = t;
         }
       }
     }
   }
   // Is particle born/dead prior to entering/exiting?
   // - sim_hit must have been assigned the mean position: (entrance+exit)/2
   // - Let's then check entrance/exit against sim_hit's position +/- path/2.
   //   When the latter is too short, it means that the particle firing the
   //  hit gets born or dies in the SUBVOLUME ("dying" taken here in the broad
   //  sense of undergoing a discrete physics process).
   // => We remove the corresponding bit in the <status> pattern.
   // - Note that we not only require that the path be long enough, but also
   //  that it matches exactly distances to entrance/exit.
   double norm = sqrt(Px * Px + Py * Py + Pz * Pz), at = path / 2 / norm;
   unsigned int statws = 0;
   for (int is = 0; is < 2; is++) {
     int s     = 1 - 2 * is;
     double Ix = s * at * Px, Iy = s * at * Py, Iz = s * at * Pz;
     for (int lu = 0; lu < 2; lu++) { // Lower/upper wall
       unsigned int statvs = lu ? 0x4 : 0x1;
       for (int io = 0; io < 2; io++) {
         statvs <<= io;
         if (status & statvs) {
           double t = io ? tOut : tIn;
           if (t * s < 0)
             continue;
           double dIx = t * Px - Ix, dIy = t * Py - Iy, dIz = t * Pz - Iz;
           double dist = sqrt(dIx * dIx + dIy * dIy + dIz * dIz);
           // RELAXED TOLERANCE for low energy stuff
           double tolerance = m_toleranceFactor(norm) * 20 * dd4hep::um;
           if (dist < tolerance)
             statws |= statvs;
         }
       }
     }
   }
   if (!(statws & 0x5)) /* No entrance */
     status &= ~0x5;
   if (!(statws & 0xa)) /* No exit */
     status &= ~0xa;
   // ***** OUTPUT POSITIONS
   Mz -= ref2Cur; // Go back to REFERENCE SUBVOLUME
   // End points:
   // Assign end points to walls, if not a secondary.
   if ((status & 0x5) && !isSecondary) {
     lpini[0] = Mx + tIn * Px;
     lpini[1] = My + tIn * Py;
     lpini[2] = Mz + tIn * Pz;
   } else {
     lpini[0] = Mx - at * Px;
     lpini[1] = My - at * Py;
     lpini[2] = Mz - at * Pz;
   }
   if ((status & 0xa) && !isSecondary) {
     lpend[0] = Mx + tOut * Px;
     lpend[1] = My + tOut * Py;
     lpend[2] = Mz + tOut * Pz;
   } else {
     lpend[0] = Mx + at * Px;
     lpend[1] = My + at * Py;
     lpend[2] = Mz + at * Pz;
   }
   // End points when on the walls:
   for (int lu = 0; lu < 2; lu++) {
     unsigned int statvs = lu ? 0x4 : 0x1;
     if (status & statvs) {
       lintos[lu][0] = Mx + tIn * Px;
       lintos[lu][1] = My + tIn * Py;
       lintos[lu][2] = Mz + tIn * Pz;
     }
     statvs <<= 1;
     if (status & statvs) {
       louts[lu][0] = Mx + tOut * Px;
       louts[lu][1] = My + tOut * Py;
       louts[lu][2] = Mz + tOut * Pz;
     }
   }
   return status;
 }
 
 // ***** EXTRAPOLATE
 bool cExtrapolate(const double* lpos, const double* lmom, // Input subHit
                   double rT,                              // Target radius
                   double* lext)                           // Extrapolated position @ <rT>
 {
   bool ok   = false;
   double Mx = lpos[0], My = lpos[1], Mz = lpos[2], M2 = Mx * Mx + My * My;
   double Px = lmom[0], Py = lmom[1], Pz = lmom[2];
   double a = Px * Px + Py * Py, b = Px * Mx + Py * My, c = M2 - rT * rT;
   double tF = 0;
   if (!c)
     ok = true;
   else if (a) { // P is not // to Z
     double det = b * b - a * c;
     if (det >= 0) {
       double sqdet = sqrt(det);
       for (int is = 0; is < 2; is++) {
         int s    = 1 - 2 * is;
         double t = (-b + s * sqdet) / a, norm = sqrt(a + Pz * Pz);
         // "t" may happen to be slightly <0, because of limited precision
         if (t * norm < -dd4hep::nm)
           continue;
         if (!ok ||
             // Two intersects: let's retain the earliest one.
             (ok && fabs(t) < fabs(tF))) {
           tF = t;
           ok = true;
         }
       }
     }
   }
   if (ok) {
     lext[0] = Mx + tF * Px;
     lext[1] = My + tF * Py;
     lext[2] = Mz + tF * Pz;
   }
   return ok;
 }
 bool bExtrapolate(const double* lpos, const double* lmom, // Input subHit
                   double zT,                              // Target Z
                   double* lext)                           // Extrapolated position @ <zT>
 {
   bool ok   = false;
   double Mx = lpos[0], My = lpos[1], Mz = lpos[2];
   double Px = lmom[0], Py = lmom[1], Pz = lmom[2], norm = sqrt(Px * Px + Py * Py + Pz * Pz);
   double tF = 0;
   if (Pz) {
     tF = (zT - Mz) / Pz;
     // "t" may happen to be slightly <0, because of limited precision
     ok = tF * norm > -dd4hep::nm;
   }
   if (ok) {
     lext[0] = Mx + tF * Px;
     lext[1] = My + tF * Py;
     lext[2] = Mz + tF * Pz;
   }
   return ok;
 }
 
 // ***** EXTENSION
 // At variance to EXTRAPOLATION, we take edges (phi,Z/X,Y) into account.
 // - Returns 0x1 if
 //   - there is an extrapolation between position <lpos> and target,
 //   - within edge limits,
 //   - along momentum <lmom>,
 //   - in direction <direction>.
 // - Else returns 0 or something in the "m_inconsistency" range.
 // - <lext> contains the position of farthest extension.
 unsigned int MPGDTrackerDigi::cExtension(double const* lpos, double const* lmom, // Input subHit
                                          double rT, int direction,               // Target radius
                                          double dZ, double startPhi,
                                          double endPhi, // Module parameters
                                          double* lext) const {
   unsigned int status = 0;
   double Mx = lpos[0], My = lpos[1], Mz = lpos[2];
   double Px = lmom[0], Py = lmom[1], Pz = lmom[2], norm = sqrt(Px * Px + Py * Py + Pz * Pz);
   // Move some distance away from <lpos>, which is expected to be sitting on
   // the wall of the SUBVOLUME to be ``extended''.
   const double margin = 10 * dd4hep::um;
   double t            = direction * margin / norm;
   Mx += t * Px;
   My += t * Py;
   Mz += t * Pz;
   double M2 = Mx * Mx + My * My, rIni = sqrt(M2), rLow, rUp;
   if (rIni < rT) {
     rLow = rIni;
     rUp  = rT;
   } else {
     rLow = rT;
     rUp  = rIni;
   }
   // Intersection w/ the edge in phi
   double tF = 0;
   for (double phi : {startPhi, endPhi}) {
     // M+t*P = 0 + t'*U. t = (My*Ux-Mx*Uy)/(Px*Uy-Py*Ux);
     double Ux = cos(phi), Uy = sin(phi);
     double D = Px * Uy - Py * Ux;
     if (D) { // If P not // to U
       double t = (My * Ux - Mx * Uy) / D;
       if (t * direction < 0)
         continue;
       double Ex = Mx + t * Px, Ey = My + t * Py, Ez = Mz + t * Pz;
       double rE = sqrt(Ex * Ex + Ey * Ey), phiE = atan2(Ey, Ex);
       // Note: have to discard the phi+pi solution.
       if (rLow < rE && rE < rUp && fabs(Ez) < dZ && fabs(phiE - phi) < 1) {
         status |= 0x1;
         tF = t;
       }
     }
   }
   // Intersection w/ the edge in Z
   double zLow = -dZ, zUp = +dZ;
   for (double Z : {zLow, zUp}) {
     // Mz+t*Pz = Z
     if (Pz) {
       double t = (Z - Mz) / Pz;
       if (t * direction < 0)
         continue;
       double Ex = Mx + t * Px, Ey = My + t * Py, rE = sqrt(Ex * Ex + Ey * Ey);
       double phi = atan2(Ey, Ex);
       if (rLow < rE && rE < rUp && startPhi < phi && phi < endPhi) {
         if (t < 0) {
           if (!status || (status && t > tF)) {
             status |= 0x1;
             tF = t;
           }
         } else if (t > 0) {
           if (!status || (status && t < tF)) {
             status |= 0x1;
             tF = t;
           }
         }
       }
     }
   }
   // Else intersection w/ target radius
   if (!status) {
     double a = Px * Px + Py * Py, b = Px * Mx + Py * My, c = M2 - rT * rT;
     if (!a) {           // P is // to Z (while it did no intersect the edge in Z)
       status |= 0x1000; // Inconsistency
     } else if (!c) {    // Hit is on target (while we've moved away from it)
       status |= 0x2000; // Inconsistency
     } else {
       double det = b * b - a * c;
       if (det >= 0) {
         double sqdet = sqrt(det);
         for (int is = 0; is < 2; is++) {
           int s    = 1 - 2 * is;
           double t = (-b + s * sqdet) / a;
           if (t * direction < 0)
             continue;
           double Ix = Mx + t * Px, Iy = My + t * Py, Iz = Mz + t * Pz, phi = atan2(Iy, Ix);
           if (fabs(Iz) > dZ || phi < startPhi || endPhi < phi)
             continue;
           if (!(status & 0x1) ||
               // Two intersects: let's retain the earliest one.
               ((status & 0x1) && fabs(t) < fabs(tF))) {
             tF = t;
             status |= 0x1;
           }
         }
       }
     }
   }
   if (status & 0x1) {
     lext[0] = Mx + tF * Px;
     lext[1] = My + tF * Py;
     lext[2] = Mz + tF * Pz;
   }
   return status;
 }
 unsigned int MPGDTrackerDigi::bExtension(const double* lpos, const double* lmom, // Input subHit
                                          double zT, int direction,               // Target Z
                                          double dX, double dY, // Module parameters
                                          double* lext) const {
   unsigned int status = 0;
   double Mx = lpos[0], My = lpos[1], Mxy[2] = {Mx, My};
   double Px = lmom[0], Py = lmom[1], Pxy[2] = {Px, Py};
   double Mz = lpos[2], Pz = lmom[2];
   double norm = sqrt(Px * Px + Py * Py + Pz * Pz);
   // Move some distance away from <lpos>, which is expected to be sitting on
   // the wall of the SUBVOLUME to be ``extended''.
   const double margin = 10 * dd4hep::um;
   double t            = direction * margin / norm;
   Mx += t * Px;
   My += t * Py;
   Mz += t * Pz;
   double &zIni = Mz, zLow, zUp;
   if (zIni < zT) {
     zLow = zIni;
     zUp  = zT;
   } else {
     zLow = zT;
     zUp  = zIni;
   }
   // Intersection w/ the edge in X,Y
   double tF       = 0;
   double xyLow[2] = {-dX, -dY}, xyUp[2] = {+dX, +dY};
   for (int xy = 0; xy < 2; xy++) {
     int yx       = 1 - xy;
     double a_Low = xyLow[xy], a_Up = xyUp[xy], Ma = Mxy[xy], Pa = Pxy[xy];
     double b_Low = xyLow[yx], b_Up = xyUp[yx], Mb = Mxy[yx], Pb = Pxy[yx];
     for (double A : {a_Low, a_Up}) {
       // Ma+t*Pa = A
       if (Pa) {
         double t = (A - Ma) / Pa;
         if (t * direction < 0)
           continue;
         double Eb = Mb + t * Pb, Ez = Mz + t * Pz;
         if (zLow < Ez && Ez < zUp && b_Low < Eb && Eb < b_Up) {
           if (!status || (status && fabs(t) < fabs(tF))) {
             status |= 0x1;
             tF = t;
           }
         }
       }
     }
   }
   // Else intersection w/ target Z
   if (!status) {
     if (Pz) {
       tF = (zT - Mz) / Pz;
       if (tF * direction > 0)
         status = 0x1;
     }
   }
   if (status) {
     lext[0] = Mx + tF * Px;
     lext[1] = My + tF * Py;
     lext[2] = Mz + tF * Pz;
   }
   return status;
 }
 
 double getRef2Cur(DetElement refVol, DetElement curVol) {
   // TGeoHMatrix: In order to avoid a "dangling-reference" warning,
   // let's take a copy of the matrix instead of a reference to it.
   const TGeoHMatrix toRefVol = refVol.nominal().worldTransformation();
   const TGeoHMatrix toCurVol = curVol.nominal().worldTransformation();
   const double* TRef         = toRefVol.GetTranslation();
   const double* TCur         = toCurVol.GetTranslation();
   // For some reason, it has to be "Ref-Cur", while I (Y.B) would have expected the opposite...
   double gdT[3];
   for (int i = 0; i < 3; i++)
     gdT[i] = TRef[i] - TCur[i];
   double ldT[3];
   toRefVol.MasterToLocalVect(gdT, ldT);
   return ldT[2];
 }
 
 std::string inconsistency(const edm4hep::EventHeader& event, unsigned int status, CellID cID,
                           const double* lpos, const double* lmom) {
   using edm4eic::unit::GeV, dd4hep::mm;
   return fmt::format("Event {}#{}, SimHit 0x{:016x} @ {:.2f},{:.2f},{:.2f} mm, P = "
                      "{:.2f},{:.2f},{:.2f} GeV inconsistency 0x{:x}",
                      event.getRunNumber(), event.getEventNumber(), cID, lpos[0] / mm, lpos[1] / mm,
                      lpos[2] / mm, lmom[0] / GeV, lmom[1] / GeV, lmom[2] / GeV, status);
 }
 std::string oddity(const edm4hep::EventHeader& event, unsigned int status, double dist, CellID cID,
                    const double* lpos, const double* lmom, CellID cJD, const double* lpoj,
                    const double* lmoj) {
   using edm4eic::unit::GeV, dd4hep::mm;
   return fmt::format("Event {}#{}, Bizarre SimHit sequence: 0x{:016x} @ {:.4f},{:.4f},{:.4f} mm, P "
                      "= {:.2f},{:.2f},{:.2f} GeV and 0x{:016x} @ {:.4f},{:.4f},{:.4f} mm, P = "
                      "{:.2f},{:.2f},{:.2f} GeV: status 0x{:x}, distance {:.4f}",
                      event.getRunNumber(), event.getEventNumber(), cID, lpos[0] / mm, lpos[1] / mm,
                      lpos[2] / mm, lmom[0] / GeV, lmom[1] / GeV, lmom[2] / GeV, cJD, lpoj[0] / mm,
                      lpoj[1] / mm, lpoj[2] / mm, lmoj[0] / GeV, lmoj[1] / GeV, lmoj[2] / GeV,
                      status, dist);
 }
 
 bool MPGDTrackerDigi::samePMO(const edm4hep::SimTrackerHit& sim_hit,
                               const edm4hep::SimTrackerHit& sim_hjt) const {
   // Status:
   // 0: Same Particle, same Module, same Origin
   // 0x1: Not same
   // Particle
   bool sameParticle = sim_hjt.getParticle() == sim_hit.getParticle();
   // Module
   CellID vID      = sim_hit.getCellID() & m_volumeBits;
   CellID refID    = vID & m_moduleBits; // => the middle slice
   CellID vJD      = sim_hjt.getCellID() & m_volumeBits;
   CellID refJD    = vJD & m_moduleBits; // => the middle slice
   bool sameModule = refJD == refID;
   // Origin
   // Note: edm4hep::SimTrackerHit possesses an "Overlay" quality. Since I don't
   // know what this is, I ignore it.
   bool isSecondary = sim_hit.isProducedBySecondary();
   bool jsSecondary = sim_hjt.isProducedBySecondary();
   bool sameOrigin  = jsSecondary == isSecondary;
   return sameParticle && sameModule && sameOrigin;
 }
 
 double outInDistance(int shape, int orientation, double lintos[][3], double louts[][3],
                      double* lmom, double* lmoj) {
   // Outgoing/incoming distance
   bool ok;
   double lExt[3];
   double lmOI[3];
   for (int i = 0; i < 3; i++)
     lmOI[i] = (lmom[i] + lmoj[i]) / 2;
   double *lOut, *lInto;
   if (orientation > 0) {
     lOut  = louts[1];
     lInto = lintos[0];
   } else if (orientation < 0) {
     lOut  = louts[0];
     lInto = lintos[1];
   } else {
     lOut  = louts[0];
     lInto = lintos[0];
   }
   if (shape == 0) { // "TGeoTubeSeg"
     double rInto = sqrt(lInto[0] * lInto[0] + lInto[1] * lInto[1]);
     ok           = cExtrapolate(lOut, lmOI, rInto, lExt);
   } else { // "TGeoBBox"
     ok = bExtrapolate(lOut, lmOI, lInto[2], lExt);
   }
   if (ok) {
     double dist2 = 0;
     for (int i = 0; i < 3; i++) {
       double d = lExt[i] - lInto[i];
       dist2 += d * d;
     }
     return sqrt(dist2);
   } else
     return -1;
 }
 
 unsigned int MPGDTrackerDigi::extendHit(CellID refID, std::vector<std::uint64_t>& cIDs,
                                         int direction, double* lpini, double* lmini, double* lpend,
                                         double* lmend) const {
   unsigned int status         = 0;
   const VolumeManager& volman = m_detector->volumeManager();
   DetElement refVol           = volman.lookupDetElement(refID);
   const auto& shape           = refVol.solid();
   double *lpoE, *lmoE; // Starting position/momentum
   if (direction < 0) {
     lpoE = lpini;
     lmoE = lmini;
   } else {
     lpoE = lpend;
     lmoE = lmend;
   }
   // Let's target successively both extreme SUBVOLUMES, disregarding those
   // already contributing to the COALESCED hit.
   // => This proscribes reEntrance extension, i.e. extending past exit point
   //   on a path reEntering into the same SUBVOLUME.
   //    That option could make sense if the reEntering path is at grazing
   //   incidence w.r.t. SUBVOLUME inner wall. But otherwise, it leads to a
   //   very large difference between initial and extended hit.
   for (int rankE : {0, 4}) {
     CellID vIDE      = refID | m_stripIDs[rankE];
     int alreadyThere = 0;
     for (int i = 0; i < (int)cIDs.size(); i++) {
       if ((cIDs[i] & m_volumeBits) == vIDE) {
         alreadyThere = 1;
         break;
       }
     }
     if (alreadyThere)
       continue;
     if (std::find(cIDs.begin(), cIDs.end(), vIDE) != cIDs.end())
       continue;
     DetElement volE = volman.lookupDetElement(vIDE);
     double lext[3];
     if (std::string_view{shape.type()} == "TGeoTubeSeg") {
       const Tube& tExt = volE.solid();
       double R         = rankE == 0 ? tExt.rMin() : tExt.rMax();
       double startPhi  = tExt.startPhi() * radian;
       startPhi -= 2 * TMath::Pi();
       double endPhi = tExt.endPhi() * radian;
       endPhi -= 2 * TMath::Pi();
       double dZ = tExt.dZ();
       status    = cExtension(lpoE, lmoE, R, direction, dZ, startPhi, endPhi, lext);
     } else if (std::string_view{shape.type()} == "TGeoBBox") {
       double ref2E    = getRef2Cur(refVol, volE);
       const Box& bExt = volE.solid();
       double Z        = rankE == 0 ? -bExt.z() : +bExt.z();
       Z -= ref2E;
       double dX = bExt.x(), dY = bExt.y();
       status = bExtension(lpoE, lmoE, Z, direction, dX, dY, lext);
     } else {
       critical(R"(Bad input data: CellID {:x} has invalid shape "{}")", refID, shape.type());
       throw std::runtime_error(R"(Inconsistency: Inappropriate SimHits fed to "MPGDTrackerDigi".)");
     }
     if (status != 0x1)
       continue;
     if (direction < 0) {
       for (int i = 0; i < 3; i++)
         lpini[i] = lext[i];
     } else {
       for (int i = 0; i < 3; i++)
         lpend[i] = lext[i];
     }
     break;
   }
   return status;
 }
 
 bool MPGDTrackerDigi::denyExtension(const edm4hep::SimTrackerHit& sim_hit, double depth) const {
   // Non COALESCED secondary: do not extend...
   //  ...if in HELPER SUBVOLUME: if it is TRAVERSING, it's probably
   //   merely because SUBVOLUME is very thin.
   CellID vID          = sim_hit.getCellID() & m_volumeBits;
   bool isHelperVolume = m_stripRank(vID) != 0 && m_stripRank(vID) != 4;
   //  ...else if path length is negligible compared to potential
   //    extension (here, we cannot avoid using a built-in: 10%).
   const double fraction = .10;
   const double edmm = edm4eic::unit::mm, ed2dd = dd4hep::mm / edmm;
   bool smallPathLength = sim_hit.getPathLength() * ed2dd < fraction * depth;
   return isHelperVolume || smallPathLength;
 }
 
 void MPGDTrackerDigi::flagUnexpected(const edm4hep::EventHeader& event, int shape, double expected,
                                      const edm4hep::SimTrackerHit& sim_hit, double* lpini,
                                      double* lpend, double* lpos, double* lmom) const {
   //  Expectations:
   // I) Primary particle: position = middle of overall sensitive volume.
   // II) Secondary particle: no diff w.r.t. initial.
   //  These expectations are naive ones. When a delta ray is created w/in a
   // sensitive SUBVOLUME, it creates one distinct SimHit and the primary itself
   // no longer spans the whole SUBVOLUME nor sits at the middle. The path of
   // the delta ray is short, typically, but may still turn out to be extendable.
   //  Therefore expectations (I) and (II) are not systematically fulfilled.
   //  The "flagUnexpected" method is mainly there as a placeholder for a
   // debugging tool that would require further development.
   double Rnew2 = 0, Znew, diff2 = 0;
   for (int i = 0; i < 3; i++) {
     double neu = (lpini[i] + lpend[i]) / 2, alt = lpos[i];
     double d = neu - alt;
     diff2 += d * d;
     if (i != 2)
       Rnew2 += neu * neu;
     if (i == 2)
       Znew = neu;
   }
   double found = shape ? Znew : sqrt(Rnew2), residual = found - expected;
   bool isSecondary = sim_hit.isProducedBySecondary();
   bool isPrimary =
       !isSecondary && sqrt(lmom[0] * lmom[0] + lmom[1] * lmom[1] + lmom[2] * lmom[2]) > .1 * GeV;
   if ((fabs(residual) > .000001 && isPrimary) || (sqrt(diff2) > .000001 && isSecondary)) {
     debug("Event {}#{}, SimHit 0x{:016x} origin {:d}: d{:c} = {:.5f} diff = {:.5f}",
           event.getRunNumber(), event.getEventNumber(), sim_hit.getCellID(), isSecondary,
           shape ? 'Z' : 'R', residual, sqrt(diff2));
   }
 }
 
 // ***** CLUSTERIZATION
 // Returns:
 // 0: OK
 // 1: input hit is beyond limits
 int MPGDTrackerDigi::get2HitCluster(CellID refID,
                                     Position& locPos,  // In DD4hep frame
                                     double surfPos[2], // In Surface frame
                                     int pn,            // 'p' or 'n' strip
                                     std::default_random_engine& generator, Cluster& cluster) const {
   //Sim2IDs sim2IDs;
   // Master CellID, from "locPos"
   CellID stripID = m_stripIDs[pn ? 3 : 1]; // 'p' is 2nd in line, 'n' is 4th.
   const Position dummy(0, 0, 0);
   CellID masterID = m_seg->cellID(locPos, dummy, refID | stripID);
   // Retrieve StripParameters (for current sensor, current strip)
   const StripParameters* pars;
   CellID sensorStripID = masterID & m_sensorStripBits;
   try {
     pars = &m_stripParameters.at(sensorStripID);
   } catch (const std::out_of_range& oor) {
     critical(R"(Error retrieving StripParameters for readout "{}", cellID 0x{:0>16x}: {}.)",
              m_cfg.readout, masterID, oor.what());
     throw std::runtime_error("Error retrieving StripParameters");
   }
   // Abscissa = coordinate along measurement axis.
   // hA = simHit   Abscissa (or abscissa of input hit)
   // mA = master   Abscissa (or abscissa of strip fired by input hit)
   // sA = smeared  Abscissa
   // nA = neighbor Abscissa
   const double& sigma = m_cfg.stripResolutions[pn];
   const double min = pars->min, max = pars->max, pitch = pars->pitch;
   double hA = surfPos[pn];
   if (hA < min - (m_truncation - .1 * dd4hep::cm) * sigma ||
       hA > max + (m_truncation - .1 * dd4hep::cm) * sigma) {
     // Exclude hits beyond limits.
     // - In standard situations, min/max are extrema extremorum.
     // - Here we allow for a little more (.1 cm), to cope width whatever
     //  non-standard case.
     // - The limit not to exceed in any case is when the excess would prevent
     //  the _truncated_ Gaussian smearing from moving us back w/in [min,max].
     return 1;
   }
   FieldID stripNum = m_id_dec->get(masterID, pars->index);
   double mA        = m_binToPosition(stripNum, pitch, pars->offset);
   std::normal_distribution<double> gaussian;
   auto stripGauss = [&](double abscissa, double sigma, double min, double max) {
     // Truncated Gaussian: +/-n*sigmas or readout extrema
     double low = abscissa - m_truncation * sigma;
     double up  = abscissa + m_truncation * sigma;
     if (min > low)
       low = min;
     if (max < up)
       up = max;
     double x;
     do {
       x = abscissa + gaussian(generator) * sigma;
     } while (x < low || up < x);
     return x;
   };
   double sA = stripGauss(hA, sigma, min, max);
   // ***** CLUSTER
   // Are we on the edge? Edge being extreme half-cell .
   if (sA < min + pitch / 2 || sA > max - pitch / 2) {
     // If indeed, single-hit cluster
     cluster.push_back({masterID, 1});
   } else {
     // Else two-hit cluster
     CellID neighID, inc = m_stripIncs[pn];
     double nA;
     // (In|de)crement neighbor w.r.t. master:
     // - Simple addition works well in most cases...
     // - ...But not when the overall coordinate field is =0 and the pCoordinate
     //  is decremented: we then get, overall, 0xffffffff, i.e. -1 (if decrement
     //  is 1, that is), nor when vice versa, pCoordinate of 0xffff (=-1) is
     //  incremented.
     // - In practice, the pCoordinate is the only one affected. Yet, to be on
     //  the safe side, let's ensure the spectator coordinate (i.e. "1-pn") is
     //  left unchanged for both (p|n)Coordinates.
     auto incrementID = [&](int pn) {
       CellID spectatorBits = pn == 0 ? ((CellID)0xffff) << 48 : ((CellID)0xffff) << 32;
       CellID iniID         = masterID & spectatorBits;
       neighID              = masterID + inc;
       neighID &= ~spectatorBits;
       neighID |= iniID;
     };
     auto decrementID = [&](int pn) {
       CellID spectatorBits = pn == 0 ? ((CellID)0xffff) << 48 : ((CellID)0xffff) << 32;
       CellID iniID         = masterID & spectatorBits;
       neighID              = masterID - inc;
       neighID &= ~spectatorBits;
       neighID |= iniID;
     };
     if (sA < mA) {
       decrementID(pn);
       nA = mA - pitch;
     } else {
       incrementID(pn);
       nA = mA + pitch;
     }
     // Amplitude Faction (Note: It can't be but >0, see "init").
     double fn = (sA - mA) / (nA - mA);
     cluster.push_back({masterID, 1 - fn});
     cluster.push_back({neighID, fn});
   }
   return 0;
 }
 
 } // namespace eicrecon