From a2f1d7ec80d8d38adc053eb0fc8cee020757aef7 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Mon, 5 Oct 2026 13:11:43 +0200 Subject: [PATCH 01/13] MatchCosmics: include the last compatible vertex in the DCA veto vtIDMax is the id of the last compatible vertex (inclusive), but the veto loop stopped before it: a seed associated to a single vertex (vtIDMin == vtIDMax) was never checked, so with --use-pv-info the DCA cut never removed collision tracks associated to one vertex. Co-Authored-By: Claude Opus 5.5 --- Detectors/GlobalTracking/src/MatchCosmics.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Detectors/GlobalTracking/src/MatchCosmics.cxx b/Detectors/GlobalTracking/src/MatchCosmics.cxx index ca6bb50e3e9a8..7aec323a17f64 100644 --- a/Detectors/GlobalTracking/src/MatchCosmics.cxx +++ b/Detectors/GlobalTracking/src/MatchCosmics.cxx @@ -68,9 +68,9 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) trc.matchID = Reject; // reject track continue; } - if (mMatchParams->dcaCutChi2[trc.origID.getSource()] > 0.f && mUsePVInfo && (std::abs(trc.getY()) < mMatchParams->fiducialRIP && std::abs(trc.getZ()) < mMatchParams->fiducialZIP)) { + if (mMatchParams->dcaCutChi2[trc.origID.getSource()] > 0.f && mUsePVInfo && trc.vtIDMin >= 0 && (std::abs(trc.getY()) < mMatchParams->fiducialRIP && std::abs(trc.getZ()) < mMatchParams->fiducialZIP)) { // do the propagation only if we are in the fiducial IP range. - for (int iv = trc.vtIDMin; iv < trc.vtIDMax; iv++) { + for (int iv = trc.vtIDMin; iv <= trc.vtIDMax; iv++) { // vtIDMax is the last compatible vertex (inclusive); vtIDMin < 0: no compatible vertex const auto& pv = data.getPrimaryVertex(iv); o2::track::TrackParCov trcatPV(trc); o2::dataformats::DCA dca; From 1472e6a7e044c207205529e2309b7d21491b2a89 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Tue, 6 Oct 2026 19:39:36 +0200 Subject: [PATCH 02/13] MatchCosmics: selection of cosmics in collision data, named presets In collision data nearly all seeds are collision tracks and TPC-only legs had their z ignored, so any two tracks back-to-back in the transverse plane were paired (PbPb: minutes per TF, tens of thousands of fake pairs). New settings, all off by default and switched on together by the preset physics-v1 (o2-cosmics-match-workflow --cosmics-preset physics-v1; --configKeyValues can still override single keys, unknown names are fatal): - Seed cuts: minSeedNClTPC, minSeedDCAxy and minSeedDCAxyNSigma reject short TPC fragments and seeds passing the beam line (cosmics through the beam line are indistinguishable from collision tracks anyway). - z test of TPC-only legs (constrainTPCOnlyZ): a one-side TPC-only track has z relative to its time0, z(t) = z + side*vD*(t-time0); a CE-crossing or non-TPC-only track has an absolute z. The legs are compared at a common time (same side, or one leg with known time); for legs on opposite sides z continuity fixes the common time, which must lie in both time brackets. That time, with its error sigma_z / (2 vD) (~0.3 mus), is used for the refit and as the time of the cosmic (before: the centre of the bracket overlap, tens of mus, which put the legs at the wrong z: cosmics run 562658 A x C chi2Match median 3.5e4 -> 43). For all other pairs the time error of the cosmic stays the half-width of the bracket overlap. - Same-half veto (vetoSameHalf): pairs whose two legs lie on the same side of the closest approach in the transverse plane (two pieces of one leg) are rejected; skipped when a leg's reference point is within 20 cm of the closest approach. - maxChi2Match: cut on the chi2 of the refitted top/bottom legs. - minPtOppositeSides: for TPC-only legs on opposite sides z continuity defines the time, so z does not reject random pairs (94 % of the PbPb candidates). Both legs (in checkPair, before the winner selection) and the refitted cosmic need pT above it (scaled with the field). Unconditionally, the seeds and the refit use the muon mass (and charge 1) for the material corrections instead of the TPC track's PID from dE/dx (a leg tagged as kaon or proton got that mass): cosmic MC 6987 -> 7003 correct pairs (0 wrong), cosmics run 562658 357 -> 358; the refitted q/pt changes by +-0.4 % (MC) and +-1 % (562658) (q16-q84, pT < 20 GeV). physics-v1: pT > 1 GeV, |DCA_xy| > 3 cm and > 10 sigma, >= 30 TPC clusters; systematic errors for the leg comparison y 0.5 cm, snp 0.02, q/pt 0.05 (true pairs had pulls 2-4x too wide) with crude chi2 < 50; open per-parameter windows except tgl (3 sigma); chi2Match < 1000; z test and same-half veto on; minPtOppositeSides 2 GeV. Validation, TPC-only legs: cosmic MC with PbPb-like distortions 7003 of 8638 cosmics correct (81.1 %, 0 wrong pairs); cosmics run 562658 358 cosmics; PbPb 2025 LHC25an 567939 (19 kHz) 0.37 cosmics per TF (3.90 without the pT cut), 568041 reconstructed from CTF 0.30 per TF (3203 in 10546 TFs; the PbPb numbers before the muon mass, which leaves 4 PbPb TFs unchanged). Co-Authored-By: Claude Opus 5.5 --- .../include/GlobalTracking/MatchCosmics.h | 13 +- .../GlobalTracking/MatchCosmicsParams.h | 14 ++ Detectors/GlobalTracking/src/MatchCosmics.cxx | 121 ++++++++++++++++-- .../GlobalTracking/src/MatchCosmicsParams.cxx | 31 +++++ .../src/cosmics-match-workflow.cxx | 8 +- 5 files changed, 172 insertions(+), 15 deletions(-) diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h index 28cf1ab846412..acfc75f0679ba 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h @@ -17,6 +17,8 @@ #define ALICEO2_MATCH_COSMICS #include +#include +#include #include #include "ReconstructionDataFormats/TrackCosmics.h" #include "ReconstructionDataFormats/GlobalTrackID.h" @@ -67,7 +69,8 @@ class MatchCosmics RejTime, RejProp, RejChi2, - RejOther + RejOther, + RejSameHalf }; using InfoAccessor = o2d::AbstractRefAccessor; // there is no unique structure, so the default return type is dummy (int) @@ -78,6 +81,8 @@ class MatchCosmics int id1 = MinusOne; ///< id of 2nd parnter float chi2 = -1.f; ///< matching chi2 int next = MinusOne; ///< index of eventual next record + float tCommon = 0.f; ///< common time [mus] fixed by z continuity of TPC-only legs on opposite TPC sides + float tCommonErr = -1.f; ///< its 1 sigma error [mus]; < 0: not fixed, the refit uses the centre of the time-bracket overlap and the cosmic's time error is the overlap's half-width }; struct TrackSeed : public o2::track::TrackParCov { @@ -86,6 +91,9 @@ class MatchCosmics int matchID = MinusOne; ///< entry (none if MinusOne) of its match in the vector of matches short vtIDMin = -1; ///< id of the 1st compatible vertex short vtIDMax = -1; ///< id of the last compatible vertex + float tRef = 0.f; ///< time [mus] the z of the parameters refers to (TPC-only: the TrackTPC time0; others: bracket centre) + int8_t tpcSide = 0; ///< TPC-only seed with clusters on one side: +1 A, -1 C (z = z(t) - side*vD*(t-tRef)); 0: z absolute + std::array xyzRef{}; ///< global position of the reference point before the propagation to the DCA (same-half veto) }; void setTPCCorrMaps(const o2::gpu::TPCFastTransformPOD* maph); void setTPCVDrift(const o2::tpc::VDriftCorrFact& v); @@ -129,7 +137,7 @@ class MatchCosmics private: void updateTimeDependentParams(); RejFlag checkPair(int i, int j); - void registerMatch(int i, int j, float chi2); + void registerMatch(int i, int j, float chi2, float tCommon = 0.f, float tCommonErr = -1.f); void suppressMatch(int partner0, int partner1); void createSeeds(const o2::globaltracking::RecoContainer& data); bool validateMatch(int partner0); @@ -154,6 +162,7 @@ class MatchCosmics bool mUseMC = true; float mITSROFrameLengthMUS = 0.; float mQ2PtCutoff = 1e9; + float mQ2PtCutoffOppositeSides = 1e9; const MatchCosmicsParams* mMatchParams = nullptr; std::vector mCosmicTracks; diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h index 4b34135d83693..9307d92ee5e3d 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h @@ -18,6 +18,7 @@ #include "CommonUtils/ConfigurableParamHelper.h" #include "DetectorsBase/Propagator.h" #include "ReconstructionDataFormats/GlobalTrackID.h" +#include namespace o2 { @@ -29,10 +30,17 @@ struct MatchCosmicsParams : public o2::conf::ConfigurableParamHelper= this [cm] (0: no cut; rejects collision tracks in physics data) + float minSeedDCAxyNSigma = 0.f; // use only tracks with |DCA_xy| >= this * sigma(DCA_xy) (0: no cut; poorly measured collision tracks) + bool constrainTPCOnlyZ = false; // TPC-only legs: test z at a common time (same side, or a leg with known time), else require the time implied by z continuity in both brackets + bool vetoSameHalf = false; // reject pairs whose two legs lie on the same side of the closest approach (two pieces of one leg) float nSigmaTError = 4.f; // number of sigmas on track time error for matching (except for TPC which provides an interval) float tpcExtraZError2 = 1.f; // extra error^2 on the TPC-only track Z coordinate float fiducialRIP = 1.0f; // consider track having |Y@x=0|< this as passing DCA cut (if requested) @@ -44,6 +52,12 @@ struct MatchCosmicsParams : public o2::conf::ConfigurableParamHelperminSeedDCAxy || trc.getY() * trc.getY() < mMatchParams->minSeedDCAxyNSigma * mMatchParams->minSeedDCAxyNSigma * trc.getSigmaY2()) { + // passes close to the beam line, absolutely or within its errors: indistinguishable from collision tracks + trc.matchID = Reject; + continue; + } if (mMatchParams->dcaCutChi2[trc.origID.getSource()] > 0.f && mUsePVInfo && trc.vtIDMin >= 0 && (std::abs(trc.getY()) < mMatchParams->fiducialRIP && std::abs(trc.getZ()) < mMatchParams->fiducialZIP)) { // do the propagation only if we are in the fiducial IP range. for (int iv = trc.vtIDMin; iv <= trc.vtIDMax; iv++) { // vtIDMax is the last compatible vertex (inclusive); vtIDMin < 0: no compatible vertex @@ -180,9 +185,16 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) for (auto winRID : mWinners) { const auto& rec = mRecords[winRID]; int poolEntryID[2] = {rec.id0, rec.id1}; - const o2::track::TrackParCov outerLegs[2] = {data.getTrackParamOut(mSeeds[rec.id0].origID), data.getTrackParamOut(mSeeds[rec.id1].origID)}; + o2::track::TrackParCov outerLegs[2] = {data.getTrackParamOut(mSeeds[rec.id0].origID), data.getTrackParamOut(mSeeds[rec.id1].origID)}; + for (auto& leg : outerLegs) { + leg.setPID(o2::track::PID::Muon, true); // as the seeds + } auto tOverlap = mSeeds[rec.id0].tBracket.getOverlap(mSeeds[rec.id1].tBracket); float t0 = tOverlap.mean(), dt = tOverlap.delta() * 0.5; + if (rec.tCommonErr >= 0.f) { // TPC-only legs on opposite sides: their z continuity fixes the time, refit both legs with it + t0 = rec.tCommon; + dt = rec.tCommonErr; + } auto pnt0 = outerLegs[0].getXYZGlo(), pnt1 = outerLegs[1].getXYZGlo(); int btm = 0, top = 1; // we fit topward from bottom @@ -304,10 +316,19 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) // calculate weighted average of 2 legs and chi2 o2::track::TrackParCov::MatrixDSym5 cov5; float chi2Match = trCosmBtm.getPredictedChi2(trCosmTop, cov5); + if (mMatchParams->maxChi2Match >= 0.f && chi2Match > mMatchParams->maxChi2Match) { + LOG(debug) << "Top/Bottom refitted legs disagree, chi2Match " << chi2Match; + continue; + } if (!trCosmBtm.update(trCosmTop, cov5)) { LOG(debug) << "Top/Bottom update failed"; continue; } + // TPC-only legs on opposite sides: the legs' pT is required in checkPair, the refitted cosmic's here + if (mSeeds[rec.id0].tpcSide * mSeeds[rec.id1].tpcSide < 0 && std::abs(trCosmBtm.getQ2Pt()) > mQ2PtCutoffOppositeSides) { + LOG(debug) << "Cosmic with legs on opposite TPC sides below minPtOppositeSides"; + continue; + } // create final track mCosmicTracks.emplace_back(mSeeds[poolEntryID[btm]].origID, mSeeds[poolEntryID[top]].origID, trCosmBtm, trCosmTop, chi2, chi2Match, nclTot, t0, dt); if (mUseMC) { @@ -424,10 +445,18 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) return (rej = RejTime); // since the brackets are sorted in tmin, all following tbj will also exceed tbi } float chi2 = 1.e9f; + float tCommon = 0.f; // time fixed by z continuity of TPC-only legs on opposite sides, used by the refit + float tCommonErr = -1.f; // its error (< 0: not fixed) // check // 1) crude check on tgl and q/pt (if B!=0). Note: back-to-back tracks will have mutually params (see TrackPar::invertParam) while (1) { + // TPC-only legs on opposite sides: their z continuity defines the time, so z does not reject random pairs of collision tracks; require + // the pT of a cosmic for both legs already here, so that such a pair cannot win against the true partner of one of its legs + if (seed0.tpcSide * seed1.tpcSide < 0 && std::max(std::abs(seed0.getQ2Pt()), std::abs(seed1.getQ2Pt())) > mQ2PtCutoffOppositeSides) { + rej = RejQ2Pt; + break; + } auto dTgl = seed0.getTgl() + seed1.getTgl(); if (dTgl * dTgl > (mMatchParams->systSigma2[o2::track::kTgl] + seed0.getSigmaTgl2() + seed1.getSigmaTgl2()) * mMatchParams->crudeNSigma2Cut[o2::track::kTgl]) { rej = RejTgl; @@ -440,6 +469,27 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) break; } } + if (mMatchParams->vetoSameHalf) { + // a cosmic has its two legs on opposite sides of its closest approach to the beam line: project the legs' reference points (before + // the propagation to the DCA) on the transverse direction of seed0 at its DCA; pieces of one leg are on the same side. Transverse + // only: the z of a TPC-only track refers to its own time0, so z differences between the legs are meaningless. Skipped for + // reference points close to the DCA (e.g. ITS-containing tracks), where the sign is undefined. + std::array pca{}; + seed0.getXYZGlo(pca); + const float phi = seed0.getAlpha() + std::asin(seed0.getSnp()); + const float dir[2] = {std::cos(phi), std::sin(phi)}; + float proj0 = 0.f; + float proj1 = 0.f; + for (int k = 0; k < 2; k++) { + proj0 += (seed0.xyzRef[k] - pca[k]) * dir[k]; + proj1 += (seed1.xyzRef[k] - pca[k]) * dir[k]; + } + constexpr float MinDist = 20.f; // cm + if (std::abs(proj0) > MinDist && std::abs(proj1) > MinDist && proj0 * proj1 > 0.f) { + rej = RejSameHalf; + break; + } + } o2::track::TrackParCov seed1Inv = seed1; seed1Inv.invert(); for (int i = 0; i < o2::track::kNParams; i++) { // add systematic error @@ -462,9 +512,47 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) break; } bool ignoreZ = seed0.origID.getSource() == o2d::GlobalTrackID::TPC || seed1.origID.getSource() == o2d::GlobalTrackID::TPC; - // RSTODO this is simplification: one should constraint the TPC-only track Z by the time of other candidate (at least their difference of both tracks are TPC only). - // If the shift is large, eventually the tracks need to be refitted. - if (!ignoreZ) { // cut on Z makes no sense for TPC only tracks + if (ignoreZ && mMatchParams->constrainTPCOnlyZ) { + // a TPC-only track with clusters on one side has z relative to its time0: z(t) = z + side * vD * (t - tRef); a CE-crossing or + // non-TPC-only track has an absolute z (side 0). Bring both legs to a common time where possible and test z; for legs on opposite + // sides, z continuity fixes the common time, which must lie in both time brackets. + const int side0 = seed0.tpcSide; + const int side1 = seed1.tpcSide; + const float sigZ2 = (seed0.getSigmaZ2() + seed1Inv.getSigmaZ2()) * mMatchParams->crudeNSigma2Cut[o2::track::kZ]; + if (side0 == 0 || side1 == 0 || side0 == side1) { + float dZ = seed0.getZ() - seed1Inv.getZ(); + float dZTimeTol = 0.f; // the time of a non-TPC absolute leg is only known within its bracket (tRef is the bracket centre) + if (side0 != 0 && side1 != 0) { // same side: the z offset is fixed by the difference of the reference times + dZ -= side0 * mTPCVDrift * (seed0.tRef - seed1.tRef); + } else if (side1 != 0) { // seed0 absolute: move seed1 to the time of seed0 + dZ -= side1 * mTPCVDrift * (seed0.tRef - seed1.tRef); + if (seed0.origID.getSource() != o2d::GlobalTrackID::TPC) { + dZTimeTol = 0.5f * mTPCVDrift * seed0.tBracket.delta(); + } + } else if (side0 != 0) { // seed1 absolute: move seed0 to the time of seed1 + dZ += side0 * mTPCVDrift * (seed1.tRef - seed0.tRef); + if (seed1.origID.getSource() != o2d::GlobalTrackID::TPC) { + dZTimeTol = 0.5f * mTPCVDrift * seed1.tBracket.delta(); + } + } + const float dZTol = std::sqrt(sigZ2) + dZTimeTol; + if (dZ * dZ > dZTol * dZTol) { + rej = RejZ; + break; + } + } else { // opposite sides + const float t = 0.5f * (side0 * (seed1Inv.getZ() - seed0.getZ()) / mTPCVDrift + seed0.tRef + seed1.tRef); + const float tTol = std::sqrt(sigZ2) / (2.f * mTPCVDrift); + if (t < std::max(seed0.tBracket.getMin(), seed1.tBracket.getMin()) - tTol || t > std::min(seed0.tBracket.getMax(), seed1.tBracket.getMax()) + tTol) { + rej = RejZ; + break; + } + tCommon = t; + tCommonErr = std::sqrt(seed0.getSigmaZ2() + seed1Inv.getSigmaZ2()) / (2.f * mTPCVDrift); + } + } + // the z of a TPC-only leg refers to its own time0: it is tested above at a common time with constrainTPCOnlyZ, otherwise ignored + if (!ignoreZ) { // both legs have an absolute z auto dZ = seed0.getZ() - seed1Inv.getZ(); if (dZ * dZ > (seed0.getSigmaZ2() + seed1Inv.getSigmaZ2()) * mMatchParams->crudeNSigma2Cut[o2::track::kZ]) { rej = RejZ; @@ -484,8 +572,8 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) break; } rej = Accept; - registerMatch(i, j, chi2); - registerMatch(j, i, chi2); // the reverse reference can be also done in a separate loop + registerMatch(i, j, chi2, tCommon, tCommonErr); + registerMatch(j, i, chi2, tCommon, tCommonErr); // the reverse reference can be also done in a separate loop LOG(debug) << "Chi2 = " << chi2 << " NMatches " << mRecords.size(); break; } @@ -497,7 +585,8 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) if (seed1I.rotate(seed0.getAlpha()) && o2::base::Propagator::Instance()->PropagateToXBxByBz(seed1I, seed0.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { int rejI = int(rej); (*mDBGOut) << "match" - << "tf=" << mTFCount << "seed0=" << seed0 << "seed1=" << seed1I << "chi2Match=" << chi2 << "rej=" << rejI << "\n"; + << "tf=" << mTFCount << "seed0=" << seed0 << "seed1=" << seed1I << "chi2Match=" << chi2 << "rej=" << rejI + << "side0=" << int(seed0.tpcSide) << "side1=" << int(seed1.tpcSide) << "tCommon=" << tCommon << "tCommonErr=" << tCommonErr << "\n"; } } #endif @@ -506,11 +595,11 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) } //________________________________________________________ -void MatchCosmics::registerMatch(int i, int j, float chi2) +void MatchCosmics::registerMatch(int i, int j, float chi2, float tCommon, float tCommonErr) { /// register track index j as a match for track index i int newRef = mRecords.size(); - auto& matchRec = mRecords.emplace_back(MatchRecord{i, j, chi2, MinusOne}); + auto& matchRec = mRecords.emplace_back(MatchRecord{i, j, chi2, MinusOne, tCommon, tCommonErr}); auto* best = &mSeeds[i].matchID; while (*best > MinusOne) { auto& oldMatchRec = mRecords[*best]; @@ -540,7 +629,7 @@ void MatchCosmics::createSeeds(const o2::globaltracking::RecoContainer& data) return true; } if constexpr (isTPCTrack()) { - if (!this->mMatchParams->allowTPCOnly) { + if (!this->mMatchParams->allowTPCOnly || _tr.getNClusters() < this->mMatchParams->minSeedNClTPC) { return true; } // unconstrained TPC track, with t0 = TrackTPC.getTime0+0.5*(DeltaFwd-DeltaBwd) and terr = 0.5*(DeltaFwd+DeltaBwd) in TimeBins @@ -554,9 +643,14 @@ void MatchCosmics::createSeeds(const o2::globaltracking::RecoContainer& data) } terr += this->mMatchParams->timeToleranceMUS; trackEntry[_origID] = mSeeds.size(); - mSeeds.emplace_back(TrackSeed{_tr, {t0 - terr, t0 + terr}, _origID, MinusOne}); + auto& seed = mSeeds.emplace_back(TrackSeed{_tr, {t0 - terr, t0 + terr}, _origID, MinusOne}); + seed.setPID(o2::track::PID::Muon, true); // muon mass and charge for the material corrections, whatever the leg's dE/dx PID + seed.getXYZGlo(seed.xyzRef); + seed.tRef = t0; if constexpr (isTPCTrack()) { - mSeeds.back().setCov(this->mMatchParams->tpcExtraZError2 + _tr.getCov()[o2::track::kSigZ2], o2::track::kSigZ2); + seed.setCov(this->mMatchParams->tpcExtraZError2 + _tr.getCov()[o2::track::kSigZ2], o2::track::kSigZ2); + seed.tRef = _tr.getTime0() * this->mTPCTBinMUS; // the z of a TPC-only track refers to its time0 + seed.tpcSide = _tr.hasASideClustersOnly() ? 1 : (_tr.hasCSideClustersOnly() ? -1 : 0); } return true; } else { @@ -603,10 +697,13 @@ void MatchCosmics::updateTimeDependentParams() mBz = o2::base::Propagator::Instance()->getNominalBz(); mFieldON = std::abs(mBz) > 0.01; mQ2PtCutoff = 1.f / std::max(0.05f, mMatchParams->minSeedPt); + mQ2PtCutoffOppositeSides = mMatchParams->minPtOppositeSides > 0.f ? 1.f / mMatchParams->minPtOppositeSides : 1e9; if (mFieldON) { mQ2PtCutoff *= 5.00668 / std::abs(mBz); + mQ2PtCutoffOppositeSides *= 5.00668 / std::abs(mBz); } else { mQ2PtCutoff = 1e9; + mQ2PtCutoffOppositeSides = 1e9; } } diff --git a/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx b/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx index f14ae04897c68..fa746111ede0f 100644 --- a/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx +++ b/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx @@ -14,4 +14,35 @@ /// \author ruben.shahoyan@cern.ch #include "GlobalTracking/MatchCosmicsParams.h" +#include "Framework/Logger.h" +#include + O2ParamImpl(o2::globaltracking::MatchCosmicsParams); + +std::string o2::globaltracking::getMatchCosmicsPreset(const std::string& name) +{ + // physics-v1, tuned on PbPb 2025 (LHC25an 567939, 19 kHz) for fakes, on the cosmics run 562658 and on a cosmic MC with PbPb-like + // distortions for the efficiency (TPC-only legs): + // - seed cuts against collision tracks: pT, transverse DCA to the beam line in cm and in sigma, TPC clusters + // - realistic systematic errors for the leg comparison (true pairs had y / snp / q/pt pulls 2-4x too wide), so that the crude pair + // chi2 cut is meaningful; per-parameter windows open except tgl (3 sigma: the main discriminant against random pairs) + // - loose cut on the chi2 of the refitted legs; z test of TPC-only legs at a common time and same-half veto switched on + // - TPC-only legs on opposite sides (time from z continuity, ~94 % of the PbPb candidates were random pairs): both legs and the + // refitted cosmic above 2 GeV (offline: PbPb candidates / 10.5, cosmic MC efficiency 81.8 -> 80.6 %) + static const std::map presets{ + {"physics-v1", + "cosmicsMatch.minSeedPt=1;cosmicsMatch.minSeedDCAxy=3;cosmicsMatch.minSeedDCAxyNSigma=10;cosmicsMatch.minSeedNClTPC=30;" + "cosmicsMatch.crudeChi2Cut=50;cosmicsMatch.systSigma2[0]=0.25;cosmicsMatch.systSigma2[2]=4e-4;cosmicsMatch.systSigma2[4]=2.5e-3;" + "cosmicsMatch.crudeNSigma2Cut[0]=144;cosmicsMatch.crudeNSigma2Cut[2]=144;cosmicsMatch.crudeNSigma2Cut[3]=9;cosmicsMatch.crudeNSigma2Cut[4]=144;" + "cosmicsMatch.maxChi2Match=1000;cosmicsMatch.constrainTPCOnlyZ=true;cosmicsMatch.vetoSameHalf=true;cosmicsMatch.minPtOppositeSides=2"}}; + auto it = presets.find(name); + if (it == presets.end()) { + std::string known; + for (const auto& [k, v] : presets) { + known += " " + k; + } + LOGP(fatal, "Unknown cosmics matching preset {}, known:{}", name, known); + return {}; + } + return it->second; +} diff --git a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx index 67e2fd6cdc4cd..85724125b1c35 100644 --- a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx @@ -26,6 +26,7 @@ #include "DetectorsCommonDataFormats/DetID.h" #include "GlobalTrackingWorkflowReaders/TrackTPCITSReaderSpec.h" #include "GlobalTrackingWorkflow/CosmicsMatchingSpec.h" +#include "GlobalTracking/MatchCosmicsParams.h" #include "GlobalTrackingWorkflow/TrackCosmicsWriterSpec.h" #include "Algorithm/RangeTokenizer.h" #include "DetectorsRaw/HBFUtilsInitializer.h" @@ -52,6 +53,7 @@ void customize(std::vector& workflowOptions) {"disable-root-input", o2::framework::VariantType::Bool, false, {"disable root-files input reader"}}, {"disable-root-output", o2::framework::VariantType::Bool, false, {"disable root-files output writer"}}, {"use-pv-info", o2::framework::VariantType::Bool, false, {"request primary vertex for relevant cuts in the collision/cosmics interleaved data"}}, + {"cosmics-preset", VariantType::String, "", {"named set of cosmicsMatch settings applied before --configKeyValues (which can override single keys): physics-v1 = cosmics in collision data"}}, {"track-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of sources to use"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}}; o2::itsmft::DPLAlpideParamInitializer::addITSConfigOption(options); @@ -82,7 +84,11 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) WorkflowSpec specs; GID::mask_t alowedSources = GID::getSourcesMask("ITS,TPC,ITS-TPC,TPC-TRD,TPC-TOF,TPC-TRD-TOF,ITS-TPC-TOF,ITS-TPC-TRD-TOF"); - // Update the (declared) parameters if changed from the command line + // Update the (declared) parameters if changed from the command line: first an eventual preset, then the explicit key=values + auto preset = configcontext.options().get("cosmics-preset"); + if (!preset.empty()) { + o2::conf::ConfigurableParam::updateFromString(o2::globaltracking::getMatchCosmicsPreset(preset)); + } o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); // write the configuration used for the workflow o2::conf::ConfigurableParam::writeINI("o2match-cosmics-workflow_configuration.ini"); From 7479fef3069a4ec22c09e7f19ba3f48b53970f39 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Tue, 6 Oct 2026 22:08:58 +0200 Subject: [PATCH 03/13] MatchCosmics: energy-loss sign along the muon's flight; optional energy-loss sign in the TPC refit The propagator applies the energy-loss correction along the direction of the track parameters (loss forward, gain backward), and refitWinners fits a cosmic from the bottom leg to the top one: after the inward refit of the bottom leg its parameters are inverted and carried up through the top leg, i.e. they describe a particle moving upward. A cosmic muon flies downward, from the top leg to the bottom leg, so only the inward refit of the bottom leg (backward along its parameters = backward along the flight) had the physical sign; the propagation of the inverted bottom leg to the top leg, the outward refit of the top leg, the inward refit of the top leg, its propagation to the bottom parameter and the propagation of the top seeds to the DCA had the opposite one. The sign is now imposed along the flight in all of them, with the propagator's signCorr argument. GPUO2InterfaceRefit / GPUTrackingRefit::RefitTrackAsTrackParCov (TrackParCov + cluster references) get the analogous optional last argument eLossSign, passed to their propagations (0 = along the parameters as before, +1 gain, -1 loss; flipped on the optional way to the reference X when that runs against the refit direction); the default keeps the behaviour for all other users. MC truth (cosmic MC, 1073 muons crossing the inner region; true parameters at the bottom TPC inner radius propagated back through the closest approach to the top one, compared to the true momentum there, q16 / q50 / q84 in MeV): true loss between the two 14.5 / 16.1 / 19.9 MeV (median for |z| at the closest approach 0-25 / 25-50 / 50-75 cm); predicted - true with the sign of the upward fit -14.5 / -15.2 / -18.6 MeV (as without material: the gain inward and the loss outward cancel), with the sign along the flight +1.2 / +2.1 / +1.0 MeV (spread q16-q84 about -2 ... +9 MeV). The matched cosmics hardly change because the two halves' errors cancel in their average: cosmic MC 7003 -> 7006 correct pairs (0 wrong), chi2Match median unchanged; cosmics run 562658 358 cosmics before and after; refitted q/pt within +-0.05 %. The halves themselves (e.g. the momentum at the TPC entrance) change by the full loss. Co-Authored-By: Claude Opus 5.5 --- Detectors/GlobalTracking/src/MatchCosmics.cxx | 38 +++++++++++++------ .../Interface/GPUO2InterfaceRefit.cxx | 2 +- .../Interface/GPUO2InterfaceRefit.h | 3 +- GPU/GPUTracking/Refit/GPUTrackingRefit.cxx | 10 +++-- GPU/GPUTracking/Refit/GPUTrackingRefit.h | 5 ++- 5 files changed, 39 insertions(+), 19 deletions(-) diff --git a/Detectors/GlobalTracking/src/MatchCosmics.cxx b/Detectors/GlobalTracking/src/MatchCosmics.cxx index 55f80dc8c9b83..2e1fa08b1c146 100644 --- a/Detectors/GlobalTracking/src/MatchCosmics.cxx +++ b/Detectors/GlobalTracking/src/MatchCosmics.cxx @@ -45,6 +45,15 @@ using namespace o2::globaltracking; using GTrackID = o2d::GlobalTrackID; using MatCorrType = o2::base::Propagator::MatCorrType; +namespace +{ +// energy-loss sign of the propagations: the propagator applies the loss along the direction of the track parameters, and the refit runs +// from the bottom leg up through the top leg (parameters of an upward-moving particle); the cosmic muon flies from the top leg to the +// bottom leg, so the sign is imposed along its flight +constexpr int ELossGain = 1; +constexpr int ELossLoss = -1; +} // namespace + //________________________________________________________ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) { @@ -64,7 +73,12 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) for (int i = 0; i < ntr; i++) { auto& trc = mSeeds[i]; if (trc.matchID != Reject) { - if (!prop->propagateToDCABxByBz(v, trc, mMatchParams->maxStep, mMatchParams->matCorr)) { + // a cosmic muon flies downward: along a leg whose outward direction points up (top leg) the inward propagation follows the flight, + // so the energy is lost; along the bottom leg it goes back in the flight, so the energy is gained + std::array momentum{}; + trc.getPxPyPzGlo(momentum); + const int eLossSign = momentum[1] > 0.f ? ELossLoss : ELossGain; + if (!prop->propagateToDCABxByBz(v, trc, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, nullptr, eLossSign)) { trc.matchID = Reject; // reject track continue; } @@ -154,7 +168,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) } } - auto refitITSTrack = [this, &data, &itsTracksROF, &itsClusters](o2::track::TrackParCov& trFit, GTrackID gidx, float& chi2, bool inward = false) { + auto refitITSTrack = [this, &data, &itsTracksROF, &itsClusters](o2::track::TrackParCov& trFit, GTrackID gidx, float& chi2, bool inward, int eLossSign) { const auto& itsTrOrig = data.getITSTrack(gidx); int nclRefit = 0, ncl = itsTrOrig.getNumberOfClusters(), rof = itsTracksROF[gidx.getIndex()]; const auto& itsTrackClusRefs = data.getITSTracksClusterRefs(); @@ -170,7 +184,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) for (int icl = from; icl != to; icl += step) { // ITS clusters are referred in layer decreasing order const auto& clus = itsClusters[itsTrackClusRefs[clEntry + icl]]; float alpha = geomITS->getSensorRefAlpha(clus.getSensorID()), x = clus.getX(); - if (!trFit.rotate(alpha) || !propagator->propagateToX(trFit, x, propagator->getNominalBz(), this->mMatchParams->maxSnp, this->mMatchParams->maxStep, this->mMatchParams->matCorr)) { + if (!trFit.rotate(alpha) || !propagator->propagateToX(trFit, x, propagator->getNominalBz(), this->mMatchParams->maxSnp, this->mMatchParams->maxStep, this->mMatchParams->matCorr, nullptr, eLossSign)) { break; } chi2 += trFit.getPredictedChi2(clus); @@ -221,7 +235,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) if (!mFieldON) { trCosm.setQ2Pt(-o2::track::kMostProbablePt); } - int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, false, false); // inward refit, reset + int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, false, false, ELossGain); // inward refit, reset if (retVal < 0) { // refit failed LOG(debug) << "Inward refit of btm TPC track failed."; continue; @@ -247,7 +261,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) } trCosm.invert(); if (!trCosm.rotate(mSeeds[poolEntryID[top]].getAlpha()) || - !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, mSeeds[poolEntryID[top]].getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { + !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, mSeeds[poolEntryID[top]].getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossGain)) { LOG(debug) << "Rotation/propagation of btm-track to top-track frame failed."; continue; } @@ -260,7 +274,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) // is there ITS sub-track? if (gidxListTop[GTrackID::ITS].isIndexSet()) { - auto nclfit = refitITSTrack(trCosm, gidxListTop[GTrackID::ITS], chi2, false); + auto nclfit = refitITSTrack(trCosm, gidxListTop[GTrackID::ITS], chi2, false, ELossGain); if (nclfit < 0) { continue; } @@ -273,13 +287,13 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) if (trCosm.getX() * trCosm.getX() + trCosm.getY() * trCosm.getY() <= o2::constants::geom::XTPCInnerRef * o2::constants::geom::XTPCInnerRef) { float xtogo = 0; if (!trCosm.getXatLabR(o2::constants::geom::XTPCInnerRef, xtogo, mBz, o2::track::DirOutward) || - !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, xtogo, mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { + !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, xtogo, mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossGain)) { LOG(debug) << "Propagation to inner TPC boundary X=" << xtogo << " failed"; continue; } } const auto& tpcTrOrig = data.getTPCTrack(gidxListTop[GTrackID::TPC]); - int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, true, false); // outward refit, no reset + int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, true, false, ELossGain); // outward refit, no reset if (retVal < 0) { // refit failed LOG(debug) << "Outward refit of top TPC track failed."; continue; @@ -293,15 +307,15 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) auto trCosmTop = outerLegs[top]; if (gidxListTop[GTrackID::TPC].isIndexSet()) { // inward refit in TPC const auto& tpcTrOrig = data.getTPCTrack(gidxListTop[GTrackID::TPC]); - int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosmTop, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2Dummy, false, true); // inward refit, reset + int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosmTop, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2Dummy, false, true, ELossLoss); // inward refit, reset if (retVal < 0) { // refit failed - LOG(debug) << "Outward refit of top TPC track failed."; + LOG(debug) << "Inward refit of top TPC track failed."; continue; } // inward refit in TPC } // is there ITS sub-track ? if (gidxListTop[GTrackID::ITS].isIndexSet()) { - auto nclfit = refitITSTrack(trCosmTop, gidxListTop[GTrackID::ITS], chi2Dummy, true); + auto nclfit = refitITSTrack(trCosmTop, gidxListTop[GTrackID::ITS], chi2Dummy, true, ELossLoss); if (nclfit < 0) { continue; } @@ -309,7 +323,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) } // ITS refit // propagate to bottom param if (!trCosmTop.rotate(trCosmBtm.getAlpha()) || - !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosmTop, trCosmBtm.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { + !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosmTop, trCosmBtm.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossLoss)) { LOG(debug) << "Rotation/propagation of top-track to bottom-track frame failed."; continue; } diff --git a/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx b/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx index cd184b3820533..aa9a036c4fb41 100644 --- a/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx +++ b/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.cxx @@ -136,7 +136,7 @@ void GPUO2InterfaceRefit::updateCalib(const TPCFastTransformPOD* trans, float bz int32_t GPUO2InterfaceRefit::RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward, bool resetCov) { return mRefit->RefitTrackAsGPU(trk, outward, resetCov); } int32_t GPUO2InterfaceRefit::RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward, bool resetCov) { return mRefit->RefitTrackAsTrackParCov(trk, outward, resetCov); } int32_t GPUO2InterfaceRefit::RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2, bool outward, bool resetCov) { return mRefit->RefitTrackAsGPU(trk, clusRef, time0, chi2, outward, resetCov); } -int32_t GPUO2InterfaceRefit::RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2, bool outward, bool resetCov) { return mRefit->RefitTrackAsTrackParCov(trk, clusRef, time0, chi2, outward, resetCov); } +int32_t GPUO2InterfaceRefit::RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2, bool outward, bool resetCov, int32_t eLossSign) { return mRefit->RefitTrackAsTrackParCov(trk, clusRef, time0, chi2, outward, resetCov, eLossSign); } void GPUO2InterfaceRefit::setIgnoreErrorsAtTrackEnds(bool v) { mRefit->mIgnoreErrorsOnTrackEnds = v; } void GPUO2InterfaceRefit::setTrackReferenceX(float v) { mParam->rec.tpc.trackReferenceX = v; } diff --git a/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h b/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h index f85a376b9185a..052f728821cac 100644 --- a/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h +++ b/GPU/GPUTracking/Interface/GPUO2InterfaceRefit.h @@ -67,7 +67,8 @@ class GPUO2InterfaceRefit int32_t RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false); int32_t RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false); int32_t RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false); - int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false); + // eLossSign: energy-loss sign of the propagations, 0 along the track parameters (loss forward, gain backward), +1 gain, -1 loss + int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false, int32_t eLossSign = 0); void setTrackReferenceX(float v); void setIgnoreErrorsAtTrackEnds(bool v); void updateCalib(const o2::gpu::TPCFastTransformPOD* trans, float bzNominalGPU); diff --git a/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx b/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx index 9e86e4627a2fb..c8d6b00d916cd 100644 --- a/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx +++ b/GPU/GPUTracking/Refit/GPUTrackingRefit.cxx @@ -221,6 +221,7 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov convertTrack::propagator>(trk, trkX, prop, &TrackParCovChi2); int32_t begin = 0, count; float tOffset; + [[maybe_unused]] int32_t eLossSign = 0; // energy-loss sign of the TrackParCov propagations (0: from the direction, outward = loss) if constexpr (std::is_same_v) { count = trkX.NClusters(); tOffset = trkX.GetParam().GetTOffset(); @@ -230,6 +231,7 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov } else if constexpr (std::is_same_v) { count = trkX.clusRef.getEntries(); tOffset = trkX.time0; + eLossSign = trkX.eLossSign; } else { static_assert("Invalid template"); } @@ -357,7 +359,7 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov IgnoreErrors(trk.getSnp()); return -1; } - if (!prop->PropagateToXBxByBz(trk, x, constants::MAX_SIN_PHI_LOW)) { + if (!prop->PropagateToXBxByBz(trk, x, constants::MAX_SIN_PHI_LOW, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT, nullptr, eLossSign)) { IgnoreErrors(trk.getSnp()); return -2; } @@ -401,11 +403,13 @@ GPUd() int32_t GPUTrackingRefit::RefitTrack(T& trkX, bool outward, bool resetCov constexpr float kDeg2Rad = M_PI / 180.f; constexpr float kSectAngle = 2 * M_PI / 18.f; if (mPparam->rec.tpc.trackReferenceX <= 500) { - if (prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX)) { + // a forced energy-loss sign holds along the refit direction: the way to the reference X can run against it + const int32_t eLossSignRef = (mPparam->rec.tpc.trackReferenceX > trk.getX()) == outward ? eLossSign : -eLossSign; + if (prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX, Propagator::MAX_SIN_PHI, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT, nullptr, eLossSignRef)) { if (CAMath::Abs(trk.getY()) > trk.getX() * CAMath::Tan(kSectAngle / 2.f)) { float newAlpha = trk.getAlpha() + CAMath::Round(CAMath::ATan2(trk.getY(), trk.getX()) / kDeg2Rad / 20.f) * kSectAngle; GPUTPCGMTrackParam::NormalizeAlpha(newAlpha); - trk.rotate(newAlpha) && prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX); + trk.rotate(newAlpha) && prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX, Propagator::MAX_SIN_PHI, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT, nullptr, eLossSignRef); } } } diff --git a/GPU/GPUTracking/Refit/GPUTrackingRefit.h b/GPU/GPUTracking/Refit/GPUTrackingRefit.h index 70c9fd47d90f6..71f8f11bee3eb 100644 --- a/GPU/GPUTracking/Refit/GPUTrackingRefit.h +++ b/GPU/GPUTracking/Refit/GPUTrackingRefit.h @@ -73,15 +73,16 @@ class GPUTrackingRefit const o2::tpc::TrackTPCClusRef& clusRef; float time0; float* chi2; + int32_t eLossSign = 0; // energy-loss sign of the TrackParCov propagations: 0 along the track parameters (loss forward), +1 gain, -1 loss }; GPUd() int32_t RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false) { TrackParCovWithArgs x{trk, clusRef, time0, chi2}; return RefitTrack(x, outward, resetCov); } - GPUd() int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false) + GPUd() int32_t RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false, int32_t eLossSign = 0) { - TrackParCovWithArgs x{trk, clusRef, time0, chi2}; + TrackParCovWithArgs x{trk, clusRef, time0, chi2, eLossSign}; return RefitTrack(x, outward, resetCov); } From 9bd67e9d9f67a18c7a485652c00f5bcd82de3f6b Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Thu, 8 Oct 2026 08:57:58 +0200 Subject: [PATCH 04/13] Cosmics: collect the raw clusters of matched cosmics and tag them with TOF (--enable-cluster-output) The async reconstruction stores neither TPC tracks nor TPC clusters, so the leg references of TrackCosmics cannot be resolved offline. New device cosmics-cluster-collector (+ cosmics-full-writer -> o2_cosmics_full.root, tree cosmicsFull, one entry per TF): - per cosmic (o2::dataformats::CosmicTrack, DataFormatsGlobalTracking): the TrackCosmics, the TPC tracks of both legs, the raw ClusterNative + sector/row of each leg's attached clusters and of all clusters in a road around the leg (flags: attached / road / used by any TPC track / found with the absolute time on the other TPC side), the ITS compact clusters (chip, anchor pixel, pattern ID, ROF BC; pattern bytes for patterns not in the dictionary), TOF clusters (raw time, ToT, channel), TRD tracklet words (+ trigger BC), MC label; per TF: TFIDInfo and vDrift / t0 of the TPC transformation. Only raw detector data are stored: TPC transformation, TOF calibration, ITS dictionary and geometry are applied offline. - TPC road: each leg's helix through all pad rows of its side, in the frame of its own time0, at the real x of the row, points beyond the leg's closest approach skipped; predicted real (y, z) mapped to nominal coordinates with the inverse correction; clusters within --corridor-width (1 cm) taken. The other TPC side is searched when the time of the cosmic is known (CE-crossing leg, z-continuity time with error < --max-abs-time-err 0.5 mus, or a TOF time). - TOF tag: all TOF clusters within the cosmic's time window and the road (--tof-road-width 5 cm) are candidates per leg; the top/bottom pair whose time difference matches the muon's flight along the helix (|t_top - t_bottom + L/c| < --tof-flight-tolerance 2 ns, z of one-side legs shifted to the pair's mean time) gets HitTOFFlight and the cosmic its TOF time (CosmicTrack::timeTOFMUS); without such a pair the closest hit per leg is kept. The same search in the impossible order (bottom hit first) only finds accidental pairs: its best score is stored next to the flagged pair's (scoreTOFPair / scoreTOFReversed) as QA of the flag's background. - The TOF time (--tof-time-error 0.1 mus) then fixes the other-side TPC corridor, the TRD road (trigger and z, --trd-road-width) and the ITS road (readout frame and z, --its-road-width); without it the TPC time is used. - Duplicates: per TF the cosmics are ranked (TOF time first, then more attached TPC clusters); one sharing >= 30 % of the smaller one's TPC clusters with a better one gets duplicateOf = that one's entry (a leg split into two TPC tracks). Nothing is removed. - Polish: cosmics with a TOF time are refitted at it (TPC-only legs: bottom leg inward to the closest approach to the beam line, top leg inward to the same point, halves combined; muon mass and energy loss along the flight as in the matcher), CosmicTrack::polished with chi2MatchPolished (< 0: not polished). The matcher's time is the bracket-overlap centre for one-side legs on the same side (tens of mus off), which puts the inner material at the wrong z. Check: polished at the matcher's time instead, the cosmic MC reproduces the matcher (chi2Match median 79.1 vs 79.2, q/pt within 7e-4). - Legs nearly parallel to the pad rows (closest approach to the beam line inside the TPC or at its inner edge): the row-by-row road cannot reach the rows, there the helix is followed in steps of path length and the rows within the road width are searched. - A TOF time can be negative (cosmic just before the TF start): the flag is CosmicTrack::hasTOFTime() (scoreTOFPair >= 0). - At most --max-cosmics-per-tf (100) cosmics per TF are written (protection against fake-dominated settings). - --debug-tree: cosmics_collector_debug.root, tree "cosmics" with one entry per cosmic (matching and timing quantities, vectors of the transformed TPC clusters, road points, ITS / TOF / TRD hits). Off by default. Validation: the road re-finds 99.8 % of the legs' own clusters in the cosmics run 562658 and 97 % in PbPb 2025 with the v10n correction map. PbPb 2025 LHC25an 567939 (600 TFs, preset physics-v1): 30 of 224 cosmics TOF-tagged (10 with legs on opposite TPC sides, 20 on one side), 0 pairs in the impossible order (flag purity >= 90 % at 95 % CL), flight-time residual centred at 0 within 0.2 ns; 10 duplicates. 568041 reconstructed from CTF (10546 TFs): 3203 cosmics, 435 TOF-tagged non-duplicates, 90 duplicates. CPU: matcher + collector ~0.6 s per PbPb TF; output ~1.5 kB per TF. Co-Authored-By: Claude Opus 5.5 --- .../Detectors/GlobalTracking/CMakeLists.txt | 1 + .../DataFormatsGlobalTracking/CosmicTrack.h | 127 ++ .../src/DataFormatsGlobalTrackingLinkDef.h | 12 + .../GlobalTrackingWorkflow/CMakeLists.txt | 1 + .../CosmicsClusterCollectorSpec.h | 31 + .../TrackCosmicsWriterSpec.h | 3 + .../src/CosmicsClusterCollectorSpec.cxx | 1641 +++++++++++++++++ .../src/TrackCosmicsWriterSpec.cxx | 17 + .../src/cosmics-match-workflow.cxx | 24 +- 9 files changed, 1856 insertions(+), 1 deletion(-) create mode 100644 DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/CosmicTrack.h create mode 100644 Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h create mode 100644 Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx diff --git a/DataFormats/Detectors/GlobalTracking/CMakeLists.txt b/DataFormats/Detectors/GlobalTracking/CMakeLists.txt index b219de73f5b47..53ce0ad335805 100644 --- a/DataFormats/Detectors/GlobalTracking/CMakeLists.txt +++ b/DataFormats/Detectors/GlobalTracking/CMakeLists.txt @@ -45,4 +45,5 @@ o2_target_root_dictionary( DataFormatsGlobalTracking HEADERS include/DataFormatsGlobalTracking/FilteredRecoTF.h include/DataFormatsGlobalTracking/TrackTuneParams.h + include/DataFormatsGlobalTracking/CosmicTrack.h ) diff --git a/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/CosmicTrack.h b/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/CosmicTrack.h new file mode 100644 index 0000000000000..ec7d85fc58793 --- /dev/null +++ b/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/CosmicTrack.h @@ -0,0 +1,127 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file CosmicTrack.h +/// \brief Matched cosmic track with the raw clusters of its legs and of the road around them, for offline refits + +#ifndef ALICEO2_COSMIC_TRACK_H +#define ALICEO2_COSMIC_TRACK_H + +#include +#include +#include +#include "ReconstructionDataFormats/TrackCosmics.h" +#include "DataFormatsTPC/TrackTPC.h" +#include "DataFormatsTPC/ClusterNative.h" +#include "SimulationDataFormat/MCCompLabel.h" + +namespace o2::dataformats +{ + +/// raw TPC cluster with its address; transformed coordinates are not stored (re-transform offline with the calibration of the TF) +struct CosmicTPCCluster { + enum Flags : uint8_t { + Attached = 0x1, ///< attached to the TPC track of this leg + Corridor = 0x2, ///< found in the road around the leg + Used = 0x4, ///< attached to some TPC track (for corridor clusters: another track, e.g. a split piece of the leg) + AbsTime = 0x8 ///< found on the other TPC side with the absolute time of the cosmic instead of the time0 of the leg: CosmicTrack::timeTOFMUS + ///< if >= 0, else the time of CosmicTrack::cosmic + }; + o2::tpc::ClusterNative cl{}; ///< raw cluster: time, pad, widths, charges, flags + uint8_t sector = 0; + uint8_t row = 0; + uint8_t flags = 0; + bool isAttached() const { return flags & Attached; } + bool isCorridor() const { return flags & Corridor; } + bool isUsed() const { return flags & Used; } + bool isAbsTime() const { return flags & AbsTime; } + ClassDefNV(CosmicTPCCluster, 1); +}; + +/// origin of the ITS / TOF / TRD hits of a cosmic +enum CosmicHitFlags : uint8_t { + HitMatched = 0x1, ///< part of the leg's matched global track + HitRoad = 0x2, ///< found by the road search around the cosmic + HitTOFFlight = 0x4 ///< TOF road hit of the top/bottom pair whose time difference matches the muon's flight (fixes the cosmic's time) +}; + +/// ITS cluster of a leg: the raw compact cluster (chip, anchor pixel, pattern ID) and its ROF; no coordinates (the dictionary and the +/// geometry are applied offline) +struct CosmicITSCluster { + uint16_t chipID = 0; + uint16_t row = 0; + uint16_t col = 0; + uint16_t pattID = 0; + int32_t pattEntry = -1; ///< start of the pattern bytes in CosmicTrack::itsPatterns (pattern not in the dictionary or group pattern), -1: none + int32_t rofBC = 0; ///< start of the cluster's ROF in BCs since the start of the TF + uint8_t leg = 0; ///< 0 bottom, 1 top + uint8_t flags = 0; ///< CosmicHitFlags + ClassDefNV(CosmicITSCluster, 1); +}; + +/// TOF cluster of a leg: matched to its track or found on the road +struct CosmicTOFCluster { + double timeRaw = 0.; ///< raw TOF time [ps] (the calibration is applied offline) + float tot = 0.f; ///< time over threshold + int32_t channel = -1; + uint8_t leg = 0; + uint8_t flags = 0; ///< CosmicHitFlags + ClassDefNV(CosmicTOFCluster, 1); +}; + +/// TRD tracklet of a leg: attached to its track or found on the road +struct CosmicTRDTracklet { + uint64_t word = 0; ///< raw Tracklet64 word + int32_t trigBC = 0; ///< BC of its trigger since the start of the TF + uint8_t layer = 0; + uint8_t leg = 0; + uint8_t flags = 0; ///< CosmicHitFlags + ClassDefNV(CosmicTRDTracklet, 1); +}; + +/// matched cosmic with everything needed for an offline refit +struct CosmicTrack { + o2::dataformats::TrackCosmics cosmic{}; ///< matcher output (time in mus; the leg references are only valid within the TF) + o2::tpc::TrackTPC tpcBottom{}; ///< TPC part of the bottom leg (default if none); z refers to its time0; its cluster + ///< references are only valid within the TF: the clusters are in clTPCBottom + o2::tpc::TrackTPC tpcTop{}; ///< TPC part of the top leg + std::vector clTPCBottom; ///< attached + corridor TPC clusters of the bottom leg + std::vector clTPCTop; ///< same for the top leg + std::vector clITS; ///< ITS clusters of the legs' matched tracks and on the road + std::vector itsPatterns; ///< pattern bytes (row span, column span, bitmap) of the ITS clusters that need them + std::vector clTOF; ///< TOF clusters of the legs' matched tracks and on the road + std::vector trdTracklets; ///< TRD tracklets of the legs' matched tracks and on the road + float timeTOFMUS = -1.f; ///< time of the cosmic [mus since the TF start, negative before it] from its HitTOFFlight pair (valid if hasTOFTime()) + float scoreTOFPair = -1.f; ///< score of the HitTOFFlight pair: road and flight-time residuals squared in units of the cuts (< 0: none) + float scoreTOFReversed = -1.f; ///< same for the best pair in the impossible order (bottom hit first) = an accidental coincidence: + ///< QA of the flag's background; comparable to scoreTOFPair = ambiguous TOF time (< 0: none) + int duplicateOf = -1; ///< entry (in the TF) of the cosmic this one duplicates: the same muon matched twice (e.g. a leg split + ///< into two TPC tracks): >= 30 % of the smaller one's TPC clusters shared with a better one (TOF time, + ///< more attached clusters) (-1: none) + o2::track::TrackParCov polished{}; ///< refit at the TOF time (TPC-only legs, muon mass, energy loss along the flight), at the closest + ///< approach to the beam line; valid if chi2MatchPolished >= 0 + float chi2MatchPolished = -1.f; ///< chi2 of the polished bottom and top halves at the closest approach (< 0: not polished) + o2::MCCompLabel label{}; ///< MC label of the cosmic (MC only) + + bool hasTOFTime() const { return scoreTOFPair >= 0.f; } ///< a HitTOFFlight pair gave the cosmic its time timeTOFMUS + ClassDefNV(CosmicTrack, 1); +}; + +/// per-TF quantities of the TPC transformation used in the reconstruction +struct CosmicsTFInfo { + float vDrift = 0.f; ///< drift velocity of the transformation [cm/time bin] + float t0 = 0.f; ///< time offset of the transformation [time bins] + ClassDefNV(CosmicsTFInfo, 1); +}; + +} // namespace o2::dataformats + +#endif diff --git a/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h b/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h index d3519bbaad607..82e483febb9f1 100644 --- a/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h +++ b/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h @@ -22,4 +22,16 @@ #pragma link C++ class o2::globaltracking::TrackTuneParams + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::globaltracking::TrackTuneParams> + ; +#pragma link C++ class o2::dataformats::CosmicTPCCluster + ; +#pragma link C++ class std::vector < o2::dataformats::CosmicTPCCluster> + ; +#pragma link C++ class o2::dataformats::CosmicITSCluster + ; +#pragma link C++ class std::vector < o2::dataformats::CosmicITSCluster> + ; +#pragma link C++ class o2::dataformats::CosmicTOFCluster + ; +#pragma link C++ class std::vector < o2::dataformats::CosmicTOFCluster> + ; +#pragma link C++ class o2::dataformats::CosmicTRDTracklet + ; +#pragma link C++ class std::vector < o2::dataformats::CosmicTRDTracklet> + ; +#pragma link C++ class o2::dataformats::CosmicTrack + ; +#pragma link C++ class std::vector < o2::dataformats::CosmicTrack> + ; +#pragma link C++ class o2::dataformats::CosmicsTFInfo + ; + #endif diff --git a/Detectors/GlobalTrackingWorkflow/CMakeLists.txt b/Detectors/GlobalTrackingWorkflow/CMakeLists.txt index 6f29ab1930f95..69f5ff9379944 100644 --- a/Detectors/GlobalTrackingWorkflow/CMakeLists.txt +++ b/Detectors/GlobalTrackingWorkflow/CMakeLists.txt @@ -26,6 +26,7 @@ o2_add_library(GlobalTrackingWorkflow src/StrangenessTrackingWriterSpec.cxx src/CosmicsMatchingSpec.cxx src/TrackCosmicsWriterSpec.cxx + src/CosmicsClusterCollectorSpec.cxx src/TOFMatcherSpec.cxx src/TOFMatchChecker.cxx src/HMPMatcherSpec.cxx diff --git a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h new file mode 100644 index 0000000000000..2db4d81c2293e --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h @@ -0,0 +1,31 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @file CosmicsClusterCollectorSpec.h +/// @brief Collects the raw clusters of matched cosmics (attached + road around the legs) for offline refits + +#ifndef O2_COSMICS_CLUSTER_COLLECTOR_SPEC_H +#define O2_COSMICS_CLUSTER_COLLECTOR_SPEC_H + +#include "Framework/DataProcessorSpec.h" +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "DetectorsCommonDataFormats/DetID.h" + +namespace o2::globaltracking +{ + +/// create a processor spec collecting the clusters of the cosmics found by the cosmics-matcher +/// roadDets: detectors (ITS, TOF, TRD) whose hits along the cosmic are collected besides the TPC road +framework::DataProcessorSpec getCosmicsClusterCollectorSpec(o2::dataformats::GlobalTrackID::mask_t src, bool useMC, bool itsStag, o2::detectors::DetID::mask_t roadDets); + +} // namespace o2::globaltracking + +#endif diff --git a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TrackCosmicsWriterSpec.h b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TrackCosmicsWriterSpec.h index 9cae3afa7b1d5..fd149b7c3c401 100644 --- a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TrackCosmicsWriterSpec.h +++ b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/TrackCosmicsWriterSpec.h @@ -27,6 +27,9 @@ namespace globaltracking /// write cosmics tracks to a root file framework::DataProcessorSpec getTrackCosmicsWriterSpec(bool useMC); +/// write cosmics with their raw clusters (cosmics-cluster-collector output) to a root file +framework::DataProcessorSpec getCosmicsFullWriterSpec(); + } // namespace globaltracking } // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx new file mode 100644 index 0000000000000..22b57fb59b0be --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx @@ -0,0 +1,1641 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @file CosmicsClusterCollectorSpec.cxx +/// @brief Collects the raw clusters of matched cosmics (attached + road around the legs) for offline refits +/// +/// The async reconstruction stores neither TPC tracks nor TPC clusters, so the leg references of TrackCosmics cannot be resolved +/// offline. For every cosmic this device stores the TPC tracks of both legs, the raw TPC clusters attached to them, all raw TPC clusters +/// in a road around each leg (gaps, split pieces, delta electrons; clusters attached to other tracks are flagged), the ITS / TOF / TRD +/// hits of the legs' matched tracks and along roads in these detectors, and per TF the quantities of the TPC transformation. +/// +/// Road: each leg's helix is propagated through all pad rows of all sectors on its TPC side, in the frame of its own time0 (the frame +/// in which its clusters are consistent). Points on the other side of the leg's closest approach to the beam line belong to the other +/// leg and are skipped. The predicted real (y, z) is mapped to nominal coordinates with the inverse correction, and clusters with +/// nominal coordinates within the corridor width of the predicted point (distance perpendicular to the track) are taken. The other TPC +/// side is searched only when the time of the cosmic is known (z-continuity time of legs on opposite sides with a small error, or a TOF +/// time); a CE-crossing leg's own frame already covers both sides. Parts of a leg nearly parallel to the pad rows (|snp| >= 0.8 in the +/// sector frame, e.g. around a closest approach to the beam line inside the TPC, where the track runs along one pad row) cannot be +/// reached row by row: there the helix is followed in steps of path length and the rows within the road width of each point are searched. +/// +/// TOF tag: all TOF clusters close to the outward extrapolation of the legs within the cosmic's time window are candidates; the top / +/// bottom pair whose time difference matches the muon's flight along the helix (HitTOFFlight) gives the cosmic its TOF time. The same +/// search in the impossible order (bottom hit first) only finds accidental pairs and is kept as QA of the flag's background. +/// +/// Roads in the other detectors (--road-detectors): TRD tracklets close to the outward extrapolation of the legs and ITS clusters close +/// to the trajectory near the beam line, within the time window of the cosmic (the TOF time if there is one, else the matcher's time, +/// which is precise for TPC-only legs on opposite sides; otherwise the time window of the legs, with a correspondingly loose z cut). +/// They are flagged as found on the road; hits of the legs' matched global tracks are flagged as matched. Cosmics sharing >= 30 % of +/// their TPC clusters with a better one (a leg split into two TPC tracks) are flagged as duplicates, nothing is removed. +/// +/// Polish: a cosmic with a TOF time and TPC-only legs is refitted at that time (as in the matcher: muon mass, energy loss along the +/// flight); one-side legs then have their real z, hence the right material, which the matcher's TPC time cannot always give. +/// +/// Besides the matcher's and the polished track, only raw detector data are stored (TPC ClusterNative, ITS compact clusters + patterns, +/// raw TOF time, TRD tracklet words); calibrations, the cluster dictionary and the geometry are applied offline. +/// +/// --debug-tree writes cosmics_collector_debug.root for test runs: tree "cosmics" with one entry per cosmic, its matching and timing +/// quantities and vectors per detector: cl* every stored TPC cluster transformed (local x, y, z in the frame used by the road, zCos in the +/// common frame of the cosmic's time, global gx, gy), road* the predicted track points of each leg (same frames; roadFrame 0 / 1 own +/// / cosmic's time frame row by row, 2 / 3 the same along the low-angle walk), its* the +/// local (with the ITS road also global) coordinates of the ITS clusters, tof* the raw and calibrated TOF time and the cluster position, +/// trd* the road tracklets and their residuals. E.g. Draw("clGy:clGx", "scoreTOF >= 0") (TOF-tagged; tTOF can be negative). + +#include +#include +#include +#include +#include +#include +#include +#include +#include "TStopwatch.h" +#include "Framework/Task.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/DeviceSpec.h" +#include "GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h" +#include "DataFormatsGlobalTracking/RecoContainer.h" +#include "DataFormatsGlobalTracking/CosmicTrack.h" +#include "ReconstructionDataFormats/TrackCosmics.h" +#include "DataFormatsTPC/TrackTPC.h" +#include "DataFormatsTPC/ClusterNative.h" +#include "DataFormatsITS/TrackITS.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DataFormatsITSMFT/TopologyDictionary.h" +#include "DataFormatsITSMFT/ClusterPattern.h" +#include "ITSMFTBase/SegmentationAlpide.h" +#include "DataFormatsTOF/Cluster.h" +#include "DataFormatsTRD/TrackTRD.h" +#include "DataFormatsTRD/Tracklet64.h" +#include "DataFormatsTRD/TriggerRecord.h" +#include "DataFormatsTRD/CalibratedTracklet.h" +#include "DataFormatsTRD/Constants.h" +#include "DataFormatsITSMFT/DPLAlpideParam.h" +#include "ITSBase/GeometryTGeo.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "CommonConstants/LHCConstants.h" +#include "CommonDataFormat/TFIDInfo.h" +#include "CommonDataFormat/InteractionRecord.h" +#include "DetectorsBase/GRPGeomHelper.h" +#include "DetectorsBase/Propagator.h" +#include "GPUO2InterfaceRefit.h" +#include "GlobalTracking/MatchCosmicsParams.h" +#include "DataFormatsTPC/WorkflowHelper.h" +#include "ReconstructionDataFormats/Vertex.h" +#include "DetectorsBase/TFIDInfoHelper.h" +#include "TPCBase/ParameterElectronics.h" +#include "MathUtils/Utils.h" +#include "MathUtils/Primitive2D.h" +#include "TPCFastTransformPOD.h" +#include "CommonUtils/TreeStreamRedirector.h" + +using namespace o2::framework; +using GTrackID = o2::dataformats::GlobalTrackID; +using TPCGeo = o2::gpu::TPCFastTransformGeoPOD; +using DetID = o2::detectors::DetID; + +namespace o2::globaltracking +{ + +namespace +{ +constexpr float TanSector = 0.17632698f; // tan(10 deg): half opening of a sector + +/// z outside the drift volume of one TPC side by more than margin +bool outsideDriftVolume(bool sideA, float z, float margin) +{ + const float zLength = TPCGeo::getTPCzLength(); + return sideA ? (z < -margin || z > zLength + margin) : (z > margin || z < -zLength - margin); +} + +/// TPC side of a track with clusters on one side only (+1 A, -1 C: its z moves with the time assumed for its clusters), 0 otherwise +int tpcSide(const o2::tpc::TrackTPC& trk) +{ + return trk.hasASideClustersOnly() ? 1 : (trk.hasCSideClustersOnly() ? -1 : 0); +} +} // namespace + +class CosmicsClusterCollectorSpec : public Task +{ + public: + CosmicsClusterCollectorSpec(std::shared_ptr dr, std::shared_ptr gr, bool useMC, DetID::mask_t roadDets) : mDataRequest(dr), mGGCCDBRequest(gr), mRoadDets(roadDets), mUseMC(useMC) {} + ~CosmicsClusterCollectorSpec() override = default; + void init(InitContext& ic) final; + void run(ProcessingContext& pc) final; + void endOfStream(EndOfStreamContext& ec) final; + void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final; + + private: + /// which side of the closest approach to the beam line (transverse) a point is on, relative to the leg's own clusters + struct LegBranch { + bool isLine = false; ///< straight line (no field / very high pT) instead of a circle + float centerX = 0.f; ///< circle centre + float centerY = 0.f; + float pcaX = 0.f; ///< point of closest approach to the beam line + float pcaY = 0.f; + float dirX = 0.f; ///< direction of the straight line + float dirY = 0.f; + int sign = 0; ///< side of the leg's clusters; 0: the leg spans both sides, accept everything + int side(float x, float y) const + { + const float orientation = isLine ? (x - pcaX) * dirX + (y - pcaY) * dirY : (pcaX - centerX) * (y - centerY) - (pcaY - centerY) * (x - centerX); + return orientation > 0.f ? 1 : -1; + } + bool accept(float x, float y) const { return sign == 0 || side(x, y) == sign; } + void init(const o2::track::TrackPar& inner, const o2::track::TrackPar& outer, float bz); + }; + + /// time of the cosmic in TPC time bins + struct CosmicTime { + bool known = false; ///< known well enough to search the other TPC side of one-side legs + float tb = 0.f; ///< time [TB] + float errTB = 0.f; ///< its error [TB] + }; + + void updateTimeDependentParams(ProcessingContext& pc); + void buildUsedMap(const RecoContainer& data, std::vector>& time0Windows); + void addTPCAttached(const RecoContainer& data, const o2::tpc::TrackTPC& trk, std::vector& out, std::unordered_set& taken) const; + /// time frame of a TPC road: the leg's z is shifted by dz and its clusters are transformed with vertexTime + struct RoadFrame { + float vertexTime; ///< vertex time used for the transformation [TB] + float dz; ///< shift of the leg's z into this frame + float zTolerance; ///< extra z tolerance from the error of the vertex time + int sectorMin; + int sectorMax; + uint8_t flag; + }; + void addTPCCorridor(const RecoContainer& data, const o2::tpc::TrackTPC& trk, const CosmicTime& cosmicTime, std::vector& out, std::unordered_set& taken) const; + void walkLowAngleRoad(const o2::tpc::ClusterNativeAccess& clusters, const o2::track::TrackPar& start, bool innerPart, float rMid, const LegBranch& branch, const RoadFrame& frame, + std::vector& out, std::unordered_set& taken) const; + void searchRow(const o2::tpc::ClusterNativeAccess& clusters, int sector, int row, float y, float z, float snp, float tgl, float vertexTime, float zTolerance, uint8_t flag, + float dxRow, std::vector& out, std::unordered_set& taken) const; + void addDebugRoadPoint(int sector, float x, float y, float z, float snp, float tgl, int row, float vertexTime, int frameCode) const; + struct ITSPattRequest { + int cosmic; ///< entry of the cosmic in the output + int cluster; ///< entry of the cluster in its clITS + int index; ///< index of the cluster in the TF + }; + void addITS(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, int icosm, std::vector& requests, std::unordered_set& matched) const; + void fillITSPatterns(const RecoContainer& data, std::vector& requests, std::vector& cosmics) const; + void addTOF(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, std::unordered_set& matched) const; + void addTRD(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, std::unordered_set& matched) const; + std::pair timeWindowMUS(const CosmicTime& cosmicTime) const; + bool predictOutward(const o2::tpc::TrackTPC& leg, const CosmicTime& cosmicTime, int sector, float x, float& y, float& z) const; + void roadTOF(const RecoContainer& data, const o2::tpc::TrackTPC* const* legs, const CosmicTime& cosmicTime, int icosm, std::vector& out, const std::unordered_set& matched, float& timeTOFMUS, float& scoreTOFPair, float& scoreTOFReversed) const; + void polish(const o2::tpc::TrackTPC* const* legs, float timeMUS, o2::gpu::GPUO2InterfaceRefit& refitter, o2::dataformats::CosmicTrack& out) const; + void roadTRD(const RecoContainer& data, const o2::tpc::TrackTPC* const* legs, const CosmicTime& cosmicTime, int icosm, std::vector& out, const std::unordered_set& matched) const; + void roadITS(const RecoContainer& data, const o2::dataformats::TrackCosmics& cosm, const CosmicTime& cosmicTime, int legsSide, int icosm, std::vector& out, std::vector& requests, + const std::unordered_set& matched) const; + void cacheITSChipCentres(); + void writeDebug(const o2::dataformats::CosmicTrack& cosm, int icosm) const; + void flagDuplicates(std::vector& cosmics) const; + void writeDebugTOF(const o2::tof::Cluster& c, int icosm, int leg, uint8_t flags) const; + + std::shared_ptr mDataRequest; + std::shared_ptr mGGCCDBRequest; + const o2::gpu::TPCFastTransformPOD* mCorrMap = nullptr; + const o2::itsmft::TopologyDictionary* mITSDict = nullptr; + std::vector mUsed; ///< TPC clusters attached to any TPC track in this TF + o2::InteractionRecord mTFStart{}; ///< first BC of the TF + float mCorridor = 1.f; ///< road half-width [cm] + float mMaxAbsTimeErr = 0.5f; ///< max. time error of the cosmic [mus] to search the other TPC side of a one-side leg + size_t mMaxCosmicsPerTF = 100; ///< cosmics processed per TF at most (protection against fake-dominated settings) + DetID::mask_t mRoadDets{}; ///< detectors searched along the road besides the TPC + float mRoadTOF = 5.f; ///< road half-width at the TOF [cm] + float mTOFFlightTol = 2.f; ///< max. deviation of the top/bottom TOF time difference from the flight time [ns] + float mTOFTimeErr = 0.1f; ///< error of the TOF time of a cosmic for the later roads [mus] (covers TPC vs TOF offsets) + float mRoadTRD = 5.f; ///< road half-width at the TRD [cm] (in z plus the pad length) + float mRoadITS = 1.5f; ///< road half-width in the ITS [cm] + std::vector> mITSChipCentres; ///< global positions of the ITS chip centres (aligned geometry) + float mTPCTBinMUS = 0.2f; ///< TPC time bin [mus] + float mBz = 0.f; + bool mUseMC = false; + mutable bool mStaggeredWarned = false; + mutable bool mNoDictWarned = false; + mutable bool mTRDMismatchWarned = false; + std::unique_ptr mDebugOut; ///< debug tree "cosmics", one entry per cosmic (--debug-tree) + int mDbgTF = 0; ///< context of the debug output: TF counter + int mDbgCosmic = 0; ///< entry of the cosmic in the TF + int mDbgLeg = 0; ///< leg (0 bottom, 1 top) + float mDbgTauC = 0.f; ///< time of the cosmic [TB]: common frame (zCos) of the debug output + /// debug quantities found while processing a cosmic (road points, TOF and TRD hits), written with its entry + struct DebugCosmic { + std::vector roadLeg; + std::vector roadFrame; + std::vector roadSector; + std::vector roadRow; + std::vector roadX; + std::vector roadY; + std::vector roadZ; + std::vector roadZCos; + std::vector roadGx; + std::vector roadGy; + std::vector roadSnp; + std::vector roadTgl; + std::vector tofLeg; + std::vector tofChannel; + std::vector tofFlags; + std::vector tofTimeRaw; + std::vector tofTime; + std::vector tofTot; + std::vector tofX; + std::vector tofY; + std::vector tofZ; + std::vector tofGx; + std::vector tofGy; + std::vector trdLeg; + std::vector trdLayer; + std::vector trdX; + std::vector trdY; + std::vector trdZ; + std::vector trdDy; + std::vector trdDz; + std::vector trdTrigMUS; + }; + mutable std::vector mDbgCosmics; ///< per cosmic of the TF + size_t mNCosmics = 0; ///< cosmics processed + size_t mNClAttached = 0; ///< attached TPC clusters stored + size_t mNClCorridor = 0; ///< road TPC clusters stored + TStopwatch mTimer; +}; + +void CosmicsClusterCollectorSpec::init(InitContext& ic) +{ + mTimer.Stop(); + mTimer.Reset(); + o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); + mCorridor = ic.options().get("corridor-width"); + mMaxAbsTimeErr = ic.options().get("max-abs-time-err"); + mMaxCosmicsPerTF = std::max(0, ic.options().get("max-cosmics-per-tf")); + mRoadTOF = ic.options().get("tof-road-width"); + mTOFFlightTol = ic.options().get("tof-flight-tolerance"); + mTOFTimeErr = ic.options().get("tof-time-error"); + mRoadTRD = ic.options().get("trd-road-width"); + mRoadITS = ic.options().get("its-road-width"); + if (ic.options().get("debug-tree")) { + const auto timesliceId = ic.services().get().inputTimesliceId; + const std::string name = timesliceId == 0 ? "cosmics_collector_debug.root" : fmt::format("cosmics_collector_debug_{}.root", timesliceId); + mDebugOut = std::make_unique(name.c_str(), "recreate"); + } +} + +void CosmicsClusterCollectorSpec::run(ProcessingContext& pc) +{ + mTimer.Start(false); + RecoContainer recoData; + recoData.collectData(pc, *mDataRequest.get()); + updateTimeDependentParams(pc); + + o2::dataformats::TFIDInfo tfID; + o2::base::TFIDInfoHelper::fillTFIDInfo(pc, tfID); + mTFStart = {0, tfID.firstTForbit}; + o2::dataformats::CosmicsTFInfo tfInfo; + tfInfo.vDrift = mCorrMap->getVDrift(); + tfInfo.t0 = mCorrMap->getT0(); + + std::vector cosmicsOut; + const auto cosmics = recoData.getCosmicTracks(); + const size_t nCosmics = std::min(cosmics.size(), mMaxCosmicsPerTF); + if (nCosmics < cosmics.size()) { + LOGP(warning, "{} cosmics in TF {}, only the first {} are written (max-cosmics-per-tf)", cosmics.size(), tfID.tfCounter, nCosmics); + } + if (nCosmics) { + // only TPC tracks whose time0 is within two drift times of a leg's time0 can share clusters with the roads + const float maxDistTB = 2.2f * TPCGeo::getTPCzLength() / mCorrMap->getVDrift(); + std::vector> time0Windows; + for (size_t ic = 0; ic < nCosmics; ic++) { + for (const auto leg : {cosmics[ic].getRefBottom(), cosmics[ic].getRefTop()}) { + const auto refs = recoData.getSingleDetectorRefs(leg); + if (refs[GTrackID::TPC].isIndexSet()) { + const float time0 = recoData.getTPCTrack(refs[GTrackID::TPC]).getTime0(); + time0Windows.emplace_back(time0 - maxDistTB, time0 + maxDistTB); + } + } + } + buildUsedMap(recoData, time0Windows); + } + std::vector pattRequests; // ITS clusters whose patterns are not in the dictionary + std::unique_ptr tpcRefitter; // polish of the cosmics with a TOF time + if (nCosmics && mRoadDets[DetID::TOF] && recoData.inputsTPCclusters) { + const auto& clusterShMap = recoData.clusterShMapTPC; + const auto& occupancyMap = recoData.occupancyMapTPC; + tpcRefitter = std::make_unique(&recoData.inputsTPCclusters->clusterIndex, mCorrMap, mBz, recoData.getTPCTracksClusterRefs().data(), 0, + clusterShMap.data(), occupancyMap.data(), occupancyMap.size(), nullptr, o2::base::Propagator::Instance()); + } + mDbgTF = tfID.tfCounter; + if (mDebugOut) { + mDbgCosmics.assign(nCosmics, DebugCosmic{}); + } + for (size_t ic = 0; ic < nCosmics; ic++) { + const auto& cosm = cosmics[ic]; + auto& out = cosmicsOut.emplace_back(); + out.cosmic = cosm; + if (mUseMC) { + out.label = recoData.getCosmicTrackMCLabel(ic); + } + const GTrackID legs[2] = {cosm.getRefBottom(), cosm.getRefTop()}; + const o2::tpc::TrackTPC* tpcLegs[2] = {nullptr, nullptr}; + std::vector* tpcCl[2] = {&out.clTPCBottom, &out.clTPCTop}; + std::unordered_set matchedITS; // hits of the legs' matched tracks, not searched again on the road + std::unordered_set matchedTOF; + std::unordered_set matchedTRD; + for (uint8_t leg = 0; leg < 2; leg++) { + auto refs = recoData.getSingleDetectorRefs(legs[leg]); + if (refs[GTrackID::TPC].isIndexSet()) { + tpcLegs[leg] = &recoData.getTPCTrack(refs[GTrackID::TPC]); + (leg == 0 ? out.tpcBottom : out.tpcTop) = *tpcLegs[leg]; + } + if (refs[GTrackID::ITS].isIndexSet()) { + addITS(recoData, refs[GTrackID::ITS], leg, out.clITS, ic, pattRequests, matchedITS); + } + if (refs[GTrackID::TOF].isIndexSet()) { + addTOF(recoData, refs[GTrackID::TOF], leg, out.clTOF, matchedTOF); + if (mDebugOut) { + writeDebugTOF(recoData.getTOFClusters()[refs[GTrackID::TOF].getIndex()], ic, leg, o2::dataformats::HitMatched); + } + } + if (refs[GTrackID::TRD].isIndexSet()) { + addTRD(recoData, refs[GTrackID::TRD], leg, out.trdTracklets, matchedTRD); + } + } + CosmicTime cosmicTime; + cosmicTime.tb = cosm.getTimeMUS().getTimeStamp() / mTPCTBinMUS; + cosmicTime.errTB = cosm.getTimeMUS().getTimeStampError() / mTPCTBinMUS; + cosmicTime.known = cosm.getTimeMUS().getTimeStampError() < mMaxAbsTimeErr; + mDbgCosmic = ic; + // TOF road first: a top/bottom hit pair matching the muon's flight gives the cosmic's time to ~ns (also for one-side legs on the same + // side, whose brackets leave it open by tens of mus); the TPC corridor of the other side and the TRD / ITS roads then use that time + if (mRoadDets[DetID::TOF]) { + roadTOF(recoData, tpcLegs, cosmicTime, ic, out.clTOF, matchedTOF, out.timeTOFMUS, out.scoreTOFPair, out.scoreTOFReversed); + } + if (out.hasTOFTime() && tpcRefitter && legs[0].getSource() == GTrackID::TPC && legs[1].getSource() == GTrackID::TPC) { + polish(tpcLegs, out.timeTOFMUS, *tpcRefitter, out); + } + CosmicTime roadTime = cosmicTime; + if (out.hasTOFTime()) { + roadTime.tb = out.timeTOFMUS / mTPCTBinMUS; + roadTime.errTB = mTOFTimeErr / mTPCTBinMUS; + roadTime.known = true; + } + mDbgTauC = roadTime.tb; + std::unordered_set taken; + for (int leg = 0; leg < 2; leg++) { // attached clusters of both legs first, so that a road never takes the other leg's clusters + if (tpcLegs[leg]) { + addTPCAttached(recoData, *tpcLegs[leg], *tpcCl[leg], taken); + mNClAttached += tpcCl[leg]->size(); + } + } + for (int leg = 0; leg < 2; leg++) { + if (tpcLegs[leg]) { + const size_t nAttached = tpcCl[leg]->size(); + mDbgLeg = leg; + addTPCCorridor(recoData, *tpcLegs[leg], roadTime, *tpcCl[leg], taken); + mNClCorridor += tpcCl[leg]->size() - nAttached; + } + } + if (mRoadDets[DetID::TRD]) { + roadTRD(recoData, tpcLegs, roadTime, ic, out.trdTracklets, matchedTRD); + } + if (mRoadDets[DetID::ITS]) { + // the refitted cosmic's z moves with the time only if both legs are TPC-only on one side (an ITS part fixes it absolutely) + int legsSide = 0; + if (tpcLegs[0] && tpcLegs[1] && legs[0].getSource() == GTrackID::TPC && legs[1].getSource() == GTrackID::TPC && tpcSide(*tpcLegs[0]) == tpcSide(*tpcLegs[1])) { + legsSide = tpcSide(*tpcLegs[0]); + } + roadITS(recoData, cosm, roadTime, legsSide, ic, out.clITS, pattRequests, matchedITS); + } + } + if (!pattRequests.empty()) { + fillITSPatterns(recoData, pattRequests, cosmicsOut); + } + flagDuplicates(cosmicsOut); + if (mDebugOut) { + for (size_t ic = 0; ic < cosmicsOut.size(); ic++) { + writeDebug(cosmicsOut[ic], ic); + } + } + mNCosmics += cosmicsOut.size(); + LOGP(info, "Collected clusters for {} cosmics in TF {}", cosmicsOut.size(), tfID.tfCounter); + pc.outputs().snapshot(Output{"GLO", "COSMFULL", 0}, cosmicsOut); + pc.outputs().snapshot(Output{"GLO", "COSMFULLTF", 0}, tfInfo); + pc.outputs().snapshot(Output{"GLO", "COSMFULLTFID", 0}, tfID); + mTimer.Stop(); +} + +void CosmicsClusterCollectorSpec::updateTimeDependentParams(ProcessingContext& pc) +{ + o2::base::GRPGeomHelper::instance().checkUpdates(pc); + mCorrMap = &o2::gpu::TPCFastTransformPOD::get(pc.inputs().get("corrMap")); + mTPCTBinMUS = o2::tpc::ParameterElectronics::Instance().ZbinWidth; + mBz = o2::base::Propagator::Instance()->getNominalBz(); + if (mRoadDets[DetID::ITS] && mITSChipCentres.empty()) { + cacheITSChipCentres(); + } +} + +void CosmicsClusterCollectorSpec::cacheITSChipCentres() +{ + auto geom = o2::its::GeometryTGeo::Instance(); + geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + mITSChipCentres.resize(geom->getNumberOfChips()); + for (int chip = 0; chip < geom->getNumberOfChips(); chip++) { + mITSChipCentres[chip] = geom->getMatrixL2G(chip) * o2::math_utils::Point3D(0.f, 0.f, 0.f); + } +} + +void CosmicsClusterCollectorSpec::buildUsedMap(const RecoContainer& data, std::vector>& time0Windows) +{ + // merge the windows, then flag the clusters of the TPC tracks whose time0 falls into one of them + std::sort(time0Windows.begin(), time0Windows.end()); + std::vector> merged; + for (const auto& window : time0Windows) { + if (!merged.empty() && window.first <= merged.back().second) { + merged.back().second = std::max(merged.back().second, window.second); + } else { + merged.push_back(window); + } + } + const auto& clusters = data.getTPCClusters(); + const auto tracks = data.getTPCTracks(); + const auto refs = data.getTPCTracksClusterRefs(); + mUsed.assign(clusters.nClustersTotal, false); + for (const auto& trk : tracks) { + const float time0 = trk.getTime0(); + const auto window = std::upper_bound(merged.begin(), merged.end(), time0, [](float t, const std::pair& w) { return t < w.first; }); + if (window == merged.begin() || time0 > std::prev(window)->second) { + continue; + } + for (int j = 0; j < trk.getNClusterReferences(); j++) { + uint8_t sector = 0; + uint8_t row = 0; + uint32_t clIdx = 0; + trk.getClusterReference(refs, j, sector, row, clIdx); + mUsed[clusters.clusterOffset[sector][row] + clIdx] = true; + } + } +} + +void CosmicsClusterCollectorSpec::addTPCAttached(const RecoContainer& data, const o2::tpc::TrackTPC& trk, std::vector& out, std::unordered_set& taken) const +{ + const auto& clusters = data.getTPCClusters(); + const auto refs = data.getTPCTracksClusterRefs(); + for (int j = 0; j < trk.getNClusterReferences(); j++) { + uint8_t sector = 0; + uint8_t row = 0; + uint32_t clIdx = 0; + trk.getClusterReference(refs, j, sector, row, clIdx); + if (!taken.insert(clusters.clusterOffset[sector][row] + clIdx).second) { + continue; + } + auto& cl = out.emplace_back(); + cl.cl = clusters.clusters[sector][row][clIdx]; + cl.sector = sector; + cl.row = row; + cl.flags = o2::dataformats::CosmicTPCCluster::Attached | o2::dataformats::CosmicTPCCluster::Used; + } +} + +void CosmicsClusterCollectorSpec::LegBranch::init(const o2::track::TrackPar& inner, const o2::track::TrackPar& outer, float bz) +{ + if (std::abs(inner.getCurvature(bz)) < 1e-5f) { // radius > 1 km: straight line + isLine = true; + const auto point = inner.getXYZGlo(); + const float phi = inner.getPhi(); + dirX = std::cos(phi); + dirY = std::sin(phi); + const float proj = point.X() * dirX + point.Y() * dirY; + pcaX = point.X() - proj * dirX; + pcaY = point.Y() - proj * dirY; + } else { + o2::math_utils::CircleXYf_t circle; + float sinAlpha = 0.f; + float cosAlpha = 0.f; + inner.getCircleParams(bz, circle, sinAlpha, cosAlpha); + centerX = circle.xC; + centerY = circle.yC; + const float centerDist = std::sqrt(centerX * centerX + centerY * centerY); + if (centerDist < 1e-3f) { // circle around the beam line: no closest approach + sign = 0; + return; + } + pcaX = centerX * (1.f - circle.rC / centerDist); + pcaY = centerY * (1.f - circle.rC / centerDist); + } + // side of the leg from its end farther from the closest approach; a leg with ends on both sides (both > 10 cm away) spans the + // closest approach and is accepted on both sides + const auto pointIn = inner.getXYZGlo(); + const auto pointOut = outer.getXYZGlo(); + const float distIn2 = (pointIn.X() - pcaX) * (pointIn.X() - pcaX) + (pointIn.Y() - pcaY) * (pointIn.Y() - pcaY); + const float distOut2 = (pointOut.X() - pcaX) * (pointOut.X() - pcaX) + (pointOut.Y() - pcaY) * (pointOut.Y() - pcaY); + const int sideIn = side(pointIn.X(), pointIn.Y()); + const int sideOut = side(pointOut.X(), pointOut.Y()); + constexpr float MinDist2 = 10.f * 10.f; + if (sideIn == sideOut) { + sign = sideIn; + } else if (std::min(distIn2, distOut2) < MinDist2) { + sign = distIn2 > distOut2 ? sideIn : sideOut; + } else { + sign = 0; + } +} + +void CosmicsClusterCollectorSpec::addTPCCorridor(const RecoContainer& data, const o2::tpc::TrackTPC& trk, const CosmicTime& cosmicTime, std::vector& out, std::unordered_set& taken) const +{ + constexpr int NSectorsA = TPCGeo::getNumberOfSectorsA(); + const float zLength = TPCGeo::getTPCzLength(); + const float vDrift = mCorrMap->getVDrift(); + const auto& clusters = data.getTPCClusters(); + const float time0Leg = trk.getTime0(); + const int legSide = tpcSide(trk); + + RoadFrame frames[2]; + int nFrames = 0; + if (legSide == 0) { // CE-crossing leg: its time0 is absolute + frames[nFrames++] = {time0Leg, 0.f, 0.f, 0, 2 * NSectorsA, 0}; + } else { + frames[nFrames++] = {time0Leg, 0.f, 0.f, legSide > 0 ? 0 : NSectorsA, legSide > 0 ? NSectorsA : 2 * NSectorsA, 0}; + if (cosmicTime.known) { // z of a one-side leg moves by side * vD * (t - time0) when its clusters are transformed with the time t + frames[nFrames++] = {cosmicTime.tb, legSide * (cosmicTime.tb - time0Leg) * vDrift, cosmicTime.errTB * vDrift, + legSide > 0 ? NSectorsA : 0, legSide > 0 ? 2 * NSectorsA : NSectorsA, o2::dataformats::CosmicTPCCluster::AbsTime}; + } + } + + const o2::track::TrackPar refPar[2] = {trk, trk.getParamOut()}; // rows below rMid: inner parameter, above: outer + LegBranch branch; + branch.init(refPar[0], refPar[1], mBz); + const auto pointIn = refPar[0].getXYZGlo(); + const auto pointOut = refPar[1].getXYZGlo(); + const float rIn = std::hypot(pointIn.X(), pointIn.Y()); + const float rOut = std::hypot(pointOut.X(), pointOut.Y()); + const float rMid = 0.5f * (rIn + rOut); + + for (int iFrame = 0; iFrame < nFrames; iFrame++) { + const auto& frame = frames[iFrame]; + for (int sector = frame.sectorMin; sector < frame.sectorMax; sector++) { + const float alpha = o2::math_utils::sector2Angle(sector % NSectorsA); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + const bool sideA = sector < NSectorsA; + for (int iPar = 0; iPar < 2; iPar++) { + auto par = refPar[iPar]; + par.setZ(par.getZ() + frame.dz); + if (!par.rotateParam(alpha)) { // the leg points away from this sector frame: same helix, opposite direction + par.invertParam(); + if (!par.rotateParam(alpha)) { + continue; + } + } + for (int row = 0; row < o2::tpc::constants::MAXGLOBALPADROW; row++) { + const float xRow = TPCGeo::getRowInfoX(row); + if ((xRow < rMid) != (iPar == 0)) { + continue; + } + auto parRow = par; + if (!parRow.propagateParamTo(xRow, mBz)) { + continue; + } + // coarse acceptance at the nominal x of the row before the more expensive move to its real x + constexpr float CoarseMargin = 10.f; // [cm] covers the shift of y and z between the nominal and the real x + if (std::abs(parRow.getY()) > xRow * TanSector + mCorridor + CoarseMargin || outsideDriftVolume(sideA, parRow.getZ(), mCorridor + frame.zTolerance + CoarseMargin)) { + continue; + } + // the corrected clusters of this row lie at the real x of the row, x + dx(y, z), not at its nominal x + float xReal = xRow; + const float zInside = sideA ? std::clamp(parRow.getZ(), 0.f, zLength) : std::clamp(parRow.getZ(), -zLength, 0.f); + mCorrMap->InverseTransformYZtoX(sector, row, parRow.getY(), zInside, xReal); + if (!parRow.propagateParamTo(xReal, mBz)) { + continue; + } + const float y = parRow.getY(); + const float z = parRow.getZ(); + if (std::abs(y) > xRow * TanSector + mCorridor) { + continue; + } + if (outsideDriftVolume(sideA, z, mCorridor + frame.zTolerance)) { + continue; + } + if (!branch.accept(xRow * cosAlpha - y * sinAlpha, xRow * sinAlpha + y * cosAlpha)) { + continue; + } + if (mDebugOut) { + addDebugRoadPoint(sector, xReal, y, z, parRow.getSnp(), parRow.getTgl(), row, frame.vertexTime, frame.flag != 0); + } + searchRow(clusters, sector, row, y, z, parRow.getSnp(), parRow.getTgl(), frame.vertexTime, frame.zTolerance, frame.flag, 0.f, out, taken); + } + } + } + // parts of the leg nearly parallel to the pad rows, which the propagation row by row cannot reach (e.g. a closest approach to the + // beam line inside the TPC: in the sector of the closest approach the track runs along one pad row); a leg whose inner end is not + // inside its outer end (all its clusters along one pad row) is walked from the inner end alone + const bool splitAtRMid = rOut > rIn; + for (int iPar = 0; iPar < (splitAtRMid ? 2 : 1); iPar++) { + auto par = refPar[iPar]; + par.setZ(par.getZ() + frame.dz); + walkLowAngleRoad(clusters, par, iPar == 0, splitAtRMid ? rMid : std::numeric_limits::max(), branch, frame, out, taken); + } + } +} + +void CosmicsClusterCollectorSpec::walkLowAngleRoad(const o2::tpc::ClusterNativeAccess& clusters, const o2::track::TrackPar& start, bool innerPart, float rMid, const LegBranch& branch, + const RoadFrame& frame, std::vector& out, std::unordered_set& taken) const +{ + constexpr int NSectorsA = TPCGeo::getNumberOfSectorsA(); + constexpr float CosSector = 0.98480775f; // cos(10 deg) + constexpr float Step = 0.5f; // [cm] path length per step, below the road width along the pad row + constexpr float MaxPath = 600.f; // [cm] protection; the walk normally ends at the TPC boundary or the closest approach + constexpr float MaxEntryPath = 20.f; // [cm] a start outside the leg's part (up to 10 cm past the closest approach, see LegBranch) is stepped into + constexpr float MinSnp = 0.8f; // below, the row crossings are well defined and found row by row + constexpr float XMargin = 10.f; // [cm] covers the shift between the nominal and the real x of a row (CoarseMargin of the row loop) + constexpr int NRows = o2::tpc::constants::MAXGLOBALPADROW; + const float zLength = TPCGeo::getTPCzLength(); + const float rLow = TPCGeo::getRowInfoX(0) - mCorridor; + const float rHigh = TPCGeo::getRowInfoX(NRows - 1) / CosSector + mCorridor; + std::array sinSector{}; + std::array cosSector{}; + for (int sector = 0; sector < NSectorsA; sector++) { + const float alpha = o2::math_utils::sector2Angle(sector); + sinSector[sector] = std::sin(alpha); + cosSector[sector] = std::cos(alpha); + } + + bool startInside = false; // the start point lies in the leg's part + for (int direction = 0; direction < 2; direction++) { + auto par = start; + if (direction == 1) { + par.invertParam(); + } + bool entered = direction == 1 && startInside; // the walk has reached the leg's part + // the start point is processed in the first direction only + for (float path = direction == 0 ? 0.f : Step; path < MaxPath; path += Step) { + // step along the helix in the frame of its direction (snp = 0), where the propagation is always defined + if (path > 0.f && (!par.rotateParam(par.getPhi()) || !par.propagateParamTo(par.getX() + Step, mBz))) { + break; + } + const auto pos = par.getXYZGlo(); + const float r = std::hypot(pos.X(), pos.Y()); + if (r > rHigh || std::abs(pos.Z()) > zLength + mCorridor + frame.zTolerance) { + break; + } + // outside the leg's part: past its closest approach to the beam line, beyond rMid or inside the inner radius; the start of the walk + // can lie there (a leg's end a few cm past the closest approach, an inner end inside the nominal inner radius) + if (r < rLow || (r < rMid) != innerPart || !branch.accept(pos.X(), pos.Y())) { + if (entered || path > MaxEntryPath) { + break; + } + continue; + } + entered = true; + startInside = startInside || path == 0.f; + const float phiDir = par.getPhi(); + const int sectorPos = o2::math_utils::angle2Sector(std::atan2(pos.Y(), pos.X())); + for (int dSector = -1; dSector <= 1; dSector++) { // the road can reach into the neighbouring sectors + const int sectorInSide = (sectorPos + dSector + NSectorsA) % NSectorsA; + for (int sector = sectorInSide; sector < 2 * NSectorsA; sector += NSectorsA) { + if (sector < frame.sectorMin || sector >= frame.sectorMax) { + continue; + } + const bool sideA = sector < NSectorsA; + if (outsideDriftVolume(sideA, pos.Z(), mCorridor + frame.zTolerance)) { + continue; + } + const float alpha = o2::math_utils::sector2Angle(sectorInSide); + const float sinAlpha = sinSector[sectorInSide]; + const float cosAlpha = cosSector[sectorInSide]; + const float x = pos.X() * cosAlpha + pos.Y() * sinAlpha; + const float y = -pos.X() * sinAlpha + pos.Y() * cosAlpha; + float snp = std::sin(phiDir - alpha); + float tgl = par.getTgl(); + if (std::cos(phiDir - alpha) < 0.f) { // the same line with its direction along +x of the sector frame, as in the row search + snp = -snp; + tgl = -tgl; + } + if (std::abs(snp) < MinSnp || std::abs(y) > x * TanSector + mCorridor) { + continue; + } + const float zInside = sideA ? std::clamp(pos.Z(), 0.f, zLength) : std::clamp(pos.Z(), -zLength, 0.f); + // rows whose real x lies within the road width of the point: the first by bisection on the nominal x, which grows with the row + int rowFirst = 0; + int rowEnd = NRows; + while (rowFirst < rowEnd) { + const int rowMid = (rowFirst + rowEnd) / 2; + if (TPCGeo::getRowInfoX(rowMid) < x - mCorridor - XMargin) { + rowFirst = rowMid + 1; + } else { + rowEnd = rowMid; + } + } + for (int row = rowFirst; row < NRows && TPCGeo::getRowInfoX(row) <= x + mCorridor + XMargin; row++) { + float xReal = TPCGeo::getRowInfoX(row); + mCorrMap->InverseTransformYZtoX(sector, row, y, zInside, xReal); + const float dxRow = xReal - x; + if (std::abs(dxRow) > mCorridor) { + continue; + } + if (mDebugOut) { + addDebugRoadPoint(sector, x, y, pos.Z(), snp, tgl, row, frame.vertexTime, 2 + (frame.flag != 0)); + } + searchRow(clusters, sector, row, y, pos.Z(), snp, tgl, frame.vertexTime, frame.zTolerance, frame.flag, dxRow, out, taken); + } + } + } + } + } +} + +void CosmicsClusterCollectorSpec::searchRow(const o2::tpc::ClusterNativeAccess& clusters, int sector, int row, float y, float z, float snp, float tgl, float vertexTime, float zTolerance, uint8_t flag, + float dxRow, std::vector& out, std::unordered_set& taken) const +{ + // (y, z) is the predicted point, dxRow the real x of the row minus the x of that point (0 when the point lies on the row) + const float zLength = TPCGeo::getTPCzLength(); + const float vDrift = mCorrMap->getVDrift(); + const float t0 = mCorrMap->getT0(); + // nominal (measured) coordinates of the predicted real point; the correction is evaluated inside the drift volume + const float zInside = sector < TPCGeo::getNumberOfSectorsA() ? std::clamp(z, 0.f, zLength) : std::clamp(z, -zLength, 0.f); + float yNominal = 0.f; + float zNominal = 0.f; + mCorrMap->InverseTransformYZtoNominalYZ(sector, row, y, zInside, yNominal, zNominal); + zNominal += z - zInside; + float padPred = 0.f; + float driftLengthPred = 0.f; + TPCGeo::convLocalToPadDriftLength(sector, row, yNominal, zNominal, padPred, driftLengthPred); + const float timePred = driftLengthPred / vDrift + t0 + vertexTime; + // the road is a cylinder around the track; its section with the pad-row plane is an ellipse with half-axes W / cos(phi) in y and + // W * sqrt(cos^2(phi) + tgl^2) / cos(phi) in z + const float cosPhi = std::sqrt((1.f - snp) * (1.f + snp)); + const float cosPhiWindow = std::max(cosPhi, 0.1f); // limits the window for tracks nearly parallel to the pad row + const float norm = 1.f / std::sqrt(1.f + tgl * tgl); + const float dirX = cosPhi * norm; // track direction + const float dirY = snp * norm; + const float dirZ = tgl * norm; + const float corridor2 = mCorridor * mCorridor; + const float windowPad = mCorridor / (cosPhiWindow * TPCGeo::getRowInfoPadWidth(row)); + const float windowTime = (mCorridor * std::sqrt(cosPhiWindow * cosPhiWindow + tgl * tgl) / cosPhiWindow + zTolerance) / vDrift; + const auto* rowClusters = clusters.clusters[sector][row]; + const uint32_t rowOffset = clusters.clusterOffset[sector][row]; + for (uint32_t k = 0; k < clusters.nClusters[sector][row]; k++) { + const auto& c = rowClusters[k]; + const float time = c.getTime(); + const float pad = c.getPad(); + if (std::abs(time - timePred) > windowTime || std::abs(pad - padPred) > windowPad) { + continue; + } + float yCluster = 0.f; + float zCluster = 0.f; + TPCGeo::convPadDriftLengthToLocal(sector, row, pad, (time - t0 - vertexTime) * vDrift, yCluster, zCluster); + const float dy = yCluster - yNominal; + float dz = zCluster - zNominal; + if (zTolerance > 0.f) { // the vertex time of this frame is known within zTolerance / vD: allow that shift along z + dz = std::abs(dz) > zTolerance ? dz - std::copysign(zTolerance, dz) : 0.f; + } + const float proj = dxRow * dirX + dy * dirY + dz * dirZ; + if (dxRow * dxRow + dy * dy + dz * dz - proj * proj > corridor2) { // distance perpendicular to the track + continue; + } + if (!taken.insert(rowOffset + k).second) { + continue; + } + auto& cl = out.emplace_back(); + cl.cl = c; + cl.sector = sector; + cl.row = row; + cl.flags = o2::dataformats::CosmicTPCCluster::Corridor | flag | (mUsed[rowOffset + k] ? o2::dataformats::CosmicTPCCluster::Used : 0); + } +} + +void CosmicsClusterCollectorSpec::addDebugRoadPoint(int sector, float x, float y, float z, float snp, float tgl, int row, float vertexTime, int frameCode) const +{ + // predicted point, z in the frame of the road and in the common frame of the cosmic's time + const float alpha = o2::math_utils::sector2Angle(sector % TPCGeo::getNumberOfSectorsA()); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + const bool sideA = sector < TPCGeo::getNumberOfSectorsA(); + auto& dbg = mDbgCosmics[mDbgCosmic]; + dbg.roadLeg.push_back(mDbgLeg); + dbg.roadFrame.push_back(frameCode); + dbg.roadSector.push_back(sector); + dbg.roadRow.push_back(row); + dbg.roadX.push_back(x); + dbg.roadY.push_back(y); + dbg.roadZ.push_back(z); + dbg.roadZCos.push_back(z + (sideA ? 1.f : -1.f) * (mDbgTauC - vertexTime) * mCorrMap->getVDrift()); + dbg.roadGx.push_back(x * cosAlpha - y * sinAlpha); + dbg.roadGy.push_back(x * sinAlpha + y * cosAlpha); + dbg.roadSnp.push_back(snp); + dbg.roadTgl.push_back(tgl); +} + +void CosmicsClusterCollectorSpec::addITS(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, int icosm, std::vector& requests, std::unordered_set& matched) const +{ + if (gid.getSource() != GTrackID::ITS) { // ITS-AB tracklets are not stored + return; + } + if (data.getITSPerLayer()) { + if (!mStaggeredWarned) { + LOGP(warning, "Staggered ITS clusters are not supported, ITS clusters of cosmics are not stored"); + mStaggeredWarned = true; + } + return; + } + const auto& trk = data.getITSTrack(gid); + const auto refs = data.getITSTracksClusterRefs(); + const auto clusters = data.getITSClusters(); + const auto rofs = data.getITSClustersROFRecords(); + for (int i = 0; i < trk.getNumberOfClusters(); i++) { + const int idx = refs[trk.getFirstClusterEntry() + i]; + const auto& c = clusters[idx]; + auto& cl = out.emplace_back(); + cl.chipID = c.getSensorID(); + cl.row = c.getRow(); + cl.col = c.getCol(); + cl.pattID = c.getPatternID(); + cl.leg = leg; + cl.flags = o2::dataformats::HitMatched; + matched.insert(idx); + if (c.getPatternID() == o2::itsmft::CompCluster::InvalidPatternID || (mITSDict && mITSDict->isGroup(c.getPatternID()))) { + requests.push_back({icosm, int(out.size()) - 1, idx}); // the pattern is in the TF's pattern stream, copied by fillITSPatterns + } + auto rof = std::upper_bound(rofs.begin(), rofs.end(), idx, [](int v, const o2::itsmft::ROFRecord& r) { return v < r.getFirstEntry(); }); + if (rof != rofs.begin()) { + cl.rofBC = std::prev(rof)->getBCData().differenceInBC(mTFStart); + } + } +} + +void CosmicsClusterCollectorSpec::addTOF(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, std::unordered_set& matched) const +{ + const auto& c = data.getTOFClusters()[gid.getIndex()]; + auto& cl = out.emplace_back(); + cl.timeRaw = c.getTimeRaw(); + cl.tot = c.getTot(); + cl.channel = c.getMainContributingChannel(); + cl.leg = leg; + cl.flags = o2::dataformats::HitMatched; + matched.insert(gid.getIndex()); +} + +std::pair CosmicsClusterCollectorSpec::timeWindowMUS(const CosmicTime& cosmicTime) const +{ + // a time fixed by z continuity is given with a 1 sigma error, otherwise the error is the half-width of the legs' time-bracket overlap + const float timeMUS = cosmicTime.tb * mTPCTBinMUS; + const float errMUS = cosmicTime.errTB * mTPCTBinMUS; + const float halfWidth = cosmicTime.known ? 5.f * errMUS + 0.2f : errMUS; + return {timeMUS - halfWidth, timeMUS + halfWidth}; +} + +bool CosmicsClusterCollectorSpec::predictOutward(const o2::tpc::TrackTPC& leg, const CosmicTime& cosmicTime, int sector, float x, float& y, float& z) const +{ + // outward continuation of a leg in the frame of a sector; z in the frame of the cosmic's time (a one-side TPC-only leg has z relative to + // its time0) + o2::track::TrackPar par = leg.getParamOut(); + if (!par.rotateParam(o2::math_utils::sector2Angle(sector % TPCGeo::getNumberOfSectorsA())) || !par.propagateParamTo(x, mBz)) { + return false; + } + y = par.getY(); + z = par.getZ() + tpcSide(leg) * (cosmicTime.tb - leg.getTime0()) * mCorrMap->getVDrift(); + return true; +} + +void CosmicsClusterCollectorSpec::roadTOF(const RecoContainer& data, const o2::tpc::TrackTPC* const* legs, const CosmicTime& cosmicTime, int icosm, std::vector& out, const std::unordered_set& matched, float& timeTOFMUS, float& scoreTOFPair, float& scoreTOFReversed) const +{ + // TOF clusters along the legs' outward continuations (the road in PbPb also contains hits of collision tracks). Preferred: the top/bottom + // pair whose time difference matches the muon's flight between them; its mean time fixes the cosmic's time, also when the legs' brackets + // leave it open by tens of mus (one-side legs on the same side), and the z of one-side legs is shifted to it. Otherwise per leg the + // cluster closest to its continuation. The same search in the impossible order (bottom hit first) only finds accidental pairs: its best + // score is kept as QA of the flag's background. + constexpr float MaxFlightMUS = 0.1f; // flight time of the muon between the TPC and the TOF, slow tails + constexpr float CmPerNS = 29.9792458f; + const auto window = timeWindowMUS(cosmicTime); + const float zTimeTol = 0.5f * (window.second - window.first) / mTPCTBinMUS * mCorrMap->getVDrift(); // z uncertainty of one-side legs + const float vDriftPerMUS = mCorrMap->getVDrift() / mTPCTBinMUS; + const float cosmicTimeMUS = cosmicTime.tb * mTPCTBinMUS; + int legSide[2] = {0, 0}; + float curvature = 0.f; // |1/R| of the muon [1/cm] for the flight path between the TOF hits + int nLegs = 0; + for (int leg = 0; leg < 2; leg++) { + if (legs[leg]) { + legSide[leg] = tpcSide(*legs[leg]); + curvature += std::abs(legs[leg]->getCurvature(mBz)); + nLegs++; + } + } + curvature = nLegs ? curvature / nLegs : 0.f; + struct Candidate { + int index; + bool matched; // hit of the leg's matched global track (already stored) + double timeNS; // since the start of the TF + float dy; + float dzAtCosmicTime; + float gx; + float gy; + float gz; + }; + std::vector candidates[2]; + const auto clusters = data.getTOFClusters(); + int best[2] = {-1, -1}; + float bestScore[2] = {1.f, 1.f}; + for (int i = 0; i < (int)clusters.size(); i++) { + const auto& c = clusters[i]; + const float timeMUS = c.getTime() * 1e-6f; // [ps] since the start of the TF + if (timeMUS < window.first - MaxFlightMUS || timeMUS > window.second + MaxFlightMUS) { + continue; + } + const bool isMatched = matched.count(i) > 0; + const float alpha = o2::math_utils::sector2Angle(c.getSector()); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + for (int leg = 0; leg < 2; leg++) { + float y = 0.f; + float z = 0.f; + if (!legs[leg] || !predictOutward(*legs[leg], cosmicTime, c.getSector(), c.getX(), y, z)) { + continue; + } + const float normY = (c.getY() - y) / mRoadTOF; + const float normZ = (c.getZ() - z) / (mRoadTOF + zTimeTol); + const float score = std::max(normY * normY, normZ * normZ); + if (score >= 1.f) { + continue; + } + if (!isMatched && score < bestScore[leg]) { + bestScore[leg] = score; + best[leg] = i; + } + candidates[leg].push_back(Candidate{i, isMatched, c.getTime() * 1e-3, c.getY() - y, c.getZ() - z, c.getX() * cosAlpha - c.getY() * sinAlpha, + c.getX() * sinAlpha + c.getY() * cosAlpha, c.getZ()}); + } + } + float bestPairScore = -1.f; + float bestReversedScore = -1.f; + int bestPair[2] = {-1, -1}; + for (const auto& c0 : candidates[0]) { + for (const auto& c1 : candidates[1]) { + if (c0.index == c1.index) { + continue; + } + // top / bottom by the hits' height: a muon goes down (for near-horizontal cosmics both orders are possible, but rare) + const auto& top = c0.gy > c1.gy ? c0 : c1; + const auto& bottom = c0.gy > c1.gy ? c1 : c0; + // flight path along the helix: arc in the transverse plane from the chord, then the dip + const float chordXY = std::hypot(top.gx - bottom.gx, top.gy - bottom.gy); + const float halfAngleSin = 0.5f * curvature * chordXY; + const float arcXY = halfAngleSin > 1e-4f && halfAngleSin < 1.f ? 2.f * std::asin(halfAngleSin) / curvature : chordXY; + const float length = std::hypot(arcXY, top.gz - bottom.gz); + const float timeDiff = float(top.timeNS - bottom.timeNS); + const float flightDev = timeDiff + length / CmPerNS; // the muon crosses the top TOF first + const float reversedDev = timeDiff - length / CmPerNS; // bottom hit first: impossible for a muon, only accidental pairs + const bool reversed = std::abs(flightDev) > mTOFFlightTol; // the two cannot both pass (L/c ~ 25 ns >> tolerance) + if (reversed && std::abs(reversedDev) > mTOFFlightTol) { + continue; + } + const double pairTimeNS = 0.5 * (c0.timeNS + c1.timeNS); + const float shiftMUS = float(pairTimeNS * 1e-3) - cosmicTimeMUS; + const float dz0 = c0.dzAtCosmicTime - legSide[0] * shiftMUS * vDriftPerMUS; + const float dz1 = c1.dzAtCosmicTime - legSide[1] * shiftMUS * vDriftPerMUS; + if (std::abs(dz0) > mRoadTOF || std::abs(dz1) > mRoadTOF) { + continue; + } + const float timeDev = reversed ? reversedDev : flightDev; + const float score = (c0.dy * c0.dy + c1.dy * c1.dy + dz0 * dz0 + dz1 * dz1) / (mRoadTOF * mRoadTOF) + timeDev * timeDev / (mTOFFlightTol * mTOFFlightTol); + if (reversed) { + if (bestReversedScore < 0.f || score < bestReversedScore) { + bestReversedScore = score; + } + continue; + } + if (bestPairScore < 0.f || score < bestPairScore) { + bestPairScore = score; + bestPair[0] = c0.index; + bestPair[1] = c1.index; + timeTOFMUS = float(pairTimeNS * 1e-3); + } + } + } + scoreTOFPair = bestPairScore; + scoreTOFReversed = bestReversedScore; + uint8_t flags = o2::dataformats::HitRoad; + if (bestPairScore >= 0.f) { + best[0] = bestPair[0]; + best[1] = bestPair[1]; + flags |= o2::dataformats::HitTOFFlight; + } else if (best[0] >= 0 && best[0] == best[1]) { // one TOF cluster in the roads of both legs: keep it for the closer one + best[bestScore[0] <= bestScore[1] ? 1 : 0] = -1; + } + for (int leg = 0; leg < 2; leg++) { + if (best[leg] < 0) { + continue; + } + const auto& c = clusters[best[leg]]; + if (matched.count(best[leg])) { // the leg's matched hit, already stored: flag it as part of the flight pair + for (auto& cl : out) { + if (cl.leg == leg && cl.channel == c.getMainContributingChannel() && cl.timeRaw == c.getTimeRaw()) { + cl.flags |= o2::dataformats::HitTOFFlight; + } + } + if (mDebugOut) { + auto& dbg = mDbgCosmics[icosm]; + for (size_t j = 0; j < dbg.tofFlags.size(); j++) { + if (dbg.tofLeg[j] == leg && dbg.tofChannel[j] == c.getMainContributingChannel() && dbg.tofTimeRaw[j] == c.getTimeRaw()) { + dbg.tofFlags[j] |= o2::dataformats::HitTOFFlight; + } + } + } + continue; + } + auto& cl = out.emplace_back(); + cl.timeRaw = c.getTimeRaw(); + cl.tot = c.getTot(); + cl.channel = c.getMainContributingChannel(); + cl.leg = leg; + cl.flags = flags; + if (mDebugOut) { + writeDebugTOF(c, icosm, leg, flags); + } + } +} + +void CosmicsClusterCollectorSpec::polish(const o2::tpc::TrackTPC* const* legs, float timeMUS, o2::gpu::GPUO2InterfaceRefit& refitter, o2::dataformats::CosmicTrack& out) const +{ + // refit of a cosmic with TPC-only legs at its TOF time, as MatchCosmics::refitWinners: the bottom leg inward, then to the closest approach + // to the beam line; the top leg inward, then to the same point; the two halves combined. Muon mass, energy loss along the muon's flight + // (top to bottom): going inward along the bottom leg and up to the closest approach is against the flight (gain), the top leg inward + // and down to the closest approach is with it (loss). With a precise time the one-side legs get their real z, hence also the right + // material, which the matcher's time (the legs' bracket overlap for one-side legs on the same side) cannot give. + constexpr int ELossGain = 1; + constexpr int ELossLoss = -1; + const auto& matchParams = o2::globaltracking::MatchCosmicsParams::Instance(); // propagation settings as in the matcher + const auto prop = o2::base::Propagator::Instance(); + const float timeTB = timeMUS / mTPCTBinMUS; + o2::track::TrackParCov bottom = legs[0]->getParamOut(); + bottom.setPID(o2::track::PID::Muon, true); + bottom.resetCovariance(); + if (std::abs(mBz) <= 0.01f) { + // B = 0: both legs carry the same conventional q/pt; as the matcher, flip the bottom one so that it is right after the inversion + bottom.setQ2Pt(-o2::track::kMostProbablePt); + } + if (refitter.RefitTrackAsTrackParCov(bottom, legs[0]->getClusterRef(), timeTB, nullptr, false, false, ELossGain) < 0) { + return; + } + bottom.invert(); + const o2::dataformats::VertexBase origin; + if (!prop->propagateToDCABxByBz(origin, bottom, matchParams.maxStep, matchParams.matCorr, nullptr, nullptr, ELossGain)) { + return; + } + o2::track::TrackParCov top = legs[1]->getParamOut(); + top.setPID(o2::track::PID::Muon, true); + if (refitter.RefitTrackAsTrackParCov(top, legs[1]->getClusterRef(), timeTB, nullptr, false, true, ELossLoss) < 0) { + return; + } + if (!top.rotate(bottom.getAlpha()) || !prop->PropagateToXBxByBz(top, bottom.getX(), matchParams.maxSnp, matchParams.maxStep, matchParams.matCorr, nullptr, ELossLoss)) { + return; + } + o2::track::TrackParCov::MatrixDSym5 cov5; + const float chi2Match = bottom.getPredictedChi2(top, cov5); + if (!bottom.update(top, cov5)) { + return; + } + out.polished = bottom; + out.chi2MatchPolished = chi2Match; +} + +void CosmicsClusterCollectorSpec::roadTRD(const RecoContainer& data, const o2::tpc::TrackTPC* const* legs, const CosmicTime& cosmicTime, int icosm, std::vector& out, const std::unordered_set& matched) const +{ + // per leg and layer the TRD tracklet closest to the leg's outward continuation, from triggers whose readout window can contain the cosmic + constexpr float ReadoutWindowMUS = 3.f; // a cosmic leaves tracklets only if it passes within the readout window after a trigger + constexpr float PadLength = 10.f; // [cm] longest TRD pads: the tracklet z is the pad-row centre + constexpr int NLayers = o2::trd::constants::NLAYER; + const auto window = timeWindowMUS(cosmicTime); + const float zTimeTol = 0.5f * (window.second - window.first) / mTPCTBinMUS * mCorrMap->getVDrift(); + const auto tracklets = data.getTRDTracklets(); + const auto calibrated = data.getTRDCalibratedTracklets(); + if (calibrated.size() != tracklets.size()) { + if (!mTRDMismatchWarned) { + LOGP(warning, "{} calibrated vs {} raw TRD tracklets: TRD road of cosmics skipped", calibrated.size(), tracklets.size()); + mTRDMismatchWarned = true; + } + return; + } + struct Candidate { + int tracklet = -1; + int trigBC = 0; + float score = 1.f; + }; + Candidate best[2][NLayers]; + for (const auto& trig : data.getTRDTriggerRecords()) { + const int trigBC = trig.getBCData().differenceInBC(mTFStart); + const float trigMUS = trigBC * o2::constants::lhc::LHCBunchSpacingMUS; + if (trigMUS < window.first - ReadoutWindowMUS || trigMUS > window.second) { + continue; + } + for (int it = trig.getFirstTracklet(); it < trig.getFirstTracklet() + trig.getNumberOfTracklets(); it++) { + if (matched.count(it)) { + continue; + } + const int detector = tracklets[it].getDetector(); + const int sector = detector / o2::trd::constants::NCHAMBERPERSEC; + const int layer = detector % NLayers; + const auto& point = calibrated[it]; + for (int leg = 0; leg < 2; leg++) { + float y = 0.f; + float z = 0.f; + if (!legs[leg] || !predictOutward(*legs[leg], cosmicTime, sector, point.getX(), y, z)) { + continue; + } + const float normY = (point.getY() - y) / mRoadTRD; + const float normZ = (point.getZ() - z) / (mRoadTRD + PadLength + zTimeTol); + const float score = std::max(normY * normY, normZ * normZ); + if (score < best[leg][layer].score) { + best[leg][layer] = {it, trigBC, score}; + } + } + } + } + for (int leg = 0; leg < 2; leg++) { + for (int layer = 0; layer < NLayers; layer++) { + const auto& cand = best[leg][layer]; + if (cand.tracklet < 0) { + continue; + } + auto& tr = out.emplace_back(); + tr.word = tracklets[cand.tracklet].getTrackletWord(); + tr.trigBC = cand.trigBC; + tr.layer = layer; + tr.leg = leg; + tr.flags = o2::dataformats::HitRoad; + if (mDebugOut) { + const auto& point = calibrated[cand.tracklet]; + const int sector = tracklets[cand.tracklet].getDetector() / o2::trd::constants::NCHAMBERPERSEC; + float y = 0.f; + float z = 0.f; + predictOutward(*legs[leg], cosmicTime, sector, point.getX(), y, z); // succeeded in the search + auto& dbg = mDbgCosmics[icosm]; + dbg.trdLeg.push_back(leg); + dbg.trdLayer.push_back(layer); + dbg.trdX.push_back(point.getX()); + dbg.trdY.push_back(point.getY()); + dbg.trdZ.push_back(point.getZ()); + dbg.trdDy.push_back(point.getY() - y); + dbg.trdDz.push_back(point.getZ() - z); + dbg.trdTrigMUS.push_back(cand.trigBC * float(o2::constants::lhc::LHCBunchSpacingMUS)); + } + } + } +} + +void CosmicsClusterCollectorSpec::roadITS(const RecoContainer& data, const o2::dataformats::TrackCosmics& cosm, const CosmicTime& cosmicTime, int legsSide, int icosm, std::vector& out, + std::vector& requests, const std::unordered_set& matched) const +{ + if (data.getITSPerLayer() || mITSChipCentres.empty()) { + return; + } + // trajectory near the beam line from the combined cosmic, sampled every cm along its direction (closest approach at local x = 0) + constexpr float MaxRadius = 45.f; // [cm] outer ITS barrel + margin + constexpr float ChipHalfDiag = 1.7f; // [cm] half diagonal of an ITS chip + o2::track::TrackPar par = cosm; + if (!par.rotateParam(par.getPhi())) { + return; + } + std::vector> points; + for (float x = -MaxRadius - 5.f; x <= MaxRadius + 5.f; x += 1.f) { + bool ok = false; + const auto point = par.getXYZGloAt(x, mBz, ok); + if (ok) { + points.push_back(point); + } + } + if (points.size() < 2) { + return; + } + // closest segment in the transverse plane: distance, z of the trajectory there + auto closest = [&points](float x, float y, float& dist2, float& zTraj) { + dist2 = 1e10f; + for (size_t i = 0; i + 1 < points.size(); i++) { + const float segX = points[i + 1].X() - points[i].X(); + const float segY = points[i + 1].Y() - points[i].Y(); + const float len2 = segX * segX + segY * segY; + const float frac = len2 > 0.f ? std::clamp(((x - points[i].X()) * segX + (y - points[i].Y()) * segY) / len2, 0.f, 1.f) : 0.f; + const float dx = x - (points[i].X() + frac * segX); + const float dy = y - (points[i].Y() + frac * segY); + if (dx * dx + dy * dy < dist2) { + dist2 = dx * dx + dy * dy; + zTraj = points[i].Z() + frac * (points[i + 1].Z() - points[i].Z()); + } + } + }; + float minR2 = 1e10f; + float pcaY = 0.f; + for (const auto& point : points) { + const float r2 = point.X() * point.X() + point.Y() * point.Y(); + if (r2 < minR2) { + minR2 = r2; + pcaY = point.Y(); + } + } + if (minR2 > MaxRadius * MaxRadius) { // the cosmic does not cross the ITS + return; + } + std::vector candidateChip(mITSChipCentres.size(), false); + bool anyChip = false; + const float chipCut2 = (mRoadITS + ChipHalfDiag) * (mRoadITS + ChipHalfDiag); + for (size_t chip = 0; chip < mITSChipCentres.size(); chip++) { + float dist2 = 0.f; + float zTraj = 0.f; + closest(mITSChipCentres[chip].X(), mITSChipCentres[chip].Y(), dist2, zTraj); + if (dist2 < chipCut2) { + candidateChip[chip] = true; + anyChip = true; + } + } + if (!anyChip) { + return; + } + const auto window = timeWindowMUS(cosmicTime); + const float zTimeTol = 0.5f * (window.second - window.first) / mTPCTBinMUS * mCorrMap->getVDrift(); // z of the cosmic is in the frame of its time + // the refitted cosmic has the z of its TPC time; with legs on one TPC side and a road time from the TOF its z moves by side * vD * dt + const float zShift = legsSide * (cosmicTime.tb - cosm.getTimeMUS().getTimeStamp() / mTPCTBinMUS) * mCorrMap->getVDrift(); + const float rofLengthMUS = o2::itsmft::DPLAlpideParam::Instance().roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingMUS; + // per half of the cosmic and layer the ITS cluster closest to the trajectory (the road near the beam line also contains collision clusters) + constexpr int NLayers = 7; + struct Candidate { + int index = -1; + int rofBC = 0; + float score = 1.f; + }; + Candidate best[2][NLayers]; + const auto clusters = data.getITSClusters(); + auto geom = o2::its::GeometryTGeo::Instance(); + for (const auto& rof : data.getITSClustersROFRecords()) { + const int rofBC = rof.getBCData().differenceInBC(mTFStart); + const float rofMUS = rofBC * o2::constants::lhc::LHCBunchSpacingMUS; + if (rofMUS + rofLengthMUS < window.first || rofMUS > window.second) { + continue; + } + for (int idx = rof.getFirstEntry(); idx < rof.getFirstEntry() + rof.getNEntries(); idx++) { + const auto& c = clusters[idx]; + if (!candidateChip[c.getSensorID()] || matched.count(idx)) { + continue; + } + o2::math_utils::Point3D local; + if (mITSDict && c.getPatternID() != o2::itsmft::CompCluster::InvalidPatternID && !mITSDict->isGroup(c.getPatternID())) { + local = mITSDict->getClusterCoordinates(c); + } else { // anchor pixel: good enough for a cm road + o2::itsmft::SegmentationAlpide::detectorToLocalUnchecked(c.getRow(), c.getCol(), local); + } + const auto global = geom->getMatrixL2G(c.getSensorID()) * local; + float dist2 = 0.f; + float zTraj = 0.f; + closest(global.X(), global.Y(), dist2, zTraj); + const float normZ = (global.Z() - zTraj - zShift) / (mRoadITS + zTimeTol); + const float score = std::max(dist2 / (mRoadITS * mRoadITS), normZ * normZ); + const int half = global.Y() < pcaY ? 0 : 1; // bottom / top half of the cosmic + const int layer = geom->getLayer(c.getSensorID()); + if (layer >= 0 && layer < NLayers && score < best[half][layer].score) { + best[half][layer] = {idx, rofBC, score}; + } + } + } + for (int half = 0; half < 2; half++) { + for (int layer = 0; layer < NLayers; layer++) { + const auto& cand = best[half][layer]; + if (cand.index < 0) { + continue; + } + const auto& c = clusters[cand.index]; + auto& cl = out.emplace_back(); + cl.chipID = c.getSensorID(); + cl.row = c.getRow(); + cl.col = c.getCol(); + cl.pattID = c.getPatternID(); + cl.rofBC = cand.rofBC; + cl.leg = half; + cl.flags = o2::dataformats::HitRoad; + if (c.getPatternID() == o2::itsmft::CompCluster::InvalidPatternID || (mITSDict && mITSDict->isGroup(c.getPatternID()))) { + requests.push_back({icosm, int(out.size()) - 1, cand.index}); + } + } + } +} + +void CosmicsClusterCollectorSpec::fillITSPatterns(const RecoContainer& data, std::vector& requests, std::vector& cosmics) const +{ + // the TF's pattern stream holds, in cluster order, the patterns of the clusters with an invalid or a group pattern ID + if (!mITSDict) { + if (!mNoDictWarned) { + LOGP(warning, "No ITS cluster dictionary: patterns of ITS clusters of cosmics are not stored"); + mNoDictWarned = true; + } + return; + } + std::sort(requests.begin(), requests.end(), [](const ITSPattRequest& a, const ITSPattRequest& b) { return a.index < b.index; }); + const auto clusters = data.getITSClusters(); + auto pattIt = data.getITSClustersPatterns().begin(); + size_t ir = 0; + for (int k = 0; k < (int)clusters.size() && ir < requests.size(); k++) { + const auto pattID = clusters[k].getPatternID(); + if (pattID != o2::itsmft::CompCluster::InvalidPatternID && !mITSDict->isGroup(pattID)) { + continue; + } + const auto start = pattIt; + o2::itsmft::ClusterPattern::skipPattern(pattIt); + for (; ir < requests.size() && requests[ir].index == k; ir++) { + auto& cosm = cosmics[requests[ir].cosmic]; + cosm.clITS[requests[ir].cluster].pattEntry = cosm.itsPatterns.size(); + cosm.itsPatterns.insert(cosm.itsPatterns.end(), start, pattIt); + } + } +} + +void CosmicsClusterCollectorSpec::addTRD(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, std::unordered_set& matched) const +{ + const auto& trk = data.getTrack(gid); // TPC-TRD or ITS-TPC-TRD track + const auto tracklets = data.getTRDTracklets(); + const auto trigs = data.getTRDTriggerRecords(); + for (int layer = 0; layer < 6; layer++) { + const int it = trk.getTrackletIndex(layer); + if (it < 0) { + continue; + } + auto& tr = out.emplace_back(); + tr.word = tracklets[it].getTrackletWord(); + tr.layer = layer; + tr.leg = leg; + tr.flags = o2::dataformats::HitMatched; + matched.insert(it); + auto trig = std::upper_bound(trigs.begin(), trigs.end(), it, [](int v, const o2::trd::TriggerRecord& t) { return v < t.getFirstTracklet(); }); + if (trig != trigs.begin()) { + tr.trigBC = std::prev(trig)->getBCData().differenceInBC(mTFStart); + } + } +} + +void CosmicsClusterCollectorSpec::flagDuplicates(std::vector& cosmics) const +{ + // the same muon can be matched twice, e.g. when a leg is split into two TPC tracks: the roads then collect largely the same clusters + // (PbPb 567939: 10 of 47 cosmic pairs in a TF share 56-99 % of the clusters of the smaller one, all others none). The best one (TOF time, + // then more attached clusters) stays, the others point to it. + constexpr float MinSharedFraction = 0.3f; + const size_t n = cosmics.size(); + if (n < 2) { + return; + } + std::vector> keys(n); // sorted cluster addresses: sector, row and the packed time and pad of the raw cluster + std::vector nAttached(n, 0); + for (size_t ic = 0; ic < n; ic++) { + for (const auto* legClusters : {&cosmics[ic].clTPCBottom, &cosmics[ic].clTPCTop}) { + for (const auto& c : *legClusters) { + keys[ic].push_back((uint64_t(c.sector) << 56) | (uint64_t(c.row) << 48) | (uint64_t(c.cl.timeFlagsPacked) << 16) | c.cl.padPacked); + nAttached[ic] += c.isAttached(); + } + } + std::sort(keys[ic].begin(), keys[ic].end()); + } + std::vector order(n); + std::iota(order.begin(), order.end(), 0); + std::stable_sort(order.begin(), order.end(), [&](int a, int b) { + const bool timedA = cosmics[a].hasTOFTime(); + const bool timedB = cosmics[b].hasTOFTime(); + return timedA != timedB ? timedA : nAttached[a] > nAttached[b]; + }); + std::vector kept; + for (int ic : order) { + for (int ik : kept) { + std::vector shared; + std::set_intersection(keys[ic].begin(), keys[ic].end(), keys[ik].begin(), keys[ik].end(), std::back_inserter(shared)); + const size_t nMin = std::min(keys[ic].size(), keys[ik].size()); + if (nMin > 0 && shared.size() >= MinSharedFraction * nMin) { + cosmics[ic].duplicateOf = ik; + break; + } + } + if (cosmics[ic].duplicateOf < 0) { + kept.push_back(ic); + } + } +} + +void CosmicsClusterCollectorSpec::writeDebug(const o2::dataformats::CosmicTrack& cosm, int icosm) const +{ + constexpr int NSectorsA = TPCGeo::getNumberOfSectorsA(); + // the time the other-side corridor used (and the common frame zCos): the TOF time if the TOF road found one, else the TPC time + const float timeCosmic = (cosm.hasTOFTime() ? cosm.timeTOFMUS : cosm.cosmic.getTimeMUS().getTimeStamp()) / mTPCTBinMUS; + const bool absTimeKnown = cosm.cosmic.getTimeMUS().getTimeStampError() < mMaxAbsTimeErr; + int nAttached[2] = {0, 0}; + int nRoad[2] = {0, 0}; + int side[2] = {-2, -2}; + std::vector clLeg; + std::vector clFlags; + std::vector clSector; + std::vector clRow; + std::vector clPad; + std::vector clTime; + std::vector clQMax; + std::vector clQTot; + std::vector clX; + std::vector clY; + std::vector clZ; + std::vector clZCos; + std::vector clGx; + std::vector clGy; + const o2::tpc::TrackTPC* legs[2] = {&cosm.tpcBottom, &cosm.tpcTop}; + const std::vector* legClusters[2] = {&cosm.clTPCBottom, &cosm.clTPCTop}; + for (int leg = 0; leg < 2; leg++) { + if (legs[leg]->getNClusters() == 0) { // no TPC part + continue; + } + side[leg] = tpcSide(*legs[leg]); + for (const auto& c : *legClusters[leg]) { + // transform with the vertex time used by the road: the leg's time0, or the cosmic's time for clusters found on the other side + const float vertexTime = c.isAbsTime() ? timeCosmic : legs[leg]->getTime0(); + float x = 0.f; + float y = 0.f; + float z = 0.f; + mCorrMap->Transform(c.sector, c.row, c.cl.getPad(), c.cl.getTime(), x, y, z, vertexTime); + float xCos = 0.f; + float yCos = 0.f; + float zCos = 0.f; + mCorrMap->Transform(c.sector, c.row, c.cl.getPad(), c.cl.getTime(), xCos, yCos, zCos, timeCosmic); + const float alpha = o2::math_utils::sector2Angle(c.sector % NSectorsA); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + clLeg.push_back(leg); + clFlags.push_back(c.flags); + clSector.push_back(c.sector); + clRow.push_back(c.row); + clPad.push_back(c.cl.getPad()); + clTime.push_back(c.cl.getTime()); + clQMax.push_back(c.cl.getQmax()); + clQTot.push_back(c.cl.getQtot()); + clX.push_back(x); + clY.push_back(y); + clZ.push_back(z); + clZCos.push_back(zCos); + clGx.push_back(x * cosAlpha - y * sinAlpha); + clGy.push_back(x * sinAlpha + y * cosAlpha); + (c.isAttached() ? nAttached : nRoad)[leg]++; + } + } + std::vector itsLeg; + std::vector itsFlags; + std::vector itsChip; + std::vector itsRow; + std::vector itsCol; + std::vector itsPattID; + std::vector itsHasPatt; + std::vector itsRofBC; + std::vector itsGx; + std::vector itsGy; + std::vector itsGz; + std::vector itsXLoc; + std::vector itsZLoc; + for (const auto& c : cosm.clITS) { // local coordinates on the chip from the dictionary or the stored pattern + const o2::itsmft::CompClusterExt compCluster(c.row, c.col, c.pattID, c.chipID); + o2::math_utils::Point3D local; + if (c.pattEntry >= 0) { + o2::itsmft::ClusterPattern pattern; + auto pattIt = cosm.itsPatterns.begin() + c.pattEntry; + pattern.acquirePattern(pattIt); + local = o2::itsmft::TopologyDictionary::getClusterCoordinates(compCluster, pattern, c.pattID != o2::itsmft::CompCluster::InvalidPatternID); + } else if (mITSDict && c.pattID != o2::itsmft::CompCluster::InvalidPatternID) { + local = mITSDict->getClusterCoordinates(compCluster); + } else { // pattern not available: pixel centre + o2::itsmft::SegmentationAlpide::detectorToLocalUnchecked(c.row, c.col, local); + } + o2::math_utils::Point3D global(0.f, 0.f, 0.f); + if (!mITSChipCentres.empty()) { // geometry loaded for the ITS road + global = o2::its::GeometryTGeo::Instance()->getMatrixL2G(c.chipID) * local; + } + itsLeg.push_back(c.leg); + itsFlags.push_back(c.flags); + itsChip.push_back(c.chipID); + itsRow.push_back(c.row); + itsCol.push_back(c.col); + itsPattID.push_back(c.pattID); + itsHasPatt.push_back(c.pattEntry >= 0); + itsRofBC.push_back(c.rofBC); + itsGx.push_back(global.X()); + itsGy.push_back(global.Y()); + itsGz.push_back(global.Z()); + itsXLoc.push_back(local.X()); + itsZLoc.push_back(local.Z()); + } + const auto& dbg = mDbgCosmics[icosm]; + const auto& time = cosm.cosmic.getTimeMUS(); + (*mDebugOut) << "cosmics" + << "tf=" << mDbgTF << "cosm=" << icosm << "t=" << time.getTimeStamp() << "tErr=" << time.getTimeStampError() << "absTime=" << int(absTimeKnown) + << "chi2Match=" << cosm.cosmic.getChi2Match() << "chi2Refit=" << cosm.cosmic.getChi2Refit() << "q2pt=" << cosm.cosmic.getQ2Pt() + << "tgl=" << cosm.cosmic.getTgl() << "side0=" << side[0] << "side1=" << side[1] << "nAtt0=" << nAttached[0] << "nAtt1=" << nAttached[1] + << "nRoad0=" << nRoad[0] << "nRoad1=" << nRoad[1] << "nITS=" << int(cosm.clITS.size()) << "nTOF=" << int(cosm.clTOF.size()) + << "nTRD=" << int(cosm.trdTracklets.size()) << "tTOF=" << cosm.timeTOFMUS << "scoreTOF=" << cosm.scoreTOFPair << "scoreTOFRev=" << cosm.scoreTOFReversed << "dupOf=" << cosm.duplicateOf << "polChi2Match=" << cosm.chi2MatchPolished << "polQ2Pt=" << cosm.polished.getQ2Pt() << "polZ=" << cosm.polished.getZ() + << "mcEvent=" << (cosm.label.isSet() ? cosm.label.getEventID() : -1) << "mcTrack=" << (cosm.label.isSet() ? cosm.label.getTrackID() : -1) + << "clLeg=" << clLeg << "clFlags=" << clFlags << "clSector=" << clSector << "clRow=" << clRow << "clPad=" << clPad << "clTime=" << clTime + << "clQMax=" << clQMax << "clQTot=" << clQTot << "clX=" << clX << "clY=" << clY << "clZ=" << clZ << "clZCos=" << clZCos << "clGx=" << clGx + << "clGy=" << clGy << "roadLeg=" << dbg.roadLeg << "roadFrame=" << dbg.roadFrame << "roadSector=" << dbg.roadSector << "roadRow=" << dbg.roadRow + << "roadX=" << dbg.roadX << "roadY=" << dbg.roadY << "roadZ=" << dbg.roadZ << "roadZCos=" << dbg.roadZCos << "roadGx=" << dbg.roadGx + << "roadGy=" << dbg.roadGy << "roadSnp=" << dbg.roadSnp << "roadTgl=" << dbg.roadTgl << "itsLeg=" << itsLeg << "itsFlags=" << itsFlags + << "itsChip=" << itsChip << "itsRow=" << itsRow << "itsCol=" << itsCol << "itsPattID=" << itsPattID << "itsHasPatt=" << itsHasPatt + << "itsRofBC=" << itsRofBC << "itsGx=" << itsGx << "itsGy=" << itsGy << "itsGz=" << itsGz << "itsXLoc=" << itsXLoc << "itsZLoc=" << itsZLoc + << "tofLeg=" << dbg.tofLeg << "tofChannel=" << dbg.tofChannel << "tofFlags=" << dbg.tofFlags << "tofTimeRaw=" << dbg.tofTimeRaw + << "tofTime=" << dbg.tofTime << "tofTot=" << dbg.tofTot << "tofX=" << dbg.tofX << "tofY=" << dbg.tofY << "tofZ=" << dbg.tofZ << "tofGx=" << dbg.tofGx + << "tofGy=" << dbg.tofGy << "trdLeg=" << dbg.trdLeg << "trdLayer=" << dbg.trdLayer << "trdX=" << dbg.trdX << "trdY=" << dbg.trdY << "trdZ=" << dbg.trdZ + << "trdDy=" << dbg.trdDy << "trdDz=" << dbg.trdDz << "trdTrigMUS=" << dbg.trdTrigMUS << "\n"; +} + +void CosmicsClusterCollectorSpec::writeDebugTOF(const o2::tof::Cluster& c, int icosm, int leg, uint8_t flags) const +{ + // the TOF cluster position is in the frame of its sector + const float alpha = o2::math_utils::sector2Angle(c.getSector()); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + auto& dbg = mDbgCosmics[icosm]; + dbg.tofLeg.push_back(leg); + dbg.tofChannel.push_back(c.getMainContributingChannel()); + dbg.tofFlags.push_back(flags); + dbg.tofTimeRaw.push_back(c.getTimeRaw()); + dbg.tofTime.push_back(c.getTime()); + dbg.tofTot.push_back(c.getTot()); + dbg.tofX.push_back(c.getX()); + dbg.tofY.push_back(c.getY()); + dbg.tofZ.push_back(c.getZ()); + dbg.tofGx.push_back(c.getX() * cosAlpha - c.getY() * sinAlpha); + dbg.tofGy.push_back(c.getX() * sinAlpha + c.getY() * cosAlpha); +} + +void CosmicsClusterCollectorSpec::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) +{ + if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + return; + } + if (matcher == ConcreteDataMatcher("ITS", "CLUSDICT", 0)) { + mITSDict = (const o2::itsmft::TopologyDictionary*)obj; + return; + } +} + +void CosmicsClusterCollectorSpec::endOfStream(EndOfStreamContext& ec) +{ + mDebugOut.reset(); + LOGP(info, "Cosmics cluster collector: {} cosmics, {} attached and {} road TPC clusters; Cpu: {:.3e} Real: {:.3e} s in {} slots", + mNCosmics, mNClAttached, mNClCorridor, mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); +} + +DataProcessorSpec getCosmicsClusterCollectorSpec(GTrackID::mask_t src, bool useMC, bool itsStag, DetID::mask_t roadDets) +{ + if (itsStag && roadDets[DetID::ITS]) { + LOGP(warning, "Staggered ITS clusters are not supported: ITS road of cosmics disabled"); + roadDets &= ~DetID::getMask(DetID::ITS); + } + auto dataRequest = std::make_shared(); + dataRequest->setITSPerLayer(itsStag); + dataRequest->requestTracks(src, false); + dataRequest->requestClusters(src, false); + if (roadDets[DetID::ITS]) { + dataRequest->requestITSClusters(false); + } + if (roadDets[DetID::TOF]) { + dataRequest->requestTOFClusters(false); + } + if (roadDets[DetID::TRD]) { + dataRequest->requestTRDTracklets(false); + } + dataRequest->requestCoscmicTracks(useMC); + auto ggRequest = std::make_shared(false, // orbitResetTime + false, // GRPECS + false, // GRPLHCIF + true, // GRPMagField + roadDets[DetID::TOF], // askMatLUT: polish refit of TOF-timed cosmics + roadDets[DetID::ITS] ? o2::base::GRPGeomRequest::Aligned : o2::base::GRPGeomRequest::None, // ITS road: chip positions + dataRequest->inputs, + true); + dataRequest->inputs.emplace_back("corrMap", o2::header::gDataOriginTPC, "TPCCORRMAP", 0, Lifetime::Timeframe); + + std::vector outputs; + outputs.emplace_back("GLO", "COSMFULL", 0, Lifetime::Timeframe); + outputs.emplace_back("GLO", "COSMFULLTF", 0, Lifetime::Timeframe); + outputs.emplace_back("GLO", "COSMFULLTFID", 0, Lifetime::Timeframe); + + return DataProcessorSpec{ + "cosmics-cluster-collector", + dataRequest->inputs, + outputs, + AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, useMC, roadDets)}, + Options{ + {"corridor-width", VariantType::Float, 1.f, {"half-width of the road around each leg [cm]"}}, + {"max-abs-time-err", VariantType::Float, 0.5f, {"max. time error of a cosmic [mus] to search the other TPC side of its one-side legs"}}, + {"max-cosmics-per-tf", VariantType::Int, 100, {"write at most this many cosmics per TF (with their clusters), the others are dropped"}}, + {"tof-road-width", VariantType::Float, 5.f, {"half-width of the road at the TOF [cm]"}}, + {"tof-flight-tolerance", VariantType::Float, 2.f, {"max. deviation of the top/bottom TOF time difference from the muon's flight time [ns]"}}, + {"tof-time-error", VariantType::Float, 0.1f, {"error of a cosmic's TOF time for the TPC other-side corridor and the TRD / ITS roads [mus]"}}, + {"trd-road-width", VariantType::Float, 5.f, {"half-width of the road at the TRD [cm] (in z plus the pad length)"}}, + {"its-road-width", VariantType::Float, 1.5f, {"half-width of the road in the ITS [cm]"}}, + {"debug-tree", VariantType::Bool, false, {"write cosmics_collector_debug.root with transformed clusters and road points (test runs)"}}}}; +} + +} // namespace o2::globaltracking diff --git a/Detectors/GlobalTrackingWorkflow/src/TrackCosmicsWriterSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/TrackCosmicsWriterSpec.cxx index 800978f7a4db3..3cc023a85f510 100644 --- a/Detectors/GlobalTrackingWorkflow/src/TrackCosmicsWriterSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/TrackCosmicsWriterSpec.cxx @@ -16,6 +16,8 @@ #include "DPLUtils/MakeRootTreeWriterSpec.h" #include "ReconstructionDataFormats/TrackCosmics.h" #include "SimulationDataFormat/MCCompLabel.h" +#include "DataFormatsGlobalTracking/CosmicTrack.h" +#include "CommonDataFormat/TFIDInfo.h" using namespace o2::framework; @@ -51,5 +53,20 @@ DataProcessorSpec getTrackCosmicsWriterSpec(bool useMC) ""})(); } +DataProcessorSpec getCosmicsFullWriterSpec() +{ + auto logger = [](std::vector const& cosmics) { + LOG(info) << "Writing " << cosmics.size() << " cosmics with clusters"; + }; + return MakeRootTreeWriterSpec("cosmics-full-writer", + "o2_cosmics_full.root", + "cosmicsFull", + -1, // do not limit number of events to store + 100, // periodically autosave + BranchDefinition>{InputSpec{"cosmics", "GLO", "COSMFULL", 0}, "cosmics", 1, logger}, + BranchDefinition{InputSpec{"tfinfo", "GLO", "COSMFULLTF", 0}, "tfInfo", 1}, + BranchDefinition{InputSpec{"tfid", "GLO", "COSMFULLTFID", 0}, "tfID", 1})(); +} + } // namespace globaltracking } // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx index 85724125b1c35..26817a71e42a1 100644 --- a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx @@ -28,6 +28,7 @@ #include "GlobalTrackingWorkflow/CosmicsMatchingSpec.h" #include "GlobalTracking/MatchCosmicsParams.h" #include "GlobalTrackingWorkflow/TrackCosmicsWriterSpec.h" +#include "GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h" #include "Algorithm/RangeTokenizer.h" #include "DetectorsRaw/HBFUtilsInitializer.h" #include "Framework/CallbacksPolicy.h" @@ -53,6 +54,8 @@ void customize(std::vector& workflowOptions) {"disable-root-input", o2::framework::VariantType::Bool, false, {"disable root-files input reader"}}, {"disable-root-output", o2::framework::VariantType::Bool, false, {"disable root-files output writer"}}, {"use-pv-info", o2::framework::VariantType::Bool, false, {"request primary vertex for relevant cuts in the collision/cosmics interleaved data"}}, + {"enable-cluster-output", o2::framework::VariantType::Bool, false, {"collect the raw clusters of the cosmics (legs + road around them) and write them with the cosmics"}}, + {"road-detectors", VariantType::String, "ITS,TOF,TRD", {"with --enable-cluster-output: detectors whose hits along the cosmic are collected besides the TPC road"}}, {"cosmics-preset", VariantType::String, "", {"named set of cosmicsMatch settings applied before --configKeyValues (which can override single keys): physics-v1 = cosmics in collision data"}}, {"track-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of sources to use"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}}; @@ -72,6 +75,9 @@ void customize(std::vector& policies) policies.push_back(o2::tpc::TPCSectorCompletionPolicy("cosmics-matcher", o2::tpc::TPCSectorCompletionPolicy::Config::RequireAll, InputSpec{"cluster", o2::framework::ConcreteDataTypeMatcher{"TPC", "CLUSTERNATIVE"}})()); + policies.push_back(o2::tpc::TPCSectorCompletionPolicy("cosmics-cluster-collector", + o2::tpc::TPCSectorCompletionPolicy::Config::RequireAll, + InputSpec{"cluster", o2::framework::ConcreteDataTypeMatcher{"TPC", "CLUSTERNATIVE"}})()); policies.push_back(CompletionPolicyHelpers::consumeWhenAllOrdered(".*cosm.*[W,w]riter.*")); } @@ -118,14 +124,30 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) } bool usePV = configcontext.options().get("use-pv-info"); specs.emplace_back(o2::globaltracking::getCosmicsMatchingSpec(src, usePV, useMC, doStag)); + bool clusterOutput = configcontext.options().get("enable-cluster-output"); + if (clusterOutput) { + if (!src[GID::TPC]) { + LOG(fatal) << "--enable-cluster-output needs TPC tracks in --track-sources"; + } + const auto roadDets = DetID::getMask(configcontext.options().get("road-detectors")) & DetID::getMask("ITS,TOF,TRD"); + for (auto det : {DetID::ITS, DetID::TOF, DetID::TRD}) { + if (roadDets[det]) { + srcCl |= GID::getSourceMask(det == DetID::ITS ? GID::ITS : (det == DetID::TOF ? GID::TOF : GID::TRD)); + } + } + specs.emplace_back(o2::globaltracking::getCosmicsClusterCollectorSpec(src, useMC, doStag, roadDets)); + } - o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, src, src, src, useMC, dummy); // clusters MC is not needed + o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCl, src, src, useMC, dummy); // clusters MC is not needed if (usePV) { o2::globaltracking::InputHelper::addInputSpecsPVertex(configcontext, specs, useMC); // P-vertex is always needed } if (!disableRootOut) { specs.emplace_back(o2::globaltracking::getTrackCosmicsWriterSpec(useMC)); + if (clusterOutput) { + specs.emplace_back(o2::globaltracking::getCosmicsFullWriterSpec()); + } } // configure dpl timer to inject correct firstTForbit: start from the 1st orbit of TF containing 1st sampled orbit From c40017072b28bd12b5ee7f51791d1391ef76b296 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Thu, 8 Oct 2026 08:58:13 +0200 Subject: [PATCH 05/13] MatchCosmics: compare same-side TPC-only legs refitted at a common time MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A TPC-only track has no time measurement: its time0 is the tracker's guess (GPUTPCGMTrackParam::ShiftZ: z = 0 at the beam line if the helix can come from it, otherwise the latest cluster at z = 0.521 * x, then clamped to the drift volume). For a cosmic, which does not come from the beam line, the two legs get independent guesses, typically tens of mus apart, so their clusters were distortion-corrected at different z and the legs disagree (y, snp, q/pt) although they are one track. For legs on opposite sides the time comes from z continuity, but legs on the same side were compared as tracked. New MatchCosmicsParams refitSameSideAtCommonTime (default false, not in the preset physics-v1, see below): after the cheap cuts (tgl, q/pt, same half), same-side TPC-only legs are refitted at the centre of their brackets' overlap (the time the refit of the winners uses) and compared there (MatchCosmics::refitSeedAtTime: TPC refit from the outer parameters, then the propagation to the DCA as for the seeds; energy loss along the muon's flight as in refitWinners). The TPC refitter is created once per TF in process() and shared with refitWinners. The matcher debug tree writes the legs compared and a commonTime flag. Cosmic MC with PbPb-like distortions (190 TFs, preset physics-v1, off -> on): correctly matched muons 6848 -> 7324 of 8639 (79.3 -> 84.8 %); same side 77.0 -> 84.2 %, its loss to the crude chi2 cut 6.9 -> 0.8 %; A×C unchanged; 0 wrong-muon pairs in both; matcher +12 ms per TF (~100 refits). With the true MC time instead of the window centre the same-side pass rate of the chi2 cut would be 99.2 % (centre 98.0 %, legs at their own time0s 84.7 %). PbPb 2025 (4 TFs): ~530 refits per TF, +15 ms per TF. Not in the preset: in PbPb (568041, 310 CTFs, 16541 TFs) it adds mostly pairs of collision tracks, whose time0 was about right and which the refit at the window centre mis-corrects: +352 same-side cosmics but only +5 TOF-tagged (2.6 % of the new ones vs 15.3 % of the common ones), and ~35 A×C cosmics lose a leg to a refitted same-side pair. Co-Authored-By: Claude Opus 5.5 --- .../include/GlobalTracking/MatchCosmics.h | 7 +- .../GlobalTracking/MatchCosmicsParams.h | 7 +- Detectors/GlobalTracking/src/MatchCosmics.cxx | 133 ++++++++++++------ 3 files changed, 103 insertions(+), 44 deletions(-) diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h index acfc75f0679ba..86efff3416d82 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h @@ -41,7 +41,8 @@ class VDriftCorrFact; namespace gpu { class TPCFastTransformPOD; -} +class GPUO2InterfaceRefit; +} // namespace gpu namespace globaltracking { @@ -137,6 +138,7 @@ class MatchCosmics private: void updateTimeDependentParams(); RejFlag checkPair(int i, int j); + bool refitSeedAtTime(const TrackSeed& seed, float timeMUS, TrackSeed& out); void registerMatch(int i, int j, float chi2, float tCommon = 0.f, float tCommonErr = -1.f); void suppressMatch(int partner0, int partner1); void createSeeds(const o2::globaltracking::RecoContainer& data); @@ -164,6 +166,9 @@ class MatchCosmics float mQ2PtCutoff = 1e9; float mQ2PtCutoffOppositeSides = 1e9; const MatchCosmicsParams* mMatchParams = nullptr; + const o2::globaltracking::RecoContainer* mRecoData = nullptr; ///< inputs of the TF being processed + o2::gpu::GPUO2InterfaceRefit* mTPCRefitter = nullptr; ///< TPC refitter of the TF being processed (owned by process()) + size_t mNRefitsCommonTime = 0; ///< seeds refitted at the common time of a same-side pair in this TF std::vector mCosmicTracks; std::vector mCosmicTracksLbl; diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h index 9307d92ee5e3d..8e51c33e96288 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h @@ -38,9 +38,10 @@ struct MatchCosmicsParams : public o2::conf::ConfigurableParamHelper= this [cm] (0: no cut; rejects collision tracks in physics data) - float minSeedDCAxyNSigma = 0.f; // use only tracks with |DCA_xy| >= this * sigma(DCA_xy) (0: no cut; poorly measured collision tracks) - bool constrainTPCOnlyZ = false; // TPC-only legs: test z at a common time (same side, or a leg with known time), else require the time implied by z continuity in both brackets - bool vetoSameHalf = false; // reject pairs whose two legs lie on the same side of the closest approach (two pieces of one leg) + float minSeedDCAxyNSigma = 0.f; // use only tracks with |DCA_xy| >= this * sigma(DCA_xy) (0: no cut; poorly measured collision tracks) + bool constrainTPCOnlyZ = false; // TPC-only legs: test z at a common time (same side, or a leg with known time), else require the time implied by z continuity in both brackets + bool vetoSameHalf = false; // reject pairs whose two legs lie on the same side of the closest approach (two pieces of one leg) + bool refitSameSideAtCommonTime = false; // TPC-only legs on the same side: compare them refitted at the centre of their brackets' overlap instead of at their own time0s float nSigmaTError = 4.f; // number of sigmas on track time error for matching (except for TPC which provides an interval) float tpcExtraZError2 = 1.f; // extra error^2 on the TPC-only track Z coordinate float fiducialRIP = 1.0f; // consider track having |Y@x=0|< this as passing DCA cut (if requested) diff --git a/Detectors/GlobalTracking/src/MatchCosmics.cxx b/Detectors/GlobalTracking/src/MatchCosmics.cxx index 2e1fa08b1c146..cfd089c09a867 100644 --- a/Detectors/GlobalTracking/src/MatchCosmics.cxx +++ b/Detectors/GlobalTracking/src/MatchCosmics.cxx @@ -111,6 +111,17 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) } } } + // TPC refitter of this TF: same-side TPC-only legs compared at a common time (checkPair) and the refit of the winners + std::unique_ptr tpcRefitter; + if (data.inputsTPCclusters) { + tpcRefitter = std::make_unique(&data.inputsTPCclusters->clusterIndex, mTPCCorrMaps, mBz, data.getTPCTracksClusterRefs().data(), 0, + data.clusterShMapTPC.data(), data.occupancyMapTPC.data(), data.occupancyMapTPC.size(), nullptr, + o2::base::Propagator::Instance()); + } + mTPCRefitter = tpcRefitter.get(); + mRecoData = &data; + mNRefitsCommonTime = 0; + // sort in time bracket lower edge, putting rejected tracks in the end std::vector sortID(ntr); std::iota(sortID.begin(), sortID.end(), 0); @@ -134,6 +145,11 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) selectWinners(); refitWinners(data); + if (mNRefitsCommonTime) { + LOGP(info, "{} seeds refitted at the common time of same-side pairs", mNRefitsCommonTime); + } + mTPCRefitter = nullptr; + mRecoData = nullptr; mTFCount++; } @@ -144,16 +160,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) LOG(info) << "Refitting " << mWinners.size() << " winner matches"; int count = 0; auto tpcTBinMUSInv = 1. / mTPCTBinMUS; - const auto& tpcClusRefs = data.getTPCTracksClusterRefs(); - const auto& tpcClusShMap = data.clusterShMapTPC; - const auto& tpcClusOccMap = data.occupancyMapTPC; - std::unique_ptr tpcRefitter; - if (data.inputsTPCclusters) { - tpcRefitter = std::make_unique(&data.inputsTPCclusters->clusterIndex, - mTPCCorrMaps, mBz, - tpcClusRefs.data(), 0, tpcClusShMap.data(), - tpcClusOccMap.data(), tpcClusOccMap.size(), nullptr, o2::base::Propagator::Instance()); - } + auto* tpcRefitter = mTPCRefitter; // created in process() const auto& itsClusters = prepareITSClusters(data); // RS FIXME: this is probably a temporary solution, since ITS tracking over boundaries will likely change the TrackITS format @@ -461,6 +468,9 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) float chi2 = 1.e9f; float tCommon = 0.f; // time fixed by z continuity of TPC-only legs on opposite sides, used by the refit float tCommonErr = -1.f; // its error (< 0: not fixed) + TrackSeed seed0Common; // same-side TPC-only legs refitted at a common time (refitSameSideAtCommonTime) + TrackSeed seed1Common; + bool commonTime = false; // check // 1) crude check on tgl and q/pt (if B!=0). Note: back-to-back tracks will have mutually params (see TrackPar::invertParam) @@ -504,49 +514,58 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) break; } } - o2::track::TrackParCov seed1Inv = seed1; + // TPC-only legs on the same side: the tracker transformed each leg's clusters with its own time0, a guess (the two legs of a cosmic + // get guesses tens of mus apart), so the legs were distortion-corrected at different z and disagree although they are one track; + // compare them refitted at one common time, the centre of their brackets' overlap (the time the refit of the winners uses) + if (mMatchParams->refitSameSideAtCommonTime && seed0.tpcSide != 0 && seed0.tpcSide == seed1.tpcSide) { // tpcSide != 0: TPC-only one-side legs + const float tPair = seed0.tBracket.getOverlap(seed1.tBracket).mean(); + commonTime = refitSeedAtTime(seed0, tPair, seed0Common) && refitSeedAtTime(seed1, tPair, seed1Common); + } + const TrackSeed& leg0 = commonTime ? seed0Common : seed0; + const TrackSeed& leg1 = commonTime ? seed1Common : seed1; + o2::track::TrackParCov seed1Inv = leg1; seed1Inv.invert(); for (int i = 0; i < o2::track::kNParams; i++) { // add systematic error seed1Inv.updateCov(mMatchParams->systSigma2[i], o2::track::DiagMap[i]); } - if (!seed1Inv.rotate(seed0.getAlpha()) || - !o2::base::Propagator::Instance()->PropagateToXBxByBz(seed1Inv, seed0.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { + if (!seed1Inv.rotate(leg0.getAlpha()) || + !o2::base::Propagator::Instance()->PropagateToXBxByBz(seed1Inv, leg0.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { rej = RejProp; break; } - auto dSnp = seed0.getSnp() - seed1Inv.getSnp(); - if (dSnp * dSnp > (seed0.getSigmaSnp2() + seed1Inv.getSigmaSnp2()) * mMatchParams->crudeNSigma2Cut[o2::track::kSnp]) { + auto dSnp = leg0.getSnp() - seed1Inv.getSnp(); + if (dSnp * dSnp > (leg0.getSigmaSnp2() + seed1Inv.getSigmaSnp2()) * mMatchParams->crudeNSigma2Cut[o2::track::kSnp]) { rej = RejSnp; break; } - auto dY = seed0.getY() - seed1Inv.getY(); - if (dY * dY > (seed0.getSigmaY2() + seed1Inv.getSigmaY2()) * mMatchParams->crudeNSigma2Cut[o2::track::kY]) { + auto dY = leg0.getY() - seed1Inv.getY(); + if (dY * dY > (leg0.getSigmaY2() + seed1Inv.getSigmaY2()) * mMatchParams->crudeNSigma2Cut[o2::track::kY]) { rej = RejY; break; } - bool ignoreZ = seed0.origID.getSource() == o2d::GlobalTrackID::TPC || seed1.origID.getSource() == o2d::GlobalTrackID::TPC; + bool ignoreZ = leg0.origID.getSource() == o2d::GlobalTrackID::TPC || leg1.origID.getSource() == o2d::GlobalTrackID::TPC; if (ignoreZ && mMatchParams->constrainTPCOnlyZ) { // a TPC-only track with clusters on one side has z relative to its time0: z(t) = z + side * vD * (t - tRef); a CE-crossing or // non-TPC-only track has an absolute z (side 0). Bring both legs to a common time where possible and test z; for legs on opposite // sides, z continuity fixes the common time, which must lie in both time brackets. - const int side0 = seed0.tpcSide; - const int side1 = seed1.tpcSide; - const float sigZ2 = (seed0.getSigmaZ2() + seed1Inv.getSigmaZ2()) * mMatchParams->crudeNSigma2Cut[o2::track::kZ]; + const int side0 = leg0.tpcSide; + const int side1 = leg1.tpcSide; + const float sigZ2 = (leg0.getSigmaZ2() + seed1Inv.getSigmaZ2()) * mMatchParams->crudeNSigma2Cut[o2::track::kZ]; if (side0 == 0 || side1 == 0 || side0 == side1) { - float dZ = seed0.getZ() - seed1Inv.getZ(); + float dZ = leg0.getZ() - seed1Inv.getZ(); float dZTimeTol = 0.f; // the time of a non-TPC absolute leg is only known within its bracket (tRef is the bracket centre) if (side0 != 0 && side1 != 0) { // same side: the z offset is fixed by the difference of the reference times - dZ -= side0 * mTPCVDrift * (seed0.tRef - seed1.tRef); - } else if (side1 != 0) { // seed0 absolute: move seed1 to the time of seed0 - dZ -= side1 * mTPCVDrift * (seed0.tRef - seed1.tRef); - if (seed0.origID.getSource() != o2d::GlobalTrackID::TPC) { - dZTimeTol = 0.5f * mTPCVDrift * seed0.tBracket.delta(); + dZ -= side0 * mTPCVDrift * (leg0.tRef - leg1.tRef); + } else if (side1 != 0) { // leg0 absolute: move leg1 to the time of leg0 + dZ -= side1 * mTPCVDrift * (leg0.tRef - leg1.tRef); + if (leg0.origID.getSource() != o2d::GlobalTrackID::TPC) { + dZTimeTol = 0.5f * mTPCVDrift * leg0.tBracket.delta(); } - } else if (side0 != 0) { // seed1 absolute: move seed0 to the time of seed1 - dZ += side0 * mTPCVDrift * (seed1.tRef - seed0.tRef); - if (seed1.origID.getSource() != o2d::GlobalTrackID::TPC) { - dZTimeTol = 0.5f * mTPCVDrift * seed1.tBracket.delta(); + } else if (side0 != 0) { // leg1 absolute: move leg0 to the time of leg1 + dZ += side0 * mTPCVDrift * (leg1.tRef - leg0.tRef); + if (leg1.origID.getSource() != o2d::GlobalTrackID::TPC) { + dZTimeTol = 0.5f * mTPCVDrift * leg1.tBracket.delta(); } } const float dZTol = std::sqrt(sigZ2) + dZTimeTol; @@ -555,20 +574,20 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) break; } } else { // opposite sides - const float t = 0.5f * (side0 * (seed1Inv.getZ() - seed0.getZ()) / mTPCVDrift + seed0.tRef + seed1.tRef); + const float t = 0.5f * (side0 * (seed1Inv.getZ() - leg0.getZ()) / mTPCVDrift + leg0.tRef + leg1.tRef); const float tTol = std::sqrt(sigZ2) / (2.f * mTPCVDrift); - if (t < std::max(seed0.tBracket.getMin(), seed1.tBracket.getMin()) - tTol || t > std::min(seed0.tBracket.getMax(), seed1.tBracket.getMax()) + tTol) { + if (t < std::max(leg0.tBracket.getMin(), leg1.tBracket.getMin()) - tTol || t > std::min(leg0.tBracket.getMax(), leg1.tBracket.getMax()) + tTol) { rej = RejZ; break; } tCommon = t; - tCommonErr = std::sqrt(seed0.getSigmaZ2() + seed1Inv.getSigmaZ2()) / (2.f * mTPCVDrift); + tCommonErr = std::sqrt(leg0.getSigmaZ2() + seed1Inv.getSigmaZ2()) / (2.f * mTPCVDrift); } } // the z of a TPC-only leg refers to its own time0: it is tested above at a common time with constrainTPCOnlyZ, otherwise ignored if (!ignoreZ) { // both legs have an absolute z - auto dZ = seed0.getZ() - seed1Inv.getZ(); - if (dZ * dZ > (seed0.getSigmaZ2() + seed1Inv.getSigmaZ2()) * mMatchParams->crudeNSigma2Cut[o2::track::kZ]) { + auto dZ = leg0.getZ() - seed1Inv.getZ(); + if (dZ * dZ > (leg0.getSigmaZ2() + seed1Inv.getSigmaZ2()) * mMatchParams->crudeNSigma2Cut[o2::track::kZ]) { rej = RejZ; break; } @@ -580,7 +599,7 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) seed1Inv.setCov(0., o2::track::CovarMap[o2::track::kZ][o2::track::kQ2Pt]); } // calculate chi2 (expensive) - chi2 = seed0.getPredictedChi2(seed1Inv); + chi2 = leg0.getPredictedChi2(seed1Inv); if (chi2 > mMatchParams->crudeChi2Cut) { rej = RejChi2; break; @@ -594,12 +613,14 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) #ifdef _ALLOW_DEBUG_TREES_ if (mDBGOut && ((rej == Accept && isDebugFlag(MatchTreeAccOnly)) || isDebugFlag(MatchTreeAll))) { - auto seed1I = seed1; + const TrackSeed& dbgLeg0 = commonTime ? seed0Common : seed0; // the legs compared (refitted at the common time if they were) + auto seed1I = commonTime ? seed1Common : seed1; seed1I.invert(); - if (seed1I.rotate(seed0.getAlpha()) && o2::base::Propagator::Instance()->PropagateToXBxByBz(seed1I, seed0.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { + if (seed1I.rotate(dbgLeg0.getAlpha()) && o2::base::Propagator::Instance()->PropagateToXBxByBz(seed1I, dbgLeg0.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr)) { int rejI = int(rej); + int commonTimeI = commonTime; (*mDBGOut) << "match" - << "tf=" << mTFCount << "seed0=" << seed0 << "seed1=" << seed1I << "chi2Match=" << chi2 << "rej=" << rejI + << "tf=" << mTFCount << "seed0=" << dbgLeg0 << "seed1=" << seed1I << "chi2Match=" << chi2 << "rej=" << rejI << "commonTime=" << commonTimeI << "side0=" << int(seed0.tpcSide) << "side1=" << int(seed1.tpcSide) << "tCommon=" << tCommon << "tCommonErr=" << tCommonErr << "\n"; } } @@ -608,6 +629,38 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) return rej; } +//________________________________________________________ +bool MatchCosmics::refitSeedAtTime(const TrackSeed& seed, float timeMUS, TrackSeed& out) +{ + // TPC-only seed refitted with its clusters transformed at the time timeMUS, then treated as the seeds in createSeeds and process() (muon, + // extra z error, propagation to the DCA to the beam line); energy loss along the muon's flight in the refit and the propagation, as in + // refitWinners. Its z then refers to timeMUS + if (!mTPCRefitter || !mRecoData) { + return false; + } + const auto& tpcTrack = mRecoData->getTPCTrack(seed.origID); + o2::track::TrackParCov trk = tpcTrack.getParamOut(); + trk.setPID(o2::track::PID::Muon, true); + trk.resetCovariance(); + // the muon flies downward: inward along a leg pointing up (top leg) is along its flight (loss), along the bottom leg against it (gain) + std::array momentum{}; + trk.getPxPyPzGlo(momentum); + const int eLossSign = momentum[1] > 0.f ? ELossLoss : ELossGain; + if (mTPCRefitter->RefitTrackAsTrackParCov(trk, tpcTrack.getClusterRef(), timeMUS / mTPCTBinMUS, nullptr, false, true, eLossSign) < 0) { + return false; + } + trk.setCov(mMatchParams->tpcExtraZError2 + trk.getSigmaZ2(), o2::track::kSigZ2); + const o2::dataformats::VertexBase v; + if (!o2::base::Propagator::Instance()->propagateToDCABxByBz(v, trk, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, nullptr, eLossSign)) { + return false; + } + out = seed; + static_cast(out) = trk; + out.tRef = timeMUS; + mNRefitsCommonTime++; + return true; +} + //________________________________________________________ void MatchCosmics::registerMatch(int i, int j, float chi2, float tCommon, float tCommonErr) { From 9a493f4665d775cbaca844624c59b9ff5ae0fa88 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Thu, 8 Oct 2026 09:02:01 +0200 Subject: [PATCH 06/13] MatchCosmics: TOF flight pair of the legs as selection criterion New MatchCosmicsParams tofFlightSelection (default off, on in the preset physics-v1; needs TOF clusters, the workflow then adds them to the matcher's inputs): every pair of TPC-only legs accepted by checkPair is tested for a top / bottom TOF cluster pair along the outward continuations of its legs whose time difference matches the muon's flight between them along the helix (tofFlightTolerance), with the mean time in the pair's time window and both legs' y and z within tofRoad at that time (the test of the cluster collector). A confirmed pair takes the TOF time as its common time (error tofTimeError), so the winner is refitted at it, and is put in front of all unconfirmed partners of both legs (chi2 order within each group): the mutual-best selection prefers it, and a random pair can no longer take a leg of a TOF-confirmed cosmic. If the refit of a confirmed winner at the TOF time fails (in practice chi2Match > maxChi2Match), it is refitted at the time it has without its TOF flight pair (z-continuity time of legs on opposite sides, else the centre of the overlap of the legs' time brackets): the winner selection has already released the other candidates of both legs, so dropping the winner loses both legs. PbPb 568041: 24 TOF-tagged same-side cosmics, all real (legs on one line far from the beam line, TOF pair in the physical order without a reversed-order alternative), mostly on the C side nearer the readout at the TOF time; A-C cosmics show no TOF-TPC time offset there. The refit loop body is a lambda refitting a winner at a given time. The TOF candidates of a leg are searched once per TF, in the time-sorted TOF clusters of its time bracket and in the sector it points to and its neighbours. Debug tree "match": tofScore. Distorted cosmic MC (19 x 10 TFs, PbPb-like distortions, physics-v1 without the common-time refit, without -> with tofFlightSelection): matched muons 6848 -> 6871 of 8639, 0 wrong pairs both, 2134 confirmed pairs, 9 fallbacks. PbPb 568041 reconstructed from CTF (310 CTFs, 16541 TFs, vs. the same CTFs without the TOF selection, older code, before the fallback): cosmics 4782 -> 4760, TOF-tagged 681 -> 697; A-C: 11 untagged pairings replaced by 12 TOF-tagged ones. Matcher CPU 0.536 -> 0.545 s per TF. Co-Authored-By: Claude Opus 5.5 --- .../include/GlobalTracking/MatchCosmics.h | 26 +- .../GlobalTracking/MatchCosmicsParams.h | 4 + Detectors/GlobalTracking/src/MatchCosmics.cxx | 235 ++++++++++++++++-- .../GlobalTracking/src/MatchCosmicsParams.cxx | 6 +- .../CosmicsMatchingSpec.h | 2 +- .../src/CosmicsMatchingSpec.cxx | 5 +- .../src/cosmics-match-workflow.cxx | 6 +- 7 files changed, 254 insertions(+), 30 deletions(-) diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h index 86efff3416d82..b945f84864584 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h @@ -83,7 +83,20 @@ class MatchCosmics float chi2 = -1.f; ///< matching chi2 int next = MinusOne; ///< index of eventual next record float tCommon = 0.f; ///< common time [mus] fixed by z continuity of TPC-only legs on opposite TPC sides - float tCommonErr = -1.f; ///< its 1 sigma error [mus]; < 0: not fixed, the refit uses the centre of the time-bracket overlap and the cosmic's time error is the overlap's half-width + float tCommonErr = -1.f; ///< its 1 sigma error [mus]; < 0: not fixed, the refit uses the centre of the time-bracket overlap and the cosmic's time error is the overlap's half-width + float tofScore = -1.f; ///< score of the top / bottom TOF hit pair matching the muon's flight (tofFlightSelection; < 0: none); a pair with one wins against pairs without + float tCommonNoTOF = 0.f; ///< tCommon before a TOF flight pair replaced it: the refit falls back to it if the refit at the TOF time fails + float tCommonErrNoTOF = -1.f; ///< tCommonErr before a TOF flight pair replaced it + }; + + struct TOFCandidate { ///< TOF cluster along the outward continuation of a TPC-only seed + int index = -1; ///< index of the TOF cluster + double timeNS = 0.; ///< its time since the start of the TF [ns] + float dy = 0.f; ///< cluster - predicted y in the frame of the cluster's sector [cm] + float dz = 0.f; ///< cluster - predicted z, the leg's z taken at its own reference time tRef [cm] + float gx = 0.f; ///< global position of the cluster [cm] + float gy = 0.f; + float gz = 0.f; }; struct TrackSeed : public o2::track::TrackParCov { @@ -139,7 +152,10 @@ class MatchCosmics void updateTimeDependentParams(); RejFlag checkPair(int i, int j); bool refitSeedAtTime(const TrackSeed& seed, float timeMUS, TrackSeed& out); - void registerMatch(int i, int j, float chi2, float tCommon = 0.f, float tCommonErr = -1.f); + void registerMatch(int i, int j, float chi2, float tCommon = 0.f, float tCommonErr = -1.f, float tofScore = -1.f, float tCommonNoTOF = 0.f, float tCommonErrNoTOF = -1.f); + void prepareTOFClusters(const o2::globaltracking::RecoContainer& data); + const std::vector& getTOFCandidates(int iseed); + float findTOFFlightPair(int i, int j, float tMinMUS, float tMaxMUS, float& tofTimeMUS); void suppressMatch(int partner0, int partner1); void createSeeds(const o2::globaltracking::RecoContainer& data); bool validateMatch(int partner0); @@ -169,6 +185,12 @@ class MatchCosmics const o2::globaltracking::RecoContainer* mRecoData = nullptr; ///< inputs of the TF being processed o2::gpu::GPUO2InterfaceRefit* mTPCRefitter = nullptr; ///< TPC refitter of the TF being processed (owned by process()) size_t mNRefitsCommonTime = 0; ///< seeds refitted at the common time of a same-side pair in this TF + std::vector mTOFClusterOrder; ///< TOF clusters of the TF sorted in time (tofFlightSelection) + std::vector mTOFClusterTimeMUS; ///< their times since the start of the TF [mus], same order + std::vector> mSeedTOFCandidates; ///< TOF candidates per seed, filled on first use + std::vector mSeedTOFDone; ///< the TOF candidates of the seed are filled + size_t mNTOFConfirmed = 0; ///< accepted pairs with a TOF flight pair in this TF + size_t mNTOFFallbacks = 0; ///< TOF-confirmed winners refitted at their time without TOF in this TF std::vector mCosmicTracks; std::vector mCosmicTracksLbl; diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h index 8e51c33e96288..33121f0b323fc 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h @@ -42,6 +42,10 @@ struct MatchCosmicsParams : public o2::conf::ConfigurableParamHelpertofFlightSelection) { + prepareTOFClusters(data); + } // sort in time bracket lower edge, putting rejected tracks in the end std::vector sortID(ntr); @@ -142,12 +148,18 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) } } } + if (mMatchParams->tofFlightSelection) { + LOGP(info, "{} accepted pairs confirmed by a TOF flight pair", mNTOFConfirmed); + } selectWinners(); refitWinners(data); if (mNRefitsCommonTime) { LOGP(info, "{} seeds refitted at the common time of same-side pairs", mNRefitsCommonTime); } + if (mNTOFFallbacks) { + LOGP(info, "{} TOF-confirmed winners failed the refit at the TOF time and were refitted at their time without TOF", mNTOFFallbacks); + } mTPCRefitter = nullptr; mRecoData = nullptr; @@ -203,7 +215,8 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) return nclRefit == ncl ? ncl : -1; }; - for (auto winRID : mWinners) { + // refit the legs of winner winRID at the time t0 [mus] (error dt) and add the cosmic; false: a refit, propagation or cut failed + auto refitWinner = [&](int winRID, float t0, float dt) { const auto& rec = mRecords[winRID]; int poolEntryID[2] = {rec.id0, rec.id1}; o2::track::TrackParCov outerLegs[2] = {data.getTrackParamOut(mSeeds[rec.id0].origID), data.getTrackParamOut(mSeeds[rec.id1].origID)}; @@ -211,11 +224,6 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) leg.setPID(o2::track::PID::Muon, true); // as the seeds } auto tOverlap = mSeeds[rec.id0].tBracket.getOverlap(mSeeds[rec.id1].tBracket); - float t0 = tOverlap.mean(), dt = tOverlap.delta() * 0.5; - if (rec.tCommonErr >= 0.f) { // TPC-only legs on opposite sides: their z continuity fixes the time, refit both legs with it - t0 = rec.tCommon; - dt = rec.tCommonErr; - } auto pnt0 = outerLegs[0].getXYZGlo(), pnt1 = outerLegs[1].getXYZGlo(); int btm = 0, top = 1; // we fit topward from bottom @@ -245,7 +253,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, false, false, ELossGain); // inward refit, reset if (retVal < 0) { // refit failed LOG(debug) << "Inward refit of btm TPC track failed."; - continue; + return false; } nclTot += retVal; LOG(debug) << "chi2 after btm TPC refit with " << retVal << " clusters : " << chi2 << " orig.chi2 was " << tpcTrOrig.getChi2(); @@ -270,7 +278,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) if (!trCosm.rotate(mSeeds[poolEntryID[top]].getAlpha()) || !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, mSeeds[poolEntryID[top]].getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossGain)) { LOG(debug) << "Rotation/propagation of btm-track to top-track frame failed."; - continue; + return false; } // save bottom parameter at merging point auto trCosmBtm = trCosm; @@ -283,7 +291,7 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) if (gidxListTop[GTrackID::ITS].isIndexSet()) { auto nclfit = refitITSTrack(trCosm, gidxListTop[GTrackID::ITS], chi2, false, ELossGain); if (nclfit < 0) { - continue; + return false; } LOG(debug) << "chi2 after top ITS refit with " << nclfit << " clusters : " << chi2 << " orig.chi2 was " << data.getITSTrack(gidxListTop[GTrackID::ITS]).getChi2(); nclTot += nclfit; @@ -296,14 +304,14 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) if (!trCosm.getXatLabR(o2::constants::geom::XTPCInnerRef, xtogo, mBz, o2::track::DirOutward) || !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosm, xtogo, mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossGain)) { LOG(debug) << "Propagation to inner TPC boundary X=" << xtogo << " failed"; - continue; + return false; } } const auto& tpcTrOrig = data.getTPCTrack(gidxListTop[GTrackID::TPC]); int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosm, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2, true, false, ELossGain); // outward refit, no reset if (retVal < 0) { // refit failed LOG(debug) << "Outward refit of top TPC track failed."; - continue; + return false; } // outward refit in TPC LOG(debug) << "chi2 after top TPC refit with " << retVal << " clusters : " << chi2 << " orig.chi2 was " << tpcTrOrig.getChi2(); nclTot += retVal; @@ -317,14 +325,14 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) int retVal = tpcRefitter->RefitTrackAsTrackParCov(trCosmTop, tpcTrOrig.getClusterRef(), t0 * tpcTBinMUSInv, &chi2Dummy, false, true, ELossLoss); // inward refit, reset if (retVal < 0) { // refit failed LOG(debug) << "Inward refit of top TPC track failed."; - continue; + return false; } // inward refit in TPC } // is there ITS sub-track ? if (gidxListTop[GTrackID::ITS].isIndexSet()) { auto nclfit = refitITSTrack(trCosmTop, gidxListTop[GTrackID::ITS], chi2Dummy, true, ELossLoss); if (nclfit < 0) { - continue; + return false; } nclTot += nclfit; } // ITS refit @@ -332,23 +340,23 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) if (!trCosmTop.rotate(trCosmBtm.getAlpha()) || !o2::base::Propagator::Instance()->PropagateToXBxByBz(trCosmTop, trCosmBtm.getX(), mMatchParams->maxSnp, mMatchParams->maxStep, mMatchParams->matCorr, nullptr, ELossLoss)) { LOG(debug) << "Rotation/propagation of top-track to bottom-track frame failed."; - continue; + return false; } // calculate weighted average of 2 legs and chi2 o2::track::TrackParCov::MatrixDSym5 cov5; float chi2Match = trCosmBtm.getPredictedChi2(trCosmTop, cov5); if (mMatchParams->maxChi2Match >= 0.f && chi2Match > mMatchParams->maxChi2Match) { LOG(debug) << "Top/Bottom refitted legs disagree, chi2Match " << chi2Match; - continue; + return false; } if (!trCosmBtm.update(trCosmTop, cov5)) { LOG(debug) << "Top/Bottom update failed"; - continue; + return false; } // TPC-only legs on opposite sides: the legs' pT is required in checkPair, the refitted cosmic's here if (mSeeds[rec.id0].tpcSide * mSeeds[rec.id1].tpcSide < 0 && std::abs(trCosmBtm.getQ2Pt()) > mQ2PtCutoffOppositeSides) { LOG(debug) << "Cosmic with legs on opposite TPC sides below minPtOppositeSides"; - continue; + return false; } // create final track mCosmicTracks.emplace_back(mSeeds[poolEntryID[btm]].origID, mSeeds[poolEntryID[top]].origID, trCosmBtm, trCosmTop, chi2, chi2Match, nclTot, t0, dt); @@ -357,6 +365,23 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) auto& tlb = mCosmicTracksLbl.emplace_back((nclBtm > nclTot - nclBtm ? lbl[0] : lbl[1])); tlb.setFakeFlag(lbl[0] != lbl[1]); } + return true; + }; + for (auto winRID : mWinners) { + const auto& rec = mRecords[winRID]; + // refit at the common time if one is fixed (z continuity of TPC-only legs on opposite sides, TOF flight pair), else at the centre of + // the overlap of the legs' time brackets + auto refitAt = [&](float tCommon, float tCommonErr) { + if (tCommonErr >= 0.f) { + return refitWinner(winRID, tCommon, tCommonErr); + } + auto tOverlap = mSeeds[rec.id0].tBracket.getOverlap(mSeeds[rec.id1].tBracket); + return refitWinner(winRID, tOverlap.mean(), tOverlap.delta() * 0.5f); + }; + // a TOF-confirmed winner whose refit at the TOF time fails is refitted at the time it has without its TOF flight pair + if (!refitAt(rec.tCommon, rec.tCommonErr) && rec.tofScore >= 0.f && refitAt(rec.tCommonNoTOF, rec.tCommonErrNoTOF)) { + mNTOFFallbacks++; + } } LOG(info) << "Validated " << mCosmicTracks.size() << " top-bottom tracks in TF# " << mTFCount; } @@ -468,6 +493,7 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) float chi2 = 1.e9f; float tCommon = 0.f; // time fixed by z continuity of TPC-only legs on opposite sides, used by the refit float tCommonErr = -1.f; // its error (< 0: not fixed) + float tofScore = -1.f; // score of the TOF flight pair of an accepted pair (tofFlightSelection; < 0: none) TrackSeed seed0Common; // same-side TPC-only legs refitted at a common time (refitSameSideAtCommonTime) TrackSeed seed1Common; bool commonTime = false; @@ -605,8 +631,23 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) break; } rej = Accept; - registerMatch(i, j, chi2, tCommon, tCommonErr); - registerMatch(j, i, chi2, tCommon, tCommonErr); // the reverse reference can be also done in a separate loop + const float tCommonNoTOF = tCommon; // the pair's time without a TOF flight pair: the fallback of the refit at the TOF time + const float tCommonErrNoTOF = tCommonErr; + if (mMatchParams->tofFlightSelection) { // a top / bottom TOF hit pair with the muon's flight time confirms the pair and gives its time + const bool timeFixed = tCommonErr >= 0.f; + const auto overlap = seed0.tBracket.getOverlap(seed1.tBracket); + const float tMin = timeFixed ? tCommon - mMatchParams->nSigmaTError * tCommonErr : overlap.getMin(); + const float tMax = timeFixed ? tCommon + mMatchParams->nSigmaTError * tCommonErr : overlap.getMax(); + float tofTimeMUS = 0.f; + tofScore = findTOFFlightPair(i, j, tMin, tMax, tofTimeMUS); + if (tofScore >= 0.f) { + tCommon = tofTimeMUS; + tCommonErr = mMatchParams->tofTimeError; + mNTOFConfirmed++; + } + } + registerMatch(i, j, chi2, tCommon, tCommonErr, tofScore, tCommonNoTOF, tCommonErrNoTOF); + registerMatch(j, i, chi2, tCommon, tCommonErr, tofScore, tCommonNoTOF, tCommonErrNoTOF); // the reverse reference can be also done in a separate loop LOG(debug) << "Chi2 = " << chi2 << " NMatches " << mRecords.size(); break; } @@ -621,7 +662,7 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) int commonTimeI = commonTime; (*mDBGOut) << "match" << "tf=" << mTFCount << "seed0=" << dbgLeg0 << "seed1=" << seed1I << "chi2Match=" << chi2 << "rej=" << rejI << "commonTime=" << commonTimeI - << "side0=" << int(seed0.tpcSide) << "side1=" << int(seed1.tpcSide) << "tCommon=" << tCommon << "tCommonErr=" << tCommonErr << "\n"; + << "side0=" << int(seed0.tpcSide) << "side1=" << int(seed1.tpcSide) << "tCommon=" << tCommon << "tCommonErr=" << tCommonErr << "tofScore=" << tofScore << "\n"; } } #endif @@ -662,15 +703,161 @@ bool MatchCosmics::refitSeedAtTime(const TrackSeed& seed, float timeMUS, TrackSe } //________________________________________________________ -void MatchCosmics::registerMatch(int i, int j, float chi2, float tCommon, float tCommonErr) +void MatchCosmics::prepareTOFClusters(const o2::globaltracking::RecoContainer& data) +{ + // TOF clusters of the TF sorted in time; the TOF candidates of a seed are searched on its first use + const auto clusters = data.getTOFClusters(); + mTOFClusterOrder.resize(clusters.size()); + std::iota(mTOFClusterOrder.begin(), mTOFClusterOrder.end(), 0); + std::sort(mTOFClusterOrder.begin(), mTOFClusterOrder.end(), [&clusters](int a, int b) { return clusters[a].getTime() < clusters[b].getTime(); }); + mTOFClusterTimeMUS.resize(clusters.size()); + for (size_t k = 0; k < clusters.size(); k++) { + mTOFClusterTimeMUS[k] = clusters[mTOFClusterOrder[k]].getTime() * 1e-6; // [ps] since the start of the TF + } + mSeedTOFCandidates.clear(); + mSeedTOFCandidates.resize(mSeeds.size()); + mSeedTOFDone.assign(mSeeds.size(), false); +} + +//________________________________________________________ +const std::vector& MatchCosmics::getTOFCandidates(int iseed) +{ + // TOF clusters along the outward continuation of a TPC-only seed (from its outer parameters, in the sector it points to and its two + // neighbours) within its time bracket: |dy| < tofRoad, and dz within tofRoad of the drift of a one-side leg's z over the bracket + auto& candidates = mSeedTOFCandidates[iseed]; + if (mSeedTOFDone[iseed]) { + return candidates; + } + mSeedTOFDone[iseed] = true; + const auto& seed = mSeeds[iseed]; + if (seed.origID.getSource() != GTrackID::TPC) { + return candidates; + } + constexpr float MaxFlightMUS = 0.1f; // flight time of the muon between the TPC and the TOF, slow tails + constexpr float RadiusTOF = 380.f; // [cm], to find the sector the leg points to + const auto& tpcTrack = mRecoData->getTPCTrack(seed.origID); + const o2::track::TrackPar& parOut = tpcTrack.getParamOut(); + std::array xyz{}; + std::array dir{}; + parOut.getXYZGlo(xyz); + parOut.getPxPyPzGlo(dir); + // straight line from the outer parameters to the TOF radius (the sector only) + const float a = dir[0] * dir[0] + dir[1] * dir[1]; + const float b = xyz[0] * dir[0] + xyz[1] * dir[1]; + const float c = xyz[0] * xyz[0] + xyz[1] * xyz[1] - RadiusTOF * RadiusTOF; + const float disc = b * b - a * c; + if (a <= 0.f || disc < 0.f) { + return candidates; + } + const float step = (-b + std::sqrt(disc)) / a; + const int sectorCentre = o2::math_utils::angle2Sector(std::atan2(xyz[1] + step * dir[1], xyz[0] + step * dir[0])); + constexpr int NSectors = 18; + o2::track::TrackPar parSector[3]; + bool okSector[3] = {false, false, false}; + for (int k = 0; k < 3; k++) { + parSector[k] = parOut; + okSector[k] = parSector[k].rotateParam(o2::math_utils::sector2Angle((sectorCentre + k - 1 + NSectors) % NSectors)); + } + const float road = mMatchParams->tofRoad; + // range of z(t) - z(tRef) = side * vD * (t - tRef) over the bracket + const float dzDrift0 = seed.tpcSide * mTPCVDrift * (seed.tBracket.getMin() - seed.tRef); + const float dzDrift1 = seed.tpcSide * mTPCVDrift * (seed.tBracket.getMax() - seed.tRef); + const float dzMin = std::min(dzDrift0, dzDrift1) - road; + const float dzMax = std::max(dzDrift0, dzDrift1) + road; + const auto clusters = mRecoData->getTOFClusters(); + auto first = std::lower_bound(mTOFClusterTimeMUS.begin(), mTOFClusterTimeMUS.end(), seed.tBracket.getMin() - MaxFlightMUS); + for (auto it = first; it != mTOFClusterTimeMUS.end() && *it <= seed.tBracket.getMax() + MaxFlightMUS; ++it) { + const int index = mTOFClusterOrder[it - mTOFClusterTimeMUS.begin()]; + const auto& cl = clusters[index]; + const int k = (cl.getSector() - sectorCentre + NSectors + 1) % NSectors; // 0, 1, 2 for the sectors before, at and after the centre + if (k > 2 || !okSector[k]) { + continue; + } + float y = 0.f; + float z = 0.f; + if (!parSector[k].getYZAt(cl.getX(), mBz, y, z)) { + continue; + } + const float dy = cl.getY() - y; + const float dz = cl.getZ() - z; + if (std::abs(dy) > road || dz < dzMin || dz > dzMax) { + continue; + } + const float alpha = o2::math_utils::sector2Angle(cl.getSector()); + const float sinAlpha = std::sin(alpha); + const float cosAlpha = std::cos(alpha); + candidates.push_back(TOFCandidate{index, cl.getTime() * 1e-3, dy, dz, cl.getX() * cosAlpha - cl.getY() * sinAlpha, cl.getX() * sinAlpha + cl.getY() * cosAlpha, cl.getZ()}); + } + return candidates; +} + +//________________________________________________________ +float MatchCosmics::findTOFFlightPair(int i, int j, float tMinMUS, float tMaxMUS, float& tofTimeMUS) +{ + // best pair of TOF candidates of seeds i and j whose time difference matches the muon's flight between them along the helix (the higher + // hit first), with its mean time within [tMinMUS, tMaxMUS] and the z of both legs in the road at that time. Returns its score, the + // squared residuals in units of their cuts (< 0: no pair); tofTimeMUS is the mean time of the two hits + constexpr float MaxFlightMUS = 0.1f; + constexpr float CmPerNS = 29.9792458f; + const auto& candidates0 = getTOFCandidates(i); + const auto& candidates1 = getTOFCandidates(j); + if (candidates0.empty() || candidates1.empty()) { + return -1.f; + } + const auto& seed0 = mSeeds[i]; + const auto& seed1 = mSeeds[j]; + const float curvature = 0.5f * (std::abs(mRecoData->getTPCTrack(seed0.origID).getCurvature(mBz)) + std::abs(mRecoData->getTPCTrack(seed1.origID).getCurvature(mBz))); + const float road = mMatchParams->tofRoad; + const float tolerance = mMatchParams->tofFlightTolerance; + float bestScore = -1.f; + for (const auto& c0 : candidates0) { + for (const auto& c1 : candidates1) { + if (c0.index == c1.index) { + continue; + } + const auto& top = c0.gy > c1.gy ? c0 : c1; + const auto& bottom = c0.gy > c1.gy ? c1 : c0; + // flight path along the helix: arc in the transverse plane from the chord, then the dip + const float chordXY = std::hypot(top.gx - bottom.gx, top.gy - bottom.gy); + const float halfAngleSin = 0.5f * curvature * chordXY; + const float arcXY = halfAngleSin > 1e-4f && halfAngleSin < 1.f ? 2.f * std::asin(halfAngleSin) / curvature : chordXY; + const float length = std::hypot(arcXY, top.gz - bottom.gz); + const float flightDev = float(top.timeNS - bottom.timeNS) + length / CmPerNS; // the muon crosses the top TOF first + if (std::abs(flightDev) > tolerance) { + continue; + } + const float pairTimeMUS = float(0.5e-3 * (c0.timeNS + c1.timeNS)); + if (pairTimeMUS < tMinMUS - MaxFlightMUS || pairTimeMUS > tMaxMUS + MaxFlightMUS) { + continue; + } + const float dz0 = c0.dz - seed0.tpcSide * mTPCVDrift * (pairTimeMUS - seed0.tRef); + const float dz1 = c1.dz - seed1.tpcSide * mTPCVDrift * (pairTimeMUS - seed1.tRef); + if (std::abs(dz0) > road || std::abs(dz1) > road) { + continue; + } + const float score = (c0.dy * c0.dy + c1.dy * c1.dy + dz0 * dz0 + dz1 * dz1) / (road * road) + flightDev * flightDev / (tolerance * tolerance); + if (bestScore < 0.f || score < bestScore) { + bestScore = score; + tofTimeMUS = pairTimeMUS; + } + } + } + return bestScore; +} + +//________________________________________________________ +void MatchCosmics::registerMatch(int i, int j, float chi2, float tCommon, float tCommonErr, float tofScore, float tCommonNoTOF, float tCommonErrNoTOF) { - /// register track index j as a match for track index i + /// register track index j as a match for track index i; the matches of i are ordered in chi2, those confirmed by a TOF flight pair + /// (tofScore >= 0) in front of all others int newRef = mRecords.size(); - auto& matchRec = mRecords.emplace_back(MatchRecord{i, j, chi2, MinusOne, tCommon, tCommonErr}); + auto& matchRec = mRecords.emplace_back(MatchRecord{i, j, chi2, MinusOne, tCommon, tCommonErr, tofScore, tCommonNoTOF, tCommonErrNoTOF}); + const bool confirmed = tofScore >= 0.f; auto* best = &mSeeds[i].matchID; while (*best > MinusOne) { auto& oldMatchRec = mRecords[*best]; - if (oldMatchRec.chi2 > chi2) { // insert new match in front of the old one + const bool oldConfirmed = oldMatchRec.tofScore >= 0.f; + if ((confirmed && !oldConfirmed) || (confirmed == oldConfirmed && oldMatchRec.chi2 > chi2)) { // insert new match in front of the old one matchRec.next = *best; // new record will refer to the one it is superseding *best = newRef; // the reference on the superseded record should now refer to new one break; diff --git a/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx b/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx index fa746111ede0f..7886558f7fab0 100644 --- a/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx +++ b/Detectors/GlobalTracking/src/MatchCosmicsParams.cxx @@ -29,12 +29,16 @@ std::string o2::globaltracking::getMatchCosmicsPreset(const std::string& name) // - loose cut on the chi2 of the refitted legs; z test of TPC-only legs at a common time and same-half veto switched on // - TPC-only legs on opposite sides (time from z continuity, ~94 % of the PbPb candidates were random pairs): both legs and the // refitted cosmic above 2 GeV (offline: PbPb candidates / 10.5, cosmic MC efficiency 81.8 -> 80.6 %) + // - TOF flight pair as selection criterion (needs TOF clusters, added to the inputs): pairs confirmed by a top / bottom TOF hit pair + // with the muon's flight time win the selection and are refitted at the TOF time (568041: fewer cosmics, more of them TOF-tagged); + // the common-time refit of same-side legs (refitSameSideAtCommonTime) is not used: in PbPb it added mostly collision-track pairs static const std::map presets{ {"physics-v1", "cosmicsMatch.minSeedPt=1;cosmicsMatch.minSeedDCAxy=3;cosmicsMatch.minSeedDCAxyNSigma=10;cosmicsMatch.minSeedNClTPC=30;" "cosmicsMatch.crudeChi2Cut=50;cosmicsMatch.systSigma2[0]=0.25;cosmicsMatch.systSigma2[2]=4e-4;cosmicsMatch.systSigma2[4]=2.5e-3;" "cosmicsMatch.crudeNSigma2Cut[0]=144;cosmicsMatch.crudeNSigma2Cut[2]=144;cosmicsMatch.crudeNSigma2Cut[3]=9;cosmicsMatch.crudeNSigma2Cut[4]=144;" - "cosmicsMatch.maxChi2Match=1000;cosmicsMatch.constrainTPCOnlyZ=true;cosmicsMatch.vetoSameHalf=true;cosmicsMatch.minPtOppositeSides=2"}}; + "cosmicsMatch.maxChi2Match=1000;cosmicsMatch.constrainTPCOnlyZ=true;cosmicsMatch.vetoSameHalf=true;cosmicsMatch.minPtOppositeSides=2;" + "cosmicsMatch.tofFlightSelection=true"}}; auto it = presets.find(name); if (it == presets.end()) { std::string known; diff --git a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h index 1b1d9c494dc6c..4610118501a37 100644 --- a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h +++ b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h @@ -24,7 +24,7 @@ namespace globaltracking { /// create a processor spec -framework::DataProcessorSpec getCosmicsMatchingSpec(o2::dataformats::GlobalTrackID::mask_t src, bool usePV, bool useMC, bool itsStag); +framework::DataProcessorSpec getCosmicsMatchingSpec(o2::dataformats::GlobalTrackID::mask_t src, bool usePV, bool useMC, bool itsStag, bool useTOFClusters = false); } // namespace globaltracking } // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx index 20941d8333e27..9a7c0185b44a5 100644 --- a/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx @@ -183,7 +183,7 @@ void CosmicsMatchingSpec::endOfStream(EndOfStreamContext& ec) mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); } -DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool useMC, bool itsStag) +DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool useMC, bool itsStag, bool useTOFClusters) { std::vector outputs; Options opts{ @@ -195,6 +195,9 @@ DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool dataRequest->requestTracks(src, useMC); dataRequest->requestClusters(src, false); // no MC labels for clusters needed for refit only + if (useTOFClusters) { + dataRequest->requestTOFClusters(false); // TOF flight pairs of the candidate pairs (MatchCosmicsParams::tofFlightSelection) + } if (usePV) { dataRequest->requestPrimaryVertices(useMC); } diff --git a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx index 26817a71e42a1..d29e234a32a7e 100644 --- a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx @@ -118,12 +118,16 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) } GID::mask_t srcCl = src; + const bool useTOFClusters = o2::globaltracking::MatchCosmicsParams::Instance().tofFlightSelection; + if (useTOFClusters) { + srcCl |= GID::getSourceMask(GID::TOF); + } GID::mask_t dummy; if (!configcontext.options().get("disable-root-input")) { specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt)); } bool usePV = configcontext.options().get("use-pv-info"); - specs.emplace_back(o2::globaltracking::getCosmicsMatchingSpec(src, usePV, useMC, doStag)); + specs.emplace_back(o2::globaltracking::getCosmicsMatchingSpec(src, usePV, useMC, doStag, useTOFClusters)); bool clusterOutput = configcontext.options().get("enable-cluster-output"); if (clusterOutput) { if (!src[GID::TPC]) { From b1a3a03b5cb8ee8321b414f37c10bbbb4e9565d3 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Thu, 8 Oct 2026 16:19:54 +0200 Subject: [PATCH 07/13] calib-workflow.sh: optional cosmics matching with raw-cluster output in async With CALIB_ASYNC_EXTRACTCOSMICS=1 (set by O2DPG from the JDL), the async calibration workflows add o2-cosmics-match-workflow on the TPC-only tracks with the preset COSMICS_PRESET (default physics-v1) and the raw-cluster output (o2_cosmics_full.root, always written, as the TPC time series). The road detectors are those of ITS, TOF, TRD reconstructed in the job; without TOF the preset's TOF flight selection is switched off. Checked in print mode with the O2DPG counterpart: the switch adds only this workflow; no workflow without TPC, in sync mode, or with the switch unset. Co-Authored-By: Claude Opus 5.5 --- prodtests/full-system-test/calib-workflow.sh | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/prodtests/full-system-test/calib-workflow.sh b/prodtests/full-system-test/calib-workflow.sh index a14ff3b620d45..b772d5730fc5a 100644 --- a/prodtests/full-system-test/calib-workflow.sh +++ b/prodtests/full-system-test/calib-workflow.sh @@ -79,6 +79,17 @@ if [[ $CALIB_ASYNC_EXTRACTTIMESERIES == 1 ]] ; then CONFIG_TPCTIMESERIES+=" --mult-max ${TPCTIMESERIES_MULT_MAX}" add_W o2-tpc-time-series-workflow "$DISABLE_ROOT_INPUT ${CONFIG_TPCTIMESERIES}" fi +if [[ $CALIB_ASYNC_EXTRACTCOSMICS == 1 ]] ; then + # cosmic muons in collision data: TPC-only legs matched with the preset's selection, raw clusters of each cosmic -> o2_cosmics_full.root + : ${COSMICS_PRESET:=physics-v1} + COSMICS_ROAD_DETECTORS= + for det in ITS TOF TRD; do + has_detector_reco $det && COSMICS_ROAD_DETECTORS+="${COSMICS_ROAD_DETECTORS:+,}$det" + done + COSMICS_CONFIG= + has_detector_reco TOF || COSMICS_CONFIG="cosmicsMatch.tofFlightSelection=false" # the TOF flight selection of the preset needs TOF clusters + add_W o2-cosmics-match-workflow "$DISABLE_ROOT_INPUT $DISABLE_MC --track-sources TPC --cosmics-preset ${COSMICS_PRESET} --enable-cluster-output --road-detectors ${COSMICS_ROAD_DETECTORS:-none}" "$COSMICS_CONFIG" +fi # output-proxy for aggregator if workflow_has_parameter CALIB_PROXIES; then From 22bdf412668ff394789a5af7137059efaf730be3 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Thu, 8 Oct 2026 17:33:57 +0200 Subject: [PATCH 08/13] DataFormatsGlobalTracking: rename CosmicTrack to TrackCosmicsExtended The per-cosmic output of the cluster collector contains the matcher's TrackCosmics plus the raw hits of its legs; the new name makes that relation explicit (review comment). Header CosmicTrack.h is now TrackCosmicsExtended.h; content and branch names are unchanged. Co-Authored-By: Claude Opus 5.5 --- .../Detectors/GlobalTracking/CMakeLists.txt | 2 +- .../{CosmicTrack.h => TrackCosmicsExtended.h} | 16 +++++++-------- .../src/DataFormatsGlobalTrackingLinkDef.h | 4 ++-- .../src/CosmicsClusterCollectorSpec.cxx | 20 +++++++++---------- .../src/TrackCosmicsWriterSpec.cxx | 6 +++--- 5 files changed, 24 insertions(+), 24 deletions(-) rename DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/{CosmicTrack.h => TrackCosmicsExtended.h} (92%) diff --git a/DataFormats/Detectors/GlobalTracking/CMakeLists.txt b/DataFormats/Detectors/GlobalTracking/CMakeLists.txt index 53ce0ad335805..6e7c7025adb18 100644 --- a/DataFormats/Detectors/GlobalTracking/CMakeLists.txt +++ b/DataFormats/Detectors/GlobalTracking/CMakeLists.txt @@ -45,5 +45,5 @@ o2_target_root_dictionary( DataFormatsGlobalTracking HEADERS include/DataFormatsGlobalTracking/FilteredRecoTF.h include/DataFormatsGlobalTracking/TrackTuneParams.h - include/DataFormatsGlobalTracking/CosmicTrack.h + include/DataFormatsGlobalTracking/TrackCosmicsExtended.h ) diff --git a/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/CosmicTrack.h b/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/TrackCosmicsExtended.h similarity index 92% rename from DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/CosmicTrack.h rename to DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/TrackCosmicsExtended.h index ec7d85fc58793..dc97cf0ea7d09 100644 --- a/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/CosmicTrack.h +++ b/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/TrackCosmicsExtended.h @@ -9,11 +9,11 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file CosmicTrack.h +/// \file TrackCosmicsExtended.h /// \brief Matched cosmic track with the raw clusters of its legs and of the road around them, for offline refits -#ifndef ALICEO2_COSMIC_TRACK_H -#define ALICEO2_COSMIC_TRACK_H +#ifndef ALICEO2_TRACK_COSMICS_EXTENDED_H +#define ALICEO2_TRACK_COSMICS_EXTENDED_H #include #include @@ -32,8 +32,8 @@ struct CosmicTPCCluster { Attached = 0x1, ///< attached to the TPC track of this leg Corridor = 0x2, ///< found in the road around the leg Used = 0x4, ///< attached to some TPC track (for corridor clusters: another track, e.g. a split piece of the leg) - AbsTime = 0x8 ///< found on the other TPC side with the absolute time of the cosmic instead of the time0 of the leg: CosmicTrack::timeTOFMUS - ///< if >= 0, else the time of CosmicTrack::cosmic + AbsTime = 0x8 ///< found on the other TPC side with the absolute time of the cosmic instead of the time0 of the leg: TrackCosmicsExtended::timeTOFMUS + ///< if >= 0, else the time of TrackCosmicsExtended::cosmic }; o2::tpc::ClusterNative cl{}; ///< raw cluster: time, pad, widths, charges, flags uint8_t sector = 0; @@ -60,7 +60,7 @@ struct CosmicITSCluster { uint16_t row = 0; uint16_t col = 0; uint16_t pattID = 0; - int32_t pattEntry = -1; ///< start of the pattern bytes in CosmicTrack::itsPatterns (pattern not in the dictionary or group pattern), -1: none + int32_t pattEntry = -1; ///< start of the pattern bytes in TrackCosmicsExtended::itsPatterns (pattern not in the dictionary or group pattern), -1: none int32_t rofBC = 0; ///< start of the cluster's ROF in BCs since the start of the TF uint8_t leg = 0; ///< 0 bottom, 1 top uint8_t flags = 0; ///< CosmicHitFlags @@ -88,7 +88,7 @@ struct CosmicTRDTracklet { }; /// matched cosmic with everything needed for an offline refit -struct CosmicTrack { +struct TrackCosmicsExtended { o2::dataformats::TrackCosmics cosmic{}; ///< matcher output (time in mus; the leg references are only valid within the TF) o2::tpc::TrackTPC tpcBottom{}; ///< TPC part of the bottom leg (default if none); z refers to its time0; its cluster ///< references are only valid within the TF: the clusters are in clTPCBottom @@ -112,7 +112,7 @@ struct CosmicTrack { o2::MCCompLabel label{}; ///< MC label of the cosmic (MC only) bool hasTOFTime() const { return scoreTOFPair >= 0.f; } ///< a HitTOFFlight pair gave the cosmic its time timeTOFMUS - ClassDefNV(CosmicTrack, 1); + ClassDefNV(TrackCosmicsExtended, 1); }; /// per-TF quantities of the TPC transformation used in the reconstruction diff --git a/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h b/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h index 82e483febb9f1..06d54368484d0 100644 --- a/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h +++ b/DataFormats/Detectors/GlobalTracking/src/DataFormatsGlobalTrackingLinkDef.h @@ -30,8 +30,8 @@ #pragma link C++ class std::vector < o2::dataformats::CosmicTOFCluster> + ; #pragma link C++ class o2::dataformats::CosmicTRDTracklet + ; #pragma link C++ class std::vector < o2::dataformats::CosmicTRDTracklet> + ; -#pragma link C++ class o2::dataformats::CosmicTrack + ; -#pragma link C++ class std::vector < o2::dataformats::CosmicTrack> + ; +#pragma link C++ class o2::dataformats::TrackCosmicsExtended + ; +#pragma link C++ class std::vector < o2::dataformats::TrackCosmicsExtended> + ; #pragma link C++ class o2::dataformats::CosmicsTFInfo + ; #endif diff --git a/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx index 22b57fb59b0be..ab0caf86c4857 100644 --- a/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx @@ -64,7 +64,7 @@ #include "Framework/DeviceSpec.h" #include "GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h" #include "DataFormatsGlobalTracking/RecoContainer.h" -#include "DataFormatsGlobalTracking/CosmicTrack.h" +#include "DataFormatsGlobalTracking/TrackCosmicsExtended.h" #include "ReconstructionDataFormats/TrackCosmics.h" #include "DataFormatsTPC/TrackTPC.h" #include "DataFormatsTPC/ClusterNative.h" @@ -186,19 +186,19 @@ class CosmicsClusterCollectorSpec : public Task int index; ///< index of the cluster in the TF }; void addITS(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, int icosm, std::vector& requests, std::unordered_set& matched) const; - void fillITSPatterns(const RecoContainer& data, std::vector& requests, std::vector& cosmics) const; + void fillITSPatterns(const RecoContainer& data, std::vector& requests, std::vector& cosmics) const; void addTOF(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, std::unordered_set& matched) const; void addTRD(const RecoContainer& data, GTrackID gid, uint8_t leg, std::vector& out, std::unordered_set& matched) const; std::pair timeWindowMUS(const CosmicTime& cosmicTime) const; bool predictOutward(const o2::tpc::TrackTPC& leg, const CosmicTime& cosmicTime, int sector, float x, float& y, float& z) const; void roadTOF(const RecoContainer& data, const o2::tpc::TrackTPC* const* legs, const CosmicTime& cosmicTime, int icosm, std::vector& out, const std::unordered_set& matched, float& timeTOFMUS, float& scoreTOFPair, float& scoreTOFReversed) const; - void polish(const o2::tpc::TrackTPC* const* legs, float timeMUS, o2::gpu::GPUO2InterfaceRefit& refitter, o2::dataformats::CosmicTrack& out) const; + void polish(const o2::tpc::TrackTPC* const* legs, float timeMUS, o2::gpu::GPUO2InterfaceRefit& refitter, o2::dataformats::TrackCosmicsExtended& out) const; void roadTRD(const RecoContainer& data, const o2::tpc::TrackTPC* const* legs, const CosmicTime& cosmicTime, int icosm, std::vector& out, const std::unordered_set& matched) const; void roadITS(const RecoContainer& data, const o2::dataformats::TrackCosmics& cosm, const CosmicTime& cosmicTime, int legsSide, int icosm, std::vector& out, std::vector& requests, const std::unordered_set& matched) const; void cacheITSChipCentres(); - void writeDebug(const o2::dataformats::CosmicTrack& cosm, int icosm) const; - void flagDuplicates(std::vector& cosmics) const; + void writeDebug(const o2::dataformats::TrackCosmicsExtended& cosm, int icosm) const; + void flagDuplicates(std::vector& cosmics) const; void writeDebugTOF(const o2::tof::Cluster& c, int icosm, int leg, uint8_t flags) const; std::shared_ptr mDataRequest; @@ -303,7 +303,7 @@ void CosmicsClusterCollectorSpec::run(ProcessingContext& pc) tfInfo.vDrift = mCorrMap->getVDrift(); tfInfo.t0 = mCorrMap->getT0(); - std::vector cosmicsOut; + std::vector cosmicsOut; const auto cosmics = recoData.getCosmicTracks(); const size_t nCosmics = std::min(cosmics.size(), mMaxCosmicsPerTF); if (nCosmics < cosmics.size()) { @@ -1062,7 +1062,7 @@ void CosmicsClusterCollectorSpec::roadTOF(const RecoContainer& data, const o2::t } } -void CosmicsClusterCollectorSpec::polish(const o2::tpc::TrackTPC* const* legs, float timeMUS, o2::gpu::GPUO2InterfaceRefit& refitter, o2::dataformats::CosmicTrack& out) const +void CosmicsClusterCollectorSpec::polish(const o2::tpc::TrackTPC* const* legs, float timeMUS, o2::gpu::GPUO2InterfaceRefit& refitter, o2::dataformats::TrackCosmicsExtended& out) const { // refit of a cosmic with TPC-only legs at its TOF time, as MatchCosmics::refitWinners: the bottom leg inward, then to the closest approach // to the beam line; the top leg inward, then to the same point; the two halves combined. Muon mass, energy loss along the muon's flight @@ -1324,7 +1324,7 @@ void CosmicsClusterCollectorSpec::roadITS(const RecoContainer& data, const o2::d } } -void CosmicsClusterCollectorSpec::fillITSPatterns(const RecoContainer& data, std::vector& requests, std::vector& cosmics) const +void CosmicsClusterCollectorSpec::fillITSPatterns(const RecoContainer& data, std::vector& requests, std::vector& cosmics) const { // the TF's pattern stream holds, in cluster order, the patterns of the clusters with an invalid or a group pattern ID if (!mITSDict) { @@ -1376,7 +1376,7 @@ void CosmicsClusterCollectorSpec::addTRD(const RecoContainer& data, GTrackID gid } } -void CosmicsClusterCollectorSpec::flagDuplicates(std::vector& cosmics) const +void CosmicsClusterCollectorSpec::flagDuplicates(std::vector& cosmics) const { // the same muon can be matched twice, e.g. when a leg is split into two TPC tracks: the roads then collect largely the same clusters // (PbPb 567939: 10 of 47 cosmic pairs in a TF share 56-99 % of the clusters of the smaller one, all others none). The best one (TOF time, @@ -1421,7 +1421,7 @@ void CosmicsClusterCollectorSpec::flagDuplicates(std::vector const& cosmics) { + auto logger = [](std::vector const& cosmics) { LOG(info) << "Writing " << cosmics.size() << " cosmics with clusters"; }; return MakeRootTreeWriterSpec("cosmics-full-writer", @@ -63,7 +63,7 @@ DataProcessorSpec getCosmicsFullWriterSpec() "cosmicsFull", -1, // do not limit number of events to store 100, // periodically autosave - BranchDefinition>{InputSpec{"cosmics", "GLO", "COSMFULL", 0}, "cosmics", 1, logger}, + BranchDefinition>{InputSpec{"cosmics", "GLO", "COSMFULL", 0}, "cosmics", 1, logger}, BranchDefinition{InputSpec{"tfinfo", "GLO", "COSMFULLTF", 0}, "tfInfo", 1}, BranchDefinition{InputSpec{"tfid", "GLO", "COSMFULLTFID", 0}, "tfID", 1})(); } From 16be6b2dcb16cc930d92e31c6adf7a371b8cd7d5 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Thu, 8 Oct 2026 20:23:20 +0200 Subject: [PATCH 09/13] MatchCosmics: legs near the beam line for TOF-confirmed pairs only The seed cuts minSeedDCAxy / minSeedDCAxyNSigma (3 cm / 10 sigma in physics-v1) reject the legs of collision tracks, but also the cosmics crossing the ITS inner barrel (L0 / L1 entirely, L2 mostly). Loosening them floods the matcher with collision-track pairs (PbPb 567939, 600 TFs: 1 cm / 3 sigma gives 38x more candidates, TOF-tagged 0.1 % vs 14 %). New MatchCosmicsParams minSeedDCAxyTOF (default -1: off) and minSeedDCAxyNSigmaTOF (3): with tofFlightSelection, TPC-only legs failing the seed DCA cuts but passing these looser ones are kept, flagged nearBeam, and used in TOF-confirmed pairs only: a pair with such a leg and no TOF flight pair is rejected (RejNoTOF) before it is registered, so it cannot compete in the winner selection or reach the cluster collector. Log line: seeds near the beam line and their confirmed pairs per TF. The cosmics carry TrackCosmics flags NearBeamBottom / NearBeamTop (new member, class version 2) for legs kept only by the looser cuts, so these cosmics can be selected or dropped downstream without redoing the seed cuts. Co-Authored-By: Claude Opus 5.5 --- .../ReconstructionDataFormats/TrackCosmics.h | 13 ++++++- .../include/GlobalTracking/MatchCosmics.h | 6 +++- .../GlobalTracking/MatchCosmicsParams.h | 2 ++ Detectors/GlobalTracking/src/MatchCosmics.cxx | 36 ++++++++++++++++--- 4 files changed, 51 insertions(+), 6 deletions(-) diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackCosmics.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackCosmics.h index 79a34dc585876..427995b3319fe 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackCosmics.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackCosmics.h @@ -30,6 +30,12 @@ class TrackCosmics : public o2::track::TrackParCov using timeEst = o2::dataformats::TimeStampWithError; public: + enum Flags : uint8_t { + NearBeamBottom = 0x1, ///< bottom leg failed the standard seed DCA cuts (close to the beam line, absolutely or within its errors) and + ///< passed only the looser minSeedDCAxyTOF ones (TOF-confirmed pairs only) + NearBeamTop = 0x2 ///< the same for the top leg + }; + TrackCosmics() = default; ~TrackCosmics() = default; TrackCosmics(const TrackCosmics& src) = default; @@ -59,6 +65,10 @@ class TrackCosmics : public o2::track::TrackParCov int getNClusters() const { return mNClusters; } void setNClusters(int n) { mNClusters = n; } + uint8_t getFlags() const { return mFlags; } + void setFlags(uint8_t flags) { mFlags = flags; } + bool isNearBeam() const { return mFlags & (NearBeamBottom | NearBeamTop); } + o2::track::TrackParCov& getParamOut() { return mParamOut; } const o2::track::TrackParCov& getParamOut() const { return mParamOut; } @@ -72,8 +82,9 @@ class TrackCosmics : public o2::track::TrackParCov int mNClusters = 0; ///< total number of fitted clusters timeEst mTimeMUS; ///< time estimate in ns o2::track::TrackParCov mParamOut; ///< refitted outer parameter + uint8_t mFlags = 0; ///< Flags: NearBeamBottom, NearBeamTop - ClassDefNV(TrackCosmics, 1); + ClassDefNV(TrackCosmics, 2); }; } // namespace dataformats } // namespace o2 diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h index b945f84864584..a5d46989c5fff 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h @@ -71,7 +71,8 @@ class MatchCosmics RejProp, RejChi2, RejOther, - RejSameHalf + RejSameHalf, + RejNoTOF ///< a leg kept by minSeedDCAxyTOF only, and no TOF flight pair }; using InfoAccessor = o2d::AbstractRefAccessor; // there is no unique structure, so the default return type is dummy (int) @@ -108,6 +109,7 @@ class MatchCosmics float tRef = 0.f; ///< time [mus] the z of the parameters refers to (TPC-only: the TrackTPC time0; others: bracket centre) int8_t tpcSide = 0; ///< TPC-only seed with clusters on one side: +1 A, -1 C (z = z(t) - side*vD*(t-tRef)); 0: z absolute std::array xyzRef{}; ///< global position of the reference point before the propagation to the DCA (same-half veto) + bool nearBeam = false; ///< kept by the looser minSeedDCAxyTOF cuts only: usable in TOF-confirmed pairs only }; void setTPCCorrMaps(const o2::gpu::TPCFastTransformPOD* maph); void setTPCVDrift(const o2::tpc::VDriftCorrFact& v); @@ -191,6 +193,8 @@ class MatchCosmics std::vector mSeedTOFDone; ///< the TOF candidates of the seed are filled size_t mNTOFConfirmed = 0; ///< accepted pairs with a TOF flight pair in this TF size_t mNTOFFallbacks = 0; ///< TOF-confirmed winners refitted at their time without TOF in this TF + size_t mNSeedsNearBeam = 0; ///< TPC-only seeds kept for TOF-confirmed pairs only (minSeedDCAxyTOF) in this TF + size_t mNNearBeamConfirmed = 0; ///< accepted pairs with such a seed, confirmed by a TOF flight pair, in this TF std::vector mCosmicTracks; std::vector mCosmicTracksLbl; diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h index 33121f0b323fc..6a418f216de85 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h @@ -39,6 +39,8 @@ struct MatchCosmicsParams : public o2::conf::ConfigurableParamHelper= this [cm] (0: no cut; rejects collision tracks in physics data) float minSeedDCAxyNSigma = 0.f; // use only tracks with |DCA_xy| >= this * sigma(DCA_xy) (0: no cut; poorly measured collision tracks) + float minSeedDCAxyTOF = -1.f; // with tofFlightSelection: TPC-only legs failing the two cuts above but with |DCA_xy| >= this [cm] are used in TOF-confirmed pairs only (< 0: off; e.g. cosmics crossing the ITS inner barrel) + float minSeedDCAxyNSigmaTOF = 3.f; // the same for these legs in units of sigma(DCA_xy) bool constrainTPCOnlyZ = false; // TPC-only legs: test z at a common time (same side, or a leg with known time), else require the time implied by z continuity in both brackets bool vetoSameHalf = false; // reject pairs whose two legs lie on the same side of the closest approach (two pieces of one leg) bool refitSameSideAtCommonTime = false; // TPC-only legs on the same side: compare them refitted at the centre of their brackets' overlap instead of at their own time0s diff --git a/Detectors/GlobalTracking/src/MatchCosmics.cxx b/Detectors/GlobalTracking/src/MatchCosmics.cxx index 9667c1a08e728..695ef95e71940 100644 --- a/Detectors/GlobalTracking/src/MatchCosmics.cxx +++ b/Detectors/GlobalTracking/src/MatchCosmics.cxx @@ -65,6 +65,7 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) mCosmicTracksLbl.clear(); createSeeds(data); + mNSeedsNearBeam = 0; int ntr = mSeeds.size(); const auto prop = o2::base::Propagator::Instance(); // propagate to DCA to origin. A VertexBase (origin, zero covariance) selects the TrackParCov overload of propagateToDCABxByBz: @@ -84,9 +85,17 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) continue; } if (std::abs(trc.getY()) < mMatchParams->minSeedDCAxy || trc.getY() * trc.getY() < mMatchParams->minSeedDCAxyNSigma * mMatchParams->minSeedDCAxyNSigma * trc.getSigmaY2()) { - // passes close to the beam line, absolutely or within its errors: indistinguishable from collision tracks - trc.matchID = Reject; - continue; + // passes close to the beam line, absolutely or within its errors: indistinguishable from collision tracks, unless a TOF flight pair + // confirms the cosmic (TPC-only legs above the looser minSeedDCAxyTOF cuts, used in TOF-confirmed pairs only) + const bool nearBeam = mMatchParams->tofFlightSelection && mMatchParams->minSeedDCAxyTOF >= 0.f && trc.origID.getSource() == GTrackID::TPC && + std::abs(trc.getY()) >= mMatchParams->minSeedDCAxyTOF && + trc.getY() * trc.getY() >= mMatchParams->minSeedDCAxyNSigmaTOF * mMatchParams->minSeedDCAxyNSigmaTOF * trc.getSigmaY2(); + if (!nearBeam) { + trc.matchID = Reject; + continue; + } + trc.nearBeam = true; + mNSeedsNearBeam++; } if (mMatchParams->dcaCutChi2[trc.origID.getSource()] > 0.f && mUsePVInfo && trc.vtIDMin >= 0 && (std::abs(trc.getY()) < mMatchParams->fiducialRIP && std::abs(trc.getZ()) < mMatchParams->fiducialZIP)) { // do the propagation only if we are in the fiducial IP range. @@ -124,6 +133,7 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) mNRefitsCommonTime = 0; mNTOFConfirmed = 0; mNTOFFallbacks = 0; + mNNearBeamConfirmed = 0; if (mMatchParams->tofFlightSelection) { prepareTOFClusters(data); } @@ -151,6 +161,9 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) if (mMatchParams->tofFlightSelection) { LOGP(info, "{} accepted pairs confirmed by a TOF flight pair", mNTOFConfirmed); } + if (mNSeedsNearBeam) { + LOGP(info, "{} seeds near the beam line used in TOF-confirmed pairs only, {} such pairs confirmed", mNSeedsNearBeam, mNNearBeamConfirmed); + } selectWinners(); refitWinners(data); @@ -359,7 +372,15 @@ void MatchCosmics::refitWinners(const o2::globaltracking::RecoContainer& data) return false; } // create final track - mCosmicTracks.emplace_back(mSeeds[poolEntryID[btm]].origID, mSeeds[poolEntryID[top]].origID, trCosmBtm, trCosmTop, chi2, chi2Match, nclTot, t0, dt); + auto& cosmic = mCosmicTracks.emplace_back(mSeeds[poolEntryID[btm]].origID, mSeeds[poolEntryID[top]].origID, trCosmBtm, trCosmTop, chi2, chi2Match, nclTot, t0, dt); + uint8_t flags = 0; + if (mSeeds[poolEntryID[btm]].nearBeam) { + flags |= o2d::TrackCosmics::NearBeamBottom; + } + if (mSeeds[poolEntryID[top]].nearBeam) { + flags |= o2d::TrackCosmics::NearBeamTop; + } + cosmic.setFlags(flags); if (mUseMC) { o2::MCCompLabel lbl[2] = {data.getTrackMCLabel(mSeeds[poolEntryID[btm]].origID), data.getTrackMCLabel(mSeeds[poolEntryID[top]].origID)}; auto& tlb = mCosmicTracksLbl.emplace_back((nclBtm > nclTot - nclBtm ? lbl[0] : lbl[1])); @@ -646,6 +667,13 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) mNTOFConfirmed++; } } + if (seed0.nearBeam || seed1.nearBeam) { // a leg near the beam line (minSeedDCAxyTOF): TOF-confirmed pairs only + if (tofScore < 0.f) { + rej = RejNoTOF; + break; + } + mNNearBeamConfirmed++; + } registerMatch(i, j, chi2, tCommon, tCommonErr, tofScore, tCommonNoTOF, tCommonErrNoTOF); registerMatch(j, i, chi2, tCommon, tCommonErr, tofScore, tCommonNoTOF, tCommonErrNoTOF); // the reverse reference can be also done in a separate loop LOG(debug) << "Chi2 = " << chi2 << " NMatches " << mRecords.size(); From 3f1d9f8f05823372756d2c1b1bd81b0b90963ab4 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Fri, 9 Oct 2026 10:51:11 +0200 Subject: [PATCH 10/13] CosmicsClusterCollector: ITS road keeps the 5 closest clusters per inner-barrel layer Near the beam line the road in the inner barrel is dense with collision clusters, and the TPC prediction there (~mm) does not single out the cosmic's hit: keeping only the closest cluster per half of the cosmic and layer often stores a random one. In layers 0-2 the 5 closest are now kept, best first (the first is the previous choice); the outer barrel keeps the closest. The insertion rejects a NaN score, as the single-hit comparison did. Only cosmics with |DCA_xy| < 3.8 cm cross the inner barrel, so the output grows only for those. PbPb 568041, 400 CTFs with near-beam legs (minSeedDCAxyTOF=0): first road hit of every (cosmic, half, layer) identical to the previous code, outer barrel unchanged, output +1.8 %. On near-beam tracks fitted in the outer barrel, 24 of the 31 inner-barrel hits on the track are not the closest road hit (the ranking is dominated by the z uncertainty). Co-Authored-By: Claude Opus 5.5 --- .../TrackCosmicsExtended.h | 2 +- .../src/CosmicsClusterCollectorSpec.cxx | 57 ++++++++++++------- 2 files changed, 37 insertions(+), 22 deletions(-) diff --git a/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/TrackCosmicsExtended.h b/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/TrackCosmicsExtended.h index dc97cf0ea7d09..a6fbbcc969052 100644 --- a/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/TrackCosmicsExtended.h +++ b/DataFormats/Detectors/GlobalTracking/include/DataFormatsGlobalTracking/TrackCosmicsExtended.h @@ -95,7 +95,7 @@ struct TrackCosmicsExtended { o2::tpc::TrackTPC tpcTop{}; ///< TPC part of the top leg std::vector clTPCBottom; ///< attached + corridor TPC clusters of the bottom leg std::vector clTPCTop; ///< same for the top leg - std::vector clITS; ///< ITS clusters of the legs' matched tracks and on the road + std::vector clITS; ///< ITS clusters of the legs' matched tracks and on the road (per half and layer the closest, in the inner barrel up to 5, best first) std::vector itsPatterns; ///< pattern bytes (row span, column span, bitmap) of the ITS clusters that need them std::vector clTOF; ///< TOF clusters of the legs' matched tracks and on the road std::vector trdTracklets; ///< TRD tracklets of the legs' matched tracks and on the road diff --git a/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx index ab0caf86c4857..cd64ee4ce35d1 100644 --- a/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx @@ -33,8 +33,10 @@ /// Roads in the other detectors (--road-detectors): TRD tracklets close to the outward extrapolation of the legs and ITS clusters close /// to the trajectory near the beam line, within the time window of the cosmic (the TOF time if there is one, else the matcher's time, /// which is precise for TPC-only legs on opposite sides; otherwise the time window of the legs, with a correspondingly loose z cut). -/// They are flagged as found on the road; hits of the legs' matched global tracks are flagged as matched. Cosmics sharing >= 30 % of -/// their TPC clusters with a better one (a leg split into two TPC tracks) are flagged as duplicates, nothing is removed. +/// They are flagged as found on the road; hits of the legs' matched global tracks are flagged as matched. In the ITS the road keeps per +/// half of the cosmic and layer the cluster closest to the trajectory, in the inner barrel (dense with collision clusters near the beam +/// line) the 5 closest, best first. Cosmics sharing >= 30 % of their TPC clusters with a better one (a leg split into two TPC tracks) are +/// flagged as duplicates, nothing is removed. /// /// Polish: a cosmic with a TOF time and TPC-only legs is refitted at that time (as in the matcher: muon mass, energy loss along the /// flight); one-side legs then have their real z, hence the right material, which the matcher's TPC time cannot always give. @@ -1262,14 +1264,17 @@ void CosmicsClusterCollectorSpec::roadITS(const RecoContainer& data, const o2::d // the refitted cosmic has the z of its TPC time; with legs on one TPC side and a road time from the TOF its z moves by side * vD * dt const float zShift = legsSide * (cosmicTime.tb - cosm.getTimeMUS().getTimeStamp() / mTPCTBinMUS) * mCorrMap->getVDrift(); const float rofLengthMUS = o2::itsmft::DPLAlpideParam::Instance().roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingMUS; - // per half of the cosmic and layer the ITS cluster closest to the trajectory (the road near the beam line also contains collision clusters) + // per half of the cosmic and layer the ITS clusters closest to the trajectory, best first: one in the outer barrel, several in the inner + // barrel, where the road near the beam line is dense with collision clusters and the TPC prediction (~mm) does not single out the hit constexpr int NLayers = 7; + constexpr int NLayersIB = 3; + constexpr int MaxHitsIB = 5; struct Candidate { int index = -1; int rofBC = 0; float score = 1.f; }; - Candidate best[2][NLayers]; + Candidate best[2][NLayers][MaxHitsIB]; const auto clusters = data.getITSClusters(); auto geom = o2::its::GeometryTGeo::Instance(); for (const auto& rof : data.getITSClustersROFRecords()) { @@ -1297,28 +1302,38 @@ void CosmicsClusterCollectorSpec::roadITS(const RecoContainer& data, const o2::d const float score = std::max(dist2 / (mRoadITS * mRoadITS), normZ * normZ); const int half = global.Y() < pcaY ? 0 : 1; // bottom / top half of the cosmic const int layer = geom->getLayer(c.getSensorID()); - if (layer >= 0 && layer < NLayers && score < best[half][layer].score) { - best[half][layer] = {idx, rofBC, score}; + if (layer < 0 || layer >= NLayers) { + continue; + } + auto* candidates = best[half][layer]; + int slot = (layer < NLayersIB ? MaxHitsIB : 1) - 1; + if (!(score < candidates[slot].score)) { // also rejects a NaN score + continue; + } + for (; slot > 0 && score < candidates[slot - 1].score; slot--) { + candidates[slot] = candidates[slot - 1]; } + candidates[slot] = {idx, rofBC, score}; } } for (int half = 0; half < 2; half++) { for (int layer = 0; layer < NLayers; layer++) { - const auto& cand = best[half][layer]; - if (cand.index < 0) { - continue; - } - const auto& c = clusters[cand.index]; - auto& cl = out.emplace_back(); - cl.chipID = c.getSensorID(); - cl.row = c.getRow(); - cl.col = c.getCol(); - cl.pattID = c.getPatternID(); - cl.rofBC = cand.rofBC; - cl.leg = half; - cl.flags = o2::dataformats::HitRoad; - if (c.getPatternID() == o2::itsmft::CompCluster::InvalidPatternID || (mITSDict && mITSDict->isGroup(c.getPatternID()))) { - requests.push_back({icosm, int(out.size()) - 1, cand.index}); + for (const auto& candidate : best[half][layer]) { + if (candidate.index < 0) { + break; + } + const auto& c = clusters[candidate.index]; + auto& cl = out.emplace_back(); + cl.chipID = c.getSensorID(); + cl.row = c.getRow(); + cl.col = c.getCol(); + cl.pattID = c.getPatternID(); + cl.rofBC = candidate.rofBC; + cl.leg = half; + cl.flags = o2::dataformats::HitRoad; + if (c.getPatternID() == o2::itsmft::CompCluster::InvalidPatternID || (mITSDict && mITSDict->isGroup(c.getPatternID()))) { + requests.push_back({icosm, int(out.size()) - 1, candidate.index}); + } } } } From ff9d75ea818202d753b3dab765e9a04b08ba384b Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Fri, 9 Oct 2026 13:44:18 +0200 Subject: [PATCH 11/13] MatchCosmics: reject legs that are part of primary-vertex contributors Review suggestion: a TPC track that is part of a global track contributing to a primary vertex comes from a collision, not from a cosmic. The existing veto (--use-pv-info with cosmicsMatch.discardPVContributors, default on) compared the contributor's own track ID with the seeds, so it never fired for TPC-only legs: the contributors are ITS-TPC and higher. New workflow option --pv-contributor-sources (implies --use-pv-info, ignored with a warning if discardPVContributors is off): the contributors of these global track sources reject every leg they contain (RecoContainer::getSingleDetectorRefs), e.g. their TPC track. The sources and the parent matches needed to resolve them are loaded only for this; the legs still come from --track-sources only. Log line: seeds rejected per TF (also those that are themselves contributors). PbPb 567939, 2 TFs with vertices (94 each), sources ITS-TPC,ITS-TPC-TRD, ITS-TPC-TOF,ITS-TPC-TRD-TOF: 54 % of the TPC seeds rejected; without seed DCA cuts the collision-track pairs drop from 183 / 285 to 58 / 88 per TF, with the near-beam option the near-beam seeds by 57 %; matcher CPU for the 4 TFs (2 of them without vertices) 3.1 -> 2.3 s; the physics-v1 cosmic is kept. Co-Authored-By: Claude Opus 5.5 --- .../include/GlobalTracking/MatchCosmics.h | 5 ++++ Detectors/GlobalTracking/src/MatchCosmics.cxx | 25 +++++++++++++++-- .../CosmicsMatchingSpec.h | 6 ++++- .../src/CosmicsMatchingSpec.cxx | 27 ++++++++++++++++--- .../src/cosmics-match-workflow.cxx | 13 +++++++-- 5 files changed, 68 insertions(+), 8 deletions(-) diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h index a5d46989c5fff..82621288afccd 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h @@ -118,6 +118,8 @@ class MatchCosmics void process(const o2::globaltracking::RecoContainer& data); void setUseMC(bool mc) { mUseMC = mc; } void setUsePVInfo(bool v) { mUsePVInfo = v; } + void setSeedSources(GTrackID::mask_t src) { mSeedSources = src; } ///< track sources used as legs (other loaded ones resolve PV contributors) + void setPVVetoSources(GTrackID::mask_t src) { mPVVetoSources = src; } ///< sources whose PV contributors veto the legs they contain void init(); void end(); @@ -180,6 +182,8 @@ class MatchCosmics bool mFieldON = true; bool mUsePVInfo = false; bool mUseMC = true; + GTrackID::mask_t mSeedSources{GTrackID::MASK_ALL}; ///< track sources used as legs + GTrackID::mask_t mPVVetoSources{GTrackID::MASK_NONE}; ///< sources whose PV contributors veto the legs they contain float mITSROFrameLengthMUS = 0.; float mQ2PtCutoff = 1e9; float mQ2PtCutoffOppositeSides = 1e9; @@ -195,6 +199,7 @@ class MatchCosmics size_t mNTOFFallbacks = 0; ///< TOF-confirmed winners refitted at their time without TOF in this TF size_t mNSeedsNearBeam = 0; ///< TPC-only seeds kept for TOF-confirmed pairs only (minSeedDCAxyTOF) in this TF size_t mNNearBeamConfirmed = 0; ///< accepted pairs with such a seed, confirmed by a TOF flight pair, in this TF + size_t mNSeedsPVContributors = 0; ///< seeds rejected as part of a primary-vertex contributor in this TF std::vector mCosmicTracks; std::vector mCosmicTracksLbl; diff --git a/Detectors/GlobalTracking/src/MatchCosmics.cxx b/Detectors/GlobalTracking/src/MatchCosmics.cxx index 695ef95e71940..ec9b7481d6453 100644 --- a/Detectors/GlobalTracking/src/MatchCosmics.cxx +++ b/Detectors/GlobalTracking/src/MatchCosmics.cxx @@ -161,6 +161,9 @@ void MatchCosmics::process(const o2::globaltracking::RecoContainer& data) if (mMatchParams->tofFlightSelection) { LOGP(info, "{} accepted pairs confirmed by a TOF flight pair", mNTOFConfirmed); } + if (mNSeedsPVContributors) { + LOGP(info, "{} seeds rejected as (part of) a primary-vertex contributor", mNSeedsPVContributors); + } if (mNSeedsNearBeam) { LOGP(info, "{} seeds near the beam line used in TOF-confirmed pairs only, {} such pairs confirmed", mNSeedsNearBeam, mNNearBeamConfirmed); } @@ -940,8 +943,9 @@ void MatchCosmics::createSeeds(const o2::globaltracking::RecoContainer& data) } }; - data.createTracksVariadic(creator); + data.createTracksVariadic(creator, mSeedSources); // other loaded sources only resolve the primary-vertex contributors below + mNSeedsPVContributors = 0; if (mUsePVInfo) { // if needed, veto with the primary vertex info auto trackIndex = data.getPrimaryVertexMatchedTracks(); // Global ID's for associated tracks auto vtxRefs = data.getPrimaryVertexMatchedTrackRefs(); // references from vertex to these track IDs @@ -954,11 +958,28 @@ void MatchCosmics::createSeeds(const o2::globaltracking::RecoContainer& data) auto tvid = trackIndex[it]; auto entry = trackEntry.find(tvid); if (entry == trackEntry.end()) { + // a contributor of a veto source (not itself a leg): the legs it contains (e.g. its TPC track) come from a collision + if (mMatchParams->discardPVContributors && tvid.isPVContributor() && mPVVetoSources[tvid.getSource()] && data.isTrackSourceLoaded(tvid.getSource())) { + const auto refs = data.getSingleDetectorRefs(tvid); + for (int src = 0; src < GTrackID::NSources; src++) { + if (!mSeedSources[src] || !refs[src].isIndexSet()) { + continue; + } + auto partEntry = trackEntry.find(refs[src]); + if (partEntry != trackEntry.end() && mSeeds[partEntry->second].matchID != Reject) { + mSeeds[partEntry->second].matchID = Reject; + mNSeedsPVContributors++; + } + } + } continue; } auto& seed = mSeeds[entry->second]; - if (seed.matchID == Reject || (mMatchParams->discardPVContributors && tvid.isPVContributor())) { + if (seed.matchID != Reject && mMatchParams->discardPVContributors && tvid.isPVContributor()) { // the leg itself is a contributor seed.matchID = Reject; + mNSeedsPVContributors++; + } + if (seed.matchID == Reject) { continue; } if (seed.vtIDMin < 0) { diff --git a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h index 4610118501a37..e2dff219749ff 100644 --- a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h +++ b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h @@ -24,7 +24,11 @@ namespace globaltracking { /// create a processor spec -framework::DataProcessorSpec getCosmicsMatchingSpec(o2::dataformats::GlobalTrackID::mask_t src, bool usePV, bool useMC, bool itsStag, bool useTOFClusters = false); +/// srcPVContributors plus the parent matches needed to resolve their single-detector parts (RecoContainer::getSingleDetectorRefs) +o2::dataformats::GlobalTrackID::mask_t addPVContributorParents(o2::dataformats::GlobalTrackID::mask_t srcPVContributors); + +framework::DataProcessorSpec getCosmicsMatchingSpec(o2::dataformats::GlobalTrackID::mask_t src, bool usePV, bool useMC, bool itsStag, bool useTOFClusters = false, + o2::dataformats::GlobalTrackID::mask_t srcPVContributors = {}); } // namespace globaltracking } // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx index 9a7c0185b44a5..00d53ae083f5d 100644 --- a/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx @@ -66,7 +66,7 @@ namespace globaltracking class CosmicsMatchingSpec : public Task { public: - CosmicsMatchingSpec(std::shared_ptr dr, std::shared_ptr gr, bool usePV, bool useMC) : mDataRequest(dr), mGGCCDBRequest(gr), mUsePVInfo(usePV), mUseMC(useMC) {} + CosmicsMatchingSpec(std::shared_ptr dr, std::shared_ptr gr, GTrackID::mask_t src, GTrackID::mask_t srcPVVeto, bool usePV, bool useMC) : mDataRequest(dr), mGGCCDBRequest(gr), mSeedSources(src), mPVVetoSources(srcPVVeto), mUsePVInfo(usePV), mUseMC(useMC) {} ~CosmicsMatchingSpec() override = default; void init(InitContext& ic) final; void run(ProcessingContext& pc) final; @@ -81,6 +81,8 @@ class CosmicsMatchingSpec : public Task o2::tpc::VDriftHelper mTPCVDriftHelper{}; const o2::gpu::TPCFastTransformPOD* mCorrMap{nullptr}; o2::globaltracking::MatchCosmics mMatching; // matching engine + GTrackID::mask_t mSeedSources; // track sources used as legs + GTrackID::mask_t mPVVetoSources; // sources whose PV contributors veto the legs they contain bool mUseMC = true; bool mUsePVInfo = false; TStopwatch mTimer; @@ -94,6 +96,8 @@ void CosmicsMatchingSpec::init(InitContext& ic) mMatching.setDebugFlag(ic.options().get("debug-tree-flags")); mMatching.setUseMC(mUseMC); mMatching.setUsePVInfo(mUsePVInfo); + mMatching.setSeedSources(mSeedSources); + mMatching.setPVVetoSources(mPVVetoSources); // } @@ -183,7 +187,22 @@ void CosmicsMatchingSpec::endOfStream(EndOfStreamContext& ec) mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); } -DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool useMC, bool itsStag, bool useTOFClusters) +GTrackID::mask_t addPVContributorParents(GTrackID::mask_t srcPVContributors) +{ + auto src = srcPVContributors; + if (src[GTrackID::ITSTPCTRDTOF]) { + src |= GTrackID::getSourceMask(GTrackID::ITSTPCTRD); + } + if (src[GTrackID::TPCTRDTOF]) { + src |= GTrackID::getSourceMask(GTrackID::TPCTRD); + } + if (src[GTrackID::ITSTPCTRD] || src[GTrackID::ITSTPCTOF]) { + src |= GTrackID::getSourceMask(GTrackID::ITSTPC); + } + return src; +} + +DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool useMC, bool itsStag, bool useTOFClusters, GTrackID::mask_t srcPVContributors) { std::vector outputs; Options opts{ @@ -200,6 +219,8 @@ DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool } if (usePV) { dataRequest->requestPrimaryVertices(useMC); + // global tracks whose primary-vertex contributors veto their TPC track as a leg (MatchCosmicsParams::discardPVContributors), not legs + dataRequest->requestTracks(addPVContributorParents(srcPVContributors) & ~src, useMC); // same MC flag: some requests are shared with the legs } outputs.emplace_back("GLO", "COSMICTRC", 0, Lifetime::Timeframe); @@ -224,7 +245,7 @@ DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool "cosmics-matcher", dataRequest->inputs, outputs, - AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, usePV, useMC)}, + AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, src, usePV ? srcPVContributors : GTrackID::mask_t{}, usePV, useMC)}, opts}; } diff --git a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx index d29e234a32a7e..99d3783df64e6 100644 --- a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx @@ -54,6 +54,7 @@ void customize(std::vector& workflowOptions) {"disable-root-input", o2::framework::VariantType::Bool, false, {"disable root-files input reader"}}, {"disable-root-output", o2::framework::VariantType::Bool, false, {"disable root-files output writer"}}, {"use-pv-info", o2::framework::VariantType::Bool, false, {"request primary vertex for relevant cuts in the collision/cosmics interleaved data"}}, + {"pv-contributor-sources", VariantType::String, "", {"global track sources (with TPC) loaded only to reject the TPC track of a primary-vertex contributor as a leg (cosmicsMatch.discardPVContributors, implies --use-pv-info), e.g. ITS-TPC,ITS-TPC-TRD,ITS-TPC-TOF,ITS-TPC-TRD-TOF"}}, {"enable-cluster-output", o2::framework::VariantType::Bool, false, {"collect the raw clusters of the cosmics (legs + road around them) and write them with the cosmics"}}, {"road-detectors", VariantType::String, "ITS,TOF,TRD", {"with --enable-cluster-output: detectors whose hits along the cosmic are collected besides the TPC road"}}, {"cosmics-preset", VariantType::String, "", {"named set of cosmicsMatch settings applied before --configKeyValues (which can override single keys): physics-v1 = cosmics in collision data"}}, @@ -127,7 +128,15 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt)); } bool usePV = configcontext.options().get("use-pv-info"); - specs.emplace_back(o2::globaltracking::getCosmicsMatchingSpec(src, usePV, useMC, doStag, useTOFClusters)); + GID::mask_t srcPV = GID::getSourcesMask("ITS-TPC,TPC-TRD,TPC-TOF,ITS-TPC-TRD,TPC-TRD-TOF,ITS-TPC-TOF,ITS-TPC-TRD-TOF") & + GID::getSourcesMask(configcontext.options().get("pv-contributor-sources")); + if (srcPV.any() && !o2::globaltracking::MatchCosmicsParams::Instance().discardPVContributors) { + LOG(warning) << "--pv-contributor-sources ignored: cosmicsMatch.discardPVContributors is off"; + srcPV.reset(); + } + usePV |= srcPV.any(); + specs.emplace_back(o2::globaltracking::getCosmicsMatchingSpec(src, usePV, useMC, doStag, useTOFClusters, srcPV)); + const auto srcPVLoaded = o2::globaltracking::addPVContributorParents(srcPV); // with the parents needed to resolve the contributors bool clusterOutput = configcontext.options().get("enable-cluster-output"); if (clusterOutput) { if (!src[GID::TPC]) { @@ -142,7 +151,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) specs.emplace_back(o2::globaltracking::getCosmicsClusterCollectorSpec(src, useMC, doStag, roadDets)); } - o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCl, src, src, useMC, dummy); // clusters MC is not needed + o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCl, src | srcPVLoaded, src | srcPVLoaded, useMC, dummy); // clusters MC is not needed if (usePV) { o2::globaltracking::InputHelper::addInputSpecsPVertex(configcontext, specs, useMC); // P-vertex is always needed } From e95db4da12a3be9e38783583c3bece094493b213 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Fri, 9 Oct 2026 18:58:28 +0200 Subject: [PATCH 12/13] calib-workflow.sh: PV-contributor veto for the cosmics in async With CALIB_ASYNC_EXTRACTCOSMICS=1 the cosmics matcher gets --pv-contributor-sources $VERTEXING_SOURCES when the primary vertexing runs, as the time series gets $TRACK_SOURCES: the TPC part of a vertex contributor is no cosmic leg. COSMICS_PV_VETO=0 switches it off, COSMICS_PV_SOURCES overrides the sources. The option ignores sources without TPC, so the vertexing sources are passed as they are (help text). Checked in print mode (veto on / off, without ITS: TPC-TOF) and on 4 PbPb TFs with the full vertexing source list: the same seeds rejected as with the 4 ITS-TPC sources that have contributors there. Co-Authored-By: Claude Opus 5.5 --- .../GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx | 2 +- prodtests/full-system-test/calib-workflow.sh | 8 +++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx index 99d3783df64e6..997d0a26ba868 100644 --- a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx @@ -54,7 +54,7 @@ void customize(std::vector& workflowOptions) {"disable-root-input", o2::framework::VariantType::Bool, false, {"disable root-files input reader"}}, {"disable-root-output", o2::framework::VariantType::Bool, false, {"disable root-files output writer"}}, {"use-pv-info", o2::framework::VariantType::Bool, false, {"request primary vertex for relevant cuts in the collision/cosmics interleaved data"}}, - {"pv-contributor-sources", VariantType::String, "", {"global track sources (with TPC) loaded only to reject the TPC track of a primary-vertex contributor as a leg (cosmicsMatch.discardPVContributors, implies --use-pv-info), e.g. ITS-TPC,ITS-TPC-TRD,ITS-TPC-TOF,ITS-TPC-TRD-TOF"}}, + {"pv-contributor-sources", VariantType::String, "", {"global track sources loaded only to reject the TPC track of a primary-vertex contributor as a leg (cosmicsMatch.discardPVContributors, implies --use-pv-info); sources without TPC are ignored, so the vertexing sources can be passed as they are"}}, {"enable-cluster-output", o2::framework::VariantType::Bool, false, {"collect the raw clusters of the cosmics (legs + road around them) and write them with the cosmics"}}, {"road-detectors", VariantType::String, "ITS,TOF,TRD", {"with --enable-cluster-output: detectors whose hits along the cosmic are collected besides the TPC road"}}, {"cosmics-preset", VariantType::String, "", {"named set of cosmicsMatch settings applied before --configKeyValues (which can override single keys): physics-v1 = cosmics in collision data"}}, diff --git a/prodtests/full-system-test/calib-workflow.sh b/prodtests/full-system-test/calib-workflow.sh index b772d5730fc5a..a8b67e6939857 100644 --- a/prodtests/full-system-test/calib-workflow.sh +++ b/prodtests/full-system-test/calib-workflow.sh @@ -88,7 +88,13 @@ if [[ $CALIB_ASYNC_EXTRACTCOSMICS == 1 ]] ; then done COSMICS_CONFIG= has_detector_reco TOF || COSMICS_CONFIG="cosmicsMatch.tofFlightSelection=false" # the TOF flight selection of the preset needs TOF clusters - add_W o2-cosmics-match-workflow "$DISABLE_ROOT_INPUT $DISABLE_MC --track-sources TPC --cosmics-preset ${COSMICS_PRESET} --enable-cluster-output --road-detectors ${COSMICS_ROAD_DETECTORS:-none}" "$COSMICS_CONFIG" + # the TPC part of a primary-vertex contributor comes from a collision: veto it as a leg, with the vertexing sources (those without TPC are ignored) + COSMICS_OPT= + : ${COSMICS_PV_SOURCES:=${VERTEXING_SOURCES:-}} + if [[ ${COSMICS_PV_VETO:-1} == 1 ]] && [[ $BEAMTYPE != "cosmic" ]] && has_detector_matching PRIMVTX && [[ -n ${VERTEXING_SOURCES:-} ]] && [[ -n $COSMICS_PV_SOURCES ]]; then + COSMICS_OPT+=" --pv-contributor-sources $COSMICS_PV_SOURCES" + fi + add_W o2-cosmics-match-workflow "$DISABLE_ROOT_INPUT $DISABLE_MC --track-sources TPC --cosmics-preset ${COSMICS_PRESET} --enable-cluster-output --road-detectors ${COSMICS_ROAD_DETECTORS:-none}$COSMICS_OPT" "$COSMICS_CONFIG" fi # output-proxy for aggregator From 7803bcf1787ec1d628e3967ce9622e8f0367b9e7 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Fri, 9 Oct 2026 19:24:10 +0200 Subject: [PATCH 13/13] Cosmics: load the ITS tracks of legs with an ITS part With --track-sources containing ITS-TPC (or another source with ITS) but not ITS, the matcher's winner refit and the cluster collector read the ITS track of such a leg (getITSTrack, its cluster references), which was neither requested nor read: segfault in refitWinner, then in the collector; without an ITS road there was also no ITS cluster reader (invalid workflow). The ITS tracks and clusters are now requested and read for these legs, without making standalone ITS tracks legs (seed sources unchanged). TPC-only legs (physics-v1 in async) are not affected. PbPb 567939, 4 TFs, --track-sources TPC,ITS-TPC: all TFs processed with and without the ITS road (crash before), no ITS-only legs (seeds 23318 vs 26345 with ITS as legs); TPC-only legs: unchanged, no ITS track reader. Co-Authored-By: Claude Opus 5.5 --- .../include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h | 3 +++ .../src/CosmicsClusterCollectorSpec.cxx | 2 ++ .../GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx | 6 ++++++ .../GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx | 6 ++++-- 4 files changed, 15 insertions(+), 2 deletions(-) diff --git a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h index e2dff219749ff..a04c2e1f81000 100644 --- a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h +++ b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/CosmicsMatchingSpec.h @@ -24,6 +24,9 @@ namespace globaltracking { /// create a processor spec +/// ITS tracks needed by legs with an ITS part (their refit and ITS clusters) if ITS tracks are no legs themselves +o2::dataformats::GlobalTrackID::mask_t getLegITSSources(o2::dataformats::GlobalTrackID::mask_t src); + /// srcPVContributors plus the parent matches needed to resolve their single-detector parts (RecoContainer::getSingleDetectorRefs) o2::dataformats::GlobalTrackID::mask_t addPVContributorParents(o2::dataformats::GlobalTrackID::mask_t srcPVContributors); diff --git a/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx index cd64ee4ce35d1..c91318243ec0d 100644 --- a/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/CosmicsClusterCollectorSpec.cxx @@ -65,6 +65,7 @@ #include "Framework/DataProcessorSpec.h" #include "Framework/DeviceSpec.h" #include "GlobalTrackingWorkflow/CosmicsClusterCollectorSpec.h" +#include "GlobalTrackingWorkflow/CosmicsMatchingSpec.h" #include "DataFormatsGlobalTracking/RecoContainer.h" #include "DataFormatsGlobalTracking/TrackCosmicsExtended.h" #include "ReconstructionDataFormats/TrackCosmics.h" @@ -1610,6 +1611,7 @@ DataProcessorSpec getCosmicsClusterCollectorSpec(GTrackID::mask_t src, bool useM auto dataRequest = std::make_shared(); dataRequest->setITSPerLayer(itsStag); dataRequest->requestTracks(src, false); + dataRequest->requestTracks(getLegITSSources(src), false); // the ITS clusters of legs with an ITS part are read through their ITS track dataRequest->requestClusters(src, false); if (roadDets[DetID::ITS]) { dataRequest->requestITSClusters(false); diff --git a/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx index 00d53ae083f5d..6f3cf0529e218 100644 --- a/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/CosmicsMatchingSpec.cxx @@ -187,6 +187,11 @@ void CosmicsMatchingSpec::endOfStream(EndOfStreamContext& ec) mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); } +GTrackID::mask_t getLegITSSources(GTrackID::mask_t src) +{ + return GTrackID::includesDet(o2::detectors::DetID::ITS, src) ? GTrackID::getSourceMask(GTrackID::ITS) & ~src : GTrackID::mask_t{}; +} + GTrackID::mask_t addPVContributorParents(GTrackID::mask_t srcPVContributors) { auto src = srcPVContributors; @@ -213,6 +218,7 @@ DataProcessorSpec getCosmicsMatchingSpec(GTrackID::mask_t src, bool usePV, bool dataRequest->setITSPerLayer(itsStag); dataRequest->requestTracks(src, useMC); + dataRequest->requestTracks(getLegITSSources(src), false); // the refit of legs with an ITS part reads their ITS track dataRequest->requestClusters(src, false); // no MC labels for clusters needed for refit only if (useTOFClusters) { dataRequest->requestTOFClusters(false); // TOF flight pairs of the candidate pairs (MatchCosmicsParams::tofFlightSelection) diff --git a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx index 997d0a26ba868..a14cad38fb5d2 100644 --- a/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/cosmics-match-workflow.cxx @@ -118,7 +118,9 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) src = src | GID::getSourcesMask("CTP"); } - GID::mask_t srcCl = src; + // legs with an ITS part need the ITS tracks and clusters for their refit, also if standalone ITS tracks are no legs (not added to src) + const GID::mask_t srcITS = o2::globaltracking::getLegITSSources(src); + GID::mask_t srcCl = src | srcITS; const bool useTOFClusters = o2::globaltracking::MatchCosmicsParams::Instance().tofFlightSelection; if (useTOFClusters) { srcCl |= GID::getSourceMask(GID::TOF); @@ -151,7 +153,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) specs.emplace_back(o2::globaltracking::getCosmicsClusterCollectorSpec(src, useMC, doStag, roadDets)); } - o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCl, src | srcPVLoaded, src | srcPVLoaded, useMC, dummy); // clusters MC is not needed + o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCl, src | srcPVLoaded | srcITS, src | srcPVLoaded | srcITS, useMC, dummy); // clusters MC is not needed if (usePV) { o2::globaltracking::InputHelper::addInputSpecsPVertex(configcontext, specs, useMC); // P-vertex is always needed }