diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITS.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITS.h index db7b3934513c3..156321ba359a9 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITS.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITS.h @@ -113,7 +113,8 @@ enum TrackRejFlag : int { RejectOnTgl, RejectOnQ2Pt, RejectOnChi2, - NSigmaShift = 10 + NSigmaShift = 10, + RejectoOnPIDCorr = 20 }; ///< TPC track parameters propagated to reference X, with time bracket and index of @@ -139,7 +140,7 @@ struct TrackLocTPC : public o2::track::TrackParCov { return constraint == Constrained ? 0. : (constraint == ASide ? dt : -dt); } - ClassDefNV(TrackLocTPC, 1); + ClassDefNV(TrackLocTPC, 2); }; ///< ITS track outward parameters propagated to reference X, with time bracket and index of @@ -151,6 +152,8 @@ struct TrackLocITS : public o2::track::TrackParCov { int sourceID = 0; ///< track origin id int roFrame = MinusOne; ///< ITS readout frame assigned to this track int matchID = MinusOne; ///< entry (non if MinusOne) of its matchCand struct in the mMatchesITS + float xrho = 0; ///< x*rho seen during propagation to reference X (as pion) + float dL = 0; ///< distance integrated during propagation to reference X (as pion) bool hasCloneBefore() const { return getUserField() & CloneBefore; } bool hasCloneAfter() const { return getUserField() & CloneAfter; } int getCloneShift() const { return hasCloneBefore() ? -1 : (hasCloneAfter() ? 1 : 0); } @@ -479,8 +482,8 @@ class MatchTPCITS int compareTPCITSTracks(const TrackLocITS& tITS, const TrackLocTPC& tTPC, float& chi2) const; float getPredictedChi2NoZ(const o2::track::TrackParCov& trITS, const o2::track::TrackParCov& trTPC) const; - bool propagateToRefX(o2::track::TrackParCov& trc); - void addLastTrackCloneForNeighbourSector(int sector); + bool propagateToRefX(o2::track::TrackParCov& trc, o2::track::TrackLTIntegral* lti = nullptr); + void addLastTrackCloneForNeighbourSector(int sector, o2::track::TrackLTIntegral* trackLTInt = nullptr); ///------------------- manipulations with matches records ---------------------- bool registerMatchRecordTPC(int iITS, int iTPC, float chi2, int candIC = MinusOne); @@ -575,7 +578,7 @@ class MatchTPCITS float YMaxAtXMatchingRef = 999.; ///< max Y in the sector at reference X - float mSectEdgeMargin2 = 0.; ///< crude check if ITS track should be matched also in neighbouring sector + float mSectEdgeMargin = 0.; ///< crude check if ITS track should be matched also in neighbouring sector ///< safety margin in TPC time bins when estimating TPC track tMin and tMax from ///< assigned time0 and its track Z position (converted from mTPCTimeEdgeZSafeMargin) diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITSParams.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITSParams.h index a9d382cdbf18f..3ec189deff54b 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITSParams.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITSParams.h @@ -41,6 +41,8 @@ struct MatchTPCITSParams : public o2::conf::ConfigurableParamHelperfillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2GRot) | o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L)); - mSectEdgeMargin2 = mParams->crudeAbsDiffCut[o2::track::kY] * mParams->crudeAbsDiffCut[o2::track::kY]; ///< precalculated ^2 + mSectEdgeMargin = mParams->crudeAbsDiffCut[o2::track::kY] / std::sqrt(Cos70I2); #ifdef _ALLOW_DEBUG_TREES_ // debug streamer @@ -446,7 +446,6 @@ void MatchTPCITS::addTPCSeed(const o2::track::TrackParCov& _tr, float t0, float if (mMCTruthON) { mTPCLblWork.emplace_back(mTPCTrkLabels[tpcID]); } - // cache work track index mTPCSectIndexCache[o2::math_utils::angle2Sector(trc.getAlpha())].push_back(mTPCWork.size() - 1); } @@ -638,6 +637,8 @@ bool MatchTPCITS::prepareITSData() } long nHBF = o2::base::GRPGeomHelper::getNHBFPerTF(); long maxBCs = nHBF * long(o2::constants::lhc::LHCMaxBunches); + o2::track::TrackLTIntegral trackLTInt; + trackLTInt.setTimeNotNeeded(); for (int irof = 0; irof < nROFs; irof++) { const auto& rofRec = mITSTrackROFRec[irof]; @@ -685,10 +686,14 @@ bool MatchTPCITS::prepareITSData() continue; } // make sure the track is at the ref. radius - if (!propagateToRefX(trc)) { + trackLTInt.clearFast(); + if (!propagateToRefX(trc, &trackLTInt)) { mITSWork.pop_back(); // discard failed track continue; // add to cache only those ITS tracks which reached ref.X and have reasonable snp } + trc.xrho = trackLTInt.getXRho(); // we collect seen x*rho and distance to the reference X for further PID correcrions + trc.dL = trackLTInt.getL(); + if (mMCTruthON) { mITSLblWork.emplace_back(mITSTrkLabels[it]); } @@ -702,19 +707,16 @@ bool MatchTPCITS::prepareITSData() // when propagated to Xr (in this neighbouring sector) and the edge will be (neglecting the curvature) // [(Xr*tg(10)-Yr)/(tgPhir+tg70)]^2 / cos(70)^2 // for the next sector // [(Xr*tg(10)+Yr)/(tgPhir-tg70)]^2 / cos(70)^2 // for the prev sector - // Distances to the sector edges in neighbourings sectors (at Xref in theit proper frames) + // Distances to the sector edges in neighbourings sectors (at Xref in their proper frames) float trcY = trc.getY(), tgp = trc.getSnp(); tgp /= std::sqrt((1.f - tgp) * (1.f + tgp)); // tan of track direction XY - // sector up - float dy2Up = (YMaxAtXMatchingRef - trcY) / (tgp + Tan70); - if ((dy2Up * dy2Up * Cos70I2) < mSectEdgeMargin2) { // need to check this track for matching in sector up - addLastTrackCloneForNeighbourSector(sector < (o2::constants::math::NSectors - 1) ? sector + 1 : 0); - } - // sector down - float dy2Dn = (YMaxAtXMatchingRef + trcY) / (tgp - Tan70); - if ((dy2Dn * dy2Dn * Cos70I2) < mSectEdgeMargin2) { // need to check this track for matching in sector down - addLastTrackCloneForNeighbourSector(sector > 1 ? sector - 1 : o2::constants::math::NSectors - 1); + float dyUpDn[2] = {std::abs((YMaxAtXMatchingRef - trcY) / (tgp + Tan70)), std::abs((YMaxAtXMatchingRef + trcY) / (tgp - Tan70))}; // sector up, down edge distances + // we do the cloning for closest edge only + int sel = dyUpDn[0] < dyUpDn[1] ? 0 : 1; + if (dyUpDn[sel] < mSectEdgeMargin) { // need to check this track for matching in sector up or down + int sectNeib = sel == 0 ? (sector < (o2::constants::math::NSectors - 1) ? sector + 1 : 0) : (sector > 1 ? sector - 1 : o2::constants::math::NSectors - 1); + addLastTrackCloneForNeighbourSector(sectNeib, &trackLTInt); } } } @@ -1014,7 +1016,7 @@ bool MatchTPCITS::registerMatchRecordTPC(int iITS, int iTPC, float chi2, int can } //______________________________________________ -void MatchTPCITS::registerMatchRecordITS(int iITS, int iTPC, float chi2, int candIC) +void MatchTPCITS::registerMatchRecordITS(const int iITS, int iTPC, float chi2, int candIC) { ///< register TPC match in ITS tracks match records, ordering them in quality auto& tITS = mITSWork[iITS]; @@ -1068,7 +1070,26 @@ int MatchTPCITS::compareTPCITSTracks(const TrackLocITS& tITS, const TrackLocTPC& if ((rejFlag = roughCheckDif(diff, mParams->crudeNSigma2Cut[o2::track::kTgl], RejectOnTgl + NSigmaShift))) { return rejFlag; } - diff = tITS.getParam(o2::track::kY) - tTPC.getParam(o2::track::kY); + // do we need to account for different PID hypotheses used for ITS and TPC tracks propagation to ref. X? + bool testOtherPID = false; + float itsParam[5] = {tITS.getY(), tITS.getZ(), tITS.getSnp(), tITS.getTgl(), tITS.getQ2Pt()}; + if (tTPC.getPID() > tITS.getPID() && tITS.dL > 0.f && tTPC.getP2() / tTPC.getPID().getMass2() < mParams->minBetaGammaForPIDDiff) { + o2::track::TrackPar tPID(mITSTracksArray[tITS.sourceID].getParamOut()); // clone original ITS track at highest update point + tPID.setPID(tTPC.getPID(), true); + if (!tPID.correctForELoss(tITS.xrho)) { + return RejectoOnPIDCorr; + } + float dCurv = (tPID.getQ2Pt() - tITS.getQ2Pt()) * mBz * o2::constants::math::B2C, dLEff = tITS.dL * mParams->ITSStepEffFraction, dCurvL = dCurv * dLEff; + itsParam[o2::track::kQ2Pt] = tPID.getQ2Pt(); + itsParam[o2::track::kSnp] += dCurvL; + if (std::abs(itsParam[o2::track::kSnp]) >= 1.) { + itsParam[o2::track::kSnp] = std::copysign(0.99, itsParam[o2::track::kSnp]); + } + itsParam[o2::track::kY] += dCurvL * dLEff * 0.5; + testOtherPID = true; + } + + diff = itsParam[o2::track::kY] - tTPC.getParam(o2::track::kY); if ((rejFlag = roughCheckDif(diff, mParams->crudeAbsDiffCut[o2::track::kY], RejectOnY))) { return rejFlag; } @@ -1078,7 +1099,7 @@ int MatchTPCITS::compareTPCITSTracks(const TrackLocITS& tITS, const TrackLocTPC& } if (tTPC.constraint == TrackLocTPC::Constrained) { // in continuous only constrained tracks can be compared in Z - diff = tITS.getParam(o2::track::kZ) - tTPC.getParam(o2::track::kZ); + diff = itsParam[o2::track::kZ] - tTPC.getParam(o2::track::kZ); if ((rejFlag = roughCheckDif(diff, mParams->crudeAbsDiffCut[o2::track::kZ], RejectOnZ))) { return rejFlag; } @@ -1088,7 +1109,7 @@ int MatchTPCITS::compareTPCITSTracks(const TrackLocITS& tITS, const TrackLocTPC& } } - diff = tITS.getParam(o2::track::kSnp) - tTPC.getParam(o2::track::kSnp); + diff = itsParam[o2::track::kSnp] - tTPC.getParam(o2::track::kSnp); if ((rejFlag = roughCheckDif(diff, mParams->crudeAbsDiffCut[o2::track::kSnp], RejectOnSnp))) { return rejFlag; } @@ -1097,7 +1118,7 @@ int MatchTPCITS::compareTPCITSTracks(const TrackLocITS& tITS, const TrackLocTPC& return rejFlag; } - diff = tITS.getParam(o2::track::kQ2Pt) - tTPC.getParam(o2::track::kQ2Pt); + diff = itsParam[o2::track::kQ2Pt] - tTPC.getParam(o2::track::kQ2Pt); if ((rejFlag = roughCheckDif(diff, mParams->crudeAbsDiffCut[o2::track::kQ2Pt], RejectOnQ2Pt))) { return rejFlag; } @@ -1105,9 +1126,18 @@ int MatchTPCITS::compareTPCITSTracks(const TrackLocITS& tITS, const TrackLocTPC& if ((rejFlag = roughCheckDif(diff, mParams->crudeNSigma2Cut[o2::track::kQ2Pt], RejectOnQ2Pt + NSigmaShift))) { return rejFlag; } - // calculate mutual chi2 excluding Z in continuos mode - chi2 = getPredictedChi2NoZ(tITS, tTPC); - if (chi2 > mParams->cutMatchingChi2 || chi2 < 0.) { // sometime due to the numerical stability the chi2 is negative, reject it. + // calculate mutual chi2 excluding Z in continuous mode + if (testOtherPID) { // temporarily substitute pion params by alternative ones + auto tITSAlt = tITS; + tITSAlt.setPID(tTPC.getPID()); + tITSAlt.setParam(itsParam[o2::track::kY], o2::track::kY); + tITSAlt.setParam(itsParam[o2::track::kSnp], o2::track::kSnp); + tITSAlt.setParam(itsParam[o2::track::kQ2Pt], o2::track::kQ2Pt); + chi2 = getPredictedChi2NoZ(tITSAlt, tTPC); + } else { + chi2 = getPredictedChi2NoZ(tITS, tTPC); + } + if (chi2 > mParams->cutMatchingChi2 || chi2 < 0.) { // sometimes due to the numerical stability the chi2 is negative, reject it. return RejectOnChi2; } @@ -1206,7 +1236,7 @@ float MatchTPCITS::getPredictedChi2NoZ(const o2::track::TrackParCov& trITS, cons } //______________________________________________ -void MatchTPCITS::addLastTrackCloneForNeighbourSector(int sector) +void MatchTPCITS::addLastTrackCloneForNeighbourSector(int sector, o2::track::TrackLTIntegral* trackLTInt) { // add clone of the src ITS track cache, propagate it to ref.X in requested sector // and register its index in the sector cache. Used for ITS tracks which are so close @@ -1214,11 +1244,15 @@ void MatchTPCITS::addLastTrackCloneForNeighbourSector(int sector) mITSWork.push_back(mITSWork.back()); // clone the last track defined in given sector auto& trc = mITSWork.back(); if (trc.rotate(o2::math_utils::sector2Angle(sector)) && - o2::base::Propagator::Instance()->PropagateToXBxByBz(trc, mParams->XMatchingRef, MaxSnp, 2., MatCorrType::USEMatCorrNONE)) { + o2::base::Propagator::Instance()->PropagateToXBxByBz(trc, mParams->XMatchingRef, MaxSnp, 2., mUseMatCorrFlag, trackLTInt)) { // TODO: use faster prop here, no 3d field, materials mITSSectIndexCache[sector].push_back(mITSWork.size() - 1); // register track CLONE // flag clone mITSWork.back().setCloneBefore(); + if (trackLTInt) { + mITSWork.back().xrho = trackLTInt->getXRho(); // we collect seen x*rho and distance to the reference X for further PID correcrions + mITSWork.back().dL = trackLTInt->getL(); + } mITSWork[mITSWork.size() - 2].setCloneAfter(); if (mMCTruthON) { mITSLblWork.emplace_back(mITSTrkLabels[trc.sourceID]); @@ -1229,7 +1263,7 @@ void MatchTPCITS::addLastTrackCloneForNeighbourSector(int sector) } //______________________________________________ -bool MatchTPCITS::propagateToRefX(o2::track::TrackParCov& trc) +bool MatchTPCITS::propagateToRefX(o2::track::TrackParCov& trc, o2::track::TrackLTIntegral* lti) { // propagate track to matching reference X, making sure its assigned alpha // is consistent with TPC sector @@ -1237,7 +1271,7 @@ bool MatchTPCITS::propagateToRefX(o2::track::TrackParCov& trc) bool refReached = false; refReached = mParams->XMatchingRef < 10.; // RS: tmp, to cover XMatchingRef~0 int trialsLeft = 2; - while (o2::base::Propagator::Instance()->PropagateToXBxByBz(trc, mParams->XMatchingRef, MaxSnp, 2., mUseMatCorrFlag)) { + while (o2::base::Propagator::Instance()->PropagateToXBxByBz(trc, mParams->XMatchingRef, MaxSnp, 2., mUseMatCorrFlag, lti)) { if (refReached) { break; } @@ -1422,10 +1456,30 @@ bool MatchTPCITS::refitTrackTPCITS(int iTPC, int& iITS, pmr::vector tITS.getPID() && tTPC.getP2() / tTPC.getPID().getMass2() < mParams->minBetaGammaForPIDDiff) { + // in case the TPC track hypothesis is not pion, we redo the outward propagation to ref.x with TPC PID + tracOut = mITSTracksArray[tITS.sourceID].getParamOut(); + tracOut.setPID(tTPC.getPID(), true); + if (!tracOut.rotate(tTPC.getAlpha()) || !o2::base::Propagator::Instance()->PropagateToXBxByBz(tracOut, mParams->XMatchingRef, MaxSnp, 2., mUseMatCorrFlag)) { + LOGP(debug, "Failed to rotate ITSouter with imposed PID to TPC alpha {} or propagate to X={}: {:s}", tTPC.getAlpha(), mParams->XMatchingRef, tracOut.asString()); + matchedTracks.pop_back(); // destroy failed track + return false; + } +#ifdef _ALLOW_DEBUG_TREES_ + if (mDBGOut) { + itsRefAltPID = tracOut; + } +#endif + } { float xtogo = 0; if (!tracOut.getXatLabR(o2::constants::geom::XTPCInnerRef, xtogo, mBz, o2::track::DirOutward) || @@ -1510,14 +1564,36 @@ bool MatchTPCITS::refitTrackTPCITS(int iTPC, int& iITS, pmr::vector tITS.getPID() && tITS.dL > 0.f && tTPC.getP2() / tTPC.getPID().getMass2() < mParams->minBetaGammaForPIDDiff) { + itsRefPIDCorr = mITSTracksArray[tITS.sourceID].getParamOut(); // clone original ITS track at highest update point + itsRefPIDCorr.setPID(tTPC.getPID(), true); + if (!itsRefPIDCorr.correctForELoss(tITS.xrho)) { + itsRefPIDCorr.setX(-10); + } else { + float q2ptPID = itsRefPIDCorr.getQ2Pt(); + float dCurv = (q2ptPID - tITS.getQ2Pt()) * mBz * o2::constants::math::B2C, dLEff = tITS.dL * mParams->ITSStepEffFraction, dCurvL = dCurv * dLEff; + itsRefPIDCorr = tITS; + itsRefPIDCorr.setPID(tTPC.getPID(), true); + itsRefPIDCorr.setQ2Pt(q2ptPID); + auto snp = tITS.getSnp() + dCurvL; + if (std::abs(snp) >= 1.) { + snp = std::copysign(0.99, snp); + } + itsRefPIDCorr.setSnp(snp); + itsRefPIDCorr.setY(tITS.getY() + dCurvL * dLEff * 0.5); + } + } (*mDBGOut) << "refit" - << "tpcOrig=" << mTPCTracksArray[tTPC.sourceID] << "itsOrig=" << itsTrOrig << "itsRef=" << tITS << "tpcRef=" << tTPC << "matchRefit=" << trfit << "timeCorr=" << timeC << "dTimeFT0=" << minDiffFT0 << "dTimes=" << dtimes; + << "tpcOrig=" << mTPCTracksArray[tTPC.sourceID] << "itsOrig=" << itsTrOrig << "itsRef=" << tITS << "tpcRef=" << tTPC << "matchRefit=" << trfit << "timeCorr=" << timeC << "dTimeFT0=" << minDiffFT0 << "dTimes=" << dtimes + << "itsRefAltPID=" << itsRefAltPID << "itsRefPIDCorr=" << itsRefPIDCorr; if (mMCTruthON) { (*mDBGOut) << "refit" << "itsLbl=" << mITSLblWork[iITS] << "tpcLbl=" << mTPCLblWork[iTPC]; } (*mDBGOut) << "refit" - << "\n"; + << "tf=" << mTFCount << "\n"; } #endif