14#include <fairlogger/Logger.h>
30#include <Math/SMatrix.h>
31#include <Math/SVector.h>
33#include <TGeoGlobalMagField.h>
61bool MatchTOF::mHasFillScheme =
false;
69 mFirstTForbit = firstTForbit;
79 updateTimeDependentParams();
81 mTimerMatchTPC.Reset();
82 mTimerMatchITSTPC.Reset();
85 mCalibInfoTOF.clear();
87 for (
int i = 0;
i < trkType::SIZEALL;
i++) {
88 mMatchedTracks[
i].clear();
89 mOutTOFLabels[
i].clear();
92 for (
int it = 0; it < trkType::SIZE; it++) {
94 mMatchedTracksIndex[sec][it].clear();
95 mTracksWork[sec][it].clear();
96 mTrackGid[sec][it].clear();
97 mLTinfos[sec][it].clear();
99 mTracksLblWork[sec][it].clear();
101 mTracksSectIndexCache[it][sec].clear();
102 mTracksSeed[it][sec].clear();
107 mSideTPC[sec].clear();
108 mExtraTPCFwdTime[sec].clear();
112 bool isPrepareTOFClusters = prepareTOFClusters();
114 LOGF(info,
"Timing prepareTOFCluster: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1);
116 if (!isPrepareTOFClusters) {
121 if (!prepareTPCData()) {
125 LOGF(info,
"Timing prepare TPC tracks: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1);
128 if (!prepareFITData()) {
132 LOGF(info,
"Timing prepare FIT data: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1);
135 std::array<uint32_t, 18> nMatches = {0};
137 mTimerMatchITSTPC.Start();
140 mMatchedTracksPairsSec[sec].clear();
145 if (mIsITSTPCused || mIsTPCTRDused || mIsITSTPCTRDused) {
147#pragma omp parallel for schedule(dynamic) num_threads(mNlanes)
153 mTimerMatchITSTPC.Stop();
155 mTimerMatchTPC.Start();
158#pragma omp parallel for schedule(dynamic) num_threads(mNlanes)
161 doMatchingForTPC(sec);
164 mTimerMatchTPC.Stop();
168 if (mStoreMatchable) {
171 for (
auto& matchingPair : mMatchedTracksPairsSec[sec]) {
172 int trkType = (
int)matchingPair.getTrackType();
173 int itrk = matchingPair.getIdLocal();
174 const auto& labelsTOF = mTOFClusLabels->
getLabels(matchingPair.getTOFClIndex());
175 const auto& labelTrack = mTracksLblWork[sec][
trkType][itrk];
179 for (
auto& lbl : labelsTOF) {
180 if (labelTrack == lbl) {
185 matchingPair.setFakeMatch();
189 for (
auto& matchingPair : mMatchedTracksPairsSec[sec]) {
190 int bct0 =
int((matchingPair.getSignal() - matchingPair.getLTIntegralOut().getTOF(0) + 5000) *
Geo::BC_TIME_INPS_INV);
195 if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(2)) < 600) {
196 }
else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(3)) < 600) {
197 }
else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(4)) < 600) {
198 }
else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(0)) < 600) {
199 }
else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(1)) < 600) {
200 }
else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(5)) < 600) {
201 }
else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(6)) < 600) {
202 }
else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(7)) < 600) {
203 }
else if (std::abs(tof - matchingPair.getLTIntegralOut().getTOF(8)) < 600) {
205 matchingPair.setFakeMatch();
211 LOG(
debug) <<
"...done. Now check the best matches";
212 nMatches[sec] = mMatchedTracksPairsSec[sec].size();
213 selectBestMatches(sec);
215 if (mStoreMatchable) {
216 for (
auto& matchingPair : mMatchedTracksPairsSec[sec]) {
217 trkType trkTypeSplitted = trkType::TPC;
218 auto sourceID = matchingPair.getTrackRef().getSource();
220 trkTypeSplitted = trkType::ITSTPC;
222 trkTypeSplitted = trkType::TPCTRD;
224 trkTypeSplitted = trkType::ITSTPCTRD;
226 matchingPair.setTrackType(trkTypeSplitted);
230 std::string nMatchesStr =
"Number of pairs matched per sector: ";
232 nMatchesStr += fmt::format(
"{} : {} ; ", sec, nMatches[sec]);
234 LOG(info) << nMatchesStr;
238 mIsITSTPCused =
false;
239 mIsTPCTRDused =
false;
240 mIsITSTPCTRDused =
false;
243 LOGF(info,
"Timing Do Matching: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1);
244 LOGF(info,
"Timing Do Matching Constrained: Cpu: %.3e s Real: %.3e s in %d slots", mTimerMatchITSTPC.CpuTime(), mTimerMatchITSTPC.RealTime(), mTimerMatchITSTPC.Counter() - 1);
245 LOGF(info,
"Timing Do Matching TPC : Cpu: %.3e s Real: %.3e s in %d slots", mTimerMatchTPC.CpuTime(), mTimerMatchTPC.RealTime(), mTimerMatchTPC.Counter() - 1);
251 mTPCVDrift =
v.refVDrift *
v.corrFact;
252 mTPCVDriftCorrFact =
v.corrFact;
253 mTPCVDriftRef =
v.refVDrift;
254 mTPCDriftTimeOffset =
v.getTimeOffset();
260 mTPCCorrMapsHelper = maph;
268 LOG(info) <<
"****** component for the matching of tracks to TOF clusters ******";
270 LOG(info) <<
"MC truth: " << (mMCTruthON ?
"on" :
"off");
271 LOG(info) <<
"Time tolerance: " << mTimeTolerance;
272 LOG(info) <<
"Space tolerance: " << mSpaceTolerance;
273 LOG(info) <<
"SigmaTimeCut: " << mSigmaTimeCut;
275 LOG(info) <<
"**********************************************************************";
283bool MatchTOF::prepareFITData()
293int MatchTOF::prepareInteractionTimes()
299void MatchTOF::printGrouping(
const std::vector<o2::dataformats::MatchInfoTOFReco>& origin,
const std::vector<std::vector<o2::dataformats::MatchInfoTOFReco>>& grouped)
301 printf(
"Original vector\n");
302 for (
const auto& seed : origin) {
303 printf(
"Pair: track=%d TOFcl=%d\n", seed.getIdLocal(), seed.getTOFClIndex());
306 printf(
"\nGroups\n");
308 for (
const auto& gr : grouped) {
310 printf(
"Group %d\n", ngroups);
311 for (
const auto& seed : gr) {
312 printf(
"Pair: track=%d TOFcl=%d\n", seed.getIdLocal(), seed.getTOFClIndex());
317void MatchTOF::groupingMatch(
const std::vector<o2::dataformats::MatchInfoTOFReco>& origin, std::vector<std::vector<o2::dataformats::MatchInfoTOFReco>>& grouped, std::vector<std::vector<int>>& firstEls, std::vector<std::vector<int>>& secondEls)
323 std::vector<o2::dataformats::MatchInfoTOFReco> dummy;
324 std::vector<int> dummyInt;
328 std::vector<int> alreadyUsed(origin.
size(), 0);
329 while (posInVector < origin.
size()) {
330 if (alreadyUsed[posInVector]) {
335 grouped.push_back(dummy);
336 auto& trkId = firstEls.emplace_back(dummyInt);
337 auto& cluId = secondEls.emplace_back(dummyInt);
340 const auto& seed = origin[posInVector];
341 trkId.push_back(seed.getIdLocal());
342 cluId.push_back(seed.getTOFClIndex());
344 grouped[
pos].push_back(seed);
345 alreadyUsed[posInVector] =
true;
350 for (
int i = posInVector;
i < origin.
size();
i++) {
351 if (alreadyUsed[
i]) {
354 const auto& seed = origin[
i];
356 int matchSecond = -1;
357 for (
const int& ind : trkId) {
358 if (seed.getIdLocal() == ind) {
363 for (
const int& ind : cluId) {
364 if (seed.getTOFClIndex() == ind) {
370 if (matchFirst >= 0 || matchSecond >= 0) {
371 if (matchFirst < 0) {
372 trkId.push_back(seed.getIdLocal());
374 if (matchSecond < 0) {
375 cluId.push_back(seed.getTOFClIndex());
377 grouped[
pos].push_back(seed);
378 alreadyUsed[
i] =
true;
389bool MatchTOF::prepareTPCData()
391 mNotPropagatedToTOF[trkType::UNCONS] = 0;
392 mNotPropagatedToTOF[trkType::CONSTR] = 0;
394 auto creator = [
this](
auto& trk,
GTrackID gid,
float time0,
float terr) {
395 const int nclustersMin = 0;
396 if constexpr (isTPCTrack<decltype(trk)>()) {
397 if (trk.getNClusters() < nclustersMin) {
401 if (std::abs(trk.getQ2Pt()) > mMaxInvPt) {
404 this->addTPCSeed(trk, gid, time0, terr);
406 if constexpr (isTPCITSTrack<decltype(trk)>()) {
410 this->addITSTPCSeed(trk, gid, time0, terr);
412 if constexpr (isTRDTrack<decltype(trk)>()) {
413 this->addTRDSeed(trk, gid, time0, terr);
420#pragma omp parallel for schedule(dynamic) num_threads(mNlanes)
424 propagateTPCTracks(sec);
426 if (mIsITSTPCused || mIsTPCTRDused || mIsITSTPCTRDused) {
427 propagateConstrTracks(sec);
432 std::array<float, 3> globalPos;
435 for (
auto& it : mTracksSeed[
trkType::UNCONS][sec]) {
436 auto& pair = mTracksWork[sec][trkType::UNCONS][it];
437 auto& trc = pair.first;
438 trc.getXYZGlo(globalPos);
441 int itnew = mTracksWork[sector][trkType::UNCONS].size();
443 mSideTPC[sector].push_back(mSideTPC[sec][it]);
444 mExtraTPCFwdTime[sector].push_back(mExtraTPCFwdTime[sec][it]);
445 mTracksWork[sector][trkType::UNCONS].emplace_back(pair);
446 mTrackGid[sector][trkType::UNCONS].emplace_back(mTrackGid[sec][trkType::UNCONS][it]);
449 mTracksLblWork[sector][trkType::UNCONS].emplace_back(mTracksLblWork[sec][trkType::UNCONS][it]);
451 mLTinfos[sector][trkType::UNCONS].emplace_back(mLTinfos[sec][trkType::UNCONS][it]);
452 mVZtpcOnly[sector].push_back(mVZtpcOnly[sec][it]);
453 mTracksSectIndexCache[trkType::UNCONS][sector].push_back(itnew);
457 for (
auto& it : mTracksSeed[
trkType::CONSTR][sec]) {
458 auto& pair = mTracksWork[sec][trkType::CONSTR][it];
459 auto& trc = pair.first;
460 trc.getXYZGlo(globalPos);
463 int itnew = mTracksWork[sector][trkType::CONSTR].size();
465 mTracksWork[sector][trkType::CONSTR].emplace_back(pair);
466 mTrackGid[sector][trkType::CONSTR].emplace_back(mTrackGid[sec][trkType::CONSTR][it]);
469 mTracksLblWork[sector][trkType::CONSTR].emplace_back(mTracksLblWork[sec][trkType::CONSTR][it]);
471 mLTinfos[sector][trkType::CONSTR].emplace_back(mLTinfos[sec][trkType::CONSTR][it]);
472 mTracksSectIndexCache[trkType::CONSTR][sector].push_back(itnew);
476 for (
int it = 0; it < trkType::SIZE; it++) {
478 mMatchedTracksIndex[sec][it].resize(mTracksWork[sec][it].
size());
479 std::fill(mMatchedTracksIndex[sec][it].
begin(), mMatchedTracksIndex[sec][it].
end(), -1);
484 LOG(
debug) <<
"Number of UNCONSTRAINED tracks that failed to be propagated to TOF = " << mNotPropagatedToTOF[trkType::UNCONS];
488#pragma omp parallel for schedule(dynamic) num_threads(mNlanes)
491 auto& indexCache = mTracksSectIndexCache[trkType::UNCONS][sec];
492 LOG(
debug) <<
"Sorting sector" << sec <<
" | " << indexCache.size() <<
" tracks";
493 if (!indexCache.size()) {
496 std::sort(indexCache.begin(), indexCache.end(), [
this, sec](
int a,
int b) {
497 auto& trcA = mTracksWork[sec][trkType::UNCONS][a].second;
498 auto& trcB = mTracksWork[sec][trkType::UNCONS][b].second;
499 return ((trcA.getTimeStamp() - trcA.getTimeStampError()) - (trcB.getTimeStamp() - trcB.getTimeStampError()) < 0.);
503 if (mIsITSTPCused || mIsTPCTRDused || mIsITSTPCTRDused) {
504 LOG(
debug) <<
"Number of CONSTRAINED tracks that failed to be propagated to TOF = " << mNotPropagatedToTOF[trkType::CONSTR];
508#pragma omp parallel for schedule(dynamic) num_threads(mNlanes)
511 auto& indexCache = mTracksSectIndexCache[trkType::CONSTR][sec];
512 LOG(
debug) <<
"Sorting sector" << sec <<
" | " << indexCache.size() <<
" tracks";
513 if (!indexCache.size()) {
516 std::sort(indexCache.begin(), indexCache.end(), [
this, sec](
int a,
int b) {
517 auto& trcA = mTracksWork[sec][trkType::CONSTR][a].second;
518 auto& trcB = mTracksWork[sec][trkType::CONSTR][b].second;
519 return ((trcA.getTimeStamp() - mSigmaTimeCut * trcA.getTimeStampError()) - (trcB.getTimeStamp() - mSigmaTimeCut * trcB.getTimeStampError()) < 0.);
536void MatchTOF::propagateTPCTracks(
int sec)
538 auto& trkWork = mTracksWork[sec][trkType::UNCONS];
540 for (
int it = 0; it < trkWork.size(); it++) {
546 if (!trackTune.sourceLevelTPC) {
547 if (trackTune.useTPCOuterCorr) {
548 trc.updateParams(trackTune.tpcParOuter);
554 if (!propagateToRefXWithoutCov(trc, mXRef, 10, mBz)) {
555 mNotPropagatedToTOF[trkType::UNCONS]++;
561 mNotPropagatedToTOF[trkType::UNCONS]++;
569 if (!propagateToRefX(trc, mXRef, 2, intLT0)) {
570 mNotPropagatedToTOF[trkType::UNCONS]++;
577 std::array<float, 3> globalPos;
578 trc.getXYZGlo(globalPos);
582 mTracksSectIndexCache[trkType::UNCONS][sector].push_back(it);
584 mTracksSeed[trkType::UNCONS][sec].push_back(it);
589void MatchTOF::propagateConstrTracks(
int sec)
591 std::array<float, 3> globalPos;
593 auto& trkWork = mTracksWork[sec][trkType::CONSTR];
595 for (
int it = 0; it < trkWork.size(); it++) {
600 if (!propagateToRefXWithoutCov(trc, mXRef, 2, mBz)) {
601 mNotPropagatedToTOF[trkType::CONSTR]++;
606 if (!propagateToRefX(trc, mXRef, 2, intLT0) || std::abs(trc.getZ()) >
Geo::MAXHZTOF) {
607 mNotPropagatedToTOF[trkType::CONSTR]++;
611 trc.getXYZGlo(globalPos);
616 mTracksSectIndexCache[trkType::CONSTR][sector].push_back(it);
618 mTracksSeed[trkType::CONSTR][sec].push_back(it);
625 mIsITSTPCused =
true;
627 auto trc = _tr.getParamOut();
631 addConstrainedSeed(trc, srcGID, intLT0, ts);
636 std::array<float, 3> globalPos;
637 trc.getXYZGlo(globalPos);
642 mTracksWork[sector][trkType::CONSTR].emplace_back(std::make_pair(trc, timeMUS));
643 mTrackGid[sector][trkType::CONSTR].emplace_back(srcGID);
644 mLTinfos[sector][trkType::CONSTR].emplace_back(intLT0);
655 mIsTPCTRDused =
true;
657 mIsITSTPCTRDused =
true;
662 auto trc = _tr.getOuterParam();
667 timeEst ts(time0, terr + mExtraTimeToleranceTRD);
669 addConstrainedSeed(trc, srcGID, intLT0, ts);
676 std::array<float, 3> globalPos;
686 float trackTime0 = _tr.getTime0() * mTPCTBinMUS;
687 timeInfo.setTimeStampError((_tr.getDeltaTBwd() + 5) * mTPCTBinMUS + extraErr);
689 timeInfo.setTimeStamp(trackTime0);
691 auto trc = _tr.getOuterParam();
692 trc.getXYZGlo(globalPos);
696 mSideTPC[sector].push_back(_tr.hasASideClustersOnly() ? 1 : (_tr.hasCSideClustersOnly() ? -1 : 0));
697 mExtraTPCFwdTime[sector].push_back((_tr.getDeltaTFwd() + 5) * mTPCTBinMUS + extraErr);
699 mTracksWork[sector][trkType::UNCONS].emplace_back(std::make_pair(trc, timeInfo));
700 mTrackGid[sector][trkType::UNCONS].emplace_back(srcGID);
703 mTracksLblWork[sector][trkType::UNCONS].emplace_back(mRecoCont->
getTPCTrackMCLabel(srcGID));
707 mLTinfos[sector][trkType::UNCONS].emplace_back(intLT0);
708 float vz0 = _tr.getZAt(0, mBz);
709 if (std::abs(vz0) > 9000) {
710 vz0 = _tr.getZ() - _tr.getX() * _tr.getTgl();
712 mVZtpcOnly[sector].push_back(vz0);
753bool MatchTOF::prepareTOFClusters()
757 mMCTruthON = mTOFClusLabels && mTOFClusLabels->
getNElements();
762 mTOFClusWork.clear();
770 mTOFClusSectIndexCache[sec].clear();
776 int nClusterInCurrentChunk = mTOFClustersArrayInp.size();
777 LOG(
debug) <<
"nClusterInCurrentChunk = " << nClusterInCurrentChunk;
778 mNumOfClusters += nClusterInCurrentChunk;
779 mTOFClusWork.reserve(mTOFClusWork.size() + mNumOfClusters);
780 for (
int it = 0; it < nClusterInCurrentChunk; it++) {
781 const Cluster& clOrig = mTOFClustersArrayInp[it];
783 mTOFClusWork.emplace_back(clOrig);
784 auto& cl = mTOFClusWork.back();
786 mTOFClusSectIndexCache[cl.getSector()].push_back(mTOFClusWork.size() - 1);
791 auto& indexCache = mTOFClusSectIndexCache[sec];
792 LOG(
debug) <<
"Sorting sector" << sec <<
" | " << indexCache.size() <<
" TOF clusters";
793 if (!indexCache.size()) {
796 std::sort(indexCache.begin(), indexCache.end(), [
this](
int a,
int b) {
797 auto& clA = mTOFClusWork[a];
798 auto& clB = mTOFClusWork[b];
799 return (clA.getTime() - clB.getTime()) < 0.;
803 if (mMatchedClustersIndex) {
804 delete[] mMatchedClustersIndex;
806 mMatchedClustersIndex =
new int[mNumOfClusters];
807 std::fill_n(mMatchedClustersIndex, mNumOfClusters, -1);
812void MatchTOF::doMatching(
int sec)
817 auto& cacheTOF = mTOFClusSectIndexCache[sec];
818 auto& cacheTrk = mTracksSectIndexCache[
type][sec];
819 int nTracks = cacheTrk.size(), nTOFCls = cacheTOF.size();
820 LOG(
debug) <<
"Matching sector " << sec <<
": number of tracks: " << nTracks <<
", number of TOF clusters: " << nTOFCls;
821 if (!nTracks || !nTOFCls) {
826 float deltaPos[2][3];
828 int nStepsInsideSameStrip[2] = {0, 0};
829 float deltaPosTemp[3];
830 std::array<float, 3>
pos;
831 std::array<float, 3> posBeforeProp;
836 LOG(
debug) <<
"Trying to match %d tracks" << cacheTrk.size();
837 for (
int itrk = 0; itrk < cacheTrk.size(); itrk++) {
838 for (
int ii = 0; ii < 2; ii++) {
840 nStepsInsideSameStrip[ii] = 0;
842 int nStripsCrossedInPropagation = 0;
843 auto& trackWork = mTracksWork[sec][
type][cacheTrk[itrk]];
844 auto& trefTrk = trackWork.first;
845 float pt = trefTrk.getPt();
846 auto& intLT = mLTinfos[sec][
type][cacheTrk[itrk]];
847 float timeShift = intLT.getL() * 33.35641;
851 float minTrkTime = (trackWork.second.getTimeStamp() - mSigmaTimeCut * trackWork.second.getTimeStampError()) * 1.E6 + timeShift;
852 float maxTrkTime = (trackWork.second.getTimeStamp() + mSigmaTimeCut * trackWork.second.getTimeStampError()) * 1.E6 + timeShift + 100E3;
853 const float sqrt12inv = 1. / sqrt(12.);
854 float resT = (trackWork.second.getTimeStampError() + 100E-3) * sqrt12inv;
868 for (
int ii = 0; ii < 2; ii++) {
869 for (
int iii = 0; iii < 5; iii++) {
874 int detIdTemp[5] = {-1, -1, -1, -1, -1};
876 double reachedPoint = mXRef + istep *
step;
878 while (propagateToRefX(trefTrk, reachedPoint, step, intLT) && nStripsCrossedInPropagation <= 2 && reachedPoint <
Geo::RMAX) {
881 trefTrk.getXYZGlo(
pos);
882 for (
int ii = 0; ii < 3; ii++) {
883 posFloat[ii] =
pos[ii];
893 for (
int idet = 0; idet < 5; idet++) {
894 detIdTemp[idet] = -1;
899 reachedPoint +=
step;
901 if (detIdTemp[2] == -1) {
923 if (nStripsCrossedInPropagation == 0 ||
924 (nStripsCrossedInPropagation >= 1 && (detId[nStripsCrossedInPropagation - 1][0] != detIdTemp[0] || detId[nStripsCrossedInPropagation - 1][1] != detIdTemp[1] || detId[nStripsCrossedInPropagation - 1][2] != detIdTemp[2]))) {
925 if (nStripsCrossedInPropagation == 0) {
926 LOG(
debug) <<
"We cross a strip for the first time";
928 if (nStripsCrossedInPropagation == 2) {
931 nStripsCrossedInPropagation++;
934 if (nStepsInsideSameStrip[nStripsCrossedInPropagation - 1] == 0) {
936 trkLTInt[nStripsCrossedInPropagation - 1] = intLT;
938 int detIdTemp2[5] = {0, 0, 0, 0, 0};
939 float deltaPosTemp2[3] = {deltaPosTemp[0], deltaPosTemp[1], deltaPosTemp[2]};
941 const int maxnstep = 50;
942 float xStart = trefTrk.getX();
943 float xStop = xStart;
944 trefTrk.getXYZGlo(
pos);
945 for (
int ii = 0; ii < 3; ii++) {
946 posFloat[ii] =
pos[ii];
948 while (deltaPosTemp2[1] < -0.05 && detIdTemp2[2] != -1 && nstep < maxnstep) {
951 propagateToRefXWithoutCov(trefTrk, xStop, 0.1, mBz, posFloat);
954 if (detIdTemp2[2] != -1) {
955 float dx = deltaPosTemp2[0] - deltaPosTemp[0];
956 float dy = deltaPosTemp2[1] - deltaPosTemp[1];
957 float dz = deltaPosTemp2[2] - deltaPosTemp[2];
958 updateTL(trkLTInt[nStripsCrossedInPropagation - 1], sqrt(dx * dx + dy * dy + dz * dz));
959 detIdTemp[0] = detIdTemp2[0];
960 detIdTemp[1] = detIdTemp2[1];
961 detIdTemp[2] = detIdTemp2[2];
962 detIdTemp[3] = detIdTemp2[3];
963 detIdTemp[4] = detIdTemp2[4];
964 deltaPosTemp[0] = deltaPosTemp2[0];
965 deltaPosTemp[1] = deltaPosTemp2[1];
966 deltaPosTemp[2] = deltaPosTemp2[2];
971 updateTL(trkLTInt[nStripsCrossedInPropagation - 1], -deltaPosTemp[1]);
973 detId[nStripsCrossedInPropagation - 1][0] = detIdTemp[0];
974 detId[nStripsCrossedInPropagation - 1][1] = detIdTemp[1];
975 detId[nStripsCrossedInPropagation - 1][2] = detIdTemp[2];
976 detId[nStripsCrossedInPropagation - 1][3] = detIdTemp[3];
977 detId[nStripsCrossedInPropagation - 1][4] = detIdTemp[4];
978 deltaPos[nStripsCrossedInPropagation - 1][0] = deltaPosTemp[0];
979 deltaPos[nStripsCrossedInPropagation - 1][1] = deltaPosTemp[1];
980 deltaPos[nStripsCrossedInPropagation - 1][2] = deltaPosTemp[2];
982 nStepsInsideSameStrip[nStripsCrossedInPropagation - 1]++;
995 for (Int_t imatch = 0; imatch < nStripsCrossedInPropagation; imatch++) {
997 deltaPos[imatch][0] /= nStepsInsideSameStrip[imatch];
998 deltaPos[imatch][1] /= nStepsInsideSameStrip[imatch];
999 deltaPos[imatch][2] /= nStepsInsideSameStrip[imatch];
1002 if (nStripsCrossedInPropagation == 0) {
1005 bool foundCluster =
false;
1006 for (
auto itof = itof0; itof < nTOFCls; itof++) {
1008 auto& trefTOF = mTOFClusWork[cacheTOF[itof]];
1012 if (trefTOF.getTime() < minTrkTime) {
1017 if (trefTOF.getTime() > maxTrkTime) {
1021 int mainChannel = trefTOF.getMainContributingChannel();
1028 float posCorr[3] = {0, 0, 0};
1055 float ndifInv = 1. / ndigits;
1057 posCorr[0] *= ndifInv;
1058 posCorr[1] *= ndifInv;
1059 posCorr[2] *= ndifInv;
1065 for (
auto iPropagation = 0; iPropagation < nStripsCrossedInPropagation; iPropagation++) {
1071 if (tof - trkLTInt[iPropagation].getTOF(6) > 2000) {
1075 const float errXinvMin = 1. / (mMatchParams->
maxResX * mMatchParams->
maxResX);
1076 const float errZinvMin = 1. / (mMatchParams->
maxResZ * mMatchParams->
maxResZ);
1077 float errXinv2 = 1. / (trefTrk.getSigmaY2());
1078 float errZinv2 = 1. / (trefTrk.getSigmaZ2() * cosangle);
1080 if (errXinv2 < errXinvMin) {
1081 errXinv2 = errXinvMin;
1083 if (errZinv2 < errZinvMin) {
1084 errZinv2 = errZinvMin;
1088 LOG(
debug) <<
"Propagated Track [" << itrk <<
"]: detId[" << iPropagation <<
"] = " << detId[iPropagation][0] <<
", " << detId[iPropagation][1] <<
", " << detId[iPropagation][2] <<
", " << detId[iPropagation][3] <<
", " << detId[iPropagation][4];
1089 float resX = deltaPos[iPropagation][0] - (
indices[4] - detId[iPropagation][4]) *
Geo::XPAD + posCorr[0];
1090 float resZ = deltaPos[iPropagation][2] - (
indices[3] - detId[iPropagation][3]) *
Geo::ZPAD + posCorr[2];
1091 float resY = deltaPos[iPropagation][1];
1092 float resXor = resX;
1093 float resZor = resZ;
1094 float res = TMath::Sqrt(resX * resX + resZ * resZ);
1098 }
else if (resX > 1.25) {
1101 resX = 1E-3 / (pt + 1E-3);
1106 }
else if (resZ > 1.75) {
1109 resZ = 1E-3 / (pt + 1E-3);
1112 LOG(
debug) <<
"resX = " << resX <<
", resZ = " << resZ <<
", res = " <<
res;
1113 if (
indices[0] != detId[iPropagation][0]) {
1116 if (
indices[1] != detId[iPropagation][1]) {
1119 if (
indices[2] != detId[iPropagation][2]) {
1122 float chi2 = 0.5 * (resX * resX * errXinv2 + resZ * resZ * errZinv2);
1124 if (
res < mSpaceTolerance && chi2 < mMatchParams->maxChi2) {
1125 LOG(
debug) <<
"MATCHING FOUND: We have a match! between track " << mTracksSectIndexCache[
type][sec][itrk] <<
" and TOF cluster " << mTOFClusSectIndexCache[
indices[0]][itof];
1126 foundCluster =
true;
1128 int eventIndexTOFCluster = mTOFClusSectIndexCache[
indices[0]][itof];
1129 mMatchedTracksPairsSec[sec].emplace_back(cacheTrk[itrk], eventIndexTOFCluster, mTOFClusWork[cacheTOF[itof]].
getTime(), chi2, trkLTInt[iPropagation], mTrackGid[sec][
type][cacheTrk[itrk]],
type, (trefTOF.getTime() - (minTrkTime + maxTrkTime - 100E3) * 0.5) * 1E-6, trefTOF.getZ(), resXor, resZor, resY);
1130 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setPt(pt);
1131 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResX(sqrt(1. / errXinv2));
1132 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResZ(sqrt(1. / errZinv2));
1133 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResT(resT);
1134 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setVz(0.0);
1135 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setChannel(mainChannel);
1143void MatchTOF::doMatchingForTPC(
int sec)
1147 int bc_grouping = 40;
1148 int bc_grouping_tolerance = bc_grouping + mTimeTolerance / 25;
1149 int bc_grouping_half = (bc_grouping + 1) / 2;
1153 auto& cacheTOF = mTOFClusSectIndexCache[sec];
1154 auto& cacheTrk = mTracksSectIndexCache[trkType::UNCONS][sec];
1155 int nTracks = cacheTrk.size(), nTOFCls = cacheTOF.size();
1156 LOG(
debug) <<
"Matching sector " << sec <<
": number of tracks: " << nTracks <<
", number of TOF clusters: " << nTOFCls;
1157 if (!nTracks || !nTOFCls) {
1161 float deltaPosTemp[3];
1162 std::array<float, 3>
pos;
1163 std::array<float, 3> posBeforeProp;
1168 std::vector<unsigned long> BCcand;
1170 std::vector<int> nStripsCrossedInPropagation;
1171 std::vector<std::array<std::array<int, 5>, 2>> detId;
1172 std::vector<std::array<o2::track::TrackLTIntegral, 2>> trkLTInt;
1173 std::vector<std::array<std::array<float, 3>, 2>> deltaPos;
1174 std::vector<std::array<int, 2>> nStepsInsideSameStrip;
1175 std::vector<std::array<float, 2>> Zshift;
1177 LOG(
debug) <<
"Trying to match %d tracks" << cacheTrk.size();
1179 for (
int itrk = 0; itrk < cacheTrk.size(); itrk++) {
1180 auto& trackWork = mTracksWork[sec][trkType::UNCONS][cacheTrk[itrk]];
1181 auto& trefTrk = trackWork.first;
1182 float pt = trefTrk.getPt();
1183 auto& intLT = mLTinfos[sec][trkType::UNCONS][cacheTrk[itrk]];
1185 float timeShift = intLT.getL() * 33.35641;
1189 nStripsCrossedInPropagation.clear();
1191 int side = mSideTPC[sec][cacheTrk[itrk]];
1193 double tpctime = trackWork.second.getTimeStamp();
1194 double minTrkTime = (tpctime - trackWork.second.getTimeStampError()) * 1.E6 + timeShift;
1195 minTrkTime =
int(minTrkTime / BCgranularity) * BCgranularity;
1196 double maxTrkTime = (tpctime + mExtraTPCFwdTime[sec][cacheTrk[itrk]]) * 1.E6 + timeShift;
1197 const float sqrt12inv = 1. / sqrt(12.);
1198 float resT = (maxTrkTime - minTrkTime) * sqrt12inv;
1201 for (
double tBC = minTrkTime; tBC < maxTrkTime; tBC += BCgranularity) {
1203 BCcand.emplace_back(ibc);
1204 nStripsCrossedInPropagation.emplace_back(0);
1208 int itofMax = nTOFCls;
1210 for (
auto itof = itof0; itof < nTOFCls; itof++) {
1211 auto& trefTOF = mTOFClusWork[cacheTOF[itof]];
1213 if (trefTOF.getTime() < minTrkTime) {
1218 if (trefTOF.getTime() > maxTrkTime) {
1223 if ((trefTOF.getZ() *
side < 0) && ((
side > 0) != (trackWork.first.getTgl() > 0))) {
1229 bc = (
bc / bc_grouping_half) * bc_grouping_half;
1231 bool isalreadyin =
false;
1233 for (
int k = 0; k < BCcand.size(); k++) {
1234 if (
bc == BCcand[k]) {
1240 BCcand.emplace_back(
bc);
1241 nStripsCrossedInPropagation.emplace_back(0);
1246 detId.reserve(BCcand.size());
1248 trkLTInt.reserve(BCcand.size());
1250 Zshift.reserve(BCcand.size());
1252 deltaPos.reserve(BCcand.size());
1253 nStepsInsideSameStrip.clear();
1254 nStepsInsideSameStrip.reserve(BCcand.size());
1269 int detIdTemp[5] = {-1, -1, -1, -1, -1};
1271 double reachedPoint = mXRef + istep *
step;
1274 for (
int ibc = 0; ibc < BCcand.size(); ibc++) {
1275 for (
int ii = 0; ii < 2; ii++) {
1276 nStepsInsideSameStrip[ibc][ii] = 0;
1277 for (
int iii = 0; iii < 5; iii++) {
1278 detId[ibc][ii][iii] = -1;
1282 while (propagateToRefX(trefTrk, reachedPoint, step, intLT) && reachedPoint <
Geo::RMAX) {
1285 trefTrk.getXYZGlo(
pos);
1286 for (
int ii = 0; ii < 3; ii++) {
1287 posFloat[ii] =
pos[ii];
1297 reachedPoint +=
step;
1300 for (
int ibc = 0; ibc < BCcand.size(); ibc++) {
1301 for (
int idet = 0; idet < 5; idet++) {
1302 detIdTemp[idet] = -1;
1306 posFloat[2] =
pos[2] - mTPCVDrift * (trackWork.second.getTimeStamp() - BCcand[ibc] *
Geo::BC_TIME_INPS * 1E-6);
1307 }
else if (
side < 0) {
1308 posFloat[2] =
pos[2] + mTPCVDrift * (trackWork.second.getTimeStamp() - BCcand[ibc] *
Geo::BC_TIME_INPS * 1E-6);
1310 posFloat[2] =
pos[2];
1313 float ZshiftCurrent = posFloat[2] -
pos[2];
1317 if (detIdTemp[2] == -1) {
1321 if (nStripsCrossedInPropagation[ibc] == 0 ||
1322 (nStripsCrossedInPropagation[ibc] >= 1 && (detId[ibc][nStripsCrossedInPropagation[ibc] - 1][0] != detIdTemp[0] || detId[ibc][nStripsCrossedInPropagation[ibc] - 1][1] != detIdTemp[1] || detId[ibc][nStripsCrossedInPropagation[ibc] - 1][2] != detIdTemp[2]))) {
1323 if (nStripsCrossedInPropagation[ibc] == 0) {
1324 LOG(
debug) <<
"We cross a strip for the first time";
1326 if (nStripsCrossedInPropagation[ibc] == 2) {
1329 nStripsCrossedInPropagation[ibc]++;
1333 if (nStepsInsideSameStrip[ibc][nStripsCrossedInPropagation[ibc] - 1] == 0) {
1334 trkLTInt[ibc][nStripsCrossedInPropagation[ibc] - 1] = intLT;
1336 int detIdTemp2[5] = {0, 0, 0, 0, 0};
1337 float deltaPosTemp2[3] = {deltaPosTemp[0], deltaPosTemp[1], deltaPosTemp[2]};
1339 const int maxnstep = 50;
1340 float xStart = trefTrk.getX();
1341 float xStop = xStart;
1342 trefTrk.getXYZGlo(
pos);
1343 for (
int ii = 0; ii < 3; ii++) {
1344 posFloat[ii] =
pos[ii];
1347 while (deltaPosTemp2[1] < -0.05 && detIdTemp2[2] != -1 && nstep < maxnstep) {
1350 propagateToRefXWithoutCov(trefTrk, xStop, 0.1, mBz, posFloat);
1352 posFloat[2] += ZshiftCurrent;
1355 if (detIdTemp2[2] != -1) {
1356 float dx = deltaPosTemp2[0] - deltaPosTemp[0];
1357 float dy = deltaPosTemp2[1] - deltaPosTemp[1];
1358 float dz = deltaPosTemp2[2] - deltaPosTemp[2];
1359 updateTL(trkLTInt[ibc][nStripsCrossedInPropagation[ibc] - 1], sqrt(dx * dx + dy * dy + dz * dz));
1360 detIdTemp[0] = detIdTemp2[0];
1361 detIdTemp[1] = detIdTemp2[1];
1362 detIdTemp[2] = detIdTemp2[2];
1363 detIdTemp[3] = detIdTemp2[3];
1364 detIdTemp[4] = detIdTemp2[4];
1365 deltaPosTemp[0] = deltaPosTemp2[0];
1366 deltaPosTemp[1] = deltaPosTemp2[1];
1367 deltaPosTemp[2] = deltaPosTemp2[2];
1372 updateTL(trkLTInt[ibc][nStripsCrossedInPropagation[ibc] - 1], -deltaPosTemp[1]);
1374 detId[ibc][nStripsCrossedInPropagation[ibc] - 1][0] = detIdTemp[0];
1375 detId[ibc][nStripsCrossedInPropagation[ibc] - 1][1] = detIdTemp[1];
1376 detId[ibc][nStripsCrossedInPropagation[ibc] - 1][2] = detIdTemp[2];
1377 detId[ibc][nStripsCrossedInPropagation[ibc] - 1][3] = detIdTemp[3];
1378 detId[ibc][nStripsCrossedInPropagation[ibc] - 1][4] = detIdTemp[4];
1379 deltaPos[ibc][nStripsCrossedInPropagation[ibc] - 1][0] = deltaPosTemp[0];
1380 deltaPos[ibc][nStripsCrossedInPropagation[ibc] - 1][1] = deltaPosTemp[1];
1381 deltaPos[ibc][nStripsCrossedInPropagation[ibc] - 1][2] = deltaPosTemp[2];
1383 Zshift[ibc][nStripsCrossedInPropagation[ibc] - 1] = ZshiftCurrent;
1385 nStepsInsideSameStrip[ibc][nStripsCrossedInPropagation[ibc] - 1]++;
1398 for (
int ibc = 0; ibc < BCcand.size(); ibc++) {
1401 for (Int_t imatch = 0; imatch < nStripsCrossedInPropagation[ibc]; imatch++) {
1403 deltaPos[ibc][imatch][0] /= nStepsInsideSameStrip[ibc][imatch];
1404 deltaPos[ibc][imatch][1] /= nStepsInsideSameStrip[ibc][imatch];
1405 deltaPos[ibc][imatch][2] /= nStepsInsideSameStrip[ibc][imatch];
1408 if (nStripsCrossedInPropagation[ibc] == 0) {
1412 bool foundCluster =
false;
1413 for (
auto itof = itof0; itof < itofMax; itof++) {
1415 auto& trefTOF = mTOFClusWork[cacheTOF[itof]];
1418 if (trefTOF.getTime() < minTime) {
1421 if (trefTOF.getTime() > maxTime) {
1425 int mainChannel = trefTOF.getMainContributingChannel();
1429 bool isInStrip =
false;
1430 for (
auto iPropagation = 0; iPropagation < nStripsCrossedInPropagation[ibc]; iPropagation++) {
1431 if (detId[ibc][iPropagation][1] ==
indices[1] && detId[ibc][iPropagation][2] ==
indices[2]) {
1445 float posCorr[3] = {0, 0, 0};
1472 float ndifInv = 1. / ndigits;
1474 posCorr[0] *= ndifInv;
1475 posCorr[1] *= ndifInv;
1476 posCorr[2] *= ndifInv;
1482 for (
auto iPropagation = 0; iPropagation < nStripsCrossedInPropagation[ibc]; iPropagation++) {
1483 if (detId[ibc][iPropagation][1] !=
indices[1] || detId[ibc][iPropagation][2] !=
indices[2]) {
1492 if (tof - trkLTInt[ibc][iPropagation].getTOF(6) > 2000) {
1497 if (std::abs(tof - trkLTInt[ibc][iPropagation].getTOF(2)) < 2000) {
1498 }
else if (std::abs(tof - trkLTInt[ibc][iPropagation].getTOF(3)) < 2000) {
1499 }
else if (std::abs(tof - trkLTInt[ibc][iPropagation].getTOF(4)) < 2000) {
1506 const float errXinvMin = 1. / (mMatchParams->
maxResX * mMatchParams->
maxResX);
1507 const float errZinvMin = 1. / (mMatchParams->
maxResZ * mMatchParams->
maxResZ);
1508 float errXinv2 = 1. / (trefTrk.getSigmaY2());
1509 float errZinv2 = 1. / (trefTrk.getSigmaZ2() * cosangle);
1511 if (errXinv2 < errXinvMin) {
1512 errXinv2 = errXinvMin;
1514 if (errZinv2 < errZinvMin) {
1515 errZinv2 = errZinvMin;
1519 LOG(
debug) <<
"Propagated Track [" << itrk <<
"]: detId[" << iPropagation <<
"] = " << detId[ibc][iPropagation][0] <<
", " << detId[ibc][iPropagation][1] <<
", " << detId[ibc][iPropagation][2] <<
", " << detId[ibc][iPropagation][3] <<
", " << detId[ibc][iPropagation][4];
1520 float resX = deltaPos[ibc][iPropagation][0] - (
indices[4] - detId[ibc][iPropagation][4]) *
Geo::XPAD + posCorr[0];
1521 float resZ = deltaPos[ibc][iPropagation][2] - (
indices[3] - detId[ibc][iPropagation][3]) *
Geo::ZPAD + posCorr[2];
1522 float resY = deltaPos[ibc][iPropagation][1];
1523 if (BCcand[ibc] > bcClus) {
1524 resZ += (BCcand[ibc] - bcClus) * vdriftInBC *
side;
1526 resZ -= (bcClus - BCcand[ibc]) * vdriftInBC *
side;
1528 float resXor = resX;
1529 float resZor = resZ;
1530 float res = TMath::Sqrt(resX * resX + resZ * resZ);
1534 }
else if (resX > 1.25) {
1537 resX = 1E-3 / (pt + 1E-3);
1542 }
else if (resZ > 1.75) {
1545 resZ = 1E-3 / (pt + 1E-3);
1548 if (
indices[0] != detId[ibc][iPropagation][0]) {
1551 if (
indices[1] != detId[ibc][iPropagation][1]) {
1554 if (
indices[2] != detId[ibc][iPropagation][2]) {
1558 LOG(
debug) <<
"resX = " << resX <<
", resZ = " << resZ <<
", res = " <<
res;
1559 float chi2 = mIsCosmics ? resX : 0.5 * (resX * resX * errXinv2 + resZ * resZ * errZinv2);
1561 if (
res < mSpaceTolerance && chi2 < mMatchParams->maxChi2) {
1562 LOG(
debug) <<
"MATCHING FOUND: We have a match! between track " << mTracksSectIndexCache[trkType::UNCONS][sec][itrk] <<
" and TOF cluster " << mTOFClusSectIndexCache[
indices[0]][itof];
1563 foundCluster =
true;
1566 int eventIndexTOFCluster = mTOFClusSectIndexCache[
indices[0]][itof];
1567 mMatchedTracksPairsSec[sec].emplace_back(cacheTrk[itrk], eventIndexTOFCluster, mTOFClusWork[cacheTOF[itof]].
getTime(), chi2, trkLTInt[ibc][iPropagation], mTrackGid[sec][trkType::UNCONS][cacheTrk[itrk]], trkType::UNCONS, trefTOF.getTime() * 1E-6 - tpctime, trefTOF.getZ(), resXor, resZor, resY);
1568 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setPt(pt);
1569 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResX(sqrt(1. / errXinv2));
1570 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResZ(sqrt(1. / errZinv2));
1571 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setResT(resT);
1572 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setVz(mVZtpcOnly[sec][itrk] + Zshift[ibc][iPropagation]);
1573 mMatchedTracksPairsSec[sec][mMatchedTracksPairsSec[sec].size() - 1].setChannel(mainChannel);
1587 mHasFillScheme =
true;
1593 if (FITRecPoints.size() == 0) {
1599 bool bestQuality =
false;
1600 const int distThr = 8;
1602 for (
unsigned int i = 0;
i < FITRecPoints.size();
i++) {
1603 const auto& ft = FITRecPoints[
i];
1604 if (!FT0Params.isSelected(ft)) {
1608 if (mHasFillScheme && !mFillScheme[
ir.
bc]) {
1611 bool quality = (std::abs(FITRecPoints[
i].getCollisionTime(0)) < 1000 && std::abs(FITRecPoints[
i].getVertex()) < 1000);
1612 if (bestQuality && !quality) {
1616 int dist =
bc - bct0;
1618 bool worseDistance = dist < 0 || dist > distThr || dist > distMax;
1619 if (worseDistance) {
1623 bestQuality = quality;
1631void MatchTOF::selectBestMatches(
int sec)
1633 if (mSetHighPurity) {
1634 BestMatchesHP(mMatchedTracksPairsSec[sec], mMatchedTracks, mMatchedTracksIndex[sec], mMatchedClustersIndex, mFITRecPoints, mTOFClusWork, mCalibInfoTOF, mTimestamp, mMCTruthON, mTOFClusLabels, mTracksLblWork[sec], mOutTOFLabels);
1637 BestMatches(mMatchedTracksPairsSec[sec], mMatchedTracks, mMatchedTracksIndex[sec], mMatchedClustersIndex, mFITRecPoints, mTOFClusWork, mTracksWork[sec], mCalibInfoTOF, mTimestamp, mMCTruthON, mTOFClusLabels, mTracksLblWork[sec], mOutTOFLabels, mMatchParams->
calibMaxChi2);
1640void MatchTOF::BestMatches(std::vector<o2::dataformats::MatchInfoTOFReco>& matchedTracksPairs, std::vector<o2::dataformats::MatchInfoTOF>* matchedTracks, std::vector<int>* matchedTracksIndex,
int* matchedClustersIndex,
const gsl::span<const o2::ft0::RecPoints>& FITRecPoints,
const std::vector<Cluster>& TOFClusWork,
const std::vector<matchTrack>* TracksWork, std::vector<o2::dataformats::CalibInfoTOF>& CalibInfoTOF,
unsigned long Timestamp,
bool MCTruthON,
const o2::dataformats::MCTruthContainer<o2::MCCompLabel>* TOFClusLabels,
const std::vector<o2::MCCompLabel>* TracksLblWork, std::vector<o2::MCCompLabel>* OutTOFLabels,
float calibMaxChi2)
1650 int trkType = (
int)matchingPair.getTrackType();
1652 int itrk = matchingPair.getIdLocal();
1654 int trkTypeSplitted =
trkType;
1655 auto sourceID = matchingPair.getTrackRef().getSource();
1657 trkTypeSplitted = (
int)trkType::TPCTRD;
1659 trkTypeSplitted = (
int)trkType::ITSTPCTRD;
1662 if (matchedClustersIndex[matchingPair.getTOFClIndex()] != -1) {
1665 if (matchedTracksIndex[
trkType][itrk] != -1) {
1668 int pairInd = matchedTracksIndex[
trkType][itrk];
1669 auto& prevMatching = matchedTracks[trkTypeSplitted][pairInd];
1671 if (pairInd >= matchedTracks[trkTypeSplitted].
size()) {
1672 LOG(error) <<
"Pair index out of range when trying to update TOF time. This should not happen";
1676 int istripOld = TOFClusWork[prevMatching.getTOFClIndex()].getMainContributingChannel() /
o2::tof::Geo::NPADS;
1677 int istripNew = TOFClusWork[matchingPair.getTOFClIndex()].getMainContributingChannel() /
o2::tof::Geo::NPADS;
1679 if (istripOld != istripNew) {
1683 const float cinv = 1.E+10 / TMath::C();
1684 float deltaT = (intLTnew.getL() - intLTold.getL()) * cinv;
1686 float timeNew = TOFClusWork[matchingPair.getTOFClIndex()].getTime() - deltaT;
1687 float timeOld = TOFClusWork[prevMatching.getTOFClIndex()].getTime();
1689 if (std::abs(timeNew - timeOld) < 200) {
1691 prevMatching.setSignal((timeNew + timeOld) * 0.5);
1692 float geanttime = (TOFClusWork[matchingPair.getTOFClIndex()].getTgeant() + TOFClusWork[prevMatching.getTOFClIndex()].getTgeant() - deltaT * 1E-3) * 0.5;
1693 double t0 = (TOFClusWork[matchingPair.getTOFClIndex()].getT0true() + TOFClusWork[prevMatching.getTOFClIndex()].getT0true()) * 0.5;
1694 prevMatching.setTgeant(geanttime);
1695 prevMatching.setT0true(
t0);
1696 prevMatching.setChi2(0);
1697 matchedClustersIndex[matchingPair.getTOFClIndex()] = matchedTracksIndex[
trkType][itrk];
1704 matchedTracksIndex[
trkType][itrk] = matchedTracks[trkTypeSplitted].size();
1705 matchedClustersIndex[matchingPair.getTOFClIndex()] = matchedTracksIndex[
trkType][itrk];
1707 matchingPair.setTgeant(TOFClusWork[matchingPair.getTOFClIndex()].getTgeant());
1708 matchingPair.setT0true(TOFClusWork[matchingPair.getTOFClIndex()].getT0true());
1711 const auto& tofcl = TOFClusWork[matchingPair.getTOFClIndex()];
1712 if (tofcl.getNumOfContributingChannels() > 1) {
1714 matchingPair.setHitPatternUpDown(tofcl.isAdditionalChannelSet(o2::tof::Cluster::kUp) ||
1715 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kUpLeft) ||
1716 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kUpRight) ||
1717 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kDown) ||
1718 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kDownLeft) ||
1719 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kDownRight));
1721 matchingPair.setHitPatternLeftRight(tofcl.isAdditionalChannelSet(o2::tof::Cluster::kLeft) ||
1722 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kDownLeft) ||
1723 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kUpLeft) ||
1724 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kRight) ||
1725 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kDownRight) ||
1726 tofcl.isAdditionalChannelSet(o2::tof::Cluster::kUpRight));
1732 float t0BestRes = 200;
1733 if (FITRecPoints.size() > 0) {
1735 if (
index > -1 && FITRecPoints[
index].isValidTime(1) && FITRecPoints[
index].isValidTime(2)) {
1736 t0Best += FITRecPoints[
index].getCollisionTime(0);
1740 matchingPair.setFT0Best(t0Best, t0BestRes);
1741 matchedTracks[trkTypeSplitted].push_back(matchingPair);
1748 if (FITRecPoints.size() > 0 && mIsFIT) {
1756 }
else if (!mIsFIT) {
1777 if (t0info > 0 && mIsFIT) {
1780 }
else if (!mIsFIT) {
1785 float pt = trc.getPt();
1799 if (TOFClusWork[matchingPair.getTOFClIndex()].getNumOfContributingChannels() != 1) {
1803 if (matchingPair.getChi2() < calibMaxChi2 && t0info > 0) {
1804 CalibInfoTOF.emplace_back(TOFClusWork[matchingPair.getTOFClIndex()].getMainContributingChannel(),
1805 Timestamp / 1000 +
int(TOFClusWork[matchingPair.getTOFClIndex()].getTimeRaw() * 1E-12),
1807 TOFClusWork[matchingPair.getTOFClIndex()].getTot(),
mask,
flags);
1811 const auto& labelsTOF = TOFClusLabels->
getLabels(matchingPair.getTOFClIndex());
1812 auto& labelTrack = TracksLblWork[
trkType][itrk];
1815 for (
auto& lbl : labelsTOF) {
1816 if (labelTrack == lbl) {
1820 OutTOFLabels[trkTypeSplitted].emplace_back(labelTrack).setFakeFlag(fake);
1826void MatchTOF::BestMatchesHP(std::vector<o2::dataformats::MatchInfoTOFReco>& matchedTracksPairs, std::vector<o2::dataformats::MatchInfoTOF>* matchedTracks, std::vector<int>* matchedTracksIndex,
int* matchedClustersIndex,
const gsl::span<const o2::ft0::RecPoints>& FITRecPoints,
const std::vector<Cluster>& TOFClusWork, std::vector<o2::dataformats::CalibInfoTOF>& CalibInfoTOF,
unsigned long Timestamp,
bool MCTruthON,
const o2::dataformats::MCTruthContainer<o2::MCCompLabel>* TOFClusLabels,
const std::vector<o2::MCCompLabel>* TracksLblWork, std::vector<o2::MCCompLabel>* OutTOFLabels)
1829 float chi2SeparationCut = 2;
1832 std::vector<o2::dataformats::MatchInfoTOFReco> tmpMatch;
1839 int trkType = (
int)matchingPair.getTrackType();
1841 int itrk = matchingPair.getIdLocal();
1843 bool discard = matchingPair.getChi2() > chi2S;
1845 if (matchedTracksIndex[
trkType][itrk] != -1) {
1846 auto winnerChi = tmpMatch[matchedTracksIndex[
trkType][itrk]].getChi2();
1847 if (winnerChi < 0) {
1850 if (matchingPair.getChi2() - winnerChi < chi2SeparationCut) {
1851 tmpMatch[matchedTracksIndex[
trkType][itrk]].setChi2(-1);
1856 if (matchedClustersIndex[matchingPair.getTOFClIndex()] != -1) {
1857 auto winnerChi = tmpMatch[matchedClustersIndex[matchingPair.getTOFClIndex()]].getChi2();
1858 if (winnerChi < 0) {
1861 if (matchingPair.getChi2() - winnerChi < chi2SeparationCut) {
1862 tmpMatch[matchedClustersIndex[matchingPair.getTOFClIndex()]].setChi2(-1);
1868 matchedTracksIndex[
trkType][itrk] = tmpMatch.size();
1869 matchedClustersIndex[matchingPair.getTOFClIndex()] = matchedTracksIndex[
trkType][itrk];
1870 tmpMatch.push_back(matchingPair);
1875 for (
auto& matchingPair : tmpMatch) {
1876 if (matchingPair.getChi2() <= 0) {
1879 int trkType = (
int)matchingPair.getTrackType();
1880 int itrk = matchingPair.getIdLocal();
1882 int trkTypeSplitted =
trkType;
1883 auto sourceID = matchingPair.getTrackRef().getSource();
1885 trkTypeSplitted = (
int)trkType::TPCTRD;
1887 trkTypeSplitted = (
int)trkType::ITSTPCTRD;
1889 matchedTracks[trkTypeSplitted].push_back(matchingPair);
1896 if (FITRecPoints.size() > 0 && mIsFIT) {
1912 CalibInfoTOF.emplace_back(TOFClusWork[matchingPair.getTOFClIndex()].getMainContributingChannel(),
1913 Timestamp / 1000 +
int(TOFClusWork[matchingPair.getTOFClIndex()].getTimeRaw() * 1E-12),
1914 TOFClusWork[matchingPair.getTOFClIndex()].getTimeRaw() - t0info - intLT.getTOF(
o2::track::PID::Pion),
1915 TOFClusWork[matchingPair.getTOFClIndex()].getTot(), 0);
1919 const auto& labelsTOF = TOFClusLabels->
getLabels(matchingPair.getTOFClIndex());
1920 auto& labelTrack = TracksLblWork[
trkType][itrk];
1923 for (
auto& lbl : labelsTOF) {
1924 if (labelTrack == lbl) {
1928 OutTOFLabels[trkTypeSplitted].emplace_back(labelTrack).setFakeFlag(fake);
1938 bool refReached =
false;
1939 float xStart = trc.getX();
1946 while (hasPropagated) {
1947 if (trc.getX() > xRef) {
1951 if (std::abs(trc.getY()) > trc.getX() * tanHalfSector) {
1955 if (!trc.rotate(alphaNew) != 0) {
1968 return refReached && std::abs(trc.getSnp()) < 0.95;
1971bool MatchTOF::propagateToRefXWithoutCov(
const o2::track::TrackParCov& trc,
float xRef,
float stepInCm,
float bzField)
1977 bool refReached =
false;
1978 float xStart = trcNoCov.getX();
1984 bool hasPropagated = trcNoCov.propagateParamTo(xStart + istep * stepInCm, bzField);
1985 while (hasPropagated) {
1986 if (trcNoCov.getX() > xRef) {
1990 if (std::abs(trcNoCov.getY()) > trcNoCov.getX() * tanHalfSector) {
1994 if (!trcNoCov.rotateParam(alphaNew) != 0) {
2002 hasPropagated = trcNoCov.propagateParamTo(xStart + istep * stepInCm, bzField);
2007 return refReached && std::abs(trcNoCov.getSnp()) < 0.95 && std::abs(trcNoCov.getZ()) <
Geo::MAXHZTOF;
2012 for (
int i = 0;
i < intLT.getNTOFs();
i++) {
2013 float betainv = intLT.getTOF(
i) / intLT.getL();
2014 intLT.setTOF(intLT.getTOF(
i) + deltal * betainv,
i);
2016 intLT.setL(intLT.getL() + deltal);
2020bool MatchTOF::propagateToRefXWithoutCov(
const o2::track::TrackParCov& trc,
float xRef,
float stepInCm,
float bzField,
float pos[3])
2026 bool refReached =
false;
2027 float xStart = trcNoCov.getX();
2033 bool hasPropagated = trcNoCov.propagateParamTo(xStart + istep * stepInCm, bzField);
2034 while (hasPropagated) {
2035 if (trcNoCov.getX() > xRef) {
2039 if (fabs(trcNoCov.getY()) > trcNoCov.getX() * tanHalfSector) {
2043 if (!trcNoCov.rotateParam(alphaNew) != 0) {
2051 hasPropagated = trcNoCov.propagateParamTo(xStart + istep * stepInCm, bzField);
2057 trcNoCov.getXYZGlo(xyz);
2062 return refReached && std::abs(trcNoCov.getSnp()) < 0.95 && TMath::Abs(trcNoCov.getZ()) <
Geo::MAXHZTOF;
2076void MatchTOF::updateTimeDependentParams()
2080 mTPCTBinMUS = elParam.ZbinWidth;
2081 mTPCTBinMUSInv = 1. / mTPCTBinMUS;
2082 mTPCBin2Z = mTPCTBinMUS * mTPCVDrift;
2085 mMaxInvPt = std::abs(mBz) > 0.1 ? 1. / (std::abs(mBz) * 0.05) : 999.;
2088 float scale = mTPCCorrMapsHelper->getInstLumiCTP();
2092 mCovDiagInner = trackTune.getCovInnerTotal(scale);
2093 mCovDiagOuter = trackTune.getCovOuterTotal(scale);
2099 auto&
match = mMatchedTracks[trkType::TPC][matchedID];
2101 const auto& tofCl = mTOFClustersArrayInp[
match.getTOFClIndex()];
2102 const auto& intLT =
match.getLTIntegralOut();
2105 auto timeTOFMUS = (tofCl.getTime() - intLT.getTOF(tpcTrOrig.getPID())) * 1e-6;
2106 auto timeTOFTB = timeTOFMUS * mTPCTBinMUSInv;
2107 auto deltaTBins = timeTOFTB - tpcTrOrig.getTime0();
2108 float timeErr = 0.010;
2109 auto dzCorr = deltaTBins * mTPCBin2Z;
2111 if (mTPCClusterIdxStruct) {
2112 trConstr.o2::track::TrackParCov::operator=(tpcTrOrig.getOuterParam());
2115 trConstr.o2::track::TrackParCov::operator=(tpcTrOrig);
2121 auto zTrack = trConstr.getZ();
2122 auto zTrackOut = trConstrOut.getZ();
2124 if (tpcTrOrig.hasASideClustersOnly()) {
2126 zTrackOut += dzCorr;
2127 }
else if (tpcTrOrig.hasCSideClustersOnly()) {
2129 zTrackOut -= dzCorr;
2133 trConstr.setZ(zTrack);
2134 trConstrOut.setZ(zTrackOut);
2138 if (mTPCClusterIdxStruct) {
2141 if (mTPCRefitter->RefitTrackAsTrackParCov(trConstr, tpcTrOrig.getClusterRef(), timeTOFTB, &chi2,
false,
true) < 0) {
2142 LOGP(
debug,
"Inward Refit failed {}", trConstr.asString());
2146 if (!trackTune.useTPCInnerCorr) {
2147 trConstr.updateParams(trackTune.tpcParInner);
2156 if (mTPCRefitter->RefitTrackAsTrackParCov(trConstrOut, tpcTrOrig.getClusterRef(), timeTOFTB, &chi2,
true,
true) < 0) {
2157 LOGP(
debug,
"Outward Refit failed {}", trConstrOut.asString());
2168 if (mTPCClusterIdxStruct) {
2169 mTPCRefitter = std::make_unique<o2::gpu::GPUO2InterfaceRefit>(mTPCClusterIdxStruct, mTPCCorrMapsHelper, mBz,
2170 mTPCTrackClusIdx.data(), 0, mTPCRefitterShMap.data(),
General auxilliary methods.
Definition of the GeometryManager class.
Header of the General Run Parameters object.
Some ALICE geometry constants of common interest.
Definition of a container to keep Monte Carlo truth external to simulation objects.
Definition of the fast magnetic field parametrization MagFieldFast.
Definition of the MagF class.
Class to perform TOF matching to global tracks.
Definition of the parameter class for the detector electronics.
Definition of the parameter class for the detector gas.
Header to collect physics constants.
Wrapper container for different reconstructed object types.
Track Length and TOF integral.
Configurable params for tracks ad hoc tuning.
calibration data from laser track calibration
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
static const MatchTOFParams & Instance()
void setTPCVDrift(const o2::tpc::VDriftCorrFact &v)
void setTPCCorrMaps(o2::gpu::CorrectionMapsHelper *maph)
void printCandidatesTOF() const
set time tolerance on track-TOF times comparison
void run(const o2::globaltracking::RecoContainer &inp, unsigned long firstTForbit=0)
< perform matching for provided input
void setDebugFlag(UInt_t flag, bool on=true)
set the name of output debug file
static void groupingMatch(const std::vector< o2::dataformats::MatchInfoTOFReco > &origin, std::vector< std::vector< o2::dataformats::MatchInfoTOFReco > > &grouped, std::vector< std::vector< int > > &firstEls, std::vector< std::vector< int > > &secondEls)
static void printGrouping(const std::vector< o2::dataformats::MatchInfoTOFReco > &origin, const std::vector< std::vector< o2::dataformats::MatchInfoTOFReco > > &grouped)
static int findFITIndex(int bc, const gsl::span< const o2::ft0::RecPoints > &FITRecPoints, unsigned long firstOrbit)
bool makeConstrainedTPCTrack(int matchedID, o2::dataformats::TrackTPCTOF &trConstr)
populate externally provided container by TOF-time-constrained TPC tracks
HMPID cluster implementation.
static constexpr Float_t RMAX
static constexpr Float_t ZPAD
static constexpr Float_t XPAD
static constexpr Double_t BC_TIME_INPS
static constexpr Double_t BC_TIME_INPS_INV
static Float_t getAngles(Int_t iplate, Int_t istrip)
static constexpr Int_t NPADS
static constexpr Int_t LATENCYWINDOW_IN_BC
static constexpr Float_t MAXHZTOF
static void getVolumeIndices(Int_t index, Int_t *detId)
static void getPadDxDyDz(const Float_t *pos, Int_t *det, Float_t *DeltaPos, int sector=-1)
static double subtractInteractionBC(double time, int &mask, bool subLatency=false)
static int getNinteractionBC()
static bool hasFillScheme()
static int getInteractionBC(int ibc)
std::array< value_t, 3 > dim3_t
bool match(const std::vector< std::string > &queries, const char *pattern)
GLboolean GLboolean GLboolean b
GLint GLint GLsizei GLint GLenum GLenum type
GLsizei GLenum const void * indices
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t0
GLboolean GLboolean GLboolean GLboolean a
constexpr float XTPCOuterRef
reference radius to propagate outer TPC track
constexpr float XTPCInnerRef
reference radius at which TPC provides the tracks
constexpr int LHCMaxBunches
constexpr float SectorSpanRad
float angle2Alpha(float phi)
int angle2Sector(float phi)
Enum< T >::Iterator begin(Enum< T >)
std::string getTime(uint64_t ts)
uint16_t bc
bunch crossing ID of interaction
o2::InteractionRecord startIR
const o2::tpc::TrackTPC & getTPCTrack(GTrackID id) const
auto getTPCITSTrackMCLabel(GTrackID id) const
void createTracksVariadic(T creator, GTrackID::mask_t srcSel=GTrackID::getSourcesMask("all")) const
auto getTPCTracksClusterRefs() const
auto getTPCTrackMCLabel(GTrackID id) const
auto getTOFClustersMCLabels() const
auto getTOFClusters() const
gsl::span< const unsigned char > clusterShMapTPC
externally set TPC clusters sharing map
auto getFT0RecPoints() const
std::unique_ptr< o2::tpc::internal::getWorkflowTPCInput_ret > inputsTPCclusters
auto getTPCTracks() const
gsl::span< const unsigned int > occupancyMapTPC
externally set TPC clusters occupancy map
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
o2::InteractionRecord ir(0, 0)