Skip to content

Commit

Permalink
account for TPC/ITS PID difference in propagation to ref.X
Browse files Browse the repository at this point in the history
  • Loading branch information
@shahor02
shahor02 committed Feb 7, 2024
1 parent 10b3d0e commit 56ff84c
Show file tree
Hide file tree
Showing 3 changed files with 114 additions and 33 deletions.
13 changes: 8 additions & 5 deletions Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITS.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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); }
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 2 additions & 0 deletions Detectors/GlobalTracking/include/GlobalTracking/MatchTPCITSParams.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,8 @@ struct MatchTPCITSParams : public o2::conf::ConfigurableParamHelper<MatchTPCITSP
float crudeNSigma2Cut[o2::track::kNParams] = {49.f, 49.f, 49.f, 49.f, 49.f};

float XMatchingRef = 70.f; ///< reference radius to propagate tracks for matching
float ITSStepEffFraction = 0.5; //< when correcting the ITS tracks for parameters difference between default PION and other PID hipothesis, use this fraction of propagated distance
float minBetaGammaForPIDDiff = 1.2; // account for difference between ITS and TPC PIDs used in propagation if TPC beta*gamma is below this

float minTPCTrackR = 50.; ///< cut on minimal TPC tracks radius to consider for matching, 666*pt_gev*B_kgaus/5
float minITSTrackR = 50.; ///< cut on minimal ITS tracks radius to consider for matching, 666*pt_gev*B_kgaus/5
Expand Down
132 changes: 104 additions & 28 deletions Detectors/GlobalTracking/src/MatchTPCITS.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,7 +220,7 @@ void MatchTPCITS::init()
// make sure T2GRot matrices are loaded into ITS geometry helper
o2::its::GeometryTGeo::Instance()->fillMatrixCache(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
Expand Down Expand Up @@ -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);
}
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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]);
}
Expand All @@ -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);
}
}
}
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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;
}
Expand All @@ -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;
}
Expand All @@ -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;
}
Expand All @@ -1097,17 +1118,26 @@ 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;
}
diff *= diff / (tITS.getDiagError2(o2::track::kQ2Pt) + tTPC.getDiagError2(o2::track::kQ2Pt));
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;
}

Expand Down Expand Up @@ -1206,19 +1236,23 @@ 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
// to their setctor edge that their matching should be checked also in the neighbouring 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]);
Expand All @@ -1229,15 +1263,15 @@ 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
constexpr float TgHalfSector = 0.17632698f;
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;
}
Expand Down Expand Up @@ -1422,10 +1456,30 @@ bool MatchTPCITS::refitTrackTPCITS(int iTPC, int& iITS, pmr::vector<o2::dataform
maxStep, MatCorrType::USEMatCorrNONE, nullptr, &trfit.getLTIntegralOut())) {
LOG(error) << "LTOF integral might be incorrect";
}
auto& tofL = trfit.getLTIntegralOut(); // this is TL integral calculated from the RefX to the DCA to the beamline, invert material integrals for outward propagation
tofL.setX2X0(-tofL.getX2X0());
tofL.setXRho(-tofL.getXRho());

// outward refit
#ifdef _ALLOW_DEBUG_TREES_
o2::track::TrackParCov itsRefAltPID;
#endif
auto& tracOut = trfit.getParamOut(); // this is a clone of ITS outward track already at the matching reference X
auto& tofL = trfit.getLTIntegralOut();
if (tTPC.getPID() > 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) ||
Expand Down Expand Up @@ -1510,14 +1564,36 @@ bool MatchTPCITS::refitTrackTPCITS(int iTPC, int& iITS, pmr::vector<o2::dataform
}
#ifdef _ALLOW_DEBUG_TREES_
if (mDBGOut) {
o2::track::TrackPar itsRefPIDCorr(tITS);
itsRefPIDCorr.setX(0);
if (tTPC.getPID() > 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

Expand Down

0 comments on commit 56ff84c

Please sign in to comment.