// 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;
  auto [raw_hits, links, associations] = output;

  // 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) {
            // create link
            auto link = links->create();
            link.setFrom(item.second);
            link.setTo(sim_hit);
            link.setWeight(1.0);
            // 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
