XRayTransport.h

Go to the documentation of this file.

// Copyright 2019-2026 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
#ifndef ALICEO2_BASE_XRAYTRANSPORT_H_
#define ALICEO2_BASE_XRAYTRANSPORT_H_
 
#include "CADSupport/O2SolidHarness.h"
 
#include "TGeoShape.h"
 
#include <algorithm>
#include <array>
#include <cmath>
#include <string>
#include <vector>
 
namespace o2
{
namespace cad
{
namespace xray
{
 
using o2::cad::harness::Point3D;
 
struct Crossing {
  double t = 0.;
  int kind = 0;
};
 
struct Robustness {
  long long rays = 0;
  long long raysWithCrossings = 0;
  long long crossings = 0;
  long long steps = 0;
  long long zeroLengthSteps = 0;    
  long long nonAdvancingSteps = 0;  
  long long unstickPushes = 0;      
  long long iterationCapHits = 0;   
  long long unterminated = 0;       
  long long oddCrossingLists = 0;   
  long long nonAlternating = 0;     
  long long duplicateCrossings = 0; 
  long long parityMismatchNearBoundary = 0;
  long long parityMismatchIntervals = 0;   
  long long originInside = 0;              
  long long boundaryWithoutTransition = 0; 
  long long originOutsideWorld = 0;
  double insideLength = 0.; 
  double seconds = 0.;
};
 
struct StepConfig {
  double push = 1.e-9;
  double zeroStep = 1.e-9;
  double unstickPush = 1.e-6;
  int maxIter = 512;
  double matchTolerance = 1.e-6;
};
 
// ------------------------------------------------------------------------------------------
// Mode (a): the direct shape-API stepping loop
// ------------------------------------------------------------------------------------------
//
// Contains() to establish the starting state, then alternating DistFromOutside/DistFromInside,
// advancing the point along the ray, until the accumulated distance leaves the raster window.
// `stepmax` is deliberately NOT used to bound the query: its semantics differ between shape
// implementations (some return the crossing, some return stepmax, some return Big), and this loop
// must be a measurement of the crossing list rather than of that convention. The window is
// enforced on the returned crossing distance instead.
 
template <typename ContainsFn, typename DistOutFn, typename DistInFn>
std::vector<Crossing> stepCrossingsWithKernels(const Point3D& origin, const Point3D& dir,
                                               double tMax, const StepConfig& cfg,
                                               Robustness& stats, ContainsFn contains,
                                               DistOutFn distFromOutside, DistInFn distFromInside)
{
  std::vector<Crossing> crossings;
  double point[3] = {origin[0], origin[1], origin[2]};
  bool inside = contains(point);
  if (inside) {
    ++stats.originInside;
  }
  double t = 0.;
  int iter = 0;
  for (; iter < cfg.maxIter; ++iter) {
    const double step = inside ? distFromInside(point, dir.data()) : distFromOutside(point, dir.data());
    ++stats.steps;
    if (!(step < TGeoShape::Big())) {
      break; // no further crossing along this ray
    }
    const double tCross = t + step;
    if (tCross > tMax) {
      break; // beyond the raster window: not this ray's business
    }
    if (step <= cfg.zeroStep) {
      ++stats.zeroLengthSteps;
    }
    crossings.push_back({tCross, inside ? -1 : +1});
    inside = !inside;
    double advance = step + cfg.push;
    if (!(advance > 0.)) {
      ++stats.nonAdvancingSteps;
      advance = cfg.unstickPush;
      ++stats.unstickPushes;
    } else if (step <= cfg.zeroStep) {
      advance = step + cfg.unstickPush;
      ++stats.unstickPushes;
    }
    t += advance;
    if (t > tMax) {
      break;
    }
    for (int k = 0; k < 3; ++k) {
      point[k] = origin[k] + t * dir[k];
    }
  }
  if (iter >= cfg.maxIter) {
    ++stats.iterationCapHits;
  }
  if (inside) {
    ++stats.unterminated;
  }
  return crossings;
}
 
inline std::vector<Crossing> stepWithShapeApi(const TGeoShape* shape, const Point3D& origin,
                                              const Point3D& dir, double tMax,
                                              const StepConfig& cfg, Robustness& stats)
{
  return stepCrossingsWithKernels(
    origin, dir, tMax, cfg, stats, [shape](const double* p) { return shape->Contains(p); },
    [shape](const double* p, const double* d) {
      return shape->DistFromOutside(p, d, 3, TGeoShape::Big(), nullptr);
    },
    [shape](const double* p, const double* d) {
      return shape->DistFromInside(p, d, 3, TGeoShape::Big(), nullptr);
    });
}
 
inline void auditCrossingList(const std::vector<Crossing>& crossings, const TGeoShape* shape,
                              const Point3D& origin, const Point3D& dir, double tMax,
                              const StepConfig& cfg, Robustness& stats)
{
  ++stats.rays;
  stats.crossings += static_cast<long long>(crossings.size());
  if (!crossings.empty()) {
    ++stats.raysWithCrossings;
  }
  if (crossings.size() % 2 != 0) {
    ++stats.oddCrossingLists;
  }
  for (size_t i = 1; i < crossings.size(); ++i) {
    if (crossings[i].kind == crossings[i - 1].kind) {
      ++stats.nonAlternating;
    }
    if (std::fabs(crossings[i].t - crossings[i - 1].t) <= cfg.matchTolerance) {
      ++stats.duplicateCrossings;
    }
  }
  // The inside-segment length: the chord integral's contribution from this ray.
  for (size_t i = 0; i + 1 < crossings.size(); i += 2) {
    if (crossings[i].kind == +1 && crossings[i + 1].kind == -1) {
      stats.insideLength += crossings[i + 1].t - crossings[i].t;
    }
  }
  // The independent check. Both stepping modes produce an alternating list *by construction*, so
  // `nonAlternating` above can never fire on them; asking the shape's own Contains() at the
  // midpoint of every interval is the only way this instrument can contradict itself.
  if (shape != nullptr) {
    std::vector<double> edges;
    edges.push_back(0.);
    for (const auto& c : crossings) {
      edges.push_back(c.t);
    }
    edges.push_back(tMax);
    bool expectInside = false;
    for (size_t i = 0; i + 1 < edges.size(); ++i) {
      const double mid = 0.5 * (edges[i] + edges[i + 1]);
      if (edges[i + 1] - edges[i] > 8. * cfg.matchTolerance) {
        double p[3];
        for (int k = 0; k < 3; ++k) {
          p[k] = origin[k] + mid * dir[k];
        }
        const bool actuallyInside = shape->Contains(p);
        if (actuallyInside != expectInside) {
          // Classify before counting. A midpoint within the match tolerance of the boundary has no
          // defined answer on either side, exactly as the sample gate's `nNoVerdict` points do;
          // counting it as a contradiction would manufacture defects out of near-tangency.
          // Safety() is only paid for on a mismatch, which is rare.
          // Safety() must be asked with the state the shape ITSELF reports; asking it with the
          // state the crossing list expects makes a plain outside point look like a boundary
          // point (TGeoBBox::Safety(p, in=true) goes negative there) and silently excuses every
          // real contradiction. That mistake made this counter read 0 on a deliberately
          // truncated list.
          if (shape->Safety(p, actuallyInside ? kTRUE : kFALSE) <= cfg.matchTolerance) {
            ++stats.parityMismatchNearBoundary;
          } else {
            ++stats.parityMismatchIntervals;
          }
        }
      }
      expectInside = !expectInside;
    }
  }
}
 
struct ListComparison {
  long long rays = 0;
  long long raysIdentical = 0;  
  long long raysStructural = 0; 
  long long matched = 0;
  long long displaced = 0; 
  long long missing = 0;   
  long long extra = 0;     
  long long kindMismatch = 0;
  double worstDeltaT = 0.; 
  Point3D worstOrigin{};
  Point3D worstDir{};
  std::string worstReason;
};
 
inline void compareLists(const std::vector<Crossing>& candidate, const std::vector<Crossing>& reference,
                         const Point3D& origin, const Point3D& dir, double tolerance, ListComparison& out)
{
  ++out.rays;
  bool sameShape = candidate.size() == reference.size();
  for (size_t i = 0; sameShape && i < candidate.size(); ++i) {
    sameShape = candidate[i].kind == reference[i].kind;
  }
  if (sameShape) {
    // Same number of crossings in the same order with the same senses: every difference is a
    // position, so report the positions and never manufacture a missing/extra pair out of one
    // displaced crossing.
    bool identical = true;
    for (size_t i = 0; i < candidate.size(); ++i) {
      const double delta = std::fabs(candidate[i].t - reference[i].t);
      ++out.matched;
      if (delta > tolerance) {
        ++out.displaced;
        identical = false;
      }
      if (delta > out.worstDeltaT) {
        out.worstDeltaT = delta;
        out.worstOrigin = origin;
        out.worstDir = dir;
        out.worstReason = delta > tolerance ? "displaced crossing" : "deltaT";
      }
    }
    out.raysIdentical += identical;
    return;
  }
 
  // Structurally different: walk both lists and attribute each unpaired crossing to the side it
  // came from. This is the branch that names a LOST wall.
  ++out.raysStructural;
  size_t i = 0;
  size_t j = 0;
  while (i < candidate.size() && j < reference.size()) {
    const double delta = candidate[i].t - reference[j].t;
    if (std::fabs(delta) <= tolerance) {
      ++out.matched;
      if (candidate[i].kind != reference[j].kind) {
        ++out.kindMismatch;
      }
      if (std::fabs(delta) > out.worstDeltaT) {
        out.worstDeltaT = std::fabs(delta);
      }
      ++i;
      ++j;
    } else if (delta < 0.) {
      ++out.extra;
      ++i;
    } else {
      ++out.missing;
      ++j;
    }
  }
  out.extra += static_cast<long long>(candidate.size() - i);
  out.missing += static_cast<long long>(reference.size() - j);
  if (out.worstReason.empty() || out.worstReason == "deltaT" ||
      out.worstReason == "displaced crossing") {
    out.worstOrigin = origin;
    out.worstDir = dir;
    out.worstReason = reference.size() > candidate.size() ? "MISSING crossing" : "EXTRA crossing";
  }
}
 
struct RayDef {
  Point3D origin{};
  Point3D dir{};
  double tMax = 0.;
  int beam = 0; 
};
 
struct Beam {
  Point3D dir{};
  Point3D u{};
  Point3D v{};
  std::string label;
};
 
struct Raster {
  int n = 0;
  std::vector<Beam> beams;
  std::vector<double> cellArea; 
  std::vector<double> windowExcess;
  std::vector<RayDef> rays;
  double transverseMargin = 0.;
  Point3D windowMin{}; 
  Point3D windowMax{};
};
 
inline double dot3(const Point3D& a, const Point3D& b)
{
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
 
inline Point3D normalize3(const Point3D& a)
{
  const double norm = std::sqrt(dot3(a, a));
  return {a[0] / norm, a[1] / norm, a[2] / norm};
}
 
inline Point3D cross3(const Point3D& a, const Point3D& b)
{
  return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]};
}
 
inline std::vector<Beam> buildFanBeams(int count)
{
  std::vector<Beam> beams;
  const double golden = 3.14159265358979323846 * (3. - std::sqrt(5.));
  for (int i = 0; i < count; ++i) {
    // Only the upper hemisphere is needed: a beam and its reverse sample the same lines.
    const double z = (count == 1) ? 1. : 1. - static_cast<double>(i) / static_cast<double>(count);
    const double radius = std::sqrt(std::max(0., 1. - z * z));
    const double theta = golden * i;
    Beam beam;
    beam.dir = normalize3({radius * std::cos(theta), radius * std::sin(theta), z});
    // A transverse frame: Gram-Schmidt off whichever axis the beam is least aligned with.
    int least = 0;
    for (int k = 1; k < 3; ++k) {
      if (std::fabs(beam.dir[k]) < std::fabs(beam.dir[least])) {
        least = k;
      }
    }
    Point3D seed{};
    seed[least] = 1.;
    const double projection = dot3(seed, beam.dir);
    beam.u = normalize3({seed[0] - projection * beam.dir[0], seed[1] - projection * beam.dir[1],
                         seed[2] - projection * beam.dir[2]});
    beam.v = cross3(beam.dir, beam.u);
    beam.label = "f" + std::to_string(i);
    beams.push_back(std::move(beam));
  }
  return beams;
}
 
inline std::vector<Beam> buildBeams(const std::string& axesSpec, double tiltDegrees)
{
  std::vector<Beam> beams;
  const double t = std::tan(tiltDegrees * 3.14159265358979323846 / 180.);
  for (const char c : axesSpec) {
    int axis = -1;
    if (c == 'x' || c == 'X') {
      axis = 0;
    } else if (c == 'y' || c == 'Y') {
      axis = 1;
    } else if (c == 'z' || c == 'Z') {
      axis = 2;
    } else {
      continue;
    }
    const int iu = (axis + 1) % 3;
    const int iv = (axis + 2) % 3;
    Point3D w{};
    Point3D u{};
    Point3D v{};
    w[axis] = 1.;
    u[iu] = 1.;
    v[iv] = 1.;
    Beam beam;
    if (t == 0.) {
      beam.dir = w;
      beam.u = u;
      beam.v = v;
      beam.label = std::string(1, "xyz"[axis]);
    } else {
      Point3D dir{w[0] + t * u[0] + 0.618 * t * v[0], w[1] + t * u[1] + 0.618 * t * v[1],
                  w[2] + t * u[2] + 0.618 * t * v[2]};
      beam.dir = normalize3(dir);
      // Gram-Schmidt the transverse frame off the original in-plane axis.
      Point3D uu{u[0] - dot3(u, beam.dir) * beam.dir[0], u[1] - dot3(u, beam.dir) * beam.dir[1],
                 u[2] - dot3(u, beam.dir) * beam.dir[2]};
      beam.u = normalize3(uu);
      beam.v = cross3(beam.dir, beam.u);
      beam.label = std::string(1, "xyz"[axis]) + "+t";
    }
    beams.push_back(std::move(beam));
  }
  return beams;
}
 
inline Raster buildRaster(const Point3D& bboxMin, const Point3D& bboxMax, int n,
                          const std::vector<Beam>& beams, double transverseMargin)
{
  Raster raster;
  raster.n = n;
  raster.beams = beams;
  raster.transverseMargin = transverseMargin;
  for (int k = 0; k < 3; ++k) {
    raster.windowMin[k] = bboxMin[k] - transverseMargin;
    raster.windowMax[k] = bboxMax[k] + transverseMargin;
  }
  for (const auto& beam : beams) {
    // Project the eight bounding-box corners into the beam frame; the window is their extent.
    double lo[3] = {1.e300, 1.e300, 1.e300};
    double hi[3] = {-1.e300, -1.e300, -1.e300};
    for (int corner = 0; corner < 8; ++corner) {
      const Point3D p{(corner & 1) ? bboxMax[0] : bboxMin[0], (corner & 2) ? bboxMax[1] : bboxMin[1],
                      (corner & 4) ? bboxMax[2] : bboxMin[2]};
      const double coordinate[3] = {dot3(p, beam.u), dot3(p, beam.v), dot3(p, beam.dir)};
      for (int k = 0; k < 3; ++k) {
        lo[k] = std::min(lo[k], coordinate[k]);
        hi[k] = std::max(hi[k], coordinate[k]);
      }
    }
    const double uLo = lo[0] - transverseMargin;
    const double vLo = lo[1] - transverseMargin;
    const double du = (hi[0] - lo[0] + 2. * transverseMargin) / n;
    const double dv = (hi[1] - lo[1] + 2. * transverseMargin) / n;
    raster.cellArea.push_back(du * dv);
    const double bboxArea = (hi[0] - lo[0]) * (hi[1] - lo[1]);
    raster.windowExcess.push_back(bboxArea > 0. ? (du * dv * n * n) / bboxArea - 1. : 0.);
    const double extent = hi[2] - lo[2];
    const double lead = 0.05 * extent + 1.e-3;
    const double wStart = lo[2] - lead;
    const int index = static_cast<int>(raster.cellArea.size()) - 1;
    for (int i = 0; i < n; ++i) {
      for (int j = 0; j < n; ++j) {
        const double uu = uLo + (i + 0.5) * du;
        const double vv = vLo + (j + 0.5) * dv;
        RayDef ray;
        ray.beam = index;
        for (int k = 0; k < 3; ++k) {
          ray.origin[k] = uu * beam.u[k] + vv * beam.v[k] + wStart * beam.dir[k];
          ray.dir[k] = beam.dir[k];
        }
        ray.tMax = extent + 2. * lead;
        raster.rays.push_back(ray);
      }
    }
  }
  return raster;
}
 
inline double chordVolume(const Raster& raster, const std::vector<double>& insideLengthPerBeam)
{
  double sum = 0.;
  size_t used = 0;
  for (size_t i = 0; i < raster.beams.size() && i < insideLengthPerBeam.size(); ++i) {
    sum += insideLengthPerBeam[i] * raster.cellArea[i];
    ++used;
  }
  return used > 0 ? sum / static_cast<double>(used) : 0.;
}
 
} // namespace xray
} // namespace cad
} // namespace o2
 
#endif