From 8ae4f15919cffca81e0665434f894608e51982ad Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Fri, 25 Sep 2026 13:14:56 +0200 Subject: [PATCH 1/6] TPC SCD: add keepClustersOnPropFail to keep tracks whose reference propagation fails TrackInterpolation drops the whole track when the propagation of the reference track to a TPC pad row (or to the TRD/TOF anchor) fails, e.g. at scdcalib.maxSnp. That propagation follows the ITS-predicted direction, so the drop selects tracks by the reference's own curvature error and biases the residuals of the surviving sample at low pt. With scdcalib.keepClustersOnPropFail=true the track is kept. TPC clusters without a reference are stored position-only: y, z = cluster position, dy = dz = 0, tgSlp = UnbinnedResid::TgSlpPositionOnly (-0x8000, which the tgSlp packing never produces; UnbinnedResid::isPositionOnly()). extrapolateTrack: all clusters after the failure. interpolateTrack: a row keeps its residual only if both the outward and the inward propagation reached it; a failure of the outer TRD/TOF anchor makes all TPC clusters position-only and suppresses the TRD/TOF residuals. Position-only clusters are excluded from validateTrack and stored in row order. The ITS and PV residuals are unaffected. ResidualsContainer::fill, staticMapCreator.C and TPCResidualReaderSpec skip position-only (and tgSlp-clamped) residuals for the binned voxel fit. Default false: output unchanged (verified on one MC TF, all fields identical). Co-Authored-By: Claude Opus 5.5 --- .../src/TPCResidualReaderSpec.cxx | 3 + .../SpacePoints/SpacePointsCalibConfParam.h | 1 + .../include/SpacePoints/TrackInterpolation.h | 6 + .../SpacePoints/macro/staticMapCreator.C | 4 +- .../SpacePoints/src/ResidualAggregator.cxx | 4 +- .../SpacePoints/src/TrackInterpolation.cxx | 195 ++++++++++++------ 6 files changed, 144 insertions(+), 69 deletions(-) diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualReaderSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualReaderSpec.cxx index b3040d99bc4f2..a1c3bc7882c3d 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualReaderSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualReaderSpec.cxx @@ -189,6 +189,9 @@ void TPCResidualReader::run(ProcessingContext& pc) } for (int i = trkInfo.idxFirstResidual; i < trkInfo.idxFirstResidual + trkInfo.nResiduals; ++i) { const auto& residIn = mUnbinnedResiduals[i]; + if (residIn.isTgSlpClamped() || residIn.isPositionOnly()) { + continue; // scdcalib.clampTgSlp / keepClustersOnPropFail: tgSlp or dy, dz not usable for the binned voxel fit + } int sec = residIn.sec; auto& residVecOut = mResidualsSector[sec]; auto& statVecOut = mVoxStatsSector[sec]; diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h index 6ef3839991d04..7bba240552215 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h @@ -49,6 +49,7 @@ struct SpacePointsCalibConfParam : public o2::conf::ConfigurableParamHelper= param::MaxTgSlp saturated (tgSlp = +-0x7fff, see UnbinnedResid::isTgSlpClamped) instead of dropping them: the cut is on the reference track's direction, so dropping selects on the reference's error + bool keepClustersOnPropFail{false}; ///< if the reference track propagation fails (maxSnp, rotation), keep the track: TPC clusters without a reference are stored position-only (y, z = cluster; dy = dz = 0; tgSlp = UnbinnedResid::TgSlpPositionOnly) instead of dropping the whole track, which selects on the reference's error float maxStep{2.f}; ///< maximum step for propagation bool debugTRDTOF{false}; ///< if true, ITS-TPC-TRD-TOF tracks and their seeding ITS-TPC-TRD track will both be interpolated and their residuals stored diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h index 91726aa9941fa..93fd816bfe26e 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h @@ -97,6 +97,10 @@ struct UnbinnedResid { /// true if tgSlp was saturated at +-param::MaxTgSlp (scdcalib.clampTgSlp): unclamped values have |tgSlp| <= 0x7fff - 1 bool isTgSlpClamped() const { return tgSlp == 0x7fff || tgSlp == -0x7fff; } + /// tgSlp marker of a position-only TPC cluster (scdcalib.keepClustersOnPropFail): no reference track at this cluster, + /// y and z are the cluster position, dy = dz = 0. Not reachable by the tgSlp packing (|tgSlp| <= 0x7fff) + static constexpr short TgSlpPositionOnly = -0x8000; + bool isPositionOnly() const { return tgSlp == TgSlpPositionOnly; } bool isTPC() const { return row < constants::MAXGLOBALPADROW; } bool isTRD() const { return row >= 160 && row < 166; } bool isTOF() const { return row == 170; } @@ -500,6 +504,8 @@ class TrackInterpolation size_t mNRejRefit = 0; size_t mNRejProp = 0; size_t mNRejLoop = 0; + size_t mNPosOnlyTracks = 0; ///< tracks kept with position-only clusters after a propagation failure (keepClustersOnPropFail) + size_t mNPosOnlyClusters = 0; ///< position-only TPC clusters stored (keepClustersOnPropFail) ClassDefNV(TrackInterpolation, 1); }; diff --git a/Detectors/TPC/calibration/SpacePoints/macro/staticMapCreator.C b/Detectors/TPC/calibration/SpacePoints/macro/staticMapCreator.C index 3bf01f21f1dce..ee5c54f32acf5 100644 --- a/Detectors/TPC/calibration/SpacePoints/macro/staticMapCreator.C +++ b/Detectors/TPC/calibration/SpacePoints/macro/staticMapCreator.C @@ -320,9 +320,9 @@ void staticMapCreator(std::string fileInput = "files.txt", if (useResidualsForVd && residualsVd.size() < 10'000'000UL) { residualsVd.push_back(residIn); } - if (residIn.isTgSlpClamped()) { + if (residIn.isTgSlpClamped() || residIn.isPositionOnly()) { // scdcalib.clampTgSlp: tgSlp saturated -- the voxel fit (dX from dY vs tan(phi)) and the map correction below - // use it, so keep this residual out of the binned residuals + // use it, so keep this residual out of the binned residuals; scdcalib.keepClustersOnPropFail: no reference (dy = dz = 0) continue; } int sec = residIn.sec; diff --git a/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx b/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx index c5594ddc40b02..7844d219deea5 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx +++ b/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx @@ -197,9 +197,9 @@ void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::sp if (!writeBinnedResid) { continue; } - if (residIn.isTgSlpClamped()) { + if (residIn.isTgSlpClamped() || residIn.isPositionOnly()) { // scdcalib.clampTgSlp: kept in the unbinned output, but its tgSlp is saturated and the voxel fit uses tgSlp (dX from - // dY vs tan(phi)), so it must not enter the binned residuals + // dY vs tan(phi)), so it must not enter the binned residuals; scdcalib.keepClustersOnPropFail: no reference, dy = dz = 0 continue; } int sec = residIn.sec; diff --git a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx index 86fca56e4f8eb..fde816ae40892 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx +++ b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx @@ -465,12 +465,26 @@ void TrackInterpolation::process() } LOGP(info, "Could process {} tracks successfully ({} rejected in refits, {} in propagation, {} as loopers), {} residuals were rejected, {} accepted", mTrackData.size(), mNRejRefit, mNRejProp, mNRejLoop, mRejectedResiduals, mClRes.size()); + if (mParams->keepClustersOnPropFail) { + LOGP(info, "keepClustersOnPropFail: {} tracks kept after a propagation failure, {} position-only TPC clusters stored", mNPosOnlyTracks, mNPosOnlyClusters); + } + mNPosOnlyTracks = 0; + mNPosOnlyClusters = 0; mRejectedResiduals = 0; mNRejRefit = 0; mNRejProp = 0; mNRejLoop = 0; } +namespace +{ +/// TPC cluster stored without a reference track (scdcalib.keepClustersOnPropFail), see UnbinnedResid::isPositionOnly +struct PositionOnlyCluster { + float y, z; + unsigned char sec, row, flags; +}; +} // namespace + void TrackInterpolation::interpolateTrack(int iSeed) { LOGP(debug, "Starting track interpolation for GID {}", mGIDs[iSeed].asString()); @@ -514,6 +528,11 @@ void TrackInterpolation::interpolateTrack(int iSeed) // store the TPC cluster positions in the cache, as well as dedx info std::array, constants::MAXGLOBALPADROW> mCacheDEDX{}; std::array multBins{}; + // keepClustersOnPropFail: a row gets a residual only if both the outward (ITS) and the inward (TRD/TOF) propagation reached + // it; the other TPC clusters are stored position-only. allLost: the outward pass or the outer anchor failed. + bool allLost = false; + std::array refOut{}, refIn{}; + std::vector posOnly; for (int iCl = trkTPC.getNClusterReferences(); iCl--;) { uint8_t sector, row; uint32_t clusterIndexInRow; @@ -545,16 +564,16 @@ void TrackInterpolation::interpolateTrack(int iSeed) if (!mCache[iRow].clAvailable) { continue; } - if (!trkWork.rotate(mCache[iRow].clAngle)) { - LOG(debug) << "Failed to rotate track during first extrapolation"; - mNRejProp++; - return; - } - if (!propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) { + if (!trkWork.rotate(mCache[iRow].clAngle) || !propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) { LOG(debug) << "Failed on first extrapolation"; - mNRejProp++; - return; + if (!mParams->keepClustersOnPropFail) { + mNRejProp++; + return; + } + allLost = true; // no outer anchor can be reached: every TPC cluster is stored position-only + break; } + refOut[iRow] = true; mCache[iRow].y[ExtOut] = trkWork.getY(); mCache[iRow].z[ExtOut] = trkWork.getZ(); mCache[iRow].sy2[ExtOut] = trkWork.getSigmaY2(); @@ -565,7 +584,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) } // start from outermost cluster with outer refit and back propagation - if (gidTable[GTrackID::TOF].isIndexSet()) { + if (!allLost && gidTable[GTrackID::TOF].isIndexSet()) { LOG(debug) << "TOF point available"; const auto& clTOF = mRecoCont->getTOFClusters()[gidTable[GTrackID::TOF]]; if (mDumpTrackPoints) { @@ -574,31 +593,23 @@ void TrackInterpolation::interpolateTrack(int iSeed) } const int clTOFSec = clTOF.getCount(); const float clTOFAlpha = o2::math_utils::sector2Angle(clTOFSec); - if (!trkWork.rotate(clTOFAlpha)) { - LOG(debug) << "Failed to rotate into TOF cluster sector frame"; - mNRejProp++; - return; - } float clTOFxyz[3] = {clTOF.getX(), clTOF.getY(), clTOF.getZ()}; if (!clTOF.isInNominalSector()) { o2::tof::Geo::alignedToNominalSector(clTOFxyz, clTOFSec); // go from the aligned to nominal sector frame } std::array clTOFYZ{clTOFxyz[1], clTOFxyz[2]}; std::array clTOFCov{mParams->sigYZ2TOF, 0.f, mParams->sigYZ2TOF}; // assume no correlation between y and z and equal cluster error sigma^2 = (3cm)^2 / 12 - if (!propagator->PropagateToXBxByBz(trkWork, clTOFxyz[0], mParams->maxSnp, mParams->maxStep, mMatCorr)) { - LOG(debug) << "Failed final propagation to TOF radius"; - mNRejProp++; - return; - } // TODO: check if reset of covariance matrix is needed here (or, in case TOF point is not available at outermost TRD layer) - if (!trkWork.update(clTOFYZ, clTOFCov)) { - LOG(debug) << "Failed to update extrapolated ITS track with TOF cluster"; - // LOGF(info, "trkWork.y=%f, cl.y=%f, trkWork.z=%f, cl.z=%f", trkWork.getY(), clTOFYZ[0], trkWork.getZ(), clTOFYZ[1]); - mNRejProp++; - return; + if (!trkWork.rotate(clTOFAlpha) || !propagator->PropagateToXBxByBz(trkWork, clTOFxyz[0], mParams->maxSnp, mParams->maxStep, mMatCorr) || !trkWork.update(clTOFYZ, clTOFCov)) { + LOG(debug) << "Failed to rotate/propagate/update the extrapolated ITS track at the TOF cluster"; + if (!mParams->keepClustersOnPropFail) { + mNRejProp++; + return; + } + allLost = true; } } - if (gidTable[GTrackID::TRD].isIndexSet()) { + if (!allLost && gidTable[GTrackID::TRD].isIndexSet()) { LOG(debug) << "TRD available"; const auto& trkTRD = mRecoCont->getITSTPCTRDTrack(gidTable[GTrackID::ITSTPCTRD]); if (mDumpTrackPoints) { @@ -611,13 +622,16 @@ void TrackInterpolation::interpolateTrack(int iSeed) if (res == -1) { // no TRD tracklet in this layer continue; } - if (res < -1) { // failed to reach this layer - return; - } - if (!trkWork.update(trkltTRDYZ, trkltTRDCov)) { - LOG(debug) << "Failed to update track at TRD layer " << iLayer; - mNRejProp++; - return; + if (res < -1 || !trkWork.update(trkltTRDYZ, trkltTRDCov)) { // failed to reach this layer or to update + LOG(debug) << "Failed to reach or update the track at TRD layer " << iLayer; + if (!mParams->keepClustersOnPropFail) { + if (res >= -1) { + mNRejProp++; // unchanged: only the update failure was counted + } + return; + } + allLost = true; + break; } } } @@ -629,7 +643,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) // go back through the TPC and store updated track positions bool outerParamStored = false; - for (int iRow = param::NPadRows; iRow--;) { + for (int iRow = param::NPadRows; !allLost && iRow--;) { if (!mCache[iRow].clAvailable) { continue; } @@ -642,17 +656,15 @@ void TrackInterpolation::interpolateTrack(int iSeed) trackData.par = trkWork; outerParamStored = true; } - if (!trkWork.rotate(mCache[iRow].clAngle)) { - LOG(debug) << "Failed to rotate track during back propagation"; - mNRejProp++; - return; - } - if (!propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) { + if (!trkWork.rotate(mCache[iRow].clAngle) || !propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) { LOG(debug) << "Failed on back propagation"; - // printf("trkX(%.2f), clX(%.2f), clY(%.2f), clZ(%.2f), alphaTOF(%.2f)\n", trkWork.getX(), param::RowX[iRow], clTOFYZ[0], clTOFYZ[1], clTOFAlpha); - mNRejProp++; - return; + if (!mParams->keepClustersOnPropFail) { + mNRejProp++; + return; + } + break; // this row and all inner ones have no inward reference: stored position-only } + refIn[iRow] = true; mCache[iRow].y[ExtIn] = trkWork.getY(); mCache[iRow].z[ExtIn] = trkWork.getZ(); mCache[iRow].sy2[ExtIn] = trkWork.getSigmaY2(); @@ -668,6 +680,11 @@ void TrackInterpolation::interpolateTrack(int iSeed) ++deltaRow; continue; } + if (!refOut[iRow] || !refIn[iRow]) { // keepClustersOnPropFail only: no reference at this row + posOnly.push_back({mCache[iRow].clY, mCache[iRow].clZ, mCache[iRow].clSec, (unsigned char)iRow, mCache[iRow].clFlags}); + ++deltaRow; + continue; + } float wTotY = 1.f / mCache[iRow].sy2[ExtOut] + 1.f / mCache[iRow].sy2[ExtIn]; float wTotZ = 1.f / mCache[iRow].sz2[ExtOut] + 1.f / mCache[iRow].sz2[ExtIn]; mCache[iRow].y[Int] = (mCache[iRow].y[ExtOut] / mCache[iRow].sy2[ExtOut] + mCache[iRow].y[ExtIn] / mCache[iRow].sy2[ExtIn]) / wTotY; @@ -713,12 +730,30 @@ void TrackInterpolation::interpolateTrack(int iSeed) mTrackValidation.clear(); // for refitted track parameters and flagging rejected clusters bool stored = false; - trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : validateTrack(trackData, mTrackValidation, clusterResiduals, true); + // keepClustersOnPropFail: a track without any interpolated residual has nothing to validate + trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : ((mParams->keepClustersOnPropFail && clusterResiduals.empty()) ? int8_t(0x1) : validateTrack(trackData, mTrackValidation, clusterResiduals, true)); if (trackData.filterFlag <= 0 || mParams->writeUnfiltered) { int nClValidated = 0; int iRow = 0; + // keepClustersOnPropFail: store the position-only clusters in row order between the residuals + size_t iPosOnly = 0; + auto flushPosOnly = [&](int rowLimit) { + for (; iPosOnly < posOnly.size() && posOnly[iPosOnly].row < rowLimit; ++iPosOnly) { + const auto& pc = posOnly[iPosOnly]; + if (std::abs(pc.y) < param::MaxY && std::abs(pc.z) < param::MaxZ) { + mClRes.emplace_back(0.f, 0.f, 0.f, pc.y, pc.z, pc.row, pc.sec, pc.flags, false); + mClRes.back().tgSlp = UnbinnedResid::TgSlpPositionOnly; + mDetInfoRes.emplace_back().setTPC(mCacheDEDX[pc.row].first, mCacheDEDX[pc.row].second); // qtot, qmax + ++nClValidated; + ++mNPosOnlyClusters; + } else { + ++mRejectedResiduals; + } + } + }; for (unsigned int iCl = 0; iCl < clusterResiduals.size(); ++iCl) { iRow += clusterResiduals[iCl].dRow; + flushPosOnly(iRow); const auto rej = trackData.filterFlag < 0 ? false : mTrackValidation.points[iCl].flagRej; if (rej && !mParams->keepRejectedResiduals) { // skip masked cluster residual continue; @@ -741,6 +776,10 @@ void TrackInterpolation::interpolateTrack(int iSeed) ++mRejectedResiduals; } } + flushPosOnly(constants::MAXGLOBALPADROW); + if (!posOnly.empty()) { + ++mNPosOnlyTracks; + } trackData.clIdx.setEntries(nClValidated); // store multiplicity info @@ -771,7 +810,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) if (!stopPropagation) { // do we have TRD residuals to add? trkWork = trkOuter; - if (gidTable[GTrackID::TRD].isIndexSet()) { + if (!allLost && gidTable[GTrackID::TRD].isIndexSet()) { // allLost: trkOuter is not a valid outer param const auto& trkTRD = mRecoCont->getITSTPCTRDTrack(gidTable[GTrackID::ITSTPCTRD]); for (int iLayer = 0; iLayer < o2::trd::constants::NLAYER; iLayer++) { std::array trkltTRDYZ{}; @@ -796,7 +835,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) } // do we have TOF residual to add? - while (gidTable[GTrackID::TOF].isIndexSet() && !stopPropagation) { + while (!allLost && gidTable[GTrackID::TOF].isIndexSet() && !stopPropagation) { const auto& clTOF = mRecoCont->getTOFClusters()[gidTable[GTrackID::TOF]]; float clTOFxyz[3] = {clTOF.getX(), clTOF.getY(), clTOF.getZ()}; if (!clTOF.isInNominalSector()) { @@ -996,6 +1035,8 @@ void TrackInterpolation::extrapolateTrack(int iSeed) uint8_t clRowPrev = constants::MAXGLOBALPADROW; // used to identify and skip split clusters on the same pad row std::array, constants::MAXGLOBALPADROW> mCacheDEDX{}; std::array multBins{}; + bool refLost = false; // keepClustersOnPropFail: the ITS extrapolation failed at an earlier cluster + std::vector posOnly; // keepClustersOnPropFail: clusters after the failure, rows ascending for (int iCl = trkTPC.getNClusterReferences(); iCl--;) { uint8_t sector, row; uint32_t clusterIndexInRow; @@ -1017,29 +1058,31 @@ void TrackInterpolation::extrapolateTrack(int iSeed) } float x = 0, y = 0, z = 0; mFastTransform->TransformIdeal(sector, row, cl.getPad(), cl.getTime(), x, y, z, clusterTimeBinOffset); - if (!trkWork.rotate(o2::math_utils::sector2Angle(sector))) { - mNRejProp++; - return; - } - if (!propagator->PropagateToXBxByBz(trkWork, x, mParams->maxSnp, mParams->maxStep, mMatCorr)) { - mNRejProp++; - return; - } - - const auto dY = y - trkWork.getY(); - const auto dZ = z - trkWork.getZ(); - const auto ty = trkWork.getY(); - const auto tz = trkWork.getZ(); - const auto snp = trkWork.getSnp(); - const auto sec = sector; unsigned char flags = cl.getFlags(); if (mTPCShClassMap[absoluteIndex] & o2::gpu::GPUTPCGMMergedTrackHit::flagShared) { flags |= o2::gpu::GPUTPCGMMergedTrackHit::flagShared; } - clusterResiduals.emplace_back(dY, dZ, ty, tz, snp, sec, row - rowPrev, flags); + if (!refLost && !(trkWork.rotate(o2::math_utils::sector2Angle(sector)) && propagator->PropagateToXBxByBz(trkWork, x, mParams->maxSnp, mParams->maxStep, mMatCorr))) { + if (!mParams->keepClustersOnPropFail) { + mNRejProp++; + return; + } + refLost = true; // scdcalib.keepClustersOnPropFail: this and all further clusters are stored position-only + } mCacheDEDX[row].first = cl.getQtot(); mCacheDEDX[row].second = cl.getQmax(); - rowPrev = row; + if (refLost) { + posOnly.push_back({y, z, sector, row, flags}); + } else { + const auto dY = y - trkWork.getY(); + const auto dZ = z - trkWork.getZ(); + const auto ty = trkWork.getY(); + const auto tz = trkWork.getZ(); + const auto snp = trkWork.getSnp(); + const auto sec = sector; + clusterResiduals.emplace_back(dY, dZ, ty, tz, snp, sec, row - rowPrev, flags); + rowPrev = row; + } int imb = int(cl.getTime() * mNTPCOccBinLengthInv); if (imb < mTPCParam->occupancyMapSize) { multBins[row] = 1 + std::max(0, imb); @@ -1065,12 +1108,30 @@ void TrackInterpolation::extrapolateTrack(int iSeed) } bool stored = false; - trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : validateTrack(trackData, mTrackValidation, clusterResiduals, false); + // keepClustersOnPropFail: a track that lost its reference before the first cluster has no residual to validate + trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : ((mParams->keepClustersOnPropFail && clusterResiduals.empty()) ? int8_t(0x1) : validateTrack(trackData, mTrackValidation, clusterResiduals, false)); if (trackData.filterFlag <= 0 || mParams->writeUnfiltered) { int nClValidated = 0, iRow = 0; unsigned int iCl = 0; + // keepClustersOnPropFail: store the position-only clusters in row order between the residuals + size_t iPosOnly = 0; + auto flushPosOnly = [&](int rowLimit) { + for (; iPosOnly < posOnly.size() && posOnly[iPosOnly].row < rowLimit; ++iPosOnly) { + const auto& pc = posOnly[iPosOnly]; + if (std::abs(pc.y) < param::MaxY && std::abs(pc.z) < param::MaxZ) { + mClRes.emplace_back(0.f, 0.f, 0.f, pc.y, pc.z, pc.row, pc.sec, pc.flags, false); + mClRes.back().tgSlp = UnbinnedResid::TgSlpPositionOnly; + mDetInfoRes.emplace_back().setTPC(mCacheDEDX[pc.row].first, mCacheDEDX[pc.row].second); // qtot, qmax + ++nClValidated; + ++mNPosOnlyClusters; + } else { + ++mRejectedResiduals; + } + } + }; for (iCl = 0; iCl < clusterResiduals.size(); ++iCl) { iRow += clusterResiduals[iCl].dRow; + flushPosOnly(iRow); if (iRow >= param::NPadRows) { // RS why do we need this? continue; } @@ -1095,6 +1156,10 @@ void TrackInterpolation::extrapolateTrack(int iSeed) ++mRejectedResiduals; } } + flushPosOnly(constants::MAXGLOBALPADROW); + if (!posOnly.empty()) { + ++mNPosOnlyTracks; + } trackData.clIdx.setEntries(nClValidated); // store multiplicity info @@ -1127,7 +1192,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed) int iSeedFull = mParentID[iSeed] == -1 ? iSeed : mParentID[iSeed]; auto gidFull = mGIDs[iSeedFull]; const auto& gidTableFull = mGIDtables[iSeedFull]; - if (gidTableFull[GTrackID::TRD].isIndexSet()) { + if (!refLost && gidTableFull[GTrackID::TRD].isIndexSet()) { // refLost: trkWork did not reach the TPC outer end const auto& trkTRD = mRecoCont->getITSTPCTRDTrack(gidTableFull[GTrackID::ITSTPCTRD]); trackData.nTrkltsTRD = trkTRD.getNtracklets(); trackData.chi2TRD = trkTRD.getChi2(); @@ -1156,7 +1221,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed) // do we have TOF residual to add? trackData.clAvailTOF = 0; - while (gidTableFull[GTrackID::TOF].isIndexSet() && !stopPropagation) { + while (!refLost && gidTableFull[GTrackID::TOF].isIndexSet() && !stopPropagation) { const auto& tofMatch = mRecoCont->getTOFMatch(gidFull); ULong64_t bclongtof = (tofMatch.getSignal() - 10000) * o2::tof::Geo::BC_TIME_INPS_INV; double t0forTOF = tofMatch.getFT0Best(); // setting t0 for TOF From 2ce8a007d9209fda13b768464ed70d680a624dbd Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Sat, 3 Oct 2026 09:55:16 +0200 Subject: [PATCH 2/6] TPC SCD: store the MC truth of the residual tracks (--enable-mc) The interpolation workflow had MC switched off ("not yet implemented"). With --enable-mc (opt-in, since the workflow also runs on data without --disable-mc) and --send-track-data, a TrackDataMC vector aligned 1:1 with the TrackData output is sent as GLO/TRKDATAMC: MC labels of the ITS-TPC part of the seeding track and of its ITS and TPC parts (flag for a fake ITS-TPC match), the truth at the ITS outer parameters (nearest ITS track reference, propagated with the workflow's material correction and the true mass to the x and alpha of TrackData::par, with the distance to that reference), the truth at the TPC entrance (first TPC track reference in time, sector frame) with the distance between it, propagated to the innermost TPC cluster of the track, and that cluster (flags loopers, wrong legs and fakes), and the PDG code. Only the ITS, TPC and ITS-TPC track labels are read; the kinematics are looked up event by event and released after each event. The residual aggregator writes the truth as branch trkMC of the trackData tree with --enable-mc. Co-Authored-By: Claude Opus 5.5 --- .../TPCInterpolationSpec.h | 7 +- .../TPCResidualAggregatorSpec.h | 21 +- .../src/TPCInterpolationSpec.cxx | 193 +++++++++++++++++- .../src/tpc-interpolation-workflow.cxx | 11 +- .../src/tpc-residual-aggregator.cxx | 11 +- .../include/SpacePoints/ResidualAggregator.h | 10 +- .../include/SpacePoints/TrackInterpolation.h | 21 ++ .../SpacePoints/src/ResidualAggregator.cxx | 18 +- .../SpacePoints/src/SpacePointCalibLinkDef.h | 2 + 9 files changed, 271 insertions(+), 23 deletions(-) diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h index bb3ebef84032c..953094b5155b7 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h @@ -23,12 +23,14 @@ #include "DetectorsBase/GRPGeomHelper.h" #include "TPCCalibration/VDriftHelper.h" #include "DataFormatsITSMFT/TopologyDictionary.h" +#include "Steer/MCKinematicsReader.h" using namespace o2::framework; namespace o2::globaltracking { struct DataRequest; +struct RecoContainer; } // namespace o2::globaltracking namespace o2 @@ -40,7 +42,7 @@ class TPCInterpolationDPL : public Task public: TPCInterpolationDPL(std::shared_ptr dr, o2::dataformats::GlobalTrackID::mask_t src, o2::dataformats::GlobalTrackID::mask_t srcMap, std::shared_ptr gr, bool useMC, bool processITSTPConly, bool sendTrackData, bool debugOutput, bool extDetResid) : mDataRequest(dr), mSources(src), mSourcesMap(srcMap), mGGCCDBRequest(gr), mUseMC(useMC), mProcessITSTPConly(processITSTPConly), mSendTrackData(sendTrackData), mDebugOutput(debugOutput), mExtDetResid(extDetResid) {} - ~TPCInterpolationDPL() override = default; + ~TPCInterpolationDPL() override; void init(InitContext& ic) final; void run(ProcessingContext& pc) final; void endOfStream(EndOfStreamContext& ec) final; @@ -48,6 +50,7 @@ class TPCInterpolationDPL : public Task private: void updateTimeDependentParams(ProcessingContext& pc); + void fillMCTruth(const o2::globaltracking::RecoContainer& recoData); o2::tpc::TrackInterpolation mInterpolation; ///< track interpolation engine std::shared_ptr mDataRequest; ///< steers the input std::shared_ptr mGGCCDBRequest; @@ -56,6 +59,8 @@ class TPCInterpolationDPL : public Task o2::dataformats::GlobalTrackID::mask_t mSources{}; ///< which input sources are configured o2::dataformats::GlobalTrackID::mask_t mSourcesMap{}; ///< possible subset of mSources specifically for map creation bool mUseMC{false}; ///< MC flag + std::unique_ptr mMCReader; ///< MC kinematics and track references (MC only) + std::vector mTrackDataMC; ///< MC truth aligned with the TrackData output (MC only) bool mProcessITSTPConly{false}; ///< should also tracks without outer point (ITS-TPC only) be processed? bool mProcessSeeds{false}; ///< process not only most complete track, but also its shorter parts bool mDebugOutput{false}; ///< add more information to the output (track points of ITS, TRD and TOF) diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualAggregatorSpec.h b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualAggregatorSpec.h index 99f20e390a09a..e17e8b49509b4 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualAggregatorSpec.h +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualAggregatorSpec.h @@ -41,7 +41,7 @@ namespace calibration class ResidualAggregatorDevice : public o2::framework::Task { public: - ResidualAggregatorDevice(std::shared_ptr req, bool trackInput, bool ctpInput, bool writeUnbinnedResiduals, bool writeBinnedResiduals, bool writeTrackData, std::shared_ptr dataRequest) : mCCDBRequest(req), mTrackInput(trackInput), mCTPInput(ctpInput), mWriteUnbinnedResiduals(writeUnbinnedResiduals), mWriteBinnedResiduals(writeBinnedResiduals), mWriteTrackData(writeTrackData), mDataRequest(dataRequest) {} + ResidualAggregatorDevice(std::shared_ptr req, bool trackInput, bool ctpInput, bool writeUnbinnedResiduals, bool writeBinnedResiduals, bool writeTrackData, bool mcInput, std::shared_ptr dataRequest) : mCCDBRequest(req), mTrackInput(trackInput), mCTPInput(ctpInput), mWriteUnbinnedResiduals(writeUnbinnedResiduals), mWriteBinnedResiduals(writeBinnedResiduals), mWriteTrackData(writeTrackData), mMCInput(mcInput), mDataRequest(dataRequest) {} void init(o2::framework::InitContext& ic) final { @@ -97,6 +97,7 @@ class ResidualAggregatorDevice : public o2::framework::Task mAggregator->setWriteBinnedResiduals(mWriteBinnedResiduals); mAggregator->setWriteUnbinnedResiduals(mWriteUnbinnedResiduals); mAggregator->setWriteTrackData(mWriteTrackData); + mAggregator->setWriteTrackDataMC(mMCInput); mAggregator->setCompression(ic.options().get("compression")); } @@ -141,6 +142,14 @@ class ResidualAggregatorDevice : public o2::framework::Task trkData.emplace(pc.inputs().get>("trkData")); trkDataPtr = &trkData.value(); } + // MC truth of the track data (optional, MC only) + const gsl::span* trkDataMCPtr = nullptr; + using trkDataMCType = std::decay_t>(""))>; + std::optional trkDataMC; + if (mMCInput) { + trkDataMC.emplace(pc.inputs().get>("trkDataMC")); + trkDataMCPtr = &trkDataMC.value(); + } // CTP lumi input (optional) const o2::ctp::LumiInfo* lumi = nullptr; using lumiDataType = std::decay_t(""))>; @@ -152,7 +161,7 @@ class ResidualAggregatorDevice : public o2::framework::Task o2::base::TFIDInfoHelper::fillTFIDInfo(pc, mAggregator->getCurrentTFInfo()); LOG(detail) << "Processing TF " << mAggregator->getCurrentTFInfo().tfCounter << " with " << trkData->size() << " tracks and " << residualsData.size() << " unbinned residuals associated to them"; - mAggregator->process(residualsData, residualsDataDet, trackRefs, trkDataPtr, lumi); + mAggregator->process(residualsData, residualsDataDet, trackRefs, trkDataPtr, trkDataMCPtr, lumi); std::chrono::duration runDuration = std::chrono::high_resolution_clock::now() - runStartTime; LOGP(debug, "Duration for run method: {} ms. From this taken for time dependent param update: {} ms", std::chrono::duration_cast(runDuration).count(), @@ -205,6 +214,7 @@ class ResidualAggregatorDevice : public o2::framework::Task bool mWriteBinnedResiduals{false}; ///< flag, whether to write binned residuals to output file bool mWriteUnbinnedResiduals{false}; ///< flag, whether to write unbinned residuals to output file bool mWriteTrackData{false}; ///< flag, whether to write track data to output file + bool mMCInput{false}; ///< flag whether to expect the MC truth of the track data as input bool mRunStopRequested{false}; ///< flag in case the run was stopped bool mInitDone{false}; ///< flag whether initialization was done for current run }; @@ -214,7 +224,7 @@ class ResidualAggregatorDevice : public o2::framework::Task namespace framework { -DataProcessorSpec getTPCResidualAggregatorSpec(bool trackInput, bool ctpInput, bool writeUnbinnedResiduals, bool writeBinnedResiduals, bool writeTrackData) +DataProcessorSpec getTPCResidualAggregatorSpec(bool trackInput, bool ctpInput, bool writeUnbinnedResiduals, bool writeBinnedResiduals, bool writeTrackData, bool mcInput = false) { std::shared_ptr dataRequest = std::make_shared(); if (ctpInput) { @@ -227,6 +237,9 @@ DataProcessorSpec getTPCResidualAggregatorSpec(bool trackInput, bool ctpInput, b inputs.emplace_back("trackRefs", "GLO", "TRKREFS"); if (trackInput) { inputs.emplace_back("trkData", "GLO", "TRKDATA"); + if (mcInput) { + inputs.emplace_back("trkDataMC", "GLO", "TRKDATAMC"); + } } auto ccdbRequest = std::make_shared(true, // orbitResetTime true, // GRPECS=true @@ -240,7 +253,7 @@ DataProcessorSpec getTPCResidualAggregatorSpec(bool trackInput, bool ctpInput, b "residual-aggregator", inputs, Outputs{}, - AlgorithmSpec{adaptFromTask(ccdbRequest, trackInput, ctpInput, writeUnbinnedResiduals, writeBinnedResiduals, writeTrackData, dataRequest)}, + AlgorithmSpec{adaptFromTask(ccdbRequest, trackInput, ctpInput, writeUnbinnedResiduals, writeBinnedResiduals, writeTrackData, trackInput && mcInput, dataRequest)}, Options{ {"sec-per-slot", VariantType::UInt32, 600u, {"number of seconds per calibration time slot (put 0 for infinite slot length)"}}, {"updateInterval", VariantType::UInt32, 6'000u, {"update interval in number of TFs (only used in case slot length is infinite)"}}, diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx index 2af349be4fd37..c95bf421b7dba 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx @@ -13,6 +13,8 @@ #include #include +#include +#include #include "DataFormatsITS/TrackITS.h" #include "ITSBase/GeometryTGeo.h" @@ -33,6 +35,9 @@ #include "Framework/ConfigParamRegistry.h" #include "Framework/ControlService.h" #include "Framework/DeviceSpec.h" +#include "Steer/MCKinematicsReader.h" +#include "SimulationDataFormat/TrackReference.h" +#include "SimulationDataFormat/O2DatabasePDG.h" using namespace o2::framework; using namespace o2::globaltracking; @@ -44,6 +49,8 @@ namespace o2 namespace tpc { +TPCInterpolationDPL::~TPCInterpolationDPL() = default; + void TPCInterpolationDPL::init(InitContext& ic) { //-------- init geometry and field --------// @@ -59,6 +66,16 @@ void TPCInterpolationDPL::init(InitContext& ic) int lane = ic.services().get().inputTimesliceId; int maxLanes = ic.services().get().maxInputTimeslices; mInterpolation.setLane(lane, maxLanes); + if (mUseMC) { + if (!mSendTrackData) { + LOG(warning) << "MC truth is stored aligned with the track data, but send-track-data is not set: no MC truth will be sent"; + } + mMCReader = std::make_unique(); + auto mcContext = ic.options().get("mc-collision-context"); + if (!mMCReader->initFromDigitContext(mcContext)) { + LOG(fatal) << "Could not initialize the MC kinematics reader from " << mcContext; + } + } } void TPCInterpolationDPL::updateTimeDependentParams(ProcessingContext& pc) @@ -149,9 +166,168 @@ void TPCInterpolationDPL::run(ProcessingContext& pc) if (mDebugOutput) { pc.outputs().snapshot(Output{"GLO", "TRKDATAEXT", 0}, mInterpolation.getTrackDataExtended()); } + if (mUseMC && mSendTrackData) { + fillMCTruth(recoData); + pc.outputs().snapshot(Output{"GLO", "TRKDATAMC", 0}, mTrackDataMC); + } mInterpolation.reset(); } +void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) +{ + // MC truth for every stored TrackData: labels of the ITS-TPC part of the seed and of its ITS and TPC parts, the truth at + // the ITS outer parameters (the ITS track reference nearest to TrackData::par, propagated to its x with the material + // correction of the workflow and the mass of the true particle) and the truth at the TPC entrance (the first TPC track + // reference in time, in the sector frame) + const auto& trkData = mInterpolation.getReferenceTracks(); + mTrackDataMC.clear(); + mTrackDataMC.resize(trkData.size()); + struct Lookup { + o2::MCCompLabel lbl; + uint32_t idx; + bool its; // ITS outer (true) or TPC entrance (false) + }; + std::vector lookups; + lookups.reserve(2 * trkData.size()); + for (size_t i = 0; i < trkData.size(); ++i) { + auto& mc = mTrackDataMC[i]; + auto gidSet = recoData.getSingleDetectorRefs(trkData[i].gid); + auto gidITS = gidSet[GTrackID::ITS].isIndexSet() ? gidSet[GTrackID::ITS] : gidSet[GTrackID::ITSAB]; + if (gidSet[GTrackID::ITSTPC].isIndexSet()) { + mc.label = recoData.getTrackMCLabel(gidSet[GTrackID::ITSTPC]); + } + if (gidITS.isIndexSet()) { + mc.labelITS = recoData.getTrackMCLabel(gidITS); + } + if (gidSet[GTrackID::TPC].isIndexSet()) { + mc.labelTPC = recoData.getTrackMCLabel(gidSet[GTrackID::TPC]); + } + if (mc.labelITS.isValid() && mc.labelTPC.isValid() && mc.labelITS.getTrackEventSourceID() != mc.labelTPC.getTrackEventSourceID()) { + mc.flags |= TrackDataMC::FakeITSTPC; + } + const auto& lblITS = mc.labelITS.isValid() ? mc.labelITS : mc.label; + const auto& lblTPC = mc.labelTPC.isValid() ? mc.labelTPC : mc.label; + if (lblITS.isValid()) { + lookups.push_back({lblITS, uint32_t(i), true}); + } + if (lblTPC.isValid()) { + lookups.push_back({lblTPC, uint32_t(i), false}); + } + } + // the reader loads the kinematics of a whole event (can be >100 MB): process event by event and release it right after + std::sort(lookups.begin(), lookups.end(), [](const Lookup& a, const Lookup& b) { + return a.lbl.getSourceID() != b.lbl.getSourceID() ? a.lbl.getSourceID() < b.lbl.getSourceID() : a.lbl.getEventID() < b.lbl.getEventID(); + }); + auto pdgToPID = [](int pdg) { + switch (std::abs(pdg)) { + case 11: + return o2::track::PID(o2::track::PID::Electron); + case 13: + return o2::track::PID(o2::track::PID::Muon); + case 321: + return o2::track::PID(o2::track::PID::Kaon); + case 2212: + return o2::track::PID(o2::track::PID::Proton); + case 1000010020: + return o2::track::PID(o2::track::PID::Deuteron); + case 1000010030: + return o2::track::PID(o2::track::PID::Triton); + case 1000020030: + return o2::track::PID(o2::track::PID::Helium3); + case 1000020040: + return o2::track::PID(o2::track::PID::Alpha); + default: + return o2::track::PID(o2::track::PID::Pion); + } + }; + auto refToPar = [&pdgToPID](const o2::TrackReference& ref, int charge, int pdg, bool sectorAlpha) { + std::array xyz{ref.X(), ref.Y(), ref.Z()}; + std::array pxyz{ref.Px(), ref.Py(), ref.Pz()}; + return o2::track::TrackPar(xyz, pxyz, charge, sectorAlpha, pdgToPID(pdg)); + }; + const auto matCorr = static_cast(mMatCorr); + auto prop = o2::base::Propagator::Instance(); + int curSrc = -1; + int curEv = -1; + for (const auto& lk : lookups) { + const auto& lbl = lk.lbl; + if (lbl.getSourceID() != curSrc || lbl.getEventID() != curEv) { + if (curSrc >= 0) { + mMCReader->releaseTracksForSourceAndEvent(curSrc, curEv); + } + curSrc = lbl.getSourceID(); + curEv = lbl.getEventID(); + } + const auto& trk = trkData[lk.idx]; + auto& mc = mTrackDataMC[lk.idx]; + const auto* mcTrk = mMCReader->getTrack(lbl); + int pdg = mcTrk ? mcTrk->GetPdgCode() : 0; + const auto* pPDG = mcTrk ? O2DatabasePDG::Instance()->GetParticle(pdg) : nullptr; + int charge = pPDG ? int(std::lround(pPDG->Charge() / 3.)) : 0; // TParticlePDG charge is in units of |e|/3 + auto refs = mMCReader->getTrackRefs(lbl.getSourceID(), lbl.getEventID(), lbl.getTrackID()); + if (lk.its) { // ITS outer: track reference of the ITS part nearest to TrackData::par + mc.pdg = pdg; + const o2::TrackReference* best = nullptr; + float bestD2 = 1e30f; + auto xyzReco = trk.par.getXYZGlo(); + for (const auto& ref : refs) { + if (ref.getDetectorId() != DetID::ITS) { + continue; + } + float dx = ref.X() - xyzReco.X(); + float dy = ref.Y() - xyzReco.Y(); + float dz = ref.Z() - xyzReco.Z(); + float d2 = dx * dx + dy * dy + dz * dz; + if (d2 < bestD2) { + bestD2 = d2; + best = &ref; + } + } + if (best && charge) { + auto par = refToPar(*best, charge, pdg, false); + if (par.rotateParam(trk.par.getAlpha()) && prop->PropagateToXBxByBz(par, trk.par.getX(), 0.999f, o2::base::Propagator::MAX_STEP, matCorr)) { + mc.parITSOut = par; + mc.distITSRef = std::sqrt(bestD2); + mc.flags |= TrackDataMC::HasITSOut; + } + } + } else { // TPC entrance: first TPC track reference in time of the TPC part + if (!mc.labelITS.isValid() && !mc.label.isValid()) { + mc.pdg = pdg; // no ITS lookup for this track + } + const o2::TrackReference* first = nullptr; + for (const auto& ref : refs) { + if (ref.getDetectorId() == DetID::TPC && (!first || ref.getTime() < first->getTime())) { + first = &ref; + } + } + if (first && charge) { + mc.parTPCIn = refToPar(*first, charge, pdg, true); + mc.flags |= TrackDataMC::HasTPCIn; + // association check (loopers, wrong leg, fake): truth at the innermost TPC cluster of the track vs that cluster + const auto& clRes = mInterpolation.getClusterResiduals(); + const UnbinnedResid* inner = nullptr; + for (int ic = trk.clIdx.getFirstEntry(); ic < trk.clIdx.getFirstEntry() + trk.clIdx.getEntries(); ++ic) { + if (clRes[ic].row < constants::MAXGLOBALPADROW && (!inner || clRes[ic].row < inner->row)) { + inner = &clRes[ic]; + } + } + if (inner) { + auto par = mc.parTPCIn; + float yCl = inner->y * param::MaxY / 0x7fff + inner->dy * param::MaxResid / 0x7fff; + float zCl = inner->z * param::MaxZ / 0x7fff + inner->dz * param::MaxResid / 0x7fff; + if (par.rotateParam(o2::math_utils::sector2Angle(inner->sec)) && prop->PropagateToXBxByBz(par, param::RowX[inner->row], 0.999f, o2::base::Propagator::MAX_STEP, matCorr)) { + mc.distTPCRef = std::hypot(par.getY() - yCl, par.getZ() - zCl); + } + } + } + } + } + if (curSrc >= 0) { + mMCReader->releaseTracksForSourceAndEvent(curSrc, curEv); + } +} + void TPCInterpolationDPL::endOfStream(EndOfStreamContext& ec) { mInterpolation.finalize(); @@ -165,14 +341,13 @@ DataProcessorSpec getTPCInterpolationSpec(GTrackID::mask_t srcCls, GTrackID::mas dataRequest->setITSPerLayer(itsStag); std::vector outputs; - if (useMC) { - LOG(fatal) << "MC usage must be disabled for this workflow, since it is not yet implemented"; + dataRequest->requestTracks(srcVtx, false); + dataRequest->requestClusters(srcCls, false); + dataRequest->requestPrimaryVertices(false); + if (useMC) { // the MC truth needs only the labels of the ITS-TPC tracks and of their ITS and TPC parts + dataRequest->requestTracks(GTrackID::getSourcesMask("ITS,TPC,ITS-TPC"), true); } - dataRequest->requestTracks(srcVtx, useMC); - dataRequest->requestClusters(srcCls, useMC); - dataRequest->requestPrimaryVertices(useMC); - auto ggRequest = std::make_shared(false, // orbitResetTime true, // GRPECS=true true, // GRPLHCIF @@ -191,6 +366,9 @@ DataProcessorSpec getTPCInterpolationSpec(GTrackID::mask_t srcCls, GTrackID::mas if (debugOutput) { outputs.emplace_back("GLO", "TRKDATAEXT", 0, Lifetime::Timeframe); } + if (useMC && sendTrackData) { + outputs.emplace_back("GLO", "TRKDATAMC", 0, Lifetime::Timeframe); + } return DataProcessorSpec{ "tpc-track-interpolation", @@ -200,7 +378,8 @@ DataProcessorSpec getTPCInterpolationSpec(GTrackID::mask_t srcCls, GTrackID::mas Options{ {"matCorrType", VariantType::Int, 2, {"material correction type (definition in Propagator.h)"}}, {"sec-per-slot", VariantType::UInt32, 300u, {"number of seconds per calibration time slot (put 0 for infinite slot length)"}}, - {"process-seeds", VariantType::Bool, false, {"do not remove duplicates, e.g. for ITS-TPC-TRD track also process its seeding ITS-TPC part"}}}}; + {"process-seeds", VariantType::Bool, false, {"do not remove duplicates, e.g. for ITS-TPC-TRD track also process its seeding ITS-TPC part"}}, + {"mc-collision-context", VariantType::String, "collisioncontext.root", {"collision context used to access the MC kinematics and track references (MC only)"}}}}; } } // namespace tpc diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx index e8b6aaa07eaba..d43c10ce9c74b 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx @@ -38,6 +38,7 @@ void customize(std::vector& workflowOptions) {"disable-root-input", VariantType::Bool, false, {"disable root-files input readers"}}, {"disable-root-output", VariantType::Bool, false, {"disable root-files output writers"}}, {"disable-mc", VariantType::Bool, false, {"disable MC propagation even if available"}}, + {"enable-mc", VariantType::Bool, false, {"store the MC truth of the track data (needs MC input and send-track-data)"}}, {"vtx-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of sources used for the vertex finding"}}, {"tracking-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of sources to use for track inter-/extrapolation"}}, {"tracking-sources-map-extraction", VariantType::String, std::string{GID::ALL}, {"can be subset of \"tracking-sources\""}}, @@ -103,8 +104,8 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); // write the configuration used for the workflow o2::conf::ConfigurableParam::writeINI("o2tpcinterpolation-workflow_configuration.ini"); - auto useMC = !configcontext.options().get("disable-mc"); - useMC = false; // force disabling MC as long as it is not implemented + // MC is opt-in: the residuals workflow is also run on data without passing disable-mc + auto useMC = configcontext.options().get("enable-mc") && !configcontext.options().get("disable-mc"); auto doStag = o2::itsmft::DPLAlpideParamInitializer::isITSStaggeringEnabled(configcontext); auto sendTrackData = configcontext.options().get("send-track-data"); auto debugOutput = configcontext.options().get("debug-output"); @@ -115,8 +116,10 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) specs.emplace_back(o2::tpc::getTPCResidualWriterSpec(sendTrackData, debugOutput)); } - o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcClusters, srcVtx, srcVtx, useMC); - o2::globaltracking::InputHelper::addInputSpecsPVertex(configcontext, specs, useMC); // P-vertex is always needed + // MC labels only for the ITS-TPC tracks and their ITS and TPC parts (see getTPCInterpolationSpec) + GID::mask_t maskTracksMC = useMC ? GID::getSourcesMask("ITS,TPC,ITS-TPC") : GID::getSourcesMask(GID::NONE); + o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcClusters, srcVtx, srcVtx, useMC, GID::getSourcesMask(GID::NONE), maskTracksMC); + o2::globaltracking::InputHelper::addInputSpecsPVertex(configcontext, specs, false); // P-vertex is always needed // configure dpl timer to inject correct firstTForbit: start from the 1st orbit of TF containing 1st sampled orbit o2::raw::HBFUtilsInitializer hbfIni(configcontext, specs); diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-residual-aggregator.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-residual-aggregator.cxx index 20e37c3bcc3b4..d3593c1c50744 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-residual-aggregator.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-residual-aggregator.cxx @@ -32,6 +32,7 @@ void customize(std::vector& workflowOptions) {"output-type", VariantType::String, "unbinnedResid,trackParams", {"Comma separated list of outputs (without spaces). Valid strings: unbinnedResid, binnedResid, trackParams"}}, {"enable-track-input", VariantType::Bool, false, {"Whether to expect track data from interpolation workflow"}}, {"enable-ctp", VariantType::Bool, false, {"Subscribe to lumi info from CTP"}}, + {"enable-mc", VariantType::Bool, false, {"Whether to expect the MC truth of the track data from interpolation workflow (requires enable-track-input)"}}, {"disable-root-input", VariantType::Bool, false, {"disable root-files input readers"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}}; o2::raw::HBFUtilsInitializer::addConfigOption(options); @@ -47,6 +48,14 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); auto trkInput = configcontext.options().get("enable-track-input"); auto ctpInput = configcontext.options().get("enable-ctp"); + auto mcInput = configcontext.options().get("enable-mc"); + if (mcInput && !trkInput) { + LOG(error) << "MC truth input requires the track input (enable-track-input), will be ignored"; + mcInput = false; + } + if (mcInput && !configcontext.options().get("disable-root-input")) { + LOG(fatal) << "MC truth input is only supported directly from the interpolation workflow (disable-root-input)"; + } bool writeUnbinnedResiduals = false; bool writeBinnedResiduals = false; @@ -78,7 +87,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) if (!configcontext.options().get("disable-root-input")) { specs.emplace_back(o2::tpc::getUnbinnedTPCResidualsReaderSpec(trkInput)); } - specs.emplace_back(getTPCResidualAggregatorSpec(trkInput, ctpInput, writeUnbinnedResiduals, writeBinnedResiduals, writeTrackData)); + specs.emplace_back(getTPCResidualAggregatorSpec(trkInput, ctpInput, writeUnbinnedResiduals, writeBinnedResiduals, writeTrackData, mcInput)); // CTP input if (ctpInput) { diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/ResidualAggregator.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/ResidualAggregator.h index 00af697da3a9b..b2c5814a45f2e 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/ResidualAggregator.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/ResidualAggregator.h @@ -45,11 +45,11 @@ struct ResidualsContainer { ResidualsContainer& operator=(const ResidualsContainer& src) = delete; ~ResidualsContainer(); - void init(const TrackResiduals* residualsEngine, std::string outputDir, bool wFile, bool wBinnedResid, bool wUnbinnedResid, bool wTrackData, int autosave, int compression, long orbitResetTime); + void init(const TrackResiduals* residualsEngine, std::string outputDir, bool wFile, bool wBinnedResid, bool wUnbinnedResid, bool wTrackData, bool wTrackDataMC, int autosave, int compression, long orbitResetTime); void fillStatisticsBranches(); uint64_t getNEntries() const { return nResidualsTotal; } - void fill(const o2::dataformats::TFIDInfo& ti, const gsl::span resid, const gsl::span detInfoRes, const gsl::span trkRefsIn, const gsl::span* trkDataIn, const o2::ctp::LumiInfo* lumiInput); + void fill(const o2::dataformats::TFIDInfo& ti, const gsl::span resid, const gsl::span detInfoRes, const gsl::span trkRefsIn, const gsl::span* trkDataIn, const gsl::span* trkDataMCIn, const o2::ctp::LumiInfo* lumiInput); void merge(ResidualsContainer* prev); void print(); void writeToFile(bool closeFileAfterwards); @@ -66,6 +66,7 @@ struct ResidualsContainer { std::vector unbinnedRes, *unbinnedResPtr{&unbinnedRes}; ///< unbinned residuals which are sent to the aggregator std::vector detInfoUnbRes, *detInfoUnbResPtr{&detInfoUnbRes}; ///< detector info associated to unbinned residuals which are sent to the aggregator std::vector trkData, *trkDataPtr{&trkData}; ///< track data and cluster ranges + std::vector trkDataMC, *trkDataMCPtr{&trkDataMC}; ///< MC truth for the track data (MC only) std::vector trackInfo, *trackInfoPtr{&trackInfo}; ///< allows to obtain track type for each unbinned residual downstream o2::ctp::LumiInfo lumiTF; ///< for each processed TF we store the lumi information in the tree of unbinned residuals uint64_t timeMS; ///< for each processed TF we store its absolute time in ms in the tree of unbinned residuals @@ -83,6 +84,7 @@ struct ResidualsContainer { bool writeBinnedResid{false}; ///< flag, whether binned residuals should be written out bool writeUnbinnedResiduals{false}; ///< flag, whether unbinned residuals should be written out bool writeTrackData{false}; ///< flag, whether full seeding track information should be written out + bool writeTrackDataMC{false}; ///< flag, whether the MC truth of the seeding tracks should be written out int autosaveInterval{0}; ///< if > 0, then the output written to file for every n-th TF // additional info @@ -94,7 +96,7 @@ struct ResidualsContainer { float TPCVDriftRef{-1.}; ///< TPC nominal drift speed in cm/microseconds float TPCDriftTimeOffsetRef{0.}; ///< TPC nominal (e.g. at the start of run) drift time bias in cm/mus - ClassDefNV(ResidualsContainer, 5); + ClassDefNV(ResidualsContainer, 6); }; class ResidualAggregator final : public o2::calibration::TimeSlotCalibration @@ -122,6 +124,7 @@ class ResidualAggregator final : public o2::calibration::TimeSlotCalibration0 then the output is written to a file for every n-th TF int mCompressionSetting{101}; ///< single integer defining the ROOT compression algorithm and level (see TFile doc for details) size_t mMinEntries; ///< the minimum number of residuals required for the map creation (per voxel) diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h index 93fd816bfe26e..efd9e8434e235 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h @@ -24,6 +24,7 @@ #include "ReconstructionDataFormats/TrackTPCITS.h" #include "ReconstructionDataFormats/MatchInfoTOF.h" #include "ReconstructionDataFormats/GlobalTrackID.h" +#include "SimulationDataFormat/MCCompLabel.h" #include "DataFormatsITSMFT/Cluster.h" #include "DataFormatsITSMFT/TrkClusRef.h" #include "DataFormatsITSMFT/TopologyDictionary.h" @@ -244,6 +245,26 @@ struct TrackData { ClassDefNV(TrackData, 12); }; +/// MC truth for a TrackData entry (stored only for MC, aligned 1:1 with the TrackData vector) +struct TrackDataMC { + enum Flags : uint8_t { HasITSOut = 0x1, ///< parITSOut is filled + HasTPCIn = 0x2, ///< parTPCIn is filled + FakeITSTPC = 0x4 }; ///< ITS and TPC parts of the track have different MC labels + o2::MCCompLabel label{}; ///< MC label of the ITS-TPC part of the seeding track + o2::MCCompLabel labelITS{}; ///< MC label of its ITS part + o2::MCCompLabel labelTPC{}; ///< MC label of its TPC part + o2::track::TrackPar parITSOut{}; ///< truth at x and alpha of TrackData::par, from the nearest ITS track reference (propagated with the material correction, true mass) + o2::track::TrackPar parTPCIn{}; ///< truth at the first TPC track reference (sector frame) + float distITSRef{-1.f}; ///< 3D distance between the ITS track reference used and TrackData::par in cm + float distTPCRef{-1.f}; ///< distance (y,z) between parTPCIn propagated to the innermost TPC cluster of the track and that cluster in cm (large: wrong leg, looper, fake) + int pdg{0}; ///< PDG code of the MC particle of the ITS part (TPC part if no ITS label) + uint8_t flags{0}; + bool hasITSOut() const { return flags & HasITSOut; } + bool hasTPCIn() const { return flags & HasTPCIn; } + bool isFakeITSTPC() const { return flags & FakeITSTPC; } + ClassDefNV(TrackDataMC, 1); +}; + /// \class TrackInterpolation /// This class is retrieving the TPC space point residuals by interpolating ITS/TRD/TOF tracks. /// The residuals are stored in the specified vectors of TPCClusterResiduals diff --git a/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx b/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx index 7844d219deea5..e8423957990c7 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx +++ b/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx @@ -77,6 +77,7 @@ ResidualsContainer::ResidualsContainer(ResidualsContainer&& rhs) unbinnedRes = std::move(rhs.unbinnedRes); trackInfo = std::move(rhs.trackInfo); trkData = std::move(rhs.trkData); + trkDataMC = std::move(rhs.trkDataMC); orbitReset = rhs.orbitReset; firstTForbit = rhs.firstTForbit; firstSeenTF = rhs.firstSeenTF; @@ -84,13 +85,14 @@ ResidualsContainer::ResidualsContainer(ResidualsContainer&& rhs) nResidualsTotal = rhs.nResidualsTotal; } -void ResidualsContainer::init(const TrackResiduals* residualsEngine, std::string outputDir, bool wFile, bool wBinnedResid, bool wUnbinnedResid, bool wTrackData, int autosave, int compression, long orbitResetTime) +void ResidualsContainer::init(const TrackResiduals* residualsEngine, std::string outputDir, bool wFile, bool wBinnedResid, bool wUnbinnedResid, bool wTrackData, bool wTrackDataMC, int autosave, int compression, long orbitResetTime) { trackResiduals = residualsEngine; writeToRootFile = wFile; writeBinnedResid = wBinnedResid; writeUnbinnedResiduals = wUnbinnedResid; writeTrackData = wTrackData; + writeTrackDataMC = wTrackData && wTrackDataMC; autosaveInterval = autosave; orbitReset = orbitResetTime; if (writeToRootFile) { @@ -129,6 +131,9 @@ void ResidualsContainer::init(const TrackResiduals* residualsEngine, std::string if (writeTrackData) { treeOutTrackData = std::make_unique("trackData", "Track information incl cluster range ref"); treeOutTrackData->Branch("trk", &trkDataPtr); + if (writeTrackDataMC) { + treeOutTrackData->Branch("trkMC", &trkDataMCPtr); + } } if (writeBinnedResid) { treeOutResiduals = std::make_unique("resid", "TPC binned residuals"); @@ -171,7 +176,7 @@ void ResidualsContainer::fillStatisticsBranches() } } -void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::span resid, const gsl::span detInfoRes, const gsl::span trkRefsIn, const gsl::span* trkDataIn, const o2::ctp::LumiInfo* lumiInput) +void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::span resid, const gsl::span detInfoRes, const gsl::span trkRefsIn, const gsl::span* trkDataIn, const gsl::span* trkDataMCIn, const o2::ctp::LumiInfo* lumiInput) { // receives large vector of unbinned residuals and fills the sector-wise vectors // with binned residuals and statistics @@ -244,8 +249,12 @@ void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::sp for (const auto& trkIn : *trkDataIn) { trkData.push_back(trkIn); } + if (writeTrackDataMC && trkDataMCIn) { + trkDataMC.assign(trkDataMCIn->begin(), trkDataMCIn->end()); + } treeOutTrackData->Fill(); trkData.clear(); + trkDataMC.clear(); } if (writeUnbinnedResiduals) { if (lumiInput) { @@ -338,6 +347,9 @@ void ResidualsContainer::merge(ResidualsContainer* prev) if (writeTrackData) { prev->treeOutTrackData->SetBranchAddress("trk", &trkDataPtr); + if (writeTrackDataMC) { + prev->treeOutTrackData->SetBranchAddress("trkMC", &trkDataMCPtr); + } for (int i = 0; i < treeOutTrackData->GetEntries(); ++i) { treeOutTrackData->GetEntry(i); prev->treeOutTrackData->Fill(); @@ -456,7 +468,7 @@ Slot& ResidualAggregator::emplaceNewSlot(bool front, TFType tStart, TFType tEnd) auto& cont = getSlots(); auto& slot = front ? cont.emplace_front(tStart, tEnd) : cont.emplace_back(tStart, tEnd); slot.setContainer(std::make_unique()); - slot.getContainer()->init(&mTrackResiduals, mOutputDir, mWriteOutput, mWriteBinnedResiduals, mWriteUnbinnedResiduals, mWriteTrackData, mAutosaveInterval, mCompressionSetting, mOrbitResetTime); + slot.getContainer()->init(&mTrackResiduals, mOutputDir, mWriteOutput, mWriteBinnedResiduals, mWriteUnbinnedResiduals, mWriteTrackData, mWriteTrackDataMC, mAutosaveInterval, mCompressionSetting, mOrbitResetTime); std::chrono::duration emplaceDuration = std::chrono::high_resolution_clock::now() - emplaceStartTime; LOGP(info, "Emplacing new calibration slot took: {} ms", std::chrono::duration_cast(emplaceDuration).count()); return slot; diff --git a/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h b/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h index e77610acb8e7e..b6db3abfecc90 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h +++ b/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h @@ -21,6 +21,8 @@ #pragma link C++ class std::vector < o2::tpc::TrackDataCompact> + ; #pragma link C++ class o2::tpc::TrackData + ; #pragma link C++ class std::vector < o2::tpc::TrackData> + ; +#pragma link C++ class o2::tpc::TrackDataMC + ; +#pragma link C++ class std::vector < o2::tpc::TrackDataMC> + ; #pragma link C++ class o2::tpc::TrackDataExtended + ; #pragma link C++ class std::vector < o2::tpc::TrackDataExtended> + ; #pragma link C++ class o2::tpc::TPCClusterResiduals + ; From e869868e6925a93cf4544cc96bd72f0facadf7c2 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Sat, 3 Oct 2026 14:41:21 +0200 Subject: [PATCH 3/6] TPC SCD: MC truth at the TRD (entrance and ideal tracklet positions) For tracks with TRD residuals (TRD seeds and ITS-TPC map seeds whose most complete track has TRD), TrackDataMC now also holds the truth at the TRD entrance (first TRD track reference in time, sector frame) and, per layer, the true y and z at the x of the tracklet in the tracklet's sector frame, i.e. the frame of the stored TRD residual (nearest TRD track reference propagated with the workflow's material correction), for the true particle of the ITS-TPC part (ideal tracklets). TrackInterpolation records which ITS-TPC-TRD track gave the TRD residuals of each stored track (getTRDGIDsSuccess, aligned with the track data). Co-Authored-By: Claude Opus 5.5 --- .../src/TPCInterpolationSpec.cxx | 72 ++++++++++++++++--- .../include/SpacePoints/TrackInterpolation.h | 13 +++- .../SpacePoints/src/TrackInterpolation.cxx | 7 ++ 3 files changed, 81 insertions(+), 11 deletions(-) diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx index c95bf421b7dba..40250254b5587 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx @@ -177,15 +177,16 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) { // MC truth for every stored TrackData: labels of the ITS-TPC part of the seed and of its ITS and TPC parts, the truth at // the ITS outer parameters (the ITS track reference nearest to TrackData::par, propagated to its x with the material - // correction of the workflow and the mass of the true particle) and the truth at the TPC entrance (the first TPC track - // reference in time, in the sector frame) + // correction of the workflow and the mass of the true particle), the truth at the TPC entrance (the first TPC track + // reference in time, in the sector frame) and, for TRD-matched seeds, the truth at the TRD entrance and the true + // positions at the x of the TRD tracklets of the track (ideal tracklets) const auto& trkData = mInterpolation.getReferenceTracks(); mTrackDataMC.clear(); mTrackDataMC.resize(trkData.size()); struct Lookup { - o2::MCCompLabel lbl; - uint32_t idx; - bool its; // ITS outer (true) or TPC entrance (false) + o2::MCCompLabel lbl{}; + uint32_t idx{0}; + uint8_t kind{0}; // 0: ITS outer, 1: TPC entrance, 2: TRD }; std::vector lookups; lookups.reserve(2 * trkData.size()); @@ -208,10 +209,13 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) const auto& lblITS = mc.labelITS.isValid() ? mc.labelITS : mc.label; const auto& lblTPC = mc.labelTPC.isValid() ? mc.labelTPC : mc.label; if (lblITS.isValid()) { - lookups.push_back({lblITS, uint32_t(i), true}); + lookups.push_back({lblITS, uint32_t(i), 0}); } if (lblTPC.isValid()) { - lookups.push_back({lblTPC, uint32_t(i), false}); + lookups.push_back({lblTPC, uint32_t(i), 1}); + } + if (mInterpolation.getTRDGIDsSuccess()[i].isIndexSet() && lblTPC.isValid()) { // track with TRD residuals: the true particle of the ITS-TPC part at the TRD + lookups.push_back({mc.label.isValid() ? mc.label : lblTPC, uint32_t(i), 2}); } } // the reader loads the kinematics of a whole event (can be >100 MB): process event by event and release it right after @@ -265,7 +269,7 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) const auto* pPDG = mcTrk ? O2DatabasePDG::Instance()->GetParticle(pdg) : nullptr; int charge = pPDG ? int(std::lround(pPDG->Charge() / 3.)) : 0; // TParticlePDG charge is in units of |e|/3 auto refs = mMCReader->getTrackRefs(lbl.getSourceID(), lbl.getEventID(), lbl.getTrackID()); - if (lk.its) { // ITS outer: track reference of the ITS part nearest to TrackData::par + if (lk.kind == 0) { // ITS outer: track reference of the ITS part nearest to TrackData::par mc.pdg = pdg; const o2::TrackReference* best = nullptr; float bestD2 = 1e30f; @@ -291,7 +295,7 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) mc.flags |= TrackDataMC::HasITSOut; } } - } else { // TPC entrance: first TPC track reference in time of the TPC part + } else if (lk.kind == 1) { // TPC entrance: first TPC track reference in time of the TPC part if (!mc.labelITS.isValid() && !mc.label.isValid()) { mc.pdg = pdg; // no ITS lookup for this track } @@ -321,6 +325,56 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) } } } + } else if (charge) { // TRD: entrance (first TRD track reference in time) and the true positions at the tracklet x of each layer + const o2::TrackReference* first = nullptr; + for (const auto& ref : refs) { + if (ref.getDetectorId() == DetID::TRD && (!first || ref.getTime() < first->getTime())) { + first = &ref; + } + } + if (!first) { + continue; + } + mc.parTRDIn = refToPar(*first, charge, pdg, true); + mc.flags |= TrackDataMC::HasTRDIn; + const auto& trkTRD = recoData.getITSTPCTRDTrack(mInterpolation.getTRDGIDsSuccess()[lk.idx]); // the TRD track of the stored TRD residuals + const auto trkltsCalib = recoData.getTRDCalibratedTracklets(); + const auto tracklets = recoData.getTRDTracklets(); + for (int iLayer = 0; iLayer < o2::trd::constants::NLAYER; iLayer++) { + int trkltIdx = trkTRD.getTrackletIndex(iLayer); + if (trkltIdx < 0) { + continue; + } + const auto& sp = trkltsCalib[trkltIdx]; // tracklet x, y, z in its sector frame, as used for the TRD residual + int sec = tracklets[trkltIdx].getDetector() / (o2::trd::constants::NLAYER * o2::trd::constants::NSTACK); + float alpha = o2::math_utils::sector2Angle(sec); + float cs = std::cos(alpha); + float sn = std::sin(alpha); + float gx = sp.getX() * cs - sp.getY() * sn; + float gy = sp.getX() * sn + sp.getY() * cs; + float gz = sp.getZ(); + const o2::TrackReference* best = nullptr; + float bestD2 = 1e30f; + for (const auto& ref : refs) { + if (ref.getDetectorId() != DetID::TRD) { + continue; + } + float dx = ref.X() - gx; + float dy = ref.Y() - gy; + float dz = ref.Z() - gz; + float d2 = dx * dx + dy * dy + dz * dz; + if (d2 < bestD2) { + bestD2 = d2; + best = &ref; + } + } + auto par = refToPar(*best, charge, pdg, false); + if (par.rotateParam(alpha) && prop->PropagateToXBxByBz(par, sp.getX(), 0.999f, o2::base::Propagator::MAX_STEP, matCorr)) { + mc.yTRD[iLayer] = par.getY(); + mc.zTRD[iLayer] = par.getZ(); + mc.trdLayerMask |= uint8_t(1) << iLayer; + } + } } } if (curSrc >= 0) { diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h index efd9e8434e235..bb26e5fd986ba 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h @@ -249,7 +249,8 @@ struct TrackData { struct TrackDataMC { enum Flags : uint8_t { HasITSOut = 0x1, ///< parITSOut is filled HasTPCIn = 0x2, ///< parTPCIn is filled - FakeITSTPC = 0x4 }; ///< ITS and TPC parts of the track have different MC labels + FakeITSTPC = 0x4, ///< ITS and TPC parts of the track have different MC labels + HasTRDIn = 0x8 }; ///< parTRDIn is filled o2::MCCompLabel label{}; ///< MC label of the ITS-TPC part of the seeding track o2::MCCompLabel labelITS{}; ///< MC label of its ITS part o2::MCCompLabel labelTPC{}; ///< MC label of its TPC part @@ -257,12 +258,17 @@ struct TrackDataMC { o2::track::TrackPar parTPCIn{}; ///< truth at the first TPC track reference (sector frame) float distITSRef{-1.f}; ///< 3D distance between the ITS track reference used and TrackData::par in cm float distTPCRef{-1.f}; ///< distance (y,z) between parTPCIn propagated to the innermost TPC cluster of the track and that cluster in cm (large: wrong leg, looper, fake) + o2::track::TrackPar parTRDIn{}; ///< truth at the first TRD track reference (sector frame), TRD-matched seeds only + float yTRD[6] = {}; ///< truth y at the x of the TRD tracklet of each layer (tracklet sector frame), see trdLayerMask + float zTRD[6] = {}; ///< truth z at the x of the TRD tracklet of each layer (tracklet sector frame), see trdLayerMask + uint8_t trdLayerMask{0}; ///< bit i set: yTRD[i], zTRD[i] filled int pdg{0}; ///< PDG code of the MC particle of the ITS part (TPC part if no ITS label) uint8_t flags{0}; bool hasITSOut() const { return flags & HasITSOut; } bool hasTPCIn() const { return flags & HasTPCIn; } bool isFakeITSTPC() const { return flags & FakeITSTPC; } - ClassDefNV(TrackDataMC, 1); + bool hasTRDIn() const { return flags & HasTRDIn; } + ClassDefNV(TrackDataMC, 2); }; /// \class TrackInterpolation @@ -445,6 +451,8 @@ class TrackInterpolation std::vector& getTrackDataCompact() { return mTrackDataCompact; } std::vector& getTrackDataExtended() { return mTrackDataExtended; } std::vector& getReferenceTracks() { return mTrackData; } + /// ITS-TPC-TRD track whose tracklets gave the TRD residuals of each stored track (not set if none), aligned with getReferenceTracks() + const std::vector& getTRDGIDsSuccess() const { return mTRDGIDsSuccess; } void setLane(int lID, int nL) { @@ -512,6 +520,7 @@ class TrackInterpolation // cache std::array mCache{{}}; ///< caching positions, covariances and angles for track extrapolations and interpolation std::vector mGIDsSuccess; ///< keep track of the GIDs which could be processed successfully + std::vector mTRDGIDsSuccess; ///< ITS-TPC-TRD track used for the TRD residuals of each stored track (not set if none) TrackValidationData mTrackValidation; diff --git a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx index fde816ae40892..17237473e85ae 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx +++ b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx @@ -807,10 +807,12 @@ void TrackInterpolation::interpolateTrack(int iSeed) } bool stopPropagation = !mExtDetResid; + GTrackID gidTRDUsed{}; if (!stopPropagation) { // do we have TRD residuals to add? trkWork = trkOuter; if (!allLost && gidTable[GTrackID::TRD].isIndexSet()) { // allLost: trkOuter is not a valid outer param + gidTRDUsed = gidTable[GTrackID::ITSTPCTRD]; const auto& trkTRD = mRecoCont->getITSTPCTRDTrack(gidTable[GTrackID::ITSTPCTRD]); for (int iLayer = 0; iLayer < o2::trd::constants::NLAYER; iLayer++) { std::array trkltTRDYZ{}; @@ -926,6 +928,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) } mGIDsSuccess.push_back(mGIDs[iSeed]); + mTRDGIDsSuccess.push_back(gidTRDUsed); mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), trackData.multStack, nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid, trackData.filterFlag); mTrackData.push_back(std::move(trackData)); stored = true; @@ -1187,12 +1190,14 @@ void TrackInterpolation::extrapolateTrack(int iSeed) } bool stopPropagation = !mExtDetResid; + GTrackID gidTRDUsed{}; if (!stopPropagation) { // do we have TRD residuals to add? int iSeedFull = mParentID[iSeed] == -1 ? iSeed : mParentID[iSeed]; auto gidFull = mGIDs[iSeedFull]; const auto& gidTableFull = mGIDtables[iSeedFull]; if (!refLost && gidTableFull[GTrackID::TRD].isIndexSet()) { // refLost: trkWork did not reach the TPC outer end + gidTRDUsed = gidTableFull[GTrackID::ITSTPCTRD]; const auto& trkTRD = mRecoCont->getITSTPCTRDTrack(gidTableFull[GTrackID::ITSTPCTRD]); trackData.nTrkltsTRD = trkTRD.getNtracklets(); trackData.chi2TRD = trkTRD.getChi2(); @@ -1319,6 +1324,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed) mTrackData.push_back(std::move(trackData)); stored = true; mGIDsSuccess.push_back(mGIDs[iSeed]); + mTRDGIDsSuccess.push_back(gidTRDUsed); mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), trackData.multStack, nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid, trackData.filterFlag); if (mDumpTrackPoints) { (*trackDataExtended).clIdx.setEntries(nClValidated); @@ -1694,6 +1700,7 @@ void TrackInterpolation::reset() mClRes.clear(); mDetInfoRes.clear(); mGIDsSuccess.clear(); + mTRDGIDsSuccess.clear(); for (auto& vec : mTrackIndices) { vec.clear(); } From d7451823c2124b0479de369b4289b082a6085d76 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Sat, 3 Oct 2026 15:13:42 +0200 Subject: [PATCH 4/6] TPC SCD: MC origin of the residual tracks (mother, production vertex, sister tracks) For the particle of the ITS-TPC part TrackDataMC now holds the mother label and PDG code (for primaries the generator-level parent), the production vertex, momentum and process, a primary flag, and sisterIdx: the index of another stored track with the same mother (cycling through all of them if more than two), e.g. to select both daughters of a K0s. pdg is now the PDG code of the particle of the ITS-TPC part, consistent with the origin and TRD fields. Co-Authored-By: Claude Opus 5.5 --- .../src/TPCInterpolationSpec.cxx | 47 ++++++++++++++++--- .../include/SpacePoints/TrackInterpolation.h | 18 +++++-- 2 files changed, 55 insertions(+), 10 deletions(-) diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx index 40250254b5587..b7e514bf76516 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx @@ -178,15 +178,16 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) // MC truth for every stored TrackData: labels of the ITS-TPC part of the seed and of its ITS and TPC parts, the truth at // the ITS outer parameters (the ITS track reference nearest to TrackData::par, propagated to its x with the material // correction of the workflow and the mass of the true particle), the truth at the TPC entrance (the first TPC track - // reference in time, in the sector frame) and, for TRD-matched seeds, the truth at the TRD entrance and the true - // positions at the x of the TRD tracklets of the track (ideal tracklets) + // reference in time, in the sector frame), for tracks with TRD residuals the truth at the TRD entrance and the true + // positions at the x of the TRD tracklets of the track (ideal tracklets), and the origin of the particle of the ITS-TPC + // part (mother, production vertex and process, other stored tracks with the same mother) const auto& trkData = mInterpolation.getReferenceTracks(); mTrackDataMC.clear(); mTrackDataMC.resize(trkData.size()); struct Lookup { o2::MCCompLabel lbl{}; uint32_t idx{0}; - uint8_t kind{0}; // 0: ITS outer, 1: TPC entrance, 2: TRD + uint8_t kind{0}; // 0: ITS outer, 1: TPC entrance, 2: TRD, 3: origin }; std::vector lookups; lookups.reserve(2 * trkData.size()); @@ -217,6 +218,9 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) if (mInterpolation.getTRDGIDsSuccess()[i].isIndexSet() && lblTPC.isValid()) { // track with TRD residuals: the true particle of the ITS-TPC part at the TRD lookups.push_back({mc.label.isValid() ? mc.label : lblTPC, uint32_t(i), 2}); } + if (lblTPC.isValid()) { // origin of the true particle of the ITS-TPC part + lookups.push_back({mc.label.isValid() ? mc.label : lblTPC, uint32_t(i), 3}); + } } // the reader loads the kinematics of a whole event (can be >100 MB): process event by event and release it right after std::sort(lookups.begin(), lookups.end(), [](const Lookup& a, const Lookup& b) { @@ -270,7 +274,6 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) int charge = pPDG ? int(std::lround(pPDG->Charge() / 3.)) : 0; // TParticlePDG charge is in units of |e|/3 auto refs = mMCReader->getTrackRefs(lbl.getSourceID(), lbl.getEventID(), lbl.getTrackID()); if (lk.kind == 0) { // ITS outer: track reference of the ITS part nearest to TrackData::par - mc.pdg = pdg; const o2::TrackReference* best = nullptr; float bestD2 = 1e30f; auto xyzReco = trk.par.getXYZGlo(); @@ -296,9 +299,6 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) } } } else if (lk.kind == 1) { // TPC entrance: first TPC track reference in time of the TPC part - if (!mc.labelITS.isValid() && !mc.label.isValid()) { - mc.pdg = pdg; // no ITS lookup for this track - } const o2::TrackReference* first = nullptr; for (const auto& ref : refs) { if (ref.getDetectorId() == DetID::TPC && (!first || ref.getTime() < first->getTime())) { @@ -325,6 +325,27 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) } } } + } else if (lk.kind == 3) { // origin: mother, production vertex and process + if (!mcTrk) { + continue; + } + mc.pdg = pdg; + if (mcTrk->isPrimary()) { + mc.flags |= TrackDataMC::IsPrimary; + } + mc.process = mcTrk->getProcess(); + mc.prodX = mcTrk->GetStartVertexCoordinatesX(); + mc.prodY = mcTrk->GetStartVertexCoordinatesY(); + mc.prodZ = mcTrk->GetStartVertexCoordinatesZ(); + mc.prodPx = mcTrk->GetStartVertexMomentumX(); + mc.prodPy = mcTrk->GetStartVertexMomentumY(); + mc.prodPz = mcTrk->GetStartVertexMomentumZ(); + int motherId = mcTrk->getMotherTrackId(); + if (motherId >= 0) { + mc.motherLabel = o2::MCCompLabel(motherId, lbl.getEventID(), lbl.getSourceID()); + const auto* mother = mMCReader->getTrack(lbl.getSourceID(), lbl.getEventID(), motherId); + mc.motherPdg = mother ? mother->GetPdgCode() : 0; + } } else if (charge) { // TRD: entrance (first TRD track reference in time) and the true positions at the tracklet x of each layer const o2::TrackReference* first = nullptr; for (const auto& ref : refs) { @@ -380,6 +401,18 @@ void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) if (curSrc >= 0) { mMCReader->releaseTracksForSourceAndEvent(curSrc, curEv); } + // stored tracks with the same mother (e.g. both daughters of a K0s): each points to the next one, cyclically + std::unordered_map> daughters; + for (uint32_t i = 0; i < mTrackDataMC.size(); ++i) { + if (mTrackDataMC[i].motherLabel.isSet()) { + daughters[mTrackDataMC[i].motherLabel.getTrackEventSourceID()].push_back(i); + } + } + for (const auto& [mother, idx] : daughters) { + for (size_t k = 0; idx.size() > 1 && k < idx.size(); ++k) { + mTrackDataMC[idx[k]].sisterIdx = idx[(k + 1) % idx.size()]; + } + } } void TPCInterpolationDPL::endOfStream(EndOfStreamContext& ec) diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h index bb26e5fd986ba..38b7fd3784462 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h @@ -250,7 +250,8 @@ struct TrackDataMC { enum Flags : uint8_t { HasITSOut = 0x1, ///< parITSOut is filled HasTPCIn = 0x2, ///< parTPCIn is filled FakeITSTPC = 0x4, ///< ITS and TPC parts of the track have different MC labels - HasTRDIn = 0x8 }; ///< parTRDIn is filled + HasTRDIn = 0x8, ///< parTRDIn is filled + IsPrimary = 0x10 }; ///< the particle of the ITS-TPC part is a primary (MCTrack::isPrimary) o2::MCCompLabel label{}; ///< MC label of the ITS-TPC part of the seeding track o2::MCCompLabel labelITS{}; ///< MC label of its ITS part o2::MCCompLabel labelTPC{}; ///< MC label of its TPC part @@ -262,13 +263,24 @@ struct TrackDataMC { float yTRD[6] = {}; ///< truth y at the x of the TRD tracklet of each layer (tracklet sector frame), see trdLayerMask float zTRD[6] = {}; ///< truth z at the x of the TRD tracklet of each layer (tracklet sector frame), see trdLayerMask uint8_t trdLayerMask{0}; ///< bit i set: yTRD[i], zTRD[i] filled - int pdg{0}; ///< PDG code of the MC particle of the ITS part (TPC part if no ITS label) + int pdg{0}; ///< PDG code of the particle of the ITS-TPC part (as for the origin and TRD fields) + o2::MCCompLabel motherLabel{}; ///< MC label of the mother of the particle of the ITS-TPC part (for primaries the generator-level parent) + int motherPdg{0}; ///< PDG code of that mother (0: none) + float prodX{0.f}; ///< production vertex x of the particle of the ITS-TPC part (global, cm) + float prodY{0.f}; ///< production vertex y (global, cm) + float prodZ{0.f}; ///< production vertex z (global, cm) + float prodPx{0.f}; ///< momentum at production x (global, GeV/c), e.g. to compare a track propagated to the vertex + float prodPy{0.f}; ///< momentum at production y (global, GeV/c) + float prodPz{0.f}; ///< momentum at production z (global, GeV/c) + int sisterIdx{-1}; ///< index in the TrackData vector of another stored track with the same mother (cycling through all of them if more than two), -1: none + uint8_t process{0}; ///< production process of the particle of the ITS-TPC part (TMCProcess) uint8_t flags{0}; bool hasITSOut() const { return flags & HasITSOut; } bool hasTPCIn() const { return flags & HasTPCIn; } bool isFakeITSTPC() const { return flags & FakeITSTPC; } bool hasTRDIn() const { return flags & HasTRDIn; } - ClassDefNV(TrackDataMC, 2); + bool isPrimary() const { return flags & IsPrimary; } + ClassDefNV(TrackDataMC, 3); }; /// \class TrackInterpolation From 333cdd578fb665217dd1a61066a6f7c4c4a7bb10 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Sat, 3 Oct 2026 10:43:29 +0200 Subject: [PATCH 5/6] MCKinematicsReader: keep the loaded tracks instead of copying them, free the baskets loadTracksForSourceAndEvent deep-copied the vector which ROOT allocated for us (a pointer to nullptr was passed, so it is ours) and deleted the original: +1x the event at peak (373 MB for a PbPb event with 6.5M MCTracks). The vector is now stored directly and the branch address reset. In addition, the decompressed baskets (~1x the event for the split MCTrack branch) are dropped once the event is the last entry of its TTree cluster, i.e. when no later entry can reuse them, which keeps small events sharing baskets as fast as before (QED, 10k events: 0.17 s either way, 0.70 s when dropping after every event). For the PbPb event: peak memory in use 1628 -> 1232 MB, held after releaseTracksForSourceAndEvent 882 -> 443 MB. Co-Authored-By: Claude Opus 5.5 --- Steer/src/MCKinematicsReader.cxx | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/Steer/src/MCKinematicsReader.cxx b/Steer/src/MCKinematicsReader.cxx index 21024dba78368..b3b4182fc24f0 100644 --- a/Steer/src/MCKinematicsReader.cxx +++ b/Steer/src/MCKinematicsReader.cxx @@ -92,9 +92,15 @@ void MCKinematicsReader::loadTracksForSourceAndEvent(int source, int event) cons std::vector* loadtracks = nullptr; br->SetAddress(&loadtracks); br->GetEntry(event); - mTracks[source][event] = new std::vector; - *mTracks[source][event] = *loadtracks; - delete loadtracks; + // ROOT allocated the vector for us and we own it (we passed a pointer to nullptr): keep it instead of copying it + mTracks[source][event] = loadtracks; + br->ResetAddress(); // the branch must not refer to the stored vector (nor to the local pointer) any more + // free the decompressed baskets (~ the size of the event) if no later entry reads them, i.e. at the end of its cluster + auto clusterIt = br->GetTree()->GetClusterIterator(event); + clusterIt.Next(); + if (event + 1 >= clusterIt.GetNextEntry()) { + br->DropBaskets("all"); + } } } } From 13be6aa43baf98d651784050f98fd12cbdd5e62e Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Sat, 3 Oct 2026 20:00:40 +0200 Subject: [PATCH 6/6] MCKinematicsReader: load and release the track references per event The track references of a source were loaded for all events at the first access and never released. They are now loaded per event on demand (with the baskets dropped at the end of a TTree cluster, as for the tracks) and releaseTracksForSourceAndEvent frees them together with the tracks. For a PbPb TF (10 signal events, 1.75M references): +304 MB held from the first access until the end before, now +135 MB at peak and +55 MB after releasing all events. Co-Authored-By: Claude Opus 5.5 --- Steer/include/Steer/MCKinematicsReader.h | 16 ++++++--- Steer/src/MCKinematicsReader.cxx | 43 +++++++++++++++++------- 2 files changed, 43 insertions(+), 16 deletions(-) diff --git a/Steer/include/Steer/MCKinematicsReader.h b/Steer/include/Steer/MCKinematicsReader.h index ae5ccf6615c56..793711c61de87 100644 --- a/Steer/include/Steer/MCKinematicsReader.h +++ b/Steer/include/Steer/MCKinematicsReader.h @@ -87,7 +87,7 @@ class MCKinematicsReader /// variant returning all tracks for source and event at once std::vector const& getTracks(int source, int event) const; - /// API to ask releasing tracks (freeing memory) for source + event + /// API to ask releasing tracks and track references (freeing memory) for source + event void releaseTracksForSourceAndEvent(int source, int event); /// variant returning all tracks for an event id (source = 0) at once @@ -129,7 +129,8 @@ class MCKinematicsReader void initTracksForSource(int source) const; void loadTracksForSourceAndEvent(int source, int eventID) const; void loadHeadersForSource(int source) const; - void loadTrackRefsForSource(int source) const; + void initTrackRefsForSource(int source) const; + void loadTrackRefsForSourceAndEvent(int source, int event) const; void initIndexedTrackRefs(std::vector& refs, o2::dataformats::MCTruthContainer& indexedrefs) const; DigitizationContext const* mDigitizationContext = nullptr; @@ -142,6 +143,7 @@ class MCKinematicsReader mutable std::vector*>> mTracks; // the in-memory track container mutable std::vector> mHeaders; // the in-memory header container mutable std::vector>> mIndexedTrackRefs; // the in-memory track ref container + mutable std::vector> mTrackRefsLoaded; // whether the track refs of a source/event are in memory bool mInitialized = false; // whether initialized }; @@ -206,11 +208,14 @@ inline gsl::span MCKinematicsReader::getTrackRefs(int source } auto& perEvent = mIndexedTrackRefs[source]; if (perEvent.size() == 0) { - loadTrackRefsForSource(source); + initTrackRefsForSource(source); } if (static_cast(event) >= perEvent.size()) { return {}; } + if (!mTrackRefsLoaded[source][event]) { + loadTrackRefsForSourceAndEvent(source, event); + } return perEvent[event].getLabels(track); } @@ -218,11 +223,14 @@ inline const std::vector& MCKinematicsReader::getTrackRefsBy { auto const& perEvent = mIndexedTrackRefs.at(source); if (perEvent.size() == 0) { - loadTrackRefsForSource(source); + initTrackRefsForSource(source); } if (static_cast(event) >= perEvent.size()) { reportMissingEvent("events of track references", source, event, perEvent.size()); } + if (!mTrackRefsLoaded[source][event]) { + loadTrackRefsForSourceAndEvent(source, event); + } return perEvent[event].getTruthArray(); } diff --git a/Steer/src/MCKinematicsReader.cxx b/Steer/src/MCKinematicsReader.cxx index b3b4182fc24f0..42cd40c90c1ae 100644 --- a/Steer/src/MCKinematicsReader.cxx +++ b/Steer/src/MCKinematicsReader.cxx @@ -111,6 +111,11 @@ void MCKinematicsReader::releaseTracksForSourceAndEvent(int source, int eventID) delete mTracks[source][eventID]; mTracks[source][eventID] = nullptr; } + // the track references of this event as well (reloaded on demand) + if (static_cast(eventID) < mTrackRefsLoaded.at(source).size() && mTrackRefsLoaded[source][eventID]) { + mIndexedTrackRefs[source][eventID] = o2::dataformats::MCTruthContainer(); + mTrackRefsLoaded[source][eventID] = false; + } } void MCKinematicsReader::loadHeadersForSource(int source) const @@ -135,31 +140,43 @@ void MCKinematicsReader::loadHeadersForSource(int source) const } } -void MCKinematicsReader::loadTrackRefsForSource(int source) const +void MCKinematicsReader::initTrackRefsForSource(int source) const { auto chain = mInputChains[source]; if (chain) { // todo: get name from NameConfig auto br = chain->GetBranch("TrackRefs"); if (br) { - std::vector* refs = nullptr; - br->SetAddress(&refs); mIndexedTrackRefs[source].resize(br->GetEntries()); - for (int event = 0; event < br->GetEntries(); ++event) { - br->GetEntry(event); - if (refs) { - // we convert the original flat vector into an indexed structure - initIndexedTrackRefs(*refs, mIndexedTrackRefs[source][event]); - delete refs; - refs = nullptr; - } - } + mTrackRefsLoaded[source].assign(br->GetEntries(), false); } else { LOG(warn) << "TrackRefs branch not found"; } } } +void MCKinematicsReader::loadTrackRefsForSourceAndEvent(int source, int event) const +{ + // todo: get name from NameConfig + auto br = mInputChains[source]->GetBranch("TrackRefs"); + std::vector* refs = nullptr; // allocated by ROOT, owned by us + br->SetAddress(&refs); + br->GetEntry(event); + if (refs) { + // we convert the original flat vector into an indexed structure + initIndexedTrackRefs(*refs, mIndexedTrackRefs[source][event]); + delete refs; + } + br->ResetAddress(); + // free the decompressed baskets if no later entry reads them, i.e. at the end of the cluster of this event + auto clusterIt = br->GetTree()->GetClusterIterator(event); + clusterIt.Next(); + if (event + 1 >= clusterIt.GetNextEntry()) { + br->DropBaskets("all"); + } + mTrackRefsLoaded[source][event] = true; +} + bool MCKinematicsReader::initFromDigitContext(o2::steer::DigitizationContext const* context) { if (mInitialized) { @@ -177,6 +194,7 @@ bool MCKinematicsReader::initFromDigitContext(o2::steer::DigitizationContext con mTracks.resize(mInputChains.size()); mHeaders.resize(mInputChains.size()); mIndexedTrackRefs.resize(mInputChains.size()); + mTrackRefsLoaded.resize(mInputChains.size()); // actual loading will be done only if someone asks // the first time for a particular source ... @@ -210,6 +228,7 @@ bool MCKinematicsReader::initFromKinematics(std::string_view name) mTracks.resize(1); mHeaders.resize(1); mIndexedTrackRefs.resize(1); + mTrackRefsLoaded.resize(1); mInitialized = true; return true;