diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmics.h index acfc75f0679ba..b945f84864584 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 { @@ -82,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 { @@ -137,7 +151,11 @@ class MatchCosmics private: void updateTimeDependentParams(); RejFlag checkPair(int i, int j); - void registerMatch(int i, int j, float chi2, float tCommon = 0.f, float tCommonErr = -1.f); + 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, 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); @@ -164,6 +182,15 @@ 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 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 9307d92ee5e3d..33121f0b323fc 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchCosmicsParams.h @@ -38,9 +38,14 @@ 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 + bool tofFlightSelection = false; // needs TOF clusters: accepted pairs of TPC-only legs pointing to a top / bottom TOF hit pair with the muon's flight time win the selection, refit at that time (if that fails, at their time without TOF) + float tofRoad = 5.f; // half-width [cm] in y and z of the road at the TOF around the outward continuation of a TPC-only leg + float tofFlightTolerance = 2.f; // max. deviation [ns] of the top / bottom TOF time difference from the flight time along the helix + float tofTimeError = 0.1f; // error [mus] of the TOF time of a confirmed cosmic, for its refit and time window (covers TPC vs TOF offsets) 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..9667c1a08e728 100644 --- a/Detectors/GlobalTracking/src/MatchCosmics.cxx +++ b/Detectors/GlobalTracking/src/MatchCosmics.cxx @@ -33,6 +33,7 @@ #include "TPCReconstruction/TPCFastTransformHelperO2.h" #include "GlobalTracking/MatchTPCITS.h" #include "CommonConstants/GeomConstants.h" +#include "MathUtils/Utils.h" #include "DataFormatsTPC/WorkflowHelper.h" #include "DataFormatsTPC/VDriftCorrFact.h" #include "TPCFastTransformPOD.h" @@ -111,6 +112,22 @@ 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; + mNTOFConfirmed = 0; + mNTOFFallbacks = 0; + if (mMatchParams->tofFlightSelection) { + prepareTOFClusters(data); + } + // sort in time bracket lower edge, putting rejected tracks in the end std::vector sortID(ntr); std::iota(sortID.begin(), sortID.end(), 0); @@ -131,9 +148,20 @@ 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; mTFCount++; } @@ -144,16 +172,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 @@ -196,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)}; @@ -204,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 @@ -238,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(); @@ -263,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; @@ -276,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; @@ -289,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; @@ -310,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 @@ -325,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); @@ -350,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; } @@ -461,6 +493,10 @@ 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; // 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 +540,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 +600,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,27 +625,44 @@ 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; } 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; } #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 - << "side0=" << int(seed0.tpcSide) << "side1=" << int(seed1.tpcSide) << "tCommon=" << tCommon << "tCommonErr=" << tCommonErr << "\n"; + << "tf=" << mTFCount << "seed0=" << dbgLeg0 << "seed1=" << seed1I << "chi2Match=" << chi2 << "rej=" << rejI << "commonTime=" << commonTimeI + << "side0=" << int(seed0.tpcSide) << "side1=" << int(seed1.tpcSide) << "tCommon=" << tCommon << "tCommonErr=" << tCommonErr << "tofScore=" << tofScore << "\n"; } } #endif @@ -609,15 +671,193 @@ MatchCosmics::RejFlag MatchCosmics::checkPair(int i, int j) } //________________________________________________________ -void MatchCosmics::registerMatch(int i, int j, float chi2, float tCommon, float tCommonErr) +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::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]) {