MatLayerCylSet.cxx

Go to the documentation of this file.

// Copyright 2019-2020 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.
 
#include "DetectorsBase/MatLayerCylSet.h"
#include "CommonConstants/MathConstants.h"
 
#ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version
#include "GPUCommonLogger.h"
#include <TFile.h>
#include <TGeoManager.h>
#include "CommonUtils/TreeStreamRedirector.h"
#include <tbb/blocked_range.h>
#include <tbb/enumerable_thread_specific.h>
#include <tbb/global_control.h>
#include <tbb/parallel_for.h>
#include <algorithm>
#include <chrono>
#include <cstdlib>
#include <vector>
//#define _DBG_LOC_ // for local debugging only
 
#endif // !GPUCA_ALIGPUCODE
#undef NDEBUG
using namespace o2::base;
 
using flatObject = o2::gpu::FlatObject;
 
#ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version
 
//________________________________________________________________________________
void MatLayerCylSet::addLayer(float rmin, float rmax, float zmax, float dz, float drphi)
{
  // add new layer checking for overlaps
  assert(mConstructionMask != Constructed);
  assert(rmin < rmax && zmax > 0 && dz > 0 && drphi > 0);
  mConstructionMask = InProgress;
  int nlr = getNLayers();
  if (!nlr) {
    // book local storage
    auto sz = sizeof(MatLayerCylSetLayout);
    o2::gpu::FlatObject::resizeArray(mFlatBufferContainer, 0, sz);
    mFlatBufferPtr = mFlatBufferContainer;
    mFlatBufferSize = sz;
    //--------------????
    get()->mRMin = 1.e99;
    get()->mRMax = 0.;
  }
 
  for (int il = 0; il < nlr; il++) {
    const auto& lr = getLayer(il);
    if (lr.getRMax() > rmin && rmax > lr.getRMin()) {
      LOG(fatal) << "new layer overlaps with layer " << il;
    }
  }
  auto* oldLayers = o2::gpu::FlatObject::resizeArray(get()->mLayers, nlr, nlr + 1);
  // dynamyc buffers of old layers were used in new ones, detach them
  for (int i = nlr; i--;) {
    oldLayers[i].clearInternalBufferPtr();
  }
  delete[] oldLayers;
  get()->mLayers[nlr].initSegmentation(rmin, rmax, zmax, dz, drphi);
  get()->mNLayers++;
  get()->mRMin = get()->mRMin > rmin ? rmin : get()->mRMin;
  get()->mRMax = get()->mRMax < rmax ? rmax : get()->mRMax;
  get()->mZMax = get()->mZMax < zmax ? zmax : get()->mZMax;
  get()->mRMin2 = get()->mRMin * get()->mRMin;
  get()->mRMax2 = get()->mRMax * get()->mRMax;
}
 
//________________________________________________________________________________
int MatLayerCylSet::getNThreadsFromEnv()
{
  const char* env = std::getenv("NTHREADS_MATBUD");
  if (!env) {
    return 1;
  }
  int n = std::atoi(env);
  if (n < 1) {
    LOG(warning) << "Ignoring invalid NTHREADS_MATBUD=" << env;
    return 1;
  }
  return n;
}
 
//________________________________________________________________________________
void MatLayerCylSet::populateFromTGeo(int ntrPerCell, int nThreads, MatbudGeomBackend backend)
{
  assert(mConstructionMask == InProgress);
 
  if (backend == MatbudGeomBackend::VECGEOM && !GeometryManager::isVecGeomAvailable()) {
    LOG(fatal) << "MatbudGeomBackend::VECGEOM requested but O2 was built without VecGeom support (TGeo2VecGeom not found at configure time)";
  }
 
  int nlr = getNLayers();
  if (!nlr) {
    LOG(error) << "The LUT is not yet initialized";
    return;
  }
  if (get()->mR2Intervals) {
    LOG(error) << "The LUT is already populated";
    return;
  }
 
  if (nThreads < 0) {
    nThreads = getNThreadsFromEnv();
  }
 
  using Clock = std::chrono::steady_clock;
  auto seconds = [](Clock::time_point a, Clock::time_point b) {
    return std::chrono::duration<double>(b - a).count();
  };
 
  if (nThreads <= 1) {
    for (int i = 0; i < nlr; i++) {
      LOG(info) << "Populating with " << ntrPerCell << " trials Lr " << i;
      get()->mLayers[i].print();
    }
    const auto tSetupStart = Clock::now();
#ifdef O2_WITH_VECGEOM
    if (backend == MatbudGeomBackend::VECGEOM) {
      // Trigger the lazy VecGeom world build/BVH init here so it counts as "setup" below,
      // not as fill time for whichever cell happens first.
      GeometryManager::vecGeomMaterialBudget(0.f, 0.f, 0.f, 0.f, 0.f, 1.f);
    }
#endif
    const auto tFillStart = Clock::now();
    for (int i = 0; i < nlr; i++) {
      get()->mLayers[i].populateFromTGeo(ntrPerCell, backend);
    }
    const auto tFillEnd = Clock::now();
    finalizeStructures();
    LOG(info) << "LUT fill: 1 thread, setup " << seconds(tSetupStart, tFillStart)
              << " s, cells " << seconds(tFillStart, tFillEnd) << " s";
    return;
  }
 
  // Cells of all layers form one flat index range so that the load is balanced across
  // threads even though layers differ a lot in cell count. layerOffsets[i] is the first
  // flat index of layer i; a binary search maps a flat index back to (layer, iz, iphi).
  std::vector<size_t> layerOffsets(nlr + 1, 0);
  for (int i = 0; i < nlr; i++) {
    LOG(info) << "Queuing " << ntrPerCell << " trials Lr " << i;
    get()->mLayers[i].print();
    const auto& lr = get()->mLayers[i];
    layerOffsets[i + 1] = layerOffsets[i] + size_t(lr.getNZBins()) * lr.getNPhiBins();
  }
  const size_t totalCells = layerOffsets[nlr];
 
  const auto tSetupStart = Clock::now();
 
  auto fillRange = [this, ntrPerCell, backend, &layerOffsets](const tbb::blocked_range<size_t>& range, TGeoNavigator* nav) {
    for (size_t idx = range.begin(); idx != range.end(); ++idx) {
      auto it = std::upper_bound(layerOffsets.begin(), layerOffsets.end(), idx);
      const int layerIdx = int(std::distance(layerOffsets.begin(), it)) - 1;
      const size_t cellInLayer = idx - layerOffsets[layerIdx];
      auto& layer = this->get()->mLayers[layerIdx];
      const int nphi = layer.getNPhiBins();
      layer.populateFromTGeo(int(cellInLayer % nphi), int(cellInLayer / nphi), ntrPerCell, nav, backend);
    }
  };
 
  Clock::time_point tFillStart, tFillEnd;
  if (backend == MatbudGeomBackend::ROOT) {
    // TGeo has to be told that several threads will navigate it, and each thread needs its own
    // navigator. SetMaxThreads() is one-way -- TGeoManager has no API to return to
    // single-threaded mode -- so we do not pretend to restore it; that is harmless because
    // meanMaterialBudget() decides whether to lock from its own argument, not from this global.
    // The navigators we book are ours, though, so those we do give back.
    gGeoManager->SetMaxThreads(nThreads);
 
    tbb::enumerable_thread_specific<TGeoNavigator*> threadNavigators(
      []() { return gGeoManager->AddNavigator(); });
 
    tFillStart = Clock::now();
    {
      tbb::global_control threadControl(tbb::global_control::max_allowed_parallelism, nThreads);
      tbb::parallel_for(tbb::blocked_range<size_t>(0, totalCells),
                        [&fillRange, &threadNavigators](const tbb::blocked_range<size_t>& range) {
                          fillRange(range, threadNavigators.local());
                        });
    }
    tFillEnd = Clock::now();
 
    for (TGeoNavigator* nav : threadNavigators) {
      gGeoManager->RemoveNavigator(nav);
    }
  } else {
    // VecGeom navigation needs no per-thread navigator bookkeeping. Trigger the lazy world
    // build/BVH init before tFillStart so it counts as "setup", not fill time.
#ifdef O2_WITH_VECGEOM
    GeometryManager::vecGeomMaterialBudget(0.f, 0.f, 0.f, 0.f, 0.f, 1.f);
#endif
    tFillStart = Clock::now();
    {
      tbb::global_control threadControl(tbb::global_control::max_allowed_parallelism, nThreads);
      tbb::parallel_for(tbb::blocked_range<size_t>(0, totalCells),
                        [&fillRange](const tbb::blocked_range<size_t>& range) {
                          fillRange(range, nullptr);
                        });
    }
    tFillEnd = Clock::now();
  }
 
  finalizeStructures();
  const auto tEnd = Clock::now();
 
  // Reported separately because only the middle term scales: the setup walks every volume
  // in the geometry (TGeoManager::SetMaxThreads) and the teardown is serial by nature.
  LOG(info) << "LUT fill: " << nThreads << " threads, setup " << seconds(tSetupStart, tFillStart)
            << " s, cells " << seconds(tFillStart, tFillEnd)
            << " s, finalize " << seconds(tFillEnd, tEnd) << " s";
}
 
//________________________________________________________________________________
void MatLayerCylSet::finalizeStructures()
{
  // build layer search structures
  assert(mConstructionMask == InProgress);
  int nlr = getNLayers();
  int nR2Int = 2 * (nlr + 1);
  o2::gpu::FlatObject::resizeArray(get()->mR2Intervals, 0, nR2Int);
  o2::gpu::FlatObject::resizeArray(get()->mInterval2LrID, 0, nR2Int);
  get()->mR2Intervals[0] = get()->mRMin2;
  get()->mR2Intervals[1] = getLayer(0).getRMax2();
  get()->mInterval2LrID[0] = 0;
  auto& nRIntervals = get()->mNRIntervals;
  nRIntervals = 1;
 
  for (int i = 1; i < nlr; i++) {
    const auto& lr = getLayer(i);
    if (o2::math_utils::sqrt(lr.getRMin2()) > o2::math_utils::sqrt(get()->mR2Intervals[nRIntervals] + Ray::Tiny)) {
      // register gap
      get()->mInterval2LrID[nRIntervals] = -1;
      get()->mR2Intervals[++nRIntervals] = lr.getRMin2();
    }
    get()->mInterval2LrID[nRIntervals] = i;
    get()->mR2Intervals[++nRIntervals] = lr.getRMax2();
  }
  delete[] o2::gpu::FlatObject::resizeArray(get()->mInterval2LrID, nR2Int, nRIntervals); // rebook with precise size
  delete[] o2::gpu::FlatObject::resizeArray(get()->mR2Intervals, nR2Int, ++nRIntervals); // rebook with precise size
  //
}
 
//________________________________________________________________________________
void MatLayerCylSet::dumpToTree(const std::string& outName) const
{
 
  o2::utils::TreeStreamRedirector dump(outName.data(), "recreate");
  for (int i = 0; i < getNLayers(); i++) {
    const auto& lr = getLayer(i);
    float r = 0.5 * (lr.getRMin() + lr.getRMax());
    // per cell dump
    int nphib = lr.getNPhiBins();
    for (int ip = 0; ip < nphib; ip++) {
      float phi = 0.5 * (lr.getPhiBinMin(ip) + lr.getPhiBinMax(ip));
      float sn, cs;
      int ips = lr.phiBin2Slice(ip);
      char merge = 0; // not mergeable
      if (ip + 1 < nphib) {
        int ips1 = lr.phiBin2Slice(ip + 1);
        merge = ips == ips1 ? -1 : lr.canMergePhiSlices(ips, ips1); // -1 for already merged
      } else {
        merge = -2; // last one
      }
      o2::math_utils::sincos(phi, sn, cs);
      float x = r * cs, y = r * sn;
      for (int iz = 0; iz < lr.getNZBins(); iz++) {
        float z = 0.5 * (lr.getZBinMin(iz) + lr.getZBinMax(iz));
        auto cell = lr.getCellPhiBin(ip, iz);
        dump << "cell"
             << "ilr=" << i << "r=" << r << "phi=" << phi << "x=" << x << "y=" << y << "z=" << z << "ip=" << ip << "ips=" << ips << "iz=" << iz
             << "mrgnxt=" << merge << "val=" << cell << "\n";
      }
    }
    //
    // statistics per layer
    MatCell mean, rms;
    lr.getMeanRMS(mean, rms);
    dump << "lay"
         << "ilr=" << i << "r=" << r << "mean=" << mean << "rms=" << rms << "\n";
  }
}
 
//________________________________________________________________________________
void MatLayerCylSet::writeToFile(const std::string& outFName)
{
 
  TFile outf(outFName.data(), "recreate");
  if (outf.IsZombie()) {
    return;
  }
  outf.WriteObjectAny(this, Class(), "ccdb_object");
  outf.Close();
}
 
void MatLayerCylSet::initLayerVoxelLU()
{
  if (mInitializedLayerVoxelLU) {
    LOG(info) << "Layer voxel already initialized; Aborting";
    return;
  }
  LOG(info) << "Initializing voxel layer lookup";
  // do some check if voxels are dimensioned correctly
  if (LayerRMax < get()->mRMax) {
    LOG(fatal) << "Cannot initialized layer voxel lookup due to dimension problem (fix constants in MatLayerCylSet.h)";
  }
  // the top bit of an entry carries the ambiguity flag, so the interval index has one bit less
  if (get()->mNRIntervals > VoxelSegmentMask) {
    LOG(fatal) << "Too many R intervals (" << get()->mNRIntervals << ") to pack into a layer voxel lookup entry";
  }
  for (int voxel = 0; voxel < NumVoxels; ++voxel) {
    // check the 2 extremes of this voxel "covering"
    const auto lowerR = voxelRMin(voxel);
    const auto upperR = voxelRMax(voxel);
    const auto lowerSegment = searchSegment(lowerR * lowerR);
    const auto upperSegment = searchSegment(upperR * upperR);
    mLayerVoxelLU[voxel] = uint16_t(lowerSegment) | (lowerSegment != upperSegment ? VoxelAmbiguousBit : uint16_t{0});
  }
  mInitializedLayerVoxelLU = true;
}
 
//________________________________________________________________________________
MatLayerCylSet* MatLayerCylSet::loadFromFile(const std::string& inpFName)
{
  TFile inpf(inpFName.data());
  if (inpf.IsZombie()) {
    LOG(error) << "Failed to open input file " << inpFName;
    return nullptr;
  }
  MatLayerCylSet* mb = reinterpret_cast<MatLayerCylSet*>(inpf.GetObjectChecked("ccdb_object", Class()));
  if (!mb && !(mb = reinterpret_cast<MatLayerCylSet*>(inpf.GetObjectChecked("MatBud", Class())))) { // for old objects
    LOG(error) << "Failed to load mat.LUT from " << inpFName;
    return nullptr;
  }
  auto rptr = rectifyPtrFromFile(mb);
  return rptr;
}
 
//________________________________________________________________________________
MatLayerCylSet* MatLayerCylSet::rectifyPtrFromFile(MatLayerCylSet* ptr)
{
  // rectify object loaded from file
  if (ptr && !ptr->get()) {
    ptr->fixPointers();
  }
  ptr->initLayerVoxelLU();
  return ptr;
}
 
//________________________________________________________________________________
void MatLayerCylSet::optimizePhiSlices(float maxRelDiff)
{
  // merge similar (whose relative budget does not differ within maxRelDiff) phi slices
  assert(mConstructionMask == InProgress);
  for (int i = getNLayers(); i--;) {
    get()->mLayers[i].optimizePhiSlices(maxRelDiff);
  }
  // flatten();  // RS: TODO
}
 
//________________________________________________________________________________
void MatLayerCylSet::print(bool data) const
{
  if (!get()) {
    printf("Not initialized yet\n");
    return;
  }
  if (mConstructionMask != Constructed) {
    LOG(warning) << "Object is not yet flattened";
  }
  for (int i = 0; i < getNLayers(); i++) {
    printf("#%3d | ", i);
    getLayer(i).print(data);
  }
  printf("%.2f < R < %.2f  %d layers with total size %.2f MB\n", getRMin(), getRMax(), getNLayers(),
         float(getFlatBufferSize()) / 1024 / 1024);
}
 
//________________________________________________________________________________
void MatLayerCylSet::scaleLayersByID(int lrFrom, int lrTo, float factor, bool _x2x0, bool _rho)
{
  lrFrom = std::max(0, std::min(lrFrom, get()->mNLayers - 1));
  lrTo = std::max(0, std::min(lrTo, get()->mNLayers - 1));
  int dir = lrFrom >= lrTo ? -1 : 1;
  lrTo += dir;
  for (int i = lrFrom; i != lrTo; i += dir) {
    get()->mLayers[i].scale(factor, _x2x0, _rho);
  }
}
 
//________________________________________________________________________________
void MatLayerCylSet::scaleLayersByR(float rFrom, float rTo, float factor, bool _x2x0, bool _rho)
{
  if (rFrom > rTo) {
    std::swap(rFrom, rTo);
  }
  Ray ray(std::max(getRMin(), rFrom), 0., 0., std::min(getRMax(), rTo), 0., 0.);
  short lmin, lmax;
  if (!getLayersRange(ray, lmin, lmax)) {
    LOGP(warn, "No layers found for {} < r < {}", rFrom, rTo);
    return;
  }
  scaleLayersByID(lmin, lmax, factor, _x2x0, _rho);
}
 
#endif 
 
#ifndef GPUCA_GPUCODE
//________________________________________________________________________________
std::size_t MatLayerCylSet::estimateFlatBufferSize() const
{
  std::size_t sz = alignSize(sizeof(MatLayerCylSetLayout), getBufferAlignmentBytes()); // hold data members
 
  sz = alignSize(sz + get()->mNLayers * sizeof(MatLayerCyl), MatLayerCyl::getClassAlignmentBytes());
  sz = alignSize(sz + (get()->mNRIntervals + 1) * sizeof(float), getBufferAlignmentBytes());
  sz = alignSize(sz + get()->mNRIntervals * sizeof(int), getBufferAlignmentBytes());
 
  for (int i = 0; i < getNLayers(); i++) {
    sz = alignSize(sz + getLayer(i).estimateFlatBufferSize(), getBufferAlignmentBytes());
  }
  return sz;
}
#endif // ! GPUCA_GPUCODE
 
//_________________________________________________________________________________________________
GPUd() MatBudget MatLayerCylSet::getMatBudget(float x0, float y0, float z0, float x1, float y1, float z1) const
{
  // get material budget traversed on the line between point0 and point1
  MatBudget rval;
  Ray ray(x0, y0, z0, x1, y1, z1);
  short lmin, lmax; // get innermost and outermost relevant layer
  if (ray.isTooShort() || !getLayersRange(ray, lmin, lmax)) {
    rval.length = ray.getDist();
    return rval;
  }
  short lrID = lmax;
  while (lrID >= lmin) { // go from outside to inside
    const auto& lr = getLayer(lrID);
    int nphiSlices = lr.getNPhiSlices();
    int nc = ray.crossLayer(lr); // determines how many crossings this ray has with this tubular layer
    for (int ic = nc; ic--;) {
      float cross1, cross2;
      ray.getCrossParams(ic, cross1, cross2); // tmax,tmin of crossing the layer
 
      auto phi0 = ray.getPhi(cross1), phi1 = ray.getPhi(cross2), dPhi = phi0 - phi1;
      auto phiID = lr.getPhiSliceID(phi0), phiIDLast = lr.getPhiSliceID(phi1);
      // account for eventual wrapping around 0
      if (dPhi > 0.f) {
        if (dPhi > o2::constants::math::PI) { // wraps around phi=0
          phiIDLast += nphiSlices;
        }
      } else {
        if (dPhi < -o2::constants::math::PI) { // wraps around phi=0
          phiID += nphiSlices;
        }
      }
 
      int stepPhiID = phiID > phiIDLast ? -1 : 1;
      bool checkMorePhi = true;
      auto tStartPhi = cross1, tEndPhi = 0.f;
      do {
        // get the path in the current phi slice
        if (phiID == phiIDLast) {
          tEndPhi = cross2;
          checkMorePhi = false;
        } else { // last phi slice still not reached
          const int boundaryPhiID = stepPhiID > 0 ? phiID + 1 : phiID;
          // phiID may be offset by one revolution to handle wrapping, but never by more.
          const int wrappedBoundaryPhiID = boundaryPhiID < nphiSlices ? boundaryPhiID : boundaryPhiID - nphiSlices;
          tEndPhi = ray.crossRadial(lr, wrappedBoundaryPhiID);
          if (tEndPhi == Ray::InvalidT) {
            break; // ray parallel to radial line, abandon check for phi bin change
          }
          const auto tMarginPhi = 1.e-6f + 1.e-5f * (cross1 - cross2);
          // if (!(tEndPhi >= cross2 - tMarginPhi) | !(tEndPhi <= cross1 + tMarginPhi)) { // use non-short-circuit | to reject eventual NANs
          if (tEndPhi < cross2 - tMarginPhi || tEndPhi > cross1 + tMarginPhi) {
            tEndPhi = cross2;
            checkMorePhi = false;
          }
        }
        auto zID = lr.getZBinID(ray.getZ(tStartPhi));
        auto zIDLast = lr.getZBinID(ray.getZ(tEndPhi));
        const int wrappedPhiID = phiID < nphiSlices ? phiID : phiID - nphiSlices;
        const auto* cellRow = lr.getCellRow(wrappedPhiID);
        // check if Zbins are crossed
 
#ifdef _DBG_LOC_
        printf("-- Zdiff (%3d : %3d) mode: t: %+e %+e\n", zID, zIDLast, tStartPhi, tEndPhi);
#endif
 
        if (zID != zIDLast) {
          auto stepZID = zID < zIDLast ? 1 : -1;
          bool checkMoreZ = true;
          auto tStartZ = tStartPhi, tEndZ = 0.f;
          do {
            if (zID == zIDLast) {
              tEndZ = tEndPhi;
              checkMoreZ = false;
            } else {
              tEndZ = ray.crossZ(lr.getZBinMin(stepZID > 0 ? zID + 1 : zID));
              if (tEndZ == Ray::InvalidT) { // track normal to Z axis, abandon Zbin change test
                break;
              }
            }
            // account materials of this step
            float step = tEndZ > tStartZ ? tEndZ - tStartZ : tStartZ - tEndZ; // the real step is ray.getDist(tEnd-tStart), will rescale all later
            const auto& cell = cellRow[zID];
            rval.meanRho += cell.meanRho * step;
            rval.meanX2X0 += cell.meanX2X0 * step;
            rval.length += step;
 
#ifdef _DBG_LOC_
            float pos0[3] = {ray.getPos(tStartZ, 0), ray.getPos(tStartZ, 1), ray.getPos(tStartZ, 2)};
            float pos1[3] = {ray.getPos(tEndZ, 0), ray.getPos(tEndZ, 1), ray.getPos(tEndZ, 2)};
            printf(
              "Lr#%3d / cross#%d : account %f<t<%f at phiSlice %d | Zbin: %3d (%3d) |[%+e %+e +%e]:[%+e %+e %+e] "
              "Step: %.3e StrpCor: %.3e\n",
              lrID, ic, tEndZ, tStartZ, wrappedPhiID, zID, zIDLast,
              pos0[0], pos0[1], pos0[2], pos1[0], pos1[1], pos1[2], step, ray.getDist(step));
#endif
 
            tStartZ = tEndZ;
            zID += stepZID;
          } while (checkMoreZ);
        } else {
          float step = tEndPhi > tStartPhi ? tEndPhi - tStartPhi : tStartPhi - tEndPhi; // the real step is |ray.getDist(tEnd-tStart)|, will rescale all later
          const auto& cell = cellRow[zID];
          rval.meanRho += cell.meanRho * step;
          rval.meanX2X0 += cell.meanX2X0 * step;
          rval.length += step;
 
#ifdef _DBG_LOC_
          float pos0[3] = {ray.getPos(tStartPhi, 0), ray.getPos(tStartPhi, 1), ray.getPos(tStartPhi, 2)};
          float pos1[3] = {ray.getPos(tEndPhi, 0), ray.getPos(tEndPhi, 1), ray.getPos(tEndPhi, 2)};
          printf(
            "Lr#%3d / cross#%d : account %f<t<%f at phiSlice %d | Zbin: %3d ----- |[%+e %+e +%e]:[%+e %+e %+e]"
            "Step: %.3e StrpCor: %.3e\n",
            lrID, ic, tEndPhi, tStartPhi, wrappedPhiID, zID,
            pos0[0], pos0[1], pos0[2], pos1[0], pos1[1], pos1[2], step, ray.getDist(step));
#endif
        }
        //
        tStartPhi = tEndPhi;
        phiID += stepPhiID;
 
      } while (checkMorePhi);
    }
    lrID--;
  } // loop over layers
 
  if (rval.length != 0.f) {
    rval.meanRho /= rval.length;                                       // average
    rval.meanX2X0 *= ray.getDist();                                    // normalize
  }
  rval.length = ray.getDist();
 
#ifdef _DBG_LOC_
  printf("<rho> = %e, x2X0 = %e  | step = %e\n", rval.meanRho, rval.meanX2X0, rval.length);
#endif
  return rval;
}
 
//_________________________________________________________________________________________________
GPUd() bool MatLayerCylSet::getLayersRange(const Ray& ray, short& lmin, short& lmax) const
{
  // get range of layers corresponding to rmin/rmax
  //
  lmin = lmax = -1;
  float rmin2, rmax2;
  ray.getMinMaxR2(rmin2, rmax2);
 
  if (rmin2 >= getRMax2() || rmax2 <= getRMin2()) {
    return false;
  }
  int lmxInt, lmnInt;
  if (!mInitializedLayerVoxelLU) {
    lmxInt = rmax2 < getRMax2() ? searchSegment(rmax2, 0) : get()->mNRIntervals - 2;
    lmnInt = rmin2 >= getRMin2() ? searchSegment(rmin2, 0, lmxInt + 1) : 0;
  } else {
    // The two lookups are independent so overlapping the pair is worth the clumsier shape.
    const bool useMax = rmax2 < getRMax2();
    const bool useMin = rmin2 >= getRMin2();
    const int ixMax = useMax ? voxelIndex(rmax2) : NumVoxels - 1;
    const int ixMin = useMin ? voxelIndex(rmin2) : 0;
    const uint16_t eMax = mLayerVoxelLU[ixMax];
    const uint16_t eMin = mLayerVoxelLU[ixMin];
    lmxInt = useMax ? resolveLayerRange(rmax2, ixMax, eMax) : get()->mNRIntervals - 2;
    lmnInt = useMin ? resolveLayerRange(rmin2, ixMin, eMin) : 0;
  }
 
  const auto* interval2LrID = get()->mInterval2LrID;
  lmax = interval2LrID[lmxInt];
  lmin = interval2LrID[lmnInt];
  // make sure lmnInt and/or lmxInt are not in the gap
  if (lmax < 0) {
    lmax = interval2LrID[lmxInt - 1]; // rmax2 is in the gap, take highest layer below rmax2
  }
  if (lmin < 0) {
    lmin = interval2LrID[lmnInt + 1]; // rmin2 is in the gap, take lowest layer above rmin2
  }
  return lmin <= lmax; // valid if both are not in the same gap
}
 
GPUd() int MatLayerCylSet::searchLayerFast(float r2, int low, int high) const
{
  // we can avoid the sqrt .. at the cost of more memory in the lookup
  const auto index = voxelIndex(r2);
  return resolveLayerRange(r2, index, mLayerVoxelLU[index]);
}
 
GPUd() int MatLayerCylSet::resolveLayerRange(float r2, int voxel, uint16_t entry) const
{
  const int layersfirst = entry & VoxelSegmentMask;
  if (entry & VoxelAmbiguousBit) {
    // Recreate the upper candidate only for the small fraction of undecided voxels
    const auto upperR = voxelRMax(voxel);
    const auto layerslast = searchSegment(upperR * upperR);
    return searchSegment(r2, layersfirst, layerslast + 1);
  }
  return layersfirst;
}
 
GPUd() int MatLayerCylSet::searchSegment(float val, int low, int high) const
{
  if (low < 0) {
    low = 0;
  }
  if (high < 0) {
    high = get()->mNRIntervals;
  }
  int mid = (low + high) >> 1;
  const auto* r2Intervals = get()->mR2Intervals;
  while (mid != low) {
    if (val < r2Intervals[mid]) {
      high = mid;
    } else {
      low = mid;
    }
    mid = (low + high) >> 1;
  }
 
  return mid;
}
 
#ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version
 
void MatLayerCylSet::flatten()
{
  // make object flat: move all content to single internally allocated buffer
  assert(mConstructionMask == InProgress);
 
  int sz = estimateFlatBufferSize();
  // create new internal buffer with total size and copy data
  delete[] o2::gpu::FlatObject::resizeArray(mFlatBufferContainer, mFlatBufferSize, sz);
  mFlatBufferPtr = mFlatBufferContainer;
  mFlatBufferSize = sz;
  int nLr = getNLayers();
 
  auto offs = alignSize(sizeof(MatLayerCylSetLayout), getBufferAlignmentBytes()); // account for the alignment
  // move array of layer pointers to the flat array
  auto* oldLayers = o2::gpu::FlatObject::resizeArray(get()->mLayers, nLr, nLr, (MatLayerCyl*)(mFlatBufferPtr + offs));
  // dynamyc buffers of old layers were used in new ones, detach them
  for (int i = nLr; i--;) {
    oldLayers[i].clearInternalBufferPtr();
  }
  delete[] oldLayers;
  offs = alignSize(offs + nLr * sizeof(MatLayerCyl), MatLayerCyl::getClassAlignmentBytes()); // account for the alignment
 
  // move array of R2 boundaries to the flat array
  const int nRBound = get()->mNRIntervals;
  delete[] o2::gpu::FlatObject::resizeArray(get()->mR2Intervals, nRBound, nRBound, (float*)(mFlatBufferPtr + offs));
  offs = alignSize(offs + nRBound * sizeof(float), getBufferAlignmentBytes()); // account for the alignment
 
  // move array of interval -> layer ID to the flat array
  delete[] o2::gpu::FlatObject::resizeArray(get()->mInterval2LrID, nRBound - 1, nRBound - 1, (int*)(mFlatBufferPtr + offs));
  offs = alignSize(offs + (nRBound - 1) * sizeof(int), getBufferAlignmentBytes()); // account for the alignment
 
  for (int il = 0; il < nLr; il++) {
    MatLayerCyl& lr = get()->mLayers[il];
    lr.flatten(mFlatBufferPtr + offs);
    offs = alignSize(offs + lr.getFlatBufferSize(), getBufferAlignmentBytes()); // account for the alignment
  }
  mConstructionMask = Constructed;
}
 
//______________________________________________
void MatLayerCylSet::moveBufferTo(char* newFlatBufferPtr)
{
  flatObject::moveBufferTo(newFlatBufferPtr);
  setActualBufferAddress(mFlatBufferPtr);
}
#endif // !GPUCA_ALIGPUCODE
 
#ifndef GPUCA_GPUCODE
//______________________________________________
void MatLayerCylSet::setFutureBufferAddress(char* futureFlatBufferPtr)
{
  fixPointers(mFlatBufferPtr, futureFlatBufferPtr, false); // flag that futureFlatBufferPtr is not valid yet
  flatObject::setFutureBufferAddress(futureFlatBufferPtr);
}
 
//______________________________________________
void MatLayerCylSet::setActualBufferAddress(char* actualFlatBufferPtr)
{
  fixPointers(actualFlatBufferPtr);
}
//______________________________________________
void MatLayerCylSet::cloneFromObject(const MatLayerCylSet& obj, char* newFlatBufferPtr)
{
  flatObject::cloneFromObject(obj, newFlatBufferPtr);
  fixPointers(mFlatBufferPtr);
  // the voxel lookup lives outside the flat buffer
  if (obj.mInitializedLayerVoxelLU) {
    std::copy(obj.mLayerVoxelLU, obj.mLayerVoxelLU + NumVoxels, mLayerVoxelLU);
    mInitializedLayerVoxelLU = true;
  }
}
 
//______________________________________________
void MatLayerCylSet::fixPointers(char* newBasePtr)
{
  // fix pointers on the internal structure of the flat buffer after retrieving it from the file
  if (newBasePtr) {
    mFlatBufferPtr = newBasePtr; // used to impose external pointer
  } else {
    mFlatBufferPtr = mFlatBufferContainer; // impose pointer after reading from file
  }
  auto offs = alignSize(sizeof(MatLayerCylSetLayout), getBufferAlignmentBytes()); // account for the alignment
  char* newPtr = mFlatBufferPtr + offs;                                           // correct pointer on MatLayerCyl*
  char* oldPtr = reinterpret_cast<char*>(get()->mLayers);                         // old pointer read from the file
  fixPointers(oldPtr, newPtr);
}
 
//______________________________________________
void MatLayerCylSet::fixPointers(char* oldPtr, char* newPtr, bool newPtrValid)
{
  // fix pointers on the internal structure of the flat buffer after retrieving it from the file
  auto* layPtr = get()->mLayers;
  get()->mLayers = flatObject::relocatePointer(oldPtr, newPtr, get()->mLayers);
  get()->mR2Intervals = flatObject::relocatePointer(oldPtr, newPtr, get()->mR2Intervals);
  get()->mInterval2LrID = flatObject::relocatePointer(oldPtr, newPtr, get()->mInterval2LrID);
  if (newPtrValid) {
    layPtr = get()->mLayers;
  }
  for (int i = 0; i < getNLayers(); i++) {
    layPtr[i].setFlatPointer(flatObject::relocatePointer(oldPtr, newPtr, layPtr[i].getFlatBufferPtr()));
    layPtr[i].fixPointers(oldPtr, newPtr);
  }
}
#endif // !GPUCA_GPUCODE
 
#ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version
 
MatLayerCylSet* MatLayerCylSet::extractCopy(float rmin, float rmax, float tolerance, const MatLayerCylSet* addTo) const
{
  // extract layers in the covering rmin-rmax range. If addTo is provided, simply substitute its layers by those from this
  if (addTo && addTo->getNLayers() != getNLayers()) {
    LOGP(fatal, "addTo has {} layers, this has {}", addTo->getNLayers(), getNLayers());
  }
  Ray ray(std::max(getRMin(), rmin), 0., 0., std::min(getRMax(), rmax), 0., 0.);
  short lmin, lmax;
  if (!getLayersRange(ray, lmin, lmax)) {
    LOGP(warn, "No layers found for {} < r < {}", rmin, rmax);
    return nullptr;
  }
  LOGP(info, "Will extract layers {}:{} (out of {} layers) for {} < r < {}", lmin, lmax, getNLayers(), rmin, rmax);
  MatLayerCylSet* copy = new MatLayerCylSet();
  int lrCount = 0, lrCounOld = 0, lrCountTot = 0;
  auto addLr = [copy, &lrCountTot](const MatLayerCyl& lr) {
    float drphi = lr.getDPhi() * (lr.getRMin() + lr.getRMax()) / 2. * 0.999;
    copy->addLayer(lr.getRMin(), lr.getRMax(), lr.getZMax(), lr.getDZ(), drphi);
    auto& lrNew = copy->getLayer(lrCountTot++);
    for (int iz = 0; iz < lrNew.getNZBins(); iz++) {
      for (int ip = 0; ip < lrNew.getNPhiBins(); ip++) {
        lrNew.getCellPhiBin(ip, iz).set(lr.getCellPhiBin(ip, iz));
      }
    }
  };
  if (addTo) {
    for (int il = 0; il < lmin; il++) {
      addLr(addTo->getLayer(il));
      lrCounOld++;
    }
  }
  for (int il = lmin; il <= lmax; il++) {
    addLr(getLayer(il));
    lrCount++;
  }
  if (addTo) {
    for (int il = lmax + 1; il < getNLayers(); il++) {
      addLr(addTo->getLayer(il));
      lrCounOld++;
    }
  }
  copy->finalizeStructures();
  copy->optimizePhiSlices(tolerance);
  copy->flatten();
  LOGP(info, "Added layers {}:{} for {}<r<{} {}", lmin, lmax, rmin, rmax, fmt::format(", {} layers were transferred from additional set", lrCounOld));
  return copy;
}
 
#endif