diff --git a/Detectors/Base/include/DetectorsBase/MatLayerCyl.h b/Detectors/Base/include/DetectorsBase/MatLayerCyl.h index 04aefed06f010..b5b0fc4db45d1 100644 --- a/Detectors/Base/include/DetectorsBase/MatLayerCyl.h +++ b/Detectors/Base/include/DetectorsBase/MatLayerCyl.h @@ -96,6 +96,7 @@ class MatLayerCyl : public o2::gpu::FlatObject // obtain material cell, cell ID must be valid GPUd() const MatCell& getCellPhiBin(int iphi, int iz) const { return mCells[getCellIDPhiBin(iphi, iz)]; } GPUd() const MatCell& getCell(int iphiSlice, int iz) const { return mCells[getCellID(iphiSlice, iz)]; } + GPUd() const MatCell* getCellRow(int iphiSlice) const { return mCells + iphiSlice * getNZBins(); } #ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version MatCell& getCellPhiBin(int iphi, int iz) diff --git a/Detectors/Base/include/DetectorsBase/MatLayerCylSet.h b/Detectors/Base/include/DetectorsBase/MatLayerCylSet.h index f2e8937ae8448..60c37eb11faeb 100644 --- a/Detectors/Base/include/DetectorsBase/MatLayerCylSet.h +++ b/Detectors/Base/include/DetectorsBase/MatLayerCylSet.h @@ -41,10 +41,10 @@ struct MatLayerCylSetLayout { float mRMin2; ///< precalculater rmin^2 float mRMax2; ///< precalculater rmax^2 int mNLayers; ///< number of layers - int mNRIntervals; ///< number of R interval boundaries (gaps are possible) + int mNRIntervals; ///< number of R interval boundaries, one more than the number of intervals (gaps are possible) MatLayerCyl* mLayers; //[mNLayers] set of cylinrical layers - float* mR2Intervals; //[mNRIntervals+1] limits of layers - int* mInterval2LrID; //[mNRIntervals] mapping from r2 interval to layer ID + float* mR2Intervals; //[mNRIntervals] limits of layers + int* mInterval2LrID; //[mNRIntervals-1] mapping from r2 interval to layer ID }; class MatLayerCylSet : public o2::gpu::FlatObject @@ -115,6 +115,16 @@ class MatLayerCylSet : public o2::gpu::FlatObject /// searches a layer based on r2 input, using a lookup table GPUd() int searchLayerFast(float r2, int low = -1, int high = -1) const; + /// resolves a layer from an already loaded lookup-table entry + GPUd() int resolveLayerRange(float r2, int voxel, uint16_t entry) const; + + /// voxel holding this radius; the caller must have checked that r2 is inside the LUT + GPUd() int voxelIndex(float r2) const { return int(o2::gpu::CAMath::Sqrt(r2) * InvVoxelRDelta); } + + /// Radial boundaries of a voxel. + GPUd() static constexpr float voxelRMin(int voxel) { return voxel * VoxelRDelta; } + GPUd() static constexpr float voxelRMax(int voxel) { return (voxel + 1) * VoxelRDelta; } + #ifndef GPUCA_GPUCODE //----------------------------------------------------------- std::size_t estimateFlatBufferSize() const; @@ -139,8 +149,10 @@ class MatLayerCylSet : public o2::gpu::FlatObject static constexpr float VoxelRDelta = 0.05; // voxel spacing for layer lookup; seems a natural choice - corresponding ~ to smallest spacing static constexpr float InvVoxelRDelta = 1.f / VoxelRDelta; static constexpr int NumVoxels = int(LayerRMax / VoxelRDelta); + static constexpr uint16_t VoxelAmbiguousBit = 0x8000u; + static constexpr uint16_t VoxelSegmentMask = 0x7fffu; - uint16_t mLayerVoxelLU[2 * NumVoxels]; //! helper structure to lookup a layer based on known radius (static dimension for easy copy to GPU) + uint16_t mLayerVoxelLU[NumVoxels]; //! first interval based on known radius, plus the ambiguity flag (static dimension for easy copy to GPU) bool mInitializedLayerVoxelLU = false; //! if the voxels have been initialized ClassDefNV(MatLayerCylSet, 1); diff --git a/Detectors/Base/src/MatLayerCylSet.cxx b/Detectors/Base/src/MatLayerCylSet.cxx index b65a7472fbcf6..2df1d694e4d3c 100644 --- a/Detectors/Base/src/MatLayerCylSet.cxx +++ b/Detectors/Base/src/MatLayerCylSet.cxx @@ -234,7 +234,7 @@ void MatLayerCylSet::finalizeStructures() o2::gpu::FlatObject::resizeArray(get()->mR2Intervals, 0, nR2Int); o2::gpu::FlatObject::resizeArray(get()->mInterval2LrID, 0, nR2Int); get()->mR2Intervals[0] = get()->mRMin2; - get()->mR2Intervals[1] = get()->mRMax2; + get()->mR2Intervals[1] = getLayer(0).getRMax2(); get()->mInterval2LrID[0] = 0; auto& nRIntervals = get()->mNRIntervals; nRIntervals = 1; @@ -319,14 +319,17 @@ void MatLayerCylSet::initLayerVoxelLU() 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 = voxel * VoxelRDelta; - const auto upperR = lowerR + VoxelRDelta; + const auto lowerR = voxelRMin(voxel); + const auto upperR = voxelRMax(voxel); const auto lowerSegment = searchSegment(lowerR * lowerR); const auto upperSegment = searchSegment(upperR * upperR); - mLayerVoxelLU[2 * voxel] = lowerSegment; - mLayerVoxelLU[2 * voxel + 1] = upperSegment; + mLayerVoxelLU[voxel] = uint16_t(lowerSegment) | (lowerSegment != upperSegment ? VoxelAmbiguousBit : uint16_t{0}); } mInitializedLayerVoxelLU = true; } @@ -477,7 +480,10 @@ GPUd() MatBudget MatLayerCylSet::getMatBudget(float x0, float y0, float z0, floa tEndPhi = cross2; checkMorePhi = false; } else { // last phi slice still not reached - tEndPhi = ray.crossRadial(lr, (stepPhiID > 0 ? phiID + 1 : phiID) % nphiSlices); + 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 } @@ -490,6 +496,8 @@ GPUd() MatBudget MatLayerCylSet::getMatBudget(float x0, float y0, float z0, floa } 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_ @@ -512,7 +520,7 @@ GPUd() MatBudget MatLayerCylSet::getMatBudget(float x0, float y0, float z0, floa } // 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 = lr.getCell(phiID % nphiSlices, zID); + const auto& cell = cellRow[zID]; rval.meanRho += cell.meanRho * step; rval.meanX2X0 += cell.meanX2X0 * step; rval.length += step; @@ -523,7 +531,7 @@ GPUd() MatBudget MatLayerCylSet::getMatBudget(float x0, float y0, float z0, floa printf( "Lr#%3d / cross#%d : account %f tStartPhi ? tEndPhi - tStartPhi : tStartPhi - tEndPhi; // the real step is |ray.getDist(tEnd-tStart)|, will rescale all later - const auto& cell = lr.getCell(phiID % nphiSlices, zID); + const auto& cell = cellRow[zID]; rval.meanRho += cell.meanRho * step; rval.meanX2X0 += cell.meanX2X0 * step; rval.length += step; @@ -543,7 +551,7 @@ GPUd() MatBudget MatLayerCylSet::getMatBudget(float x0, float y0, float z0, floa printf( "Lr#%3d / cross#%d : account %fmNRIntervals - 2; lmnInt = rmin2 >= getRMin2() ? searchSegment(rmin2, 0, lmxInt + 1) : 0; } else { - lmxInt = rmax2 < getRMax2() ? searchLayerFast(rmax2, 0) : get()->mNRIntervals - 2; - lmnInt = rmin2 >= getRMin2() ? searchLayerFast(rmin2, 0, lmxInt + 1) : 0; + // 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; @@ -605,11 +620,17 @@ GPUd() bool MatLayerCylSet::getLayersRange(const Ray& ray, short& lmin, short& l 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 = 2 * int(o2::gpu::CAMath::Sqrt(r2) * InvVoxelRDelta); - const auto layersfirst = mLayerVoxelLU[index]; - const auto layerslast = mLayerVoxelLU[index + 1]; - if (layersfirst != layerslast) { - // this means the voxel is undecided and we revert to search + 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; @@ -663,12 +684,13 @@ void MatLayerCylSet::flatten() offs = alignSize(offs + nLr * sizeof(MatLayerCyl), MatLayerCyl::getClassAlignmentBytes()); // account for the alignment // move array of R2 boundaries to the flat array - delete[] o2::gpu::FlatObject::resizeArray(get()->mR2Intervals, nLr + 1, nLr + 1, (float*)(mFlatBufferPtr + offs)); - offs = alignSize(offs + (nLr + 1) * sizeof(float), getBufferAlignmentBytes()); // account for the alignment + 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 R2 boundaries to the flat array - delete[] o2::gpu::FlatObject::resizeArray(get()->mInterval2LrID, nLr, nLr, (int*)(mFlatBufferPtr + offs)); - offs = alignSize(offs + nLr * sizeof(int), 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]; @@ -710,6 +732,11 @@ void MatLayerCylSet::cloneFromObject(const MatLayerCylSet& obj, char* newFlatBuf /// Initializes from another object, copies data to newBufferPtr 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; + } } //______________________________________________ diff --git a/Detectors/Base/test/buildMatBudLUT.C b/Detectors/Base/test/buildMatBudLUT.C index 2b371b90effa3..1ff94e5e3c8c9 100644 --- a/Detectors/Base/test/buildMatBudLUT.C +++ b/Detectors/Base/test/buildMatBudLUT.C @@ -21,6 +21,7 @@ #include #include #include +#include #endif using MatbudGeomBackend = o2::base::MatbudGeomBackend; @@ -30,6 +31,11 @@ o2::base::MatLayerCylSet mbLUT; bool testMBLUT(const std::string& lutFile = "matbud.root"); MatbudGeomBackend parseBackend(const std::string& s); +/// mR2Intervals must be non-decreasing +bool testMBLUTIntervalsSorted(const o2::base::MatLayerCylSet* lut); +/// getLayersRange() must agree with and without the voxel lookup +bool testMBLUTVoxelConsistency(o2::base::MatLayerCylSet* lut, int nRays = 5000); + /// Build the material budget LUT. nThreads < 0 takes the thread count from NTHREADS_MATBUD. /// geomBackend is "ROOT" (default) or "VECGEOM" (requires O2 built against TGeo2VecGeom). bool buildMatBudLUT(int nTst = 60, int maxLr = -1, const std::string& outFile = "matbud.root", @@ -186,6 +192,61 @@ bool testMBLUT(const std::string& lutFile) return true; } +//_______________________________________________________________________ +bool testMBLUTIntervalsSorted(const o2::base::MatLayerCylSet* lut) +{ + // searchSegment() is a binary search over mR2Intervals, enfore order + const auto* layout = lut->get(); + for (int i = 1; i < layout->mNRIntervals; i++) { // mNRIntervals counts boundaries, last index is mNRIntervals-1 + if (layout->mR2Intervals[i] < layout->mR2Intervals[i - 1]) { + LOGP(error, "mR2Intervals not monotonic at {}: {} > {}", i, layout->mR2Intervals[i - 1], layout->mR2Intervals[i]); + return false; + } + } + return true; +} + +//_______________________________________________________________________ +bool testMBLUTVoxelConsistency(o2::base::MatLayerCylSet* lut, int nRays) +{ + // The voxel lookup is only a shortcut into searchSegment(), so it must not change the answer. + if (!lut->mInitializedLayerVoxelLU) { + LOG(error) << "voxel lookup is not initialized, nothing to compare against"; + return false; + } + const float rMax = lut->getRMax(), zMax = lut->getZMax(); + TRandom rnd(20260825); + int nBad = 0, nInside = 0; + for (int i = 0; i < nRays; i++) { + float x0 = rnd.Uniform(-rMax, rMax), y0 = rnd.Uniform(-rMax, rMax), z0 = rnd.Uniform(-zMax, zMax); + float x1 = rnd.Uniform(-rMax, rMax), y1 = rnd.Uniform(-rMax, rMax), z1 = rnd.Uniform(-zMax, zMax); + o2::base::Ray ray(x0, y0, z0, x1, y1, z1); + short lmin = -1, lmax = -1, lminRef = -1, lmaxRef = -1; + const bool ok = lut->getLayersRange(ray, lmin, lmax); + lut->mInitializedLayerVoxelLU = false; // force the plain binary search + const bool okRef = lut->getLayersRange(ray, lminRef, lmaxRef); + lut->mInitializedLayerVoxelLU = true; + if (ok) { + nInside++; + } + if (ok != okRef || (ok && (lmin != lminRef || lmax != lmaxRef))) { + if (++nBad < 10) { + LOGP(error, "ray {} ({:.3f},{:.3f},{:.3f})->({:.3f},{:.3f},{:.3f}): voxel LU gives {} [{},{}], search gives {} [{},{}]", + i, x0, y0, z0, x1, y1, z1, ok, lmin, lmax, okRef, lminRef, lmaxRef); + } + } + } + if (nInside < nRays / 10) { + LOGP(error, "only {} of {} test rays crossed the LUT, the comparison is not meaningful", nInside, nRays); + return false; + } + if (nBad) { + LOGP(error, "{} of {} rays disagree between the voxel lookup and searchSegment()", nBad, nRays); + return false; + } + return true; +} + //_______________________________________________________________________ void configLayers() { diff --git a/Detectors/Base/test/testMatBudLUT.cxx b/Detectors/Base/test/testMatBudLUT.cxx index 33c3498995c90..3199333f0016b 100644 --- a/Detectors/Base/test/testMatBudLUT.cxx +++ b/Detectors/Base/test/testMatBudLUT.cxx @@ -28,5 +28,10 @@ BOOST_AUTO_TEST_CASE(MatBudLUT) matBudFile += std::to_string(getpid()) + ".root"; BOOST_CHECK(buildMatBudLUT(2, 20, matBudFile, geomPrefix + std::to_string(getpid()), "align-geom.mDetectors=none")); // generate LUT BOOST_CHECK(testMBLUT(matBudFile)); // test LUT manipulations + + o2::base::MatLayerCylSet* lut = o2::base::MatLayerCylSet::loadFromFile(matBudFile); + BOOST_REQUIRE(lut != nullptr); + BOOST_CHECK(testMBLUTIntervalsSorted(lut)); // mR2Intervals is monotonic + BOOST_CHECK(testMBLUTVoxelConsistency(lut)); // voxel lookup agrees with the plain search } } // namespace o2