14#include "GPUParam.inc"
23#include <fairlogger/Logger.h>
40#include <Math/SMatrix.h>
41#include <Math/SVector.h>
43#include <TGeoGlobalMagField.h>
67constexpr float MatchTPCITS::Tan70, MatchTPCITS::Cos70I2, MatchTPCITS::MaxSnp, MatchTPCITS::MaxTgp;
90 LOG(fatal) <<
"init() was not done yet";
95 updateTimeDependentParams();
97 mTimer[SWTot].Start(
false);
100 if (!prepareITSData() || !prepareTPCData() || !prepareFITData()) {
103 if (mVDriftCalibOn) {
104 calib.emplace_back(mTPCVDrift, mTPCDrift.
refVDrift, -999.);
105 calib.emplace_back(mTPCDriftTimeOffset, mTPCDrift.
refTimeOffset, -999.);
108 mTimer[SWDoMatching].Start(
false);
112 mTimer[SWDoMatching].Stop();
113 if constexpr (
false) {
114 mTimer[SWTot].Stop();
117 mTimer[SWTot].Start(
false);
122 bool fullMatchRefitDone =
false;
124 fullMatchRefitDone =
runAfterBurner(matchedTracks, matchLabels, ABTrackletLabels, ABTrackletClusterIDs, ABTrackletRefs, calib);
126 if (!fullMatchRefitDone) {
130 reportSizes(matchedTracks, ABTrackletRefs, ABTrackletClusterIDs, matchLabels, ABTrackletLabels, calib);
132#ifdef _ALLOW_DEBUG_TREES_
139 mTimer[SWTot].Stop();
150 for (
int i = 0;
i < NStopWatches;
i++) {
151 LOGF(info,
"Timing for %15s: Cpu: %.3e Real: %.3e s in %d slots of TF#%d", TimerName[
i], mTimer[
i].CpuTime(), mTimer[
i].RealTime(), mTimer[
i].Counter() - 1, mTFCount);
158#ifdef _ALLOW_DEBUG_TREES_
167 mMatchRecordsTPC.clear();
168 mMatchRecordsITS.clear();
169 mWinnerChi2Refit.clear();
172 mInteractions.clear();
173 mITSROFTimes.clear();
174 mITSTrackROFContMapping.clear();
175 mITSClustersArray.clear();
176 mITSClusterSizes.clear();
178 mTPCABIndexCache.clear();
179 mABWinnersIDs.clear();
180 mABClusterLinkIndex.clear();
181 mNMatchesControl = 0;
184 mITSSectIndexCache[sec].clear();
185 mITSTimeStart[sec].clear();
186 mTPCSectIndexCache[sec].clear();
187 mTPCTimeStart[sec].clear();
194 for (
int i = 0;
i < mNThreads;
i++) {
203 mTPCVDrift =
v.getVDrift();
204 mTPCDriftTimeOffset =
v.getTimeOffset();
210 mTPCCorrMapsHelper = maph;
218 LOG(error) <<
"Initialization was already done";
221 for (
int i = NStopWatches;
i--;) {
226 YMaxAtXMatchingRef = mParams->
XMatchingRef * 0.17632698;
231 if (!prop->getMatLUT() && mParams->
matCorr == o2::base::Propagator::MatCorrType::USEMatCorrLUT) {
232 LOG(warning) <<
"Requested material LUT is not loaded, switching to TGeo usage";
241#ifdef _ALLOW_DEBUG_TREES_
244 mDBGOut = std::make_unique<o2::utils::TreeStreamRedirector>(mDebugTreeFileName.data(),
"recreate");
256 if (fair::Logger::Logging(fair::Severity::info)) {
262void MatchTPCITS::updateTimeDependentParams()
267 mTPCTBinMUS = elParam.ZbinWidth;
268 mTPCTBinNS = mTPCTBinMUS * 1e3;
269 mTPCZMax = detParam.TPClength;
270 mTPCTBinMUSInv = 1. / mTPCTBinMUS;
271 assert(mITSROFrameLengthMUS > 0.0f);
272 mTPCBin2Z = mTPCTBinMUS * mTPCVDrift;
273 mZ2TPCBin = 1. / mTPCBin2Z;
274 mTPCVDriftInv = 1. / mTPCVDrift;
275 mNTPCBinsFullDrift = mTPCZMax * mZ2TPCBin;
279 mFieldON = std::abs(mBz) > 0.01;
286 mTPCmeanX0Inv = matbd.meanX2X0 / matbd.length;
289 float scale = mTPCCorrMapsHelper->getInstLumiCTP();
293 mCovDiagInner = trackTune.getCovInnerTotal(scale);
294 mCovDiagOuter = trackTune.getCovOuterTotal(scale);
298void MatchTPCITS::selectBestMatches()
301 mTimer[SWSelectBest].Start(
false);
302 int nValidated = 0, iter = 0;
307 int ntpc = mTPCWork.size(), nremaining = 0;
308 for (
int it = 0; it < ntpc; it++) {
309 auto& tTPC = mTPCWork[it];
310 if (isDisabledTPC(tTPC) || isValidatedTPC(tTPC)) {
314 if (validateTPCMatch(it)) {
316 if (mVDriftCalibOn && (!mFieldON || std::abs(tTPC.getQ2Pt()) < mParams->
maxVDriftTrackQ2Pt)) {
323 LOGP(info,
"iter {}: Validated {} of {} remaining matches", iter, nValidated, nremaining);
326 mNMatches += nValidated;
327 }
while (nValidated);
329 mTimer[SWSelectBest].Stop();
330 LOGP(info,
"Validated {} matches out of {} for {} TPC and {} ITS tracks in {} iterations", mNMatches, mNMatchesControl, mTPCWork.size(), mITSWork.size(), iter);
334bool MatchTPCITS::validateTPCMatch(
int iTPC)
336 auto& tTPC = mTPCWork[iTPC];
337 auto& rcTPC = mMatchRecordsTPC[tTPC.matchID];
339 auto& tITS = mITSWork[rcTPC.partnerID];
340 auto& rcITS = mMatchRecordsITS[tITS.matchID];
344 if (rcITS.partnerID == iTPC) {
345 int cloneID = tITS.getCloneShift();
347 cloneID += rcTPC.partnerID;
348 auto& tITSClone = mITSWork[cloneID];
349 if (isDisabledITS(tITSClone)) {
352 int nextITSCloneMatchID = tITSClone.matchID;
353 if (rcITS.isBetter(mMatchRecordsITS[nextITSCloneMatchID])) {
354 LOGP(
debug,
"Suppressing clone cloneID={} of winner clone {} of source ITS {}", cloneID, rcTPC.partnerID, tITS.sourceID);
355 while (nextITSCloneMatchID >
MinusOne) {
356 auto& rcITSClone = mMatchRecordsITS[nextITSCloneMatchID];
357 removeITSfromTPC(cloneID, rcITSClone.partnerID);
358 nextITSCloneMatchID = rcITSClone.nextRecID;
366 int nextTPC = rcTPC.nextRecID;
368 auto& rcTPCrem = mMatchRecordsTPC[nextTPC];
369 removeTPCfromITS(iTPC, rcTPCrem.partnerID);
370 nextTPC = rcTPCrem.nextRecID;
373 int itsWinID = rcTPC.partnerID;
376 int nextITS = rcITS.nextRecID;
378 auto& rcITSrem = mMatchRecordsITS[nextITS];
379 removeITSfromTPC(itsWinID, rcITSrem.partnerID);
380 nextITS = rcITSrem.nextRecID;
390int MatchTPCITS::getNMatchRecordsTPC(
const TrackLocTPC& tTPC)
const
395 recID = mMatchRecordsTPC[recID].nextRecID;
402int MatchTPCITS::getNMatchRecordsITS(
const TrackLocITS& tTPC)
const
407 auto& itsRecord = mMatchRecordsITS[recID];
408 recID = itsRecord.nextRecID;
418 const float SQRT12DInv = 2. / sqrt(12.);
422 const auto& tpcOrig = mTPCTracksArray[tpcID];
429 tpcOrig.getClusterReference(mTPCTrackClusIdx, tpcOrig.getNClusterReferences() - 1, clSect, clRow, clIdx);
433 const auto& clus = mTPCClusterIdxStruct->
clusters[clSect][clRow][clIdx];
435 if (padFromEdge > TPCGeometry.NPads(clRow) / 2) {
436 padFromEdge = TPCGeometry.NPads(clRow) - 1 - padFromEdge;
441 if (extConstrained) {
444 terr += tpcTimeBin2MUS(tpcOrig.hasBothSidesClusters() ? mParams->
safeMarginTPCITSTimeBin : mTPCTimeEdgeTSafeMargin);
446 auto& trc = mTPCWork.emplace_back(
447 TrackLocTPC{_tr, {
t0 - terr,
t0 + terr}, extConstrained ?
t0 : tpcTimeBin2MUS(tpcOrig.getTime0()),
449 terr * (extConstrained ? mTPCExtConstrainedNSigmaInv : SQRT12DInv),
460 if (srcGID.getSource() ==
GTrackID::TPC && !trackTune.sourceLevelTPC) {
461 if (trackTune.useTPCInnerCorr) {
462 trc.updateParams(trackTune.tpcParInner);
468 if (!propagateToRefX(trc)) {
473 mTPCLblWork.emplace_back(mTPCTrkLabels[tpcID]);
481bool MatchTPCITS::prepareTPCData()
484 mTimer[SWPrepTPC].Start(
false);
485 const auto& inp = *mRecoCont;
487 mTPCTracksArray = inp.getTPCTracks();
488 mTPCTrackClusIdx = inp.getTPCTracksClusterRefs();
489 mTPCClusterIdxStruct = &inp.inputsTPCclusters->clusterIndex;
490 mTPCRefitterShMap = inp.clusterShMapTPC;
491 mTPCRefitterOccMap = inp.occupancyMapTPC;
494 mTPCTrkLabels = inp.getTPCTracksMCLabels();
497 int ntr = mTPCTracksArray.size(), ntrW = 0.7 * ntr;
499 mTPCWork.reserve(ntrW);
501 mTPCLblWork.reserve(ntrW);
507 mTPCRefitter = std::make_unique<o2::gpu::GPUO2InterfaceRefit>(mTPCClusterIdxStruct, mTPCCorrMapsHelper, mBz, mTPCTrackClusIdx.data(), 0, mTPCRefitterShMap.data(), mTPCRefitterOccMap.data(), mTPCRefitterOccMap.size(),
nullptr,
o2::base::Propagator::Instance());
508 mTPCRefitter->setTrackReferenceX(900);
509 mNTPCOccBinLength = mTPCRefitter->getParam()->rec.tpc.occupancyMapTimeBins;
511 if (mNTPCOccBinLength > 1 && mTPCRefitterOccMap.size()) {
512 mNTPCOccBinLengthInv = 1. / mNTPCOccBinLength;
515 mTBinClOcc.resize(nTPCOccBins);
516 std::vector<float> mltHistTB(nTPCOccBins);
517 float sm = 0., tb = 0.5 * mNTPCOccBinLength;
518 for (
int i = 0;
i < nTPCOccBins;
i++) {
519 mltHistTB[
i] = mTPCRefitter->getParam()->GetUnscaledMult(tb);
520 tb += mNTPCOccBinLength;
522 for (
int i = nTPCOccBins;
i--;) {
524 if (
i + sumBins < nTPCOccBins) {
525 sm -= mltHistTB[
i + sumBins];
530 mTBinClOcc.resize(1);
533 auto creator = [
this](
auto& trk,
GTrackID gid,
float time0,
float terr) {
534 if constexpr (isITSTrack<decltype(trk)>()) {
539 }
else if (std::abs(trk.getQ2Pt()) > mMinTPCTrackPtInv) {
544 if constexpr (isTPCTrack<decltype(trk)>()) {
546 if (!this->mSkipTPCOnly && trk.getNClusters() > 0) {
547 resAdd = this->addTPCSeed(trk, this->tpcTimeBin2MUS(time0), this->tpcTimeBin2MUS(terr), gid, (tpcIndex = gid.getIndex()));
550 if constexpr (isTPCTOFTrack<decltype(trk)>()) {
552 resAdd = this->addTPCSeed(trk, time0, terr, gid, (tpcIndex = this->mRecoCont->
getTPCContributorGID(gid)));
554 if constexpr (isTRDTrack<decltype(trk)>()) {
556 resAdd = this->addTPCSeed(trk, time0, terr, gid, (tpcIndex = this->mRecoCont->
getTPCContributorGID(gid)));
558#ifdef _ALLOW_DEBUG_TREES_
569 int nITSROFs = mITSROFTimes.size();
572 auto& indexCache = mTPCSectIndexCache[sec];
574 LOG(info) <<
"Sorting sector" << sec <<
" | " << indexCache.size() <<
" TPC tracks";
576 if (!indexCache.size()) {
579 std::sort(indexCache.begin(), indexCache.end(), [
this](
int a,
int b) {
580 auto& trcA = mTPCWork[a];
581 auto& trcB = mTPCWork[b];
582 return (trcA.tBracket.getMax() - trcB.tBracket.getMax()) < 0.;
588 float tmax = mTPCWork[indexCache.back()].tBracket.getMax();
589 if (maxTime < tmax) {
592 int nbins = 1 + (mITSTriggered ? time2ITSROFrameTrig(tmax, 0) : time2ITSROFrameCont(tmax));
593 auto& timeStart = mTPCTimeStart[sec];
594 timeStart.resize(nbins, -1);
598 for (
int itr = 0; itr < (
int)indexCache.size(); itr++) {
599 auto& trc = mTPCWork[indexCache[itr]];
600 while (itsROF < nITSROFs && !(trc.tBracket < mITSROFTimes[itsROF])) {
603 int itsROFMatch = itsROF;
604 if (itsROFMatch && timeStart[--itsROFMatch] == -1) {
605 timeStart[itsROFMatch] = itr;
608 for (
int i = 1;
i < nbins;
i++) {
609 if (timeStart[
i] == -1) {
610 timeStart[
i] = timeStart[
i - 1];
631 mInteractionMUSLUT.clear();
633 mTimer[SWPrepTPC].Stop();
634 return mTPCWork.size() > 0;
638bool MatchTPCITS::prepareITSData()
640 static size_t errCount = 0;
641 constexpr size_t MaxErrors2Report = 10;
643 mTimer[SWPrepITS].Start(
false);
644 const auto& inp = *mRecoCont;
647 mITSClusterROFRec = inp.getITSClustersROFRecords();
648 const auto clusITS = inp.getITSClusters();
649 if (mITSClusterROFRec.empty() || clusITS.empty()) {
650 LOG(info) <<
"No ITS clusters";
653 const auto patterns = inp.getITSClustersPatterns();
654 auto pattIt = patterns.begin();
655 mITSClustersArray.reserve(clusITS.size());
656#ifdef ENABLE_UPGRADES
668 mITSClusterSizes.reserve(clusITS.size());
669 auto pattIt2 = patterns.begin();
670 for (
auto& clus : clusITS) {
671 auto pattID = clus.getPatternID();
673#ifdef ENABLE_UPGRADES
683#ifdef ENABLE_UPGRADES
685 npix = mIT3Dict->getNpixels(pattID, ib);
693 mITSClusterSizes.push_back(std::clamp(npix, 0u, 255u));
697 mITSClsLabels = inp.mcITSClusters.get();
701 mITSTracksArray = inp.getITSTracks();
702 mITSTrackClusIdx = inp.getITSTracksClusterRefs();
703 mITSTrackROFRec = inp.getITSTracksROFRecords();
705 mITSTrkLabels = inp.getITSTracksMCLabels();
707 int nROFs = mITSTrackROFRec.size();
708 mITSWork.reserve(mITSTracksArray.size());
711 const auto& lastClROF = mITSClusterROFRec[nROFs - 1];
712 int nITSClus = lastClROF.getFirstEntry() + lastClROF.getNEntries();
713 mABClusterLinkIndex.resize(nITSClus,
MinusOne);
715 mITSTimeStart[sec].resize(nROFs, -1);
720 trackLTInt.setTimeNotNeeded();
722 for (
int irof = 0; irof < nROFs; irof++) {
723 const auto& rofRec = mITSTrackROFRec[irof];
724 long nBC = rofRec.getBCData().differenceInBC(mStartIR);
725 if (nBC > maxBCs || nBC < 0) {
726 if (++errCount < MaxErrors2Report) {
727 LOGP(alarm,
"ITS ROF#{} start {} is not compatible with TF 1st orbit {} or TF length of {} HBFs",
728 irof, rofRec.getBCData().asString(), mStartIR.
asString(), nHBF);
734 if (!mITSTriggered) {
735 size_t irofCont = nBC / mITSROFrameLengthInBC;
736 if (mITSTrackROFContMapping.size() <= irofCont) {
737 mITSTrackROFContMapping.resize((1 + irofCont / 128) * 128, 0);
739 mITSTrackROFContMapping[irofCont] = irof;
742 mITSROFTimes.emplace_back(tMin, tMax);
745 mITSTimeStart[sec][irof] = mITSSectIndexCache[sec].size();
748 int trlim = rofRec.getFirstEntry() + rofRec.getNEntries();
749 for (
int it = rofRec.getFirstEntry(); it < trlim; it++) {
750 const auto& trcOrig = mITSTracksArray[it];
752 flagUsedITSClusters(trcOrig);
754 if (trcOrig.getParamOut().getX() < 1.) {
757 if (std::abs(trcOrig.getQ2Pt()) > mMinITSTrackPtInv) {
760 int nWorkTracks = mITSWork.size();
762 auto& trc = mITSWork.emplace_back(
TrackLocITS{trcOrig.getParamOut(), {tMin, tMax}, it, irof,
MinusOne});
768 trackLTInt.clearFast();
769 if (!propagateToRefX(trc, &trackLTInt)) {
773 trc.xrho = trackLTInt.getXRho();
774 trc.dL = trackLTInt.getL();
777 mITSLblWork.emplace_back(mITSTrkLabels[it]);
782 mITSSectIndexCache[sector].push_back(nWorkTracks);
790 float trcY = trc.getY(), tgp = trc.getSnp();
791 tgp /= std::sqrt((1.f - tgp) * (1.f + tgp));
793 float dyUpDn[2] = {std::abs((YMaxAtXMatchingRef - trcY) / (tgp + Tan70)), std::abs((YMaxAtXMatchingRef + trcY) / (tgp - Tan70))};
795 int sel = dyUpDn[0] < dyUpDn[1] ? 0 : 1;
796 if (dyUpDn[sel] < mSectEdgeMargin) {
798 addLastTrackCloneForNeighbourSector(sectNeib, &trackLTInt);
803 if (!mITSTriggered) {
804 int nr = mITSTrackROFContMapping.size();
805 for (
int i = 1;
i < nr;
i++) {
806 if (mITSTrackROFContMapping[
i] < mITSTrackROFContMapping[
i - 1]) {
807 mITSTrackROFContMapping[
i] = mITSTrackROFContMapping[
i - 1];
815 auto& indexCache = mITSSectIndexCache[sec];
817 LOG(info) <<
"Sorting sector" << sec <<
" | " << indexCache.size() <<
" ITS tracks";
819 if (!indexCache.size()) {
822 std::sort(indexCache.begin(), indexCache.end(), [
this](
int a,
int b) {
823 auto& trackA = mITSWork[a];
824 auto& trackB = mITSWork[b];
825 if (trackA.tBracket.getMin() < trackB.tBracket.getMin()) {
827 }
else if (trackA.tBracket.getMin() > trackB.tBracket.getMin()) {
830 return trackA.getTgl() < trackB.getTgl();
834 mTimer[SWPrepITS].Stop();
840bool MatchTPCITS::prepareFITData()
845 prepareInteractionTimes();
851void MatchTPCITS::doMatching(
int sec)
854 auto& cacheITS = mITSSectIndexCache[sec];
855 auto& cacheTPC = mTPCSectIndexCache[sec];
856 auto& timeStartTPC = mTPCTimeStart[sec];
857 auto& timeStartITS = mITSTimeStart[sec];
858 int nTracksTPC = cacheTPC.size(), nTracksITS = cacheITS.size();
859 if (!nTracksTPC || !nTracksITS) {
861 LOG(info) <<
"Matchng sector " << sec <<
" : N tracks TPC:" << nTracksTPC <<
" ITS:" << nTracksITS <<
" in sector " << sec;
867 float maxTDriftSafe = tpcTimeBin2MUS(mNTPCBinsFullDrift + mParams->
safeMarginTPCITSTimeBin + mTPCTimeEdgeTSafeMargin);
871 auto minROFITS = mITSWork[cacheITS.front()].roFrame;
873 if (minROFITS >=
int(timeStartTPC.size())) {
874 LOG(info) <<
"ITS min ROFrame " << minROFITS <<
" exceeds all cached TPC track ROF eqiuvalent " << cacheTPC.size() - 1;
878 int nCheckTPCControl = 0, nCheckITSControl = 0, nMatchesControl = 0;
879 int idxMinTPC = timeStartTPC[minROFITS];
884 for (
int itpc = idxMinTPC; itpc < nTracksTPC; itpc++) {
885 auto& trefTPC = mTPCWork[cacheTPC[itpc]];
888 auto tmn = trefTPC.tBracket.getMax() - maxTDriftSafe;
889 itsROBin = mITSTriggered ? time2ITSROFrameTrig(tmn, itsROBin) : time2ITSROFrameCont(tmn);
891 if (itsROBin >=
int(timeStartITS.size())) {
894 int iits0 = timeStartITS[itsROBin];
896 for (
auto iits = iits0; iits < nTracksITS; iits++) {
897 auto& trefITS = mITSWork[cacheITS[iits]];
899 LOG(
debug) <<
"TPC bracket: " << trefTPC.tBracket.asString() <<
" ITS bracket: " << trefITS.tBracket.asString() <<
" TPCtgl: " << trefTPC.getTgl() <<
" ITStgl: " << trefITS.getTgl();
900 if (trefTPC.tBracket < trefITS.tBracket) {
903 if (trefTPC.tBracket > trefITS.tBracket) {
908 auto deltaT = (trefITS.getZ() - trefTPC.getZ()) * mTPCVDriftInv;
909 auto timeCorr = trefTPC.getCorrectedTime(deltaT);
910 auto timeCorrErr = std::sqrt(trefITS.getSigmaZ2() + trefTPC.getSigmaZ2()) * t2nbs + mParams->
safeMarginTimeCorrErr;
911 if (mVDriftCalibOn) {
912 timeCorrErr += vdErrT * (250. - abs(trefITS.getZ()));
915 LOG(
debug) <<
"TPC range: " << trange.asString() <<
" ITS bracket: " << trefITS.tBracket.asString() <<
" DZ: " << (trefITS.getZ() - trefTPC.getZ()) <<
" DT: " << timeCorr;
917 auto cmpITSBracket = trefITS.tBracket.isOutside(timeCorr);
919 if (trefITS.tBracket.isOutside(trange)) {
924 timeCorr = trefITS.tBracket.getMin();
925 trange.setMin(timeCorr);
927 timeCorr = trefITS.tBracket.getMax();
928 trange.setMax(timeCorr);
937 int rejFlag = compareTPCITSTracks(trefITS, trefTPC, chi2);
939#ifdef _ALLOW_DEBUG_TREES_
972 auto irBracket = tBracket2IRBracket(trange);
973 if (irBracket.isInvalid()) {
977 if (checkInteractionCandidates && mInteractions.size()) {
979 int tmus = trange.getMin();
983 auto entStart = tmus <
int(mInteractionMUSLUT.size()) ? mInteractionMUSLUT[tmus] : (mInteractionMUSLUT.size() ? mInteractionMUSLUT.back() : 0);
984 for (
int ent = entStart; ent < (
int)mInteractions.size(); ent++) {
985 auto cmp = mInteractions[ent].tBracket.isOutside(trange);
999 registerMatchRecordTPC(cacheITS[iits], cacheTPC[itpc], chi2, matchedIC);
1004 LOG(info) <<
"Match sector " << sec <<
" N tracks TPC:" << nTracksTPC <<
" ITS:" << nTracksITS
1005 <<
" N TPC tracks checked: " << nCheckTPCControl <<
" (starting from " << idxMinTPC
1006 <<
"), checks: " << nCheckITSControl <<
", matches:" << nMatchesControl;
1008 mNMatchesControl += nMatchesControl;
1012void MatchTPCITS::suppressMatchRecordITS(
int itsID,
int tpcID)
1015 auto& tITS = mITSWork[itsID];
1016 int topID =
MinusOne, recordID = tITS.matchID;
1018 if (mMatchRecordsITS[recordID].partnerID == tpcID) {
1021 tITS.matchID = mMatchRecordsITS[recordID].nextRecID;
1023 mMatchRecordsITS[topID].nextRecID = mMatchRecordsITS[recordID].nextRecID;
1028 recordID = mMatchRecordsITS[recordID].nextRecID;
1033bool MatchTPCITS::registerMatchRecordTPC(
int iITS,
int iTPC,
float chi2,
int candIC)
1037 auto& tTPC = mTPCWork[iTPC];
1039 registerMatchRecordITS(iITS, iTPC, chi2, candIC);
1040 tTPC.
matchID = mMatchRecordsTPC.size();
1041 mMatchRecordsTPC.emplace_back(iITS, chi2,
MinusOne, candIC);
1047 auto& nextMatchRec = mMatchRecordsTPC[nextID];
1049 if (!nextMatchRec.isBetter(chi2, candIC)) {
1050 if (count < mParams->maxMatchCandidates) {
1053 suppressMatchRecordITS(nextMatchRec.partnerID, iTPC);
1054 registerMatchRecordITS(iITS, iTPC, chi2, candIC);
1056 nextMatchRec.chi2 = chi2;
1057 nextMatchRec.matchedIC = candIC;
1058 nextMatchRec.partnerID = iITS;
1063 nextID = nextMatchRec.nextRecID;
1069 if (count < mParams->maxMatchCandidates) {
1071 topID = tTPC.
matchID = mMatchRecordsTPC.size();
1073 topID = mMatchRecordsTPC[topID].nextRecID = mMatchRecordsTPC.size();
1076 registerMatchRecordITS(iITS, iTPC, chi2, candIC);
1077 mMatchRecordsTPC.emplace_back(iITS, chi2, nextID, candIC);
1082 suppressMatchRecordITS(mMatchRecordsTPC[nextID].partnerID, iTPC);
1084 nextID = mMatchRecordsTPC[topID].nextRecID = mMatchRecordsTPC[nextID].nextRecID;
1089 nextID = mMatchRecordsTPC[nextID].nextRecID;
1098void MatchTPCITS::registerMatchRecordITS(
const int iITS,
int iTPC,
float chi2,
int candIC)
1101 auto& tITS = mITSWork[iITS];
1102 int idnew = mMatchRecordsITS.size();
1103 auto& newRecord = mMatchRecordsITS.emplace_back(iTPC, chi2,
MinusOne, candIC);
1104 if (tITS.matchID < 0) {
1105 tITS.matchID = idnew;
1110 int topID =
MinusOne, nextRecord = tITS.matchID;
1112 auto& nextMatchRec = mMatchRecordsITS[nextRecord];
1113 if (!nextMatchRec.isBetter(chi2, candIC)) {
1114 newRecord.nextRecID = nextRecord;
1116 tITS.matchID = idnew;
1118 mMatchRecordsITS[topID].nextRecID = idnew;
1123 nextRecord = mMatchRecordsITS[nextRecord].nextRecID;
1127 mMatchRecordsITS[topID].nextRecID = idnew;
1144 if (mVDriftCalibOn) {
1146 err2Tgl += addErr * addErr;
1148 diff *= diff / err2Tgl;
1153 bool testOtherPID =
false;
1154 float itsParam[5] = {tITS.getY(), tITS.getZ(), tITS.getSnp(), tITS.getTgl(), tITS.getQ2Pt()};
1155 if (tTPC.getPID() > tITS.getPID() && tITS.
dL > 0.f && tTPC.getP2() / tTPC.getPID().getMass2() < mParams->
minBetaGammaForPIDDiff) {
1157 tPID.setPID(tTPC.getPID(),
true);
1158 if (!tPID.correctForELoss(tITS.
xrho)) {
1168 testOtherPID =
true;
1210 auto tITSAlt = tITS;
1211 tITSAlt.setPID(tTPC.getPID());
1215 chi2 = getPredictedChi2NoZ(tITSAlt, tTPC);
1217 chi2 = getPredictedChi2NoZ(tITS, tTPC);
1230 int ntpc = mTPCWork.size();
1231 printf(
"\n\nPrinting all TPC -> ITS matches for %d TPC tracks\n", ntpc);
1232 for (
int i = 0;
i < ntpc;
i++) {
1233 const auto& tTPC = mTPCWork[
i];
1234 int nm = getNMatchRecordsTPC(tTPC);
1235 printf(
"*** trackTPC#%6d %6d : Ncand = %d Best = %d\n",
i, tTPC.
sourceID, nm, tTPC.
matchID);
1238 const auto& rcTPC = mMatchRecordsTPC[recID];
1239 const auto& tITS = mITSWork[rcTPC.partnerID];
1240 printf(
" * cand %2d : ITS track %6d(src:%6d) Chi2: %.2f\n",
count, rcTPC.partnerID, tITS.
sourceID, rcTPC.chi2);
1242 recID = rcTPC.nextRecID;
1251 int nits = mITSWork.size();
1252 printf(
"\n\nPrinting all ITS -> TPC matches for %d ITS tracks\n", nits);
1254 for (
int i = 0;
i < nits;
i++) {
1255 const auto& tITS = mITSWork[
i];
1256 printf(
"*** trackITS#%6d %6d : Ncand = %d Best = %d\n",
i, tITS.
sourceID, getNMatchRecordsITS(tITS), tITS.
matchID);
1259 const auto& rcITS = mMatchRecordsITS[recID];
1260 const auto& tTPC = mTPCWork[rcITS.partnerID];
1261 printf(
" * cand %2d : TPC track %6d(src:%6d) Chi2: %.2f\n",
count, rcITS.partnerID, tTPC.
sourceID, rcITS.chi2);
1263 recID = rcITS.nextRecID;
1285 covMat(0, 0) =
static_cast<double>(trITS.getSigmaY2()) +
static_cast<double>(trTPC.getSigmaY2());
1286 covMat(1, 0) =
static_cast<double>(trITS.getSigmaSnpY()) +
static_cast<double>(trTPC.getSigmaSnpY());
1287 covMat(1, 1) =
static_cast<double>(trITS.getSigmaSnp2()) +
static_cast<double>(trTPC.getSigmaSnp2());
1288 covMat(2, 0) =
static_cast<double>(trITS.getSigmaTglY()) +
static_cast<double>(trTPC.getSigmaTglY());
1289 covMat(2, 1) =
static_cast<double>(trITS.getSigmaTglSnp()) +
static_cast<double>(trTPC.getSigmaTglSnp());
1290 covMat(2, 2) =
static_cast<double>(trITS.getSigmaTgl2()) +
static_cast<double>(trTPC.getSigmaTgl2());
1291 if (mVDriftCalibOn) {
1293 covMat(2, 2) += addErr * addErr;
1295 covMat(3, 0) =
static_cast<double>(trITS.getSigma1PtY()) +
static_cast<double>(trTPC.getSigma1PtY());
1296 covMat(3, 1) =
static_cast<double>(trITS.getSigma1PtSnp()) +
static_cast<double>(trTPC.getSigma1PtSnp());
1297 covMat(3, 2) =
static_cast<double>(trITS.getSigma1PtTgl()) +
static_cast<double>(trTPC.getSigma1PtTgl());
1298 covMat(3, 3) =
static_cast<double>(trITS.getSigma1Pt2()) +
static_cast<double>(trTPC.getSigma1Pt2());
1299 if (!covMat.Invert()) {
1300 LOG(error) <<
"Cov.matrix inversion failed: " << covMat;
1303 double chi2diag = 0., chi2ndiag = 0.,
1309 chi2diag += diff[
i] * diff[
i] * covMat(
i,
i);
1310 for (
int j =
i;
j--;) {
1311 chi2ndiag += diff[
i] * diff[
j] * covMat(
i,
j);
1314 return chi2diag + 2. * chi2ndiag;
1323 mITSWork.push_back(mITSWork.back());
1324 auto& trc = mITSWork.back();
1328 mITSSectIndexCache[sector].push_back(mITSWork.size() - 1);
1330 mITSWork.back().setCloneBefore();
1332 mITSWork.back().xrho = trackLTInt->getXRho();
1333 mITSWork.back().dL = trackLTInt->getL();
1335 mITSWork[mITSWork.size() - 2].setCloneAfter();
1337 mITSLblWork.emplace_back(mITSTrkLabels[trc.sourceID]);
1340 mITSWork.pop_back();
1349 constexpr float TgHalfSector = 0.17632698f;
1350 bool refReached =
false;
1358 if (fabs(trc.getY()) < mParams->
XMatchingRef * TgHalfSector) {
1362 if (!trialsLeft--) {
1366 if (!trc.rotate(alphaNew) != 0) {
1370 return refReached && std::abs(trc.getSnp()) < MaxSnp;
1377 printf(
"\n******************** TPC-ITS matching component ********************\n");
1379 printf(
"init is not done yet\n");
1383 printf(
"MC truth: %s\n", mMCTruthON ?
"on" :
"off");
1384 printf(
"Matching reference X: %.3f\n", mParams->
XMatchingRef);
1385 printf(
"Account Z dimension: %s\n", mCompareTracksDZ ?
"on" :
"off");
1387 printf(
"Max number ITS candidates per TPC track: %d\n", mParams->
maxMatchCandidates);
1388 printf(
"Crude cut on track params: ");
1394 printf(
"NSigma^2 cut on track params: ");
1402#ifdef _ALLOW_DEBUG_TREES_
1404 printf(
"Output debug tree (%s) file: %s\n", mDBGFlags ?
"on" :
"off", mDebugTreeFileName.data());
1406 printf(
"Debug stream flags:\n");
1411 printf(
"* winner matches\n");
1416 printf(
"**********************************************************************\n");
1423 mTimer[SWRefit].Start(
false);
1424 matchedTracks.reserve(mNMatches + mABWinnersIDs.size());
1425 matchedTracks.resize(mNMatches);
1427 matchLabels.reserve(mNMatches + mABWinnersIDs.size());
1428 matchLabels.resize(mNMatches);
1430 if (mVDriftCalibOn) {
1431 calib.reserve(mNCalibPrelim * 1.2 + 1);
1433 std::vector<int> tpcToFit;
1434 tpcToFit.reserve(mNMatches);
1435 for (
int iTPC = 0; iTPC < (
int)mTPCWork.size(); iTPC++) {
1436 const auto& tTPC = mTPCWork[iTPC];
1437 if (!isDisabledTPC(tTPC) && tTPC.
gid.testBit(0)) {
1438 tpcToFit.push_back(iTPC);
1441 LOG(
debug) <<
"Refitting winner matches";
1442 mWinnerChi2Refit.resize(mITSWork.size(), -1.f);
1443 int nToFit = (
int)tpcToFit.size();
1444 unsigned int nFailedRefit{0};
1447#pragma omp parallel for schedule(dynamic) num_threads(mNThreads) \
1448 reduction(+ : nFailedRefit)
1450 for (
int ifit = 0; ifit < nToFit; ifit++) {
1451 int iTPC = tpcToFit[ifit], iITS;
1452 const auto& tTPC = mTPCWork[iTPC];
1453 if (
refitTrackTPCITS(ifit, iTPC, iITS, matchedTracks, matchLabels, calib)) {
1454 mWinnerChi2Refit[iITS] = matchedTracks.back().getChi2Refit();
1459 LOGP(info,
"Failed {} TPC-ITS refits out of {}", nFailedRefit, nToFit);
1464 for (
int ifit = 0; ifit < nToFit; ifit++) {
1465 int itpc = tpcToFit[ifit];
1466 if (!matchedTracks[ifit].
isValid()) {
1467 while (--last > ifit && !matchedTracks[last].
isValid()) {
1470 matchedTracks[ifit] = matchedTracks[last];
1471 matchedTracks[last].invalidate();
1472 itpc = tpcToFit[last];
1474 matchLabels[ifit] = matchLabels[last];
1480 if (mDBGOut || mVDriftCalibOn) {
1486 matchedTracks.resize(mNMatches);
1488 matchLabels.resize(mNMatches);
1490 mTimer[SWRefit].Stop();
1496 const auto& tTPC = mTPCWork[iTPC];
1497 int iITS = mMatchRecordsTPC[tTPC.
matchID].partnerID;
1498 const auto& tITS = mITSWork[iITS];
1499 float minDiffFT0 = -999.,
timeC = 0.f;
1500 std::vector<float> dtimes;
1502 if (fillVDCalib || mDBGOut) {
1505 if (mInteractions.size()) {
1506 int timeC0 =
timeC - minDiffA;
1510 auto entStart = timeC0 <
int(mInteractionMUSLUT.size()) ? mInteractionMUSLUT[timeC0] : (mInteractionMUSLUT.size() ? mInteractionMUSLUT.back() : 0);
1511 for (
int ent = entStart; ent < (
int)mInteractions.size(); ent++) {
1512 float diff = mInteractions[ent].tBracket.mean() -
timeC;
1513 if (diff > minDiffA) {
1516 if (diff < -minDiffA) {
1519 dtimes.push_back(diff);
1521 minDiffA = std::abs(minDiffFT0);
1526 calib.emplace_back(tITS.getTgl(), tTPC.getTgl(), minDiffFT0);
1528#ifdef _ALLOW_DEBUG_TREES_
1532 itsRefPIDCorr.setX(0);
1533 if (tTPC.getPID() > tITS.getPID() && tITS.
dL > 0.f && tTPC.getP2() / tTPC.getPID().getMass2() < mParams->
minBetaGammaForPIDDiff) {
1534 itsRefAltPID = mITSTracksArray[tITS.
sourceID].getParamOut();
1535 itsRefPIDCorr.setPID(tTPC.getPID(),
true);
1536 itsRefPIDCorr = itsRefAltPID;
1538 if (!itsRefPIDCorr.correctForELoss(tITS.
xrho)) {
1539 itsRefPIDCorr.setX(-10);
1541 float q2ptPID = itsRefPIDCorr.getQ2Pt();
1543 itsRefPIDCorr = tITS;
1544 itsRefPIDCorr.setPID(tTPC.getPID(),
true);
1545 itsRefPIDCorr.setQ2Pt(q2ptPID);
1546 auto snp = tITS.getSnp() + dCurvL;
1547 if (std::abs(snp) >= 1.) {
1548 snp = std::copysign(0.99, snp);
1550 itsRefPIDCorr.setSnp(snp);
1551 itsRefPIDCorr.setY(tITS.getY() + dCurvL * dLEff * 0.5);
1555 itsRefAltPID.setX(-10);
1558 int tb = mTPCTracksArray[tTPC.
sourceID].getTime0() * mNTPCOccBinLengthInv;
1559 float mltTPC = tb < 0 ? mTBinClOcc[0] : (tb >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[tb]);
1560 (*mDBGOut) <<
"refit"
1561 <<
"tpcOrig=" << mTPCTracksArray[tTPC.
sourceID] <<
"itsOrig=" << mITSTracksArray[tITS.
sourceID] <<
"itsRef=" << tITS <<
"tpcRef=" << tTPC <<
"matchRefit=" <<
match
1562 <<
"timeCorr=" <<
timeC <<
"dTimeFT0=" << minDiffFT0 <<
"dTimes=" << dtimes
1563 <<
"itsRefAltPID=" << itsRefAltPID <<
"itsRefPIDCorr=" << itsRefPIDCorr;
1565 (*mDBGOut) <<
"refit"
1566 <<
"itsLbl=" << mITSLblWork[iITS] <<
"tpcLbl=" << mTPCLblWork[iTPC];
1568 (*mDBGOut) <<
"refit"
1569 <<
"multTPC=" << mltTPC
1570 <<
"multITSTr=" << mITSTrackROFRec[tITS.
roFrame].getNEntries()
1571 <<
"multITSCl=" << mITSClusterROFRec[tITS.
roFrame].getNEntries()
1572 <<
"tf=" << mTFCount <<
"\n";
1582 const float maxStep = 2.f;
1583 const auto& tTPC = mTPCWork[iTPC];
1584 const auto& tpcMatchRec = mMatchRecordsTPC[tTPC.
matchID];
1585 iITS = tpcMatchRec.partnerID;
1586 const auto& tITS = mITSWork[iITS];
1587 const auto& itsTrOrig = mITSTracksArray[tITS.
sourceID];
1588 auto& trfit = matchedTracks[slot];
1591 trfit.getParamOut().setUserField(0);
1592 trfit.setPID(tTPC.getPID(),
true);
1593 trfit.getParamOut().setPID(tTPC.getPID(),
true);
1596 if (!mCompareTracksDZ) {
1597 trfit.setZ(tITS.getZ());
1599 float deltaT = (trfit.getZ() - tTPC.getZ()) * mTPCVDriftInv;
1602 timeErr = mITSTimeResMUS;
1625 int nclRefit = 0, ncl = itsTrOrig.getNumberOfClusters();
1629 int clEntry = itsTrOrig.getFirstClusterEntry();
1633 if (mVDriftCalibOn) {
1639 for (
int icl = 0; icl < ncl; icl++) {
1640 const auto& clus = mITSClustersArray[mITSTrackClusIdx[clEntry++]];
1641 float alpha = geom->getSensorRefAlpha(clus.getSensorID()),
x = clus.getX();
1642 if (!trfit.rotate(
alpha) ||
1646 !propagator->propagateToX(trfit,
x, propagator->getNominalBz(),
1647 MaxSnp, maxStep, mUseMatCorrFlag, &trfit.getLTIntegralOut())) {
1650 chi2 += trfit.getPredictedChi2(clus);
1651 if (!trfit.update(clus)) {
1656 if (nclRefit != ncl) {
1657 LOGP(
debug,
"Refit in ITS failed after ncl={}, match between TPC track #{} and ITS track #{}", nclRefit, tTPC.
sourceID, tITS.
sourceID);
1658 LOGP(
debug,
"{:s}", trfit.asString());
1668 if (!propagator->propagateToDCA(vtxDummy.getXYZ(), trpar, propagator->getNominalBz(),
1669 maxStep, MatCorrType::USEMatCorrNONE,
nullptr, &trfit.getLTIntegralOut())) {
1670 LOG(error) <<
"LTOF integral might be incorrect";
1672 auto& tofL = trfit.getLTIntegralOut();
1673 tofL.setX2X0(-tofL.getX2X0());
1674 tofL.setXRho(-tofL.getXRho());
1677 auto& tracOut = trfit.getParamOut();
1678 if (tTPC.getPID() > tITS.getPID() && tTPC.getP2() / tTPC.getPID().getMass2() < mParams->
minBetaGammaForPIDDiff) {
1680 tracOut = mITSTracksArray[tITS.
sourceID].getParamOut();
1681 tracOut.setPID(tTPC.getPID(),
true);
1683 LOGP(
debug,
"Failed to rotate ITSouter with imposed PID to TPC alpha {} or propagate to X={}: {:s}", tTPC.getAlpha(), mParams->
XMatchingRef, tracOut.asString());
1691 !propagator->PropagateToXBxByBz(tracOut, xtogo, MaxSnp, 10., mUseMatCorrFlag, &tofL)) {
1692 LOG(
debug) <<
"Propagation to inner TPC boundary X=" << xtogo <<
" failed, Xtr=" << tracOut.getX() <<
" snp=" << tracOut.getSnp();
1696 if (mVDriftCalibOn) {
1700 auto posStart = tracOut.getXYZGlo();
1701 auto tImposed =
timeC * mTPCTBinMUSInv;
1702 if (std::abs(tImposed - mTPCTracksArray[tTPC.
sourceID].getTime0()) > 550) {
1703 LOGP(alarm,
"Impossible imposed timebin {} for TPC track time0:{}, dBwd:{} dFwd:{} TB | ZShift:{}, TShift:{}", tImposed, mTPCTracksArray[tTPC.
sourceID].getTime0(),
1704 mTPCTracksArray[tTPC.
sourceID].getDeltaTBwd(), mTPCTracksArray[tTPC.
sourceID].getDeltaTFwd(), trfit.getZ() - tTPC.getZ(), deltaT);
1705 LOGP(info,
"Trc: {}", mTPCTracksArray[tTPC.
sourceID].asString());
1709 int retVal = mTPCRefitter->RefitTrackAsTrackParCov(tracOut, mTPCTracksArray[tTPC.
sourceID].getClusterRef(), tImposed, &chi2Out,
true,
false);
1715 auto posEnd = tracOut.getXYZGlo();
1716 auto lInt = propagator->estimateLTIncrement(tracOut, posStart, posEnd);
1717 tofL.addStep(lInt, tracOut.getQ2P2());
1718 tofL.addX2X0(lInt * mTPCmeanX0Inv);
1722 if (trackTune.useTPCOuterCorr) {
1723 tracOut.updateParams(trackTune.tpcParOuter);
1729 trfit.setChi2Match(tpcMatchRec.chi2);
1730 trfit.setChi2Refit(chi2);
1731 trfit.setTimeMUS(
timeC, timeErr);
1736 matchLabels[slot] = mTPCLblWork[iTPC];
1737 matchLabels[slot].setFakeFlag(mITSLblWork[iITS] != mTPCLblWork[iTPC]);
1748 const float maxStep = 2.f;
1749 const auto& tTPC = mTPCWork[seed.
tpcWID];
1751 auto& newtr = matchedTracks.emplace_back(winLink, winLink);
1752 auto& tracOut = newtr.getParamOut();
1753 auto& tofL = newtr.getLTIntegralOut();
1756 tracOut.resetCovariance();
1757 propagator->estimateLTFast(tofL, winLink);
1759 const auto& itsClRefs = ABTrackletRefs[iITSAB];
1760 int nclRefit = 0, ncl = itsClRefs.getNClusters();
1764 for (
int icl = itsClRefs.getFirstEntry(); icl < itsClRefs.getEntriesBound(); icl++) {
1765 const auto& clus = mITSClustersArray[ABTrackletClusterIDs[icl]];
1766 float alpha = geom->getSensorRefAlpha(clus.getSensorID()),
x = clus.getX();
1767 if (!tracOut.rotate(
alpha) ||
1771 !propagator->propagateToX(tracOut,
x, propagator->getNominalBz(), MaxSnp, maxStep, mUseMatCorrFlag, &tofL)) {
1774 chi2 += tracOut.getPredictedChi2(clus);
1775 if (!tracOut.update(clus)) {
1780 if (nclRefit != ncl) {
1781 LOGP(
debug,
"AfterBurner refit in ITS failed after ncl={}, match between TPC track #{} and ITS tracklet #{}", nclRefit, tTPC.
sourceID, iITSAB);
1782 LOGP(
debug,
"{:s}", tracOut.asString());
1783 matchedTracks.pop_back();
1787 float timeC = mInteractions[seed.
ICCanID].tBracket.mean();
1788 float timeErr = mInteractions[seed.
ICCanID].tBracket.delta();
1792 !propagator->PropagateToXBxByBz(tracOut, xtogo, MaxSnp, 10., mUseMatCorrFlag, &tofL)) {
1793 LOG(
debug) <<
"Propagation to inner TPC boundary X=" << xtogo <<
" failed, Xtr=" << tracOut.getX() <<
" snp=" << tracOut.getSnp();
1794 matchedTracks.pop_back();
1798 auto posStart = tracOut.getXYZGlo();
1799 int retVal = mTPCRefitter->RefitTrackAsTrackParCov(tracOut, mTPCTracksArray[tTPC.
sourceID].getClusterRef(),
timeC * mTPCTBinMUSInv, &chi2Out,
true,
false);
1802 matchedTracks.pop_back();
1805 auto posEnd = tracOut.getXYZGlo();
1806 auto lInt = propagator->estimateLTIncrement(tracOut, posStart, posEnd);
1807 tofL.addStep(lInt, tracOut.getQ2P2());
1808 tofL.addX2X0(lInt * mTPCmeanX0Inv);
1811 if (trackTune.useTPCOuterCorr) {
1812 tracOut.updateParams(trackTune.tpcParOuter);
1819 newtr.setChi2Match(winLink.chi2Norm());
1820 newtr.setChi2Refit(chi2);
1821 newtr.setTimeMUS(
timeC, timeErr);
1829bool MatchTPCITS::refitTPCInward(
o2::track::TrackParCov& trcIn,
float& chi2,
float xTgt,
int trcID,
float timeTB)
const
1832 const auto& tpcTrOrig = mTPCTracksArray[trcID];
1834 trcIn = tpcTrOrig.getOuterParam();
1838 int retVal = mTPCRefitter->RefitTrackAsTrackParCov(trcIn, tpcTrOrig.getClusterRef(), timeTB, &chi2,
false,
true);
1840 LOG(warning) <<
"Refit failed";
1841 LOG(warning) << trcIn.asString();
1847 if (!propagator->PropagateToXBxByBz(trcIn, xTgt, MaxSnp, 2., mUseMatCorrFlag)) {
1848 LOG(
debug) <<
"Propagation to target X=" << xTgt <<
" failed, Xtr=" << trcIn.getX() <<
" snp=" << trcIn.getSnp() <<
" pT=" << trcIn.getPt();
1857int MatchTPCITS::prepareABSeeds()
1860 const auto& outerLr = mRGHelper.
layers.back();
1862 const float ROuter = outerLr.rRange.getMax() + 0.5f;
1865 for (
int iTPC = 0; iTPC < (
int)mTPCWork.size(); iTPC++) {
1866 auto& tTPC = mTPCWork[iTPC];
1867 if (isDisabledTPC(tTPC)) {
1873 !propagator->PropagateToXBxByBz(tTPC, xTgt, MaxSnp, 2., mUseMatCorrFlag)) {
1876 mTPCABIndexCache.push_back(iTPC);
1880 std::sort(mTPCABIndexCache.begin(), mTPCABIndexCache.end(), [
this](
int a,
int b) {
1881 auto& trcA = mTPCWork[a];
1882 auto& trcB = mTPCWork[b];
1883 return (trcA.tBracket.getMin() - trcB.tBracket.getMin()) < 0.;
1886 float maxTDriftSafe = tpcTimeBin2MUS(mNTPCBinsFullDrift + mParams->
safeMarginTPCITSTimeBin + mTPCTimeEdgeTSafeMargin);
1887 int nIntCand = mInteractions.size(), nTPCCand = mTPCABIndexCache.size();
1889 for (
int ic = 0; ic < nIntCand; ic++) {
1890 int icFirstSeed = mTPCABSeeds.size();
1891 auto& intCand = mInteractions[ic];
1892 auto tic = intCand.tBracket.mean();
1893 for (
int it = tpcStart; it < nTPCCand; it++) {
1894 auto& trc = mTPCWork[mTPCABIndexCache[it]];
1895 auto cmp = trc.tBracket.isOutside(intCand.tBracket);
1900 if (trc.tBracket.getMin() + maxTDriftSafe < intCand.tBracket.getMin()) {
1906 float dt = trc.getSignedDT(tic - trc.time0);
1907 float dz = dt * mTPCVDrift,
z = trc.getZ() + dz;
1908 if (outerLr.zRange.isOutside(
z, std::sqrt(trc.getSigmaZ2()) + 2.)) {
1912 auto& seed = mTPCABSeeds.emplace_back(mTPCABIndexCache[it], ic, trc);
1915 intCand.seedsRef.set(icFirstSeed, mTPCABSeeds.size() - icFirstSeed);
1917 return mTPCABIndexCache.size();
1921int MatchTPCITS::prepareInteractionTimes()
1924 const float ft0Uncertainty = 0.5e-3;
1925 int nITSROFs = mITSROFTimes.size();
1926 if (mFITInfo.size()) {
1928 for (
const auto& ft : mFITInfo) {
1932 auto fitTime = ft.getInteractionRecord().differenceInBCMUS(mStartIR);
1936 if (
size_t(fitTime) >= mInteractionMUSLUT.size()) {
1937 mInteractionMUSLUT.resize(
size_t(fitTime) + 1, -1);
1939 if (mInteractionMUSLUT[fitTime] < 0) {
1940 mInteractionMUSLUT[fitTime] = mInteractions.size();
1942 for (; rof < nITSROFs; rof++) {
1943 if (mITSROFTimes[rof] < fitTime) {
1948 if (rof >= nITSROFs) {
1955 for (
auto&
val : mInteractionMUSLUT) {
1962 return mInteractions.size();
1969 mTimer[SWABSeeds].Start(
false);
1972 int nIntCand = mInteractions.size(), nABSeeds = mTPCABSeeds.size();
1973 LOGP(info,
"AfterBurner will check {} seeds from {} TPC tracks and {} interaction candidates with {} threads", nABSeeds, mTPCABIndexCache.size(), nIntCand, mNThreads);
1974 mTimer[SWABSeeds].Stop();
1975 if (!nIntCand || !mTPCABSeeds.size()) {
1978 mTimer[SWABMatch].Start(
false);
1980 std::vector<ITSChipClustersRefs> itsChipClRefsBuff(mNThreads);
1981#ifdef ENABLE_UPGRADES
1984 std::generate(itsChipClRefsBuff.begin(), itsChipClRefsBuff.end(), []() {
1985 return ITSChipClustersRefs(o2::its::GeometryTGeo::Instance()->getNumberOfChips());
1990#pragma omp parallel for schedule(dynamic) num_threads(mNThreads)
1992 for (
int ic = 0; ic < nIntCand; ic++) {
1993 const auto& intCand = mInteractions[ic];
1994 LOGP(
debug,
"cand T {} Entries: {} : {} : {} | ITS ROF: {}", intCand.tBracket.mean(), intCand.seedsRef.getEntries(), intCand.seedsRef.getFirstEntry(), intCand.seedsRef.getEntriesBound(), intCand.rofITS);
1995 if (!intCand.seedsRef.getEntries()) {
1999 uint8_t tid = (uint8_t)omp_get_thread_num();
2003 fillClustersForAfterBurner(intCand.rofITS, 1, itsChipClRefsBuff[tid]);
2004 for (
int is = intCand.seedsRef.getFirstEntry(); is < intCand.seedsRef.getEntriesBound(); is++) {
2005 processABSeed(is, itsChipClRefsBuff[tid], tid);
2008 mTimer[SWABMatch].Stop();
2009 mTimer[SWABWinners].Start(
false);
2016 std::vector<SID> candAB;
2017 candAB.reserve(nABSeeds);
2018 mABWinnersIDs.reserve(mTPCABIndexCache.size());
2020 for (
int i = 0;
i < nABSeeds;
i++) {
2021 auto& ABSeed = mTPCABSeeds[
i];
2022 if (ABSeed.isDisabled()) {
2025 candAB.emplace_back(SID{
i, ABSeed.getLink(ABSeed.getBestLinkID()).chi2Norm()});
2027 std::sort(candAB.begin(), candAB.end(), [](SID
a, SID
b) { return a.chi2 < b.chi2; });
2028 for (
int i = 0;
i < (
int)candAB.size();
i++) {
2029 auto& ABSeed = mTPCABSeeds[candAB[
i].seedID];
2030 if (ABSeed.isDisabled()) {
2034 auto& tTPC = mTPCWork[ABSeed.tpcWID];
2040 auto bestID = ABSeed.getBestLinkID();
2041 if (ABSeed.checkLinkHasUsedClusters(bestID, mABClusterLinkIndex)) {
2042 ABSeed.setNeedAlternative();
2046 ABSeed.validate(bestID);
2047 ABSeed.flagLinkUsedClusters(bestID, mABClusterLinkIndex);
2048 mABWinnersIDs.push_back(tTPC.
matchID = candAB[
i].seedID);
2049 mNABRefsClus += ABSeed.getNLayers();
2053 mTimer[SWABWinners].Stop();
2054 mTimer[SWABRefit].Start(
false);
2055 refitABWinners(matchedTracks, matchLabels, ABTrackletLabels, ABTrackletClusterIDs, ABTrackletRefs, calib);
2056 mTimer[SWABRefit].Stop();
2067 ABTrackletClusterIDs.reserve(mNABRefsClus);
2068 ABTrackletRefs.reserve(mABWinnersIDs.size());
2070 ABTrackletLabels.reserve(mABWinnersIDs.size());
2072 if (matchedTracks.capacity() < mABWinnersIDs.size() + matchedTracks.size()) {
2073 LOGP(warn,
"need to expand matched tracks container from {} to {}", matchedTracks.capacity(), mABWinnersIDs.size() + matchedTracks.size());
2074 matchedTracks.reserve(mABWinnersIDs.size() + matchedTracks.size());
2076 matchLabels.reserve(mABWinnersIDs.size() + matchedTracks.size());
2080 std::map<o2::MCCompLabel, int> labelOccurence;
2081 auto accountClusterLabel = [&labelOccurence, itsClLabs = mITSClsLabels](
int clID) {
2082 auto labels = itsClLabs->getLabels(clID);
2083 for (
auto lab : labels) {
2085 labelOccurence[lab]++;
2090 for (
auto wid : mABWinnersIDs) {
2091 const auto& ABSeed = mTPCABSeeds[wid];
2092 int start = ABTrackletClusterIDs.size();
2093 int lID = ABSeed.winLinkID, ncl = 0;
2094 auto& clref = ABTrackletRefs.emplace_back(
start, ncl);
2096 const auto& winL = ABSeed.getLink(lID);
2098 ABTrackletClusterIDs.push_back(winL.clID);
2100 clref.pattern |= 0x1 << winL.layerID;
2101 clref.setClusterSize(winL.layerID, mITSClusterSizes[winL.clID]);
2103 accountClusterLabel(winL.clID);
2106 lID = winL.parentID;
2108 clref.setEntries(ncl);
2109 if (!
refitABTrack(ABTrackletRefs.size() - 1, ABSeed, matchedTracks, ABTrackletClusterIDs, ABTrackletRefs)) {
2110 ABTrackletRefs.pop_back();
2111 ABTrackletClusterIDs.resize(
start);
2113 labelOccurence.clear();
2129 labelOccurence.clear();
2130 ABTrackletLabels.push_back(lab);
2131 auto& lblGlo = matchLabels.emplace_back(mTPCLblWork[ABSeed.tpcWID]);
2132 lblGlo.setFakeFlag(lab != lblGlo);
2133 LOG(
debug) <<
"ABWinner ncl=" << ncl <<
" mcLBAB " << lab <<
" mcLBGlo " << lblGlo <<
" chi2=" << ABSeed.getLink(ABSeed.winLinkID).chi2Norm() <<
" pT = " << ABSeed.track.getPt();
2137 LOG(info) <<
"AfterBurner validated " << ABTrackletRefs.size() <<
" tracks";
2141void MatchTPCITS::processABSeed(
int sid,
const ITSChipClustersRefs& itsChipClRefs, uint8_t tID)
2144 auto& ABSeed = mTPCABSeeds[sid];
2145 ABSeed.threadID = tID;
2146 ABSeed.linksEntry = mABLinksPool.
threadPool[tID].size();
2149 int nextLinkID = ABSeed.firstInLr[ilr];
2150 if (nextLinkID < 0) {
2154 const auto& seedLink = ABSeed.getLink(nextLinkID);
2155 if (seedLink.isDisabled()) {
2156 nextLinkID = seedLink.nextOnLr;
2159 int next2nextLinkID = seedLink.nextOnLr;
2160 followABSeed(seedLink, itsChipClRefs, nextLinkID, ilr - 1, ABSeed);
2161 nextLinkID = next2nextLinkID;
2165 auto candID = ABSeed.getBestLinkID();
2166 if (ABSeed.isDisabled() ||
2171 mABLinksPool.
threadPool[tID].resize(
size_t(ABSeed.linksEntry));
2192 const auto& lr = mRGHelper.
layers[lrID];
2194 if (!seedC.getXatLabR(lr.rRange.getMax(), xTgt, propagator->getNominalBz(),
o2::track::DirInward) ||
2195 !propagator->propagateToX(seedC, xTgt, propagator->getNominalBz(), MaxSnp, 2., mUseMatCorrFlag)) {
2199 float zDRStep = -seedC.getTgl() * lr.rRange.delta();
2200 float errZ = std::sqrt(seedC.getSigmaZ2() + mParams->
err2ABExtraZ);
2201 if (lr.zRange.isOutside(seedC.getZ(), mParams->
nABSigmaZ * errZ + std::abs(zDRStep))) {
2204 std::vector<int> chipSelClusters;
2211 seedC.getCircleParams(propagator->getNominalBz(), trcCircle, sna, csa);
2213 seedC.getLineParams(trcLinPar, sna, csa);
2217 float phi = std::atan2(yCurr, xCurr);
2219 int nLad2Check = 0, ladIDguess = lr.getLadderID(phi), chipIDguess = lr.getChipID(seedC.getZ() + 0.5 * zDRStep);
2220 std::array<int, MaxLadderCand> lad2Check;
2221 nLad2Check = mFieldON ? findLaddersToCheckBOn(lrID, ladIDguess, trcCircle, errYFrac, lad2Check) : findLaddersToCheckBOff(lrID, ladIDguess, trcLinPar, errYFrac, lad2Check);
2223 for (
int ilad = nLad2Check; ilad--;) {
2224 int ladID = lad2Check[ilad];
2225 const auto& lad = lr.ladders[ladID];
2229 float t = 1e9, xCross, yCross;
2230 const auto& chipC = lad.chips[chipIDguess];
2232 chipC.xyEdges.circleCrossParam(trcCircle, t);
2234 chipC.xyEdges.lineCrossParam(trcLinPar, t);
2236 chipC.xyEdges.eval(t, xCross, yCross);
2237 float dx = xCross - xCurr, dy = yCross - yCurr, dst2 = dx * dx + dy * dy,
dst = sqrtf(dst2);
2239 float zCross = seedC.getZ() + seedC.getTgl() * (dst2 < 2 * (dx * xCurr + dy * yCurr) ?
dst : -
dst);
2241 for (
int ich = -1; ich < 2; ich++) {
2242 int chipID = chipIDguess + ich;
2244 if (chipID < 0 || chipID >=
static_cast<int>(lad.chips.size())) {
2247 if (lad.chips[chipID].zRange.isOutside(zCross, mParams->
nABSigmaZ * errZ)) {
2250 const auto& clRange = itsChipClRefs.
chipRefs[lad.chips[chipID].id];
2251 if (!clRange.getEntries()) {
2252 LOG(
debug) <<
"No clusters in chip range";
2256 float errYcalp = errY * (csa * chipC.csAlp + sna * chipC.snAlp);
2258 float yTrack = -xCross * chipC.snAlp + yCross * chipC.csAlp;
2259 if (!preselectChipClusters(chipSelClusters, clRange, itsChipClRefs, yTrack, zCross, tolerY, tolerZ)) {
2260 LOG(
debug) <<
"No compatible clusters found";
2265 if (!trcLC.rotate(chipC.alp) || !trcLC.propagateTo(chipC.xRef, propagator->getNominalBz())) {
2266 LOG(
debug) <<
" failed to rotate to alpha=" << chipC.alp <<
" or prop to X=" << chipC.xRef;
2271 for (
auto clID : chipSelClusters) {
2272 const auto& cls = mITSClustersArray[clID];
2273 auto chi2 = trcLC.getPredictedChi2(cls);
2277 int lnkID = registerABTrackLink(ABSeed, trcLC, clID, seedID, lrID, ladID, chi2);
2279 auto& link = ABSeed.
getLink(lnkID);
2295void MatchTPCITS::accountForOverlapsAB(
int lrSeed)
2298 LOG(warning) <<
"TODO";
2303 std::array<int, MatchTPCITS::MaxLadderCand>& lad2Check)
const
2307 const auto& lr = mRGHelper.
layers[ilr];
2308 int nacc = 0, jmp = 0;
2309 if (lr.ladders[lad0].xyEdges.seenByCircle(circle, errYFrac)) {
2310 lad2Check[nacc++] = lad0;
2312 bool doUp =
true, doDn =
true;
2315 int ldID = (lad0 + jmp) % lr.nLadders;
2316 if (lr.ladders[ldID].xyEdges.seenByCircle(circle, errYFrac)) {
2317 lad2Check[nacc++] = ldID;
2323 int ldID = lad0 - jmp;
2325 ldID += lr.nLadders;
2327 if (lr.ladders[ldID].xyEdges.seenByCircle(circle, errYFrac)) {
2328 lad2Check[nacc++] = ldID;
2339 std::array<int, MatchTPCITS::MaxLadderCand>& lad2Check)
const
2343 const auto& lr = mRGHelper.
layers[ilr];
2344 int nacc = 0, jmp = 0;
2345 if (lr.ladders[lad0].xyEdges.seenByLine(trcLinPar, errYFrac)) {
2346 lad2Check[nacc++] = lad0;
2348 bool doUp =
true, doDn =
true;
2351 int ldID = (lad0 + jmp) % lr.nLadders;
2352 if (lr.ladders[ldID].xyEdges.seenByLine(trcLinPar, errYFrac)) {
2353 lad2Check[nacc++] = ldID;
2359 int ldID = lad0 - jmp;
2361 ldID += lr.nLadders;
2363 if (lr.ladders[ldID].xyEdges.seenByLine(trcLinPar, errYFrac)) {
2364 lad2Check[nacc++] = ldID;
2391 auto& nextLink = ABSeed.
getLink(nextID);
2393 bool newIsBetter = parentID <=
MinusOne ? isBetter(chi2, nextLink.chi2) : isBetter(ABSeed.getLink(parentID).chi2NormPredict(chi2Cl), nextLink.chi2Norm());
2395 if (count < mParams->maxABLinksOnLayer) {
2396 ABSeed.
addLink(trc, clID, parentID, nextID, lr, nc, laddID, chi2);
2409 nextID = nextLink.nextOnLr;
2412 if (count < mParams->maxABLinksOnLayer) {
2432 auto tpcTrOrig = mTPCTracksArray[tTPC.
sourceID];
2455 float dDrift = (timeIC - tTPC.
time0) * mTPCBin2Z;
2472void MatchTPCITS::fillClustersForAfterBurner(
int rofStart,
int nROFs,
ITSChipClustersRefs& itsChipClRefs)
2477 int first = mITSClusterROFRec[rofStart].getFirstEntry(),
last =
first;
2478 for (
int ir = nROFs;
ir--;) {
2479 last += mITSClusterROFRec[rofStart +
ir].getNEntries();
2481 itsChipClRefs.
clear();
2482 auto& idxSort = itsChipClRefs.
clusterID;
2483 for (
int icl =
first; icl <
last; icl++) {
2484 if (mABClusterLinkIndex[icl] !=
MinusTen) {
2485 idxSort.push_back(icl);
2489 const auto& clusArr = mITSClustersArray;
2490 std::sort(idxSort.begin(), idxSort.end(), [&clusArr](
int i,
int j) {
2491 const auto &clI = clusArr[i], &clJ = clusArr[j];
2492 if (clI.getSensorID() < clJ.getSensorID()) {
2495 if (clI.getSensorID() == clJ.getSensorID()) {
2496 return clI.getZ() < clJ.getZ();
2501 int ncl = idxSort.size();
2502 int lastSens = -1, nClInSens = 0;
2504 for (
int icl = 0; icl < ncl; icl++) {
2505 const auto& clus = mITSClustersArray[idxSort[icl]];
2506 int sens = clus.getSensorID();
2507 if (sens != lastSens) {
2509 chipClRefs->setEntries(nClInSens);
2512 chipClRefs = &itsChipClRefs.
chipRefs[(lastSens = sens)];
2513 chipClRefs->setFirstEntry(icl);
2518 chipClRefs->setEntries(nClInSens);
2525 mITSTimeBiasInBC =
n;
2532 mITSROFrameLengthMUS = fums;
2533 mITSTimeResMUS = mITSROFrameLengthMUS / std::sqrt(12.f);
2534 mITSROFrameLengthMUSInv = 1. / mITSROFrameLengthMUS;
2541 mITSROFrameLengthInBC = nbc;
2543 mITSTimeResMUS = mITSROFrameLengthMUS / std::sqrt(12.f);
2544 mITSROFrameLengthMUSInv = 1. / mITSROFrameLengthMUS;
2553 if (minBC < 0 && mUseBCFilling) {
2554 if (mUseBCFilling) {
2555 mUseBCFilling =
false;
2556 LOG(warning) <<
"Disabling match validation by BunchFilling as no interacting bunches found";
2560 int bcAbove = minBC;
2565 mClosestBunchAbove[
i] = bcAbove;
2567 int bcBelow = maxBC;
2572 mClosestBunchBelow[
i] = bcBelow;
2581 if (tbrange.
getMin() > 0) {
2586 int bc = mClosestBunchAbove[irMin.bc];
2587 if (
bc < irMin.bc) {
2591 bc = mClosestBunchBelow[irMax.bc];
2592 if (
bc > irMax.bc) {
2593 if (irMax.orbit == 0) {
2600 return {irMin, irMax};
2604void MatchTPCITS::removeTPCfromITS(
int tpcID,
int itsID)
2607 auto& tITS = mITSWork[itsID];
2608 if (isValidatedITS(tITS)) {
2613 auto& rcITS = mMatchRecordsITS[next];
2614 if (rcITS.partnerID == tpcID) {
2616 tITS.
matchID = rcITS.nextRecID;
2618 mMatchRecordsITS[topID].nextRecID = rcITS.nextRecID;
2623 next = rcITS.nextRecID;
2628void MatchTPCITS::removeITSfromTPC(
int itsID,
int tpcID)
2631 auto& tTPC = mTPCWork[tpcID];
2632 if (isValidatedTPC(tTPC)) {
2637 auto& rcTPC = mMatchRecordsTPC[next];
2638 if (rcTPC.partnerID == itsID) {
2640 tTPC.
matchID = rcTPC.nextRecID;
2642 mMatchRecordsTPC[topID].nextRecID = rcTPC.nextRecID;
2647 next = rcTPC.nextRecID;
2656 for (
int icl = track.getNumberOfClusters(); icl--;) {
2657 mABClusterLinkIndex[mITSTrackClusIdx[clEntry++]] =
MinusTen;
2662int MatchTPCITS::preselectChipClusters(std::vector<int>& clVecOut,
const ClusRange& clRange,
const ITSChipClustersRefs& itsChipClRefs,
2663 float trackY,
float trackZ,
float tolerY,
float tolerZ)
const
2666 int icID = clRange.getFirstEntry();
2667 for (
int icl = clRange.getEntries(); icl--;) {
2668 int clID = itsChipClRefs.
clusterID[icID++];
2669 const auto& cls = mITSClustersArray[clID];
2670 float dz = cls.getZ() - trackZ;
2671 LOG(
debug) <<
"cl" << icl <<
'/' << clID <<
" "
2672 <<
" dZ: " << dz <<
" [" << tolerZ <<
"| dY: " << trackY - cls.getY() <<
" [" << tolerY <<
"]";
2674 float clsZ = cls.getZ();
2675 LOG(
debug) <<
"Skip the rest since " << trackZ <<
" < " << clsZ <<
"\n";
2677 }
else if (dz < -tolerZ) {
2678 LOG(
debug) <<
"Skip cluster dz=" << dz <<
" Ztr=" << trackZ <<
" zCl=" << cls.getZ();
2681 if (fabs(trackY - cls.getY()) > tolerY) {
2682 LOG(
debug) <<
"Skip cluster dy= " << trackY - cls.getY() <<
" Ytr=" << trackY <<
" yCl=" << cls.getY();
2685 clVecOut.push_back(clID);
2687 return clVecOut.size();
2698 size_t sizTotShm = 0, capTotShm = 0, sizTot = 0, capTot = 0, siz = 0, cap = 0, cnt = 0, cntCap = 0;
2704 LOGP(info,
"Size SHM, matchedTracks : size {:9} cap {:9}", siz, cap);
2710 LOGP(info,
"Size SHM, ABTrackletRefs : size {:9} cap {:9}", siz, cap);
2712 siz = ABTrackletClusterIDs.size() *
sizeof(
int);
2713 cap = ABTrackletClusterIDs.capacity() *
sizeof(
int);
2716 LOGP(info,
"Size SHM, ABTrackletClusterIDs : size {:9} cap {:9}", siz, cap);
2722 LOGP(info,
"Size SHM, matchLabels : size {:9} cap {:9}", siz, cap);
2728 LOGP(info,
"Size SHM, ABTrackletLabels : size {:9} cap {:9}", siz, cap);
2734 LOGP(info,
"Size SHM, calib : size {:9} cap {:9}", siz, cap);
2737 siz = mITSClustersArray.size() *
sizeof(
ITSCluster);
2738 cap = mITSClustersArray.capacity() *
sizeof(
ITSCluster);
2741 LOGP(info,
"Size RSS, mITSClustersArray : size {:9} cap {:9}", siz, cap);
2743 siz = mMatchRecordsTPC.size() *
sizeof(
MatchRecord);
2744 cap = mMatchRecordsTPC.capacity() *
sizeof(
MatchRecord);
2747 LOGP(info,
"Size RSS, mMatchRecordsTPC : size {:9} cap {:9}", siz, cap);
2749 siz = mMatchRecordsITS.size() *
sizeof(
MatchRecord);
2750 cap = mMatchRecordsITS.capacity() *
sizeof(
MatchRecord);
2753 LOGP(info,
"Size RSS, mMatchRecordsITS : size {:9} cap {:9}", siz, cap);
2755 siz = mITSROFTimes.size() *
sizeof(
BracketF);
2756 cap = mITSROFTimes.capacity() *
sizeof(
BracketF);
2759 LOGP(info,
"Size RSS, mITSROFTimes : size {:9} cap {:9}", siz, cap);
2765 LOGP(info,
"Size RSS, mTPCWork : size {:9} cap {:9}", siz, cap);
2771 LOGP(info,
"Size RSS, mITSWork : size {:9} cap {:9}", siz, cap);
2773 siz = mWinnerChi2Refit.size() *
sizeof(float);
2774 cap = mWinnerChi2Refit.capacity() *
sizeof(float);
2777 LOGP(info,
"Size RSS, mWinnerChi2Refit : size {:9} cap {:9}", siz, cap);
2779 siz = mTPCABSeeds.size() *
sizeof(float);
2780 cap = mTPCABSeeds.capacity() *
sizeof(float);
2783 for (
const auto&
a : mTPCABSeeds) {
2784 siz +=
a.sizeInternal();
2785 cap +=
a.capInternal();
2786 cnt +=
a.getNLinks();
2787 cntCap +=
a.getNLinks();
2791 LOGP(info,
"Size RSS, mTPCABSeeds : size {:9} cap {:9} | internals size:{}/capacity:{} for {} elements", siz, cap, cnt, cntCap, mTPCABSeeds.size());
2793 siz = mTPCABIndexCache.size() *
sizeof(
int);
2794 cap = mTPCABIndexCache.capacity() *
sizeof(
int);
2797 LOGP(info,
"Size RSS, mTPCABIndexCache : size {:9} cap {:9}", siz, cap);
2799 siz = mABWinnersIDs.size() *
sizeof(
int);
2800 cap = mABWinnersIDs.capacity() *
sizeof(
int);
2803 LOGP(info,
"Size RSS, mABWinnersIDs : size {:9} cap {:9}", siz, cap);
2805 siz = mABClusterLinkIndex.size() *
sizeof(
int);
2806 cap = mABClusterLinkIndex.capacity() *
sizeof(
int);
2809 LOGP(info,
"Size RSS, mABClusterLinkIndex : size {:9} cap {:9}", siz, cap);
2812 siz += mTPCSectIndexCache[is].size() *
sizeof(
int);
2813 cap += mTPCSectIndexCache[is].capacity() *
sizeof(
int);
2817 LOGP(info,
"Size RSS, mTPCSectIndexCache : size {:9} cap {:9}", siz, cap);
2820 siz += mITSSectIndexCache[is].size() *
sizeof(
int);
2821 cap += mITSSectIndexCache[is].capacity() *
sizeof(
int);
2825 LOGP(info,
"Size RSS, mITSSectIndexCache : size {:9} cap {:9}", siz, cap);
2828 siz += mTPCTimeStart[is].size() *
sizeof(
int);
2829 cap += mTPCTimeStart[is].capacity() *
sizeof(
int);
2833 LOGP(info,
"Size RSS, mTPCTimeStart : size {:9} cap {:9}", siz, cap);
2836 siz += mITSTimeStart[is].size() *
sizeof(
int);
2837 cap += mITSTimeStart[is].capacity() *
sizeof(
int);
2841 LOGP(info,
"Size RSS, mITSTimeStart : size {:9} cap {:9}", siz, cap);
2843 siz = mITSTrackROFContMapping.size() *
sizeof(
int);
2844 cap = mITSTrackROFContMapping.capacity() *
sizeof(
int);
2847 LOGP(info,
"Size RSS, ITSTrackROFContMapping: size {:9} cap {:9}", siz, cap);
2849 LOGP(info,
"TotalSizes/Capacities: SHM: {}/{} Heap: {}/{}", sizTotShm, capTotShm, sizTot, capTot);
2856 mNThreads =
n > 0 ?
n : 1;
2858 LOG(warning) <<
"Multithreading is not supported, imposing single thread";
2867#ifdef _ALLOW_DEBUG_TREES_
2883 mTimer[SWDBG].Start(
false);
2884 const auto& tpcOrig = mTPCTracksArray[tpcIndex];
2885 uint8_t clSect = 0, clRow = 0, prevRow = 0xff, padFromEdge = -1;
2888 std::array<bool, 152> shMap{};
2889 bool prevRawShared =
false;
2890 for (
int i = 0;
i < tpcOrig.getNClusterReferences();
i++) {
2891 tpcOrig.getClusterReference(mTPCTrackClusIdx,
i, clSect, clRow, clIdx);
2892 unsigned int absoluteIndex = mTPCClusterIdxStruct->
clusterOffset[clSect][clRow] + clIdx;
2894 if (!(prevRow == clRow && prevRawShared)) {
2898 prevRawShared =
true;
2901 const auto& clus = mTPCClusterIdxStruct->
clusters[clSect][clRow][clIdx];
2902 padFromEdge =
uint8_t(clus.getPad());
2903 if (padFromEdge > TPCGeometry.NPads(clRow) / 2) {
2904 padFromEdge = TPCGeometry.NPads(clRow) - 1 - padFromEdge;
2906 int tb = tpcOrig.getTime0() * mNTPCOccBinLengthInv;
2907 float mltTPC = tb < 0 ? mTBinClOcc[0] : (tb >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[tb]);
2908 (*mDBGOut) <<
"tpcOrig"
2909 <<
"tf=" << mTFCount
2910 <<
"index=" << tpcIndex
2912 <<
"chi2TPC=" << tpcOrig.getChi2()
2913 <<
"nClus=" << tpcOrig.getNClusters()
2914 <<
"nShared=" << nshared
2915 <<
"time0=" << tpcOrig.getTime0()
2917 <<
"minRow=" << clRow
2918 <<
"padFromEdge=" << padFromEdge
2919 <<
"multTPC=" << mltTPC;
2921 (*mDBGOut) <<
"tpcOrig"
2922 <<
"tpcLbl=" << mTPCTrkLabels[tpcIndex];
2924 (*mDBGOut) <<
"tpcOrig"
2926 mTimer[SWDBG].Stop();
2934 mTimer[SWDBG].Start(
false);
2936 auto& trackITS = mITSWork[itsID];
2937 auto& trackTPC = mTPCWork[tpcID];
2939 chi2 = getPredictedChi2NoZ(trackITS, trackTPC);
2941 (*mDBGOut) <<
"match"
2942 <<
"tf=" << mTFCount <<
"chi2Match=" << chi2 <<
"its=" << trackITS <<
"tpc=" << trackTPC <<
"tcorr=" << tCorr;
2944 (*mDBGOut) <<
"match"
2945 <<
"itsLbl=" << mITSLblWork[itsID] <<
"tpcLbl=" << mTPCLblWork[tpcID];
2947 int tb = mTPCTracksArray[trackTPC.sourceID].getTime0() * mNTPCOccBinLengthInv;
2948 float mltTPC = tb < 0 ? mTBinClOcc[0] : (tb >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[tb]);
2949 (*mDBGOut) <<
"match"
2950 <<
"rejFlag=" << rejFlag
2951 <<
"multTPC=" << mltTPC
2952 <<
"multITSTr=" << mITSTrackROFRec[trackITS.roFrame].getNEntries()
2953 <<
"multITSCl=" << mITSClusterROFRec[trackITS.roFrame].getNEntries()
2956 mTimer[SWDBG].Stop();
2964 mTimer[SWDBG].Start(
false);
2966 LOG(info) <<
"Dumping debug tree for winner matches";
2967 for (
int iits = 0; iits <
int(mITSWork.size()); iits++) {
2968 auto& tITS = mITSWork[iits];
2969 if (isDisabledITS(tITS)) {
2972 auto& itsMatchRec = mMatchRecordsITS[tITS.
matchID];
2973 int itpc = itsMatchRec.partnerID;
2974 auto& tTPC = mTPCWork[itpc];
2976 (*mDBGOut) <<
"matchWin"
2977 <<
"tf=" << mTFCount <<
"chi2Match=" << itsMatchRec.chi2 <<
"chi2Refit=" << mWinnerChi2Refit[iits] <<
"its=" << tITS <<
"tpc=" << tTPC;
2980 (*mDBGOut) <<
"matchWin"
2981 <<
"itsLbl=" << mITSLblWork[iits] <<
"tpcLbl=" << mTPCLblWork[itpc];
2983 int tb = mTPCTracksArray[tTPC.sourceID].getTime0() * mNTPCOccBinLengthInv;
2984 float mltTPC = tb < 0 ? mTBinClOcc[0] : (tb >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[tb]);
2985 (*mDBGOut) <<
"matchWin"
2986 <<
"multTPC=" << mltTPC
2987 <<
"multITSTr=" << mITSTrackROFRec[tITS.
roFrame].getNEntries()
2988 <<
"multITSCl=" << mITSClusterROFRec[tITS.
roFrame].getNEntries()
2991 mTimer[SWDBG].Stop();
General auxilliary methods.
Definition of the GeometryManager class.
Helper for geometry and GRP related CCDB requests.
Some ALICE geometry constants of common interest.
constexpr int p1()
constexpr to accelerate the coordinates changing
Definition of the GeometryTGeo class.
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 TPC ITS matching.
Definition of the Names Generator class.
Definition of the parameter class for the detector.
Definition of the parameter class for the detector electronics.
Header to collect physics constants.
Wrapper container for different reconstructed object types.
Configurable params for tracks ad hoc tuning.
Helper class to obtain TPC clusters / digits / labels from DPL.
int getLastFilledBC(int dir=-1) const
bool testBC(int bcID, int dir=-1) const
int getFirstFilledBC(int dir=-1) const
void setFakeFlag(bool v=true)
static int getNHBFPerTF()
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
void printKeyValues(bool showProv=true, bool useLogger=false, bool withPadding=true, bool showHash=true) const final
static const MatchTPCITSParams & Instance()
o2::math_utils::Bracket< float > BracketF
void setUseMatCorrFlag(MatCorrType f)
void fillCalibDebug(int ifit, int iTPC, const o2::dataformats::TrackTPCITS &match, pmr::vector< o2::dataformats::Triplet< float, float, float > > &calib)
o2::BaseCluster< float > ITSCluster
void reportTiming()
clear results of previous event reco
void setITSROFrameLengthMUS(float fums)
set ITS ROFrame duration in BC (continuous mode only)
@ MatchTreeAccOnly
fill the matching candidates tree only once the cut is passed
@ TPCOrigTree
original TPC tracks with some aux info
@ MatchTreeAll
produce matching candidates tree for all candidates
@ WinnerMatchesTree
separate debug tree for winner matches
void clear()
set Bunch filling and init helpers for validation by BCs
void refitWinners(pmr::vector< o2::dataformats::TrackTPCITS > &matchedTracks, pmr::vector< o2::MCCompLabel > &matchLabels, pmr::vector< o2::dataformats::Triplet< float, float, float > > &calib)
void refitABWinners(pmr::vector< o2::dataformats::TrackTPCITS > &matchedTracks, pmr::vector< o2::MCCompLabel > &matchLabels, pmr::vector< o2::MCCompLabel > &ABTrackletLabels, pmr::vector< int > &ABTrackletClusterIDs, pmr::vector< o2::itsmft::TrkClusRef > &ABTrackletRefs, pmr::vector< o2::dataformats::Triplet< float, float, float > > &calib)
bool refitABTrack(int iITSAB, const TPCABSeed &seed, pmr::vector< o2::dataformats::TrackTPCITS > &matchedTracks, pmr::vector< int > &ABTrackletClusterIDs, pmr::vector< o2::itsmft::TrkClusRef > &ABTrackletRefs)
void setTPCVDrift(const o2::tpc::VDriftCorrFact &v)
o2::dataformats::RangeReference< int, int > ClusRange
bool isDebugFlag(UInt_t flags) const
get debug trees flags
void printCandidatesITS() const
bool runAfterBurner(pmr::vector< o2::dataformats::TrackTPCITS > &matchedTracks, pmr::vector< o2::MCCompLabel > &matchLabels, pmr::vector< o2::MCCompLabel > &ABTrackletLabels, pmr::vector< int > &ABTrackletClusterIDs, pmr::vector< o2::itsmft::TrkClusRef > &ABTrackletRefs, pmr::vector< o2::dataformats::Triplet< float, float, float > > &calib)
void setTPCCorrMaps(o2::gpu::CorrectionMapsHelper *maph)
print settings
void reportSizes(pmr::vector< o2::dataformats::TrackTPCITS > &matchedTracks, pmr::vector< o2::itsmft::TrkClusRef > &ABTrackletRefs, pmr::vector< int > &ABTrackletClusterIDs, pmr::vector< o2::MCCompLabel > &matchLabels, pmr::vector< o2::MCCompLabel > &ABTrackletLabels, pmr::vector< o2::dataformats::Triplet< float, float, float > > &calib)
void printCandidatesTPC() const
void fillTPCITSmatchTree(int itsID, int tpcID, int rejFlag, float chi2=-1., float tCorr=0.)
bool refitTrackTPCITS(int slot, int iTPC, int &iITS, pmr::vector< o2::dataformats::TrackTPCITS > &matchedTracks, pmr::vector< o2::MCCompLabel > &matchLabels, pmr::vector< o2::dataformats::Triplet< float, float, float > > &calib)
static constexpr int MaxUpDnLadders
void run(const o2::globaltracking::RecoContainer &inp, pmr::vector< o2::dataformats::TrackTPCITS > &matchedTracks, pmr::vector< o2::itsmft::TrkClusRef > &ABTrackletRefs, pmr::vector< int > &ABTrackletClusterIDs, pmr::vector< o2::MCCompLabel > &matchLabels, pmr::vector< o2::MCCompLabel > &ABTrackletLabels, pmr::vector< o2::dataformats::Triplet< float, float, float > > &calib)
void setBunchFilling(const o2::BunchFilling &bf)
void dumpTPCOrig(bool acc, int tpcIndex)
void setDebugFlag(UInt_t flag, bool on=true)
set the name of output debug file
void setITSROFrameLengthInBC(int nbc)
void setITSTimeBiasInBC(int n)
static constexpr int NITSLayers
perform matching for provided input
UInt_t getDebugFlags() const
set or unset debug stream flag
void setNThreads(int n)
perform all initializations
static GeometryTGeo * Instance()
void fillMatrixCache(int mask) override
int getFirstClusterEntry() const
void acquirePattern(iterator &pattIt)
int getNPixels() const
Returns the number of fired pixels.
static constexpr unsigned short InvalidPatternID
int getNpixels(int n) const
Returns the number of fired pixels of the n_th element.
bool isGroup(int n) const
Returns true if the element corresponds to a group of rare topologies.
Relation isOutside(const Bracket< T > &t) const
bool match(const std::vector< std::string > &queries, const char *pattern)
GLfloat GLfloat GLfloat alpha
GLboolean GLboolean GLboolean b
GLuint GLsizei const GLchar * label
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t0
GLboolean GLboolean GLboolean GLboolean a
GLdouble GLdouble GLdouble z
uint8_t itsSharedClusterMap uint8_t
constexpr float XTPCOuterRef
reference radius to propagate outer TPC track
constexpr float XTPCInnerRef
reference radius at which TPC provides the tracks
constexpr double LHCBunchSpacingMUS
constexpr int LHCMaxBunches
constexpr double LHCOrbitMUS
constexpr double LHCBunchSpacingNS
constexpr int Validated
flags to tell the status of TPC-ITS tracks comparison
void convertCompactClusters(gsl::span< const itsmft::CompClusterExt > clusters, gsl::span< const unsigned char >::iterator &pattIt, std::vector< o2::BaseCluster< float > > &output, const its3::TopologyDictionary *dict)
convert compact clusters to 3D spacepoints
void convertCompactClusters(gsl::span< const itsmft::CompClusterExt > clusters, gsl::span< const unsigned char >::iterator &pattIt, std::vector< o2::BaseCluster< float > > &output, const itsmft::TopologyDictionary *dict)
convert compact clusters to 3D spacepoints
if(!okForPhiMin(phi0, phi1))
float angle2Alpha(float phi)
int angle2Sector(float phi)
std::tuple< float, float > rotateZ(float xL, float yL, float snAlp, float csAlp)
float sector2Angle(int sect)
std::vector< T, o2::pmr::polymorphic_allocator< T > > vector
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
bool isValid(std::string alias)
std::string asString() const
bool isSelected(const RecPoints &rp) const
float chi2
chi2 after update
uint8_t nDaughters
number of daughter links on lower layers
int8_t nContLayers
number of contributing layers
int nextOnLr
ID of the next (in quality) link on the same layer.
std::array< ClusRange, o2::its::RecoGeomHelper::getNChips()> chipRefs
std::vector< int > clusterID
o2::math_utils::Bracketf_t tBracket
std::vector< std::deque< ABTrackLink > > threadPool
bool runAfterBurner
run afterburner for TPCtrack-ITScluster matching
float nABSigmaZ
nSigma cut on afterburner track-cluster Z distance
float err2ABExtraZ
extra "systematic" error on Z
float tpcTimeICMatchingNSigma
nsigma for matching TPC corrected time and InteractionCandidate from FT0
float XMatchingRef
reference radius to propagate tracks for matching
float crudeNSigma2Cut[o2::track::kNParams]
float tfEdgeTimeToleranceMUS
corrected TPC time allowed to go out from the TF time edges by this amount
float maxVDriftUncertainty
max assumed VDrift relative uncertainty, used only in VDrift calibration mode
float globalTimeBiasMUS
global time shift to apply to assigned time, brute force way to eliminate bias wrt FIT
int requireToReachLayerAB
AB tracks should reach at least this layer from above.
int maxMatchCandidates
max allowed matching candidates per TPC track
float minTPCTrackR
cut on minimal TPC tracks radius to consider for matching, 666*pt_gev*B_kgaus/5
float minBetaGammaForPIDDiff
float safeMarginTPCITSTimeBin
safety margin (in TPC time bins) for ITS-TPC tracks time (in TPC time bins!) comparison
TimeOutliersPolicy ITSTimeOutliersPolicy
int minTPCClusters
minimum number of clusters to consider
ValidateMatchByFIT validateMatchByFIT
when comparing ITS-TPC matches, prefer those which have time of Interaction Candidate
int lowestLayerAB
lowest layer to reach in AfterBurner
float tpcExtConstrainedNSigma
nsigma to apply to externally (TRD,TOF) time-constrained TPC tracks time error
float maxVDritTimeOffset
max possible TDrift offset to calibrate
o2::base::Propagator::MatCorrType matCorr
int minContributingLayersAB
AB tracks must have at least this amount on contributing layers.
float crudeAbsDiffCut[o2::track::kNParams]
float cutMatchingChi2
cut on matching chi2
float safeMarginTimeCorrErr
safety marging (in \mus) for TPC track time corrected by ITS constraint
float cutABTrack2ClChi2
cut on AfterBurner track-cluster chi2
float maxVDriftTrackQ2Pt
use only tracks below this q/pt (with field only)
float nABSigmaY
nSigma cut on afterburner track-cluster Y distance
int verbosity
verbosit level
float minITSTrackR
cut on minimal ITS tracks radius to consider for matching, 666*pt_gev*B_kgaus/5
float err2ABExtraY
extra "systematic" error on Y
float globalTimeExtraErrorMUS
extra error to add to global time estimate
float safeMarginTPCTimeEdge
safety margin in cm when estimating TPC track tMin and tMax from assigned time0 and its track Z posit...
o2::InteractionRecord startIR
void createTracksVariadic(T creator, GTrackID::mask_t srcSel=GTrackID::getSourcesMask("all")) const
GTrackID getTPCContributorGID(GTrackID source) const
auto getFT0RecPoints() const
std::array< int, o2::its::RecoGeomHelper::getNLayers()> firstInLr
entry of 1st (best) hypothesis on each layer
int winLinkID
ID of the validated link.
void addLink(const o2::track::TrackParCov &trc, int clID, int parentID, int nextID, int lr, int nc, int laddID, float chi2)
int ICCanID
interaction candidate ID (they are sorted in increasing time)
ABTrackLink & getLink(int i)
static LinksPoolMT * gLinksPool
pool of links per thread
int8_t lowestLayer
lowest layer reached
int roFrame
ITS readout frame assigned to this track.
o2::math_utils::Bracketf_t tBracket
bracketing time in \mus
float dL
distance integrated during propagation to reference X (as pion)
int matchID
entry (non if MinusOne) of its matchCand struct in the mMatchesITS
float xrho
x*rho seen during propagation to reference X (as pion)
int sourceID
track origin id
float time0
nominal time in \mus since start of TF (time0 for bare TPC tracks, constrained time for TRD/TOF const...
int sourceID
TPC track origin in.
int matchID
entry (non if MinusOne) of its matchTPC struct in the mMatchesTPC
float getCorrectedTime(float dt) const
float timeErr
time sigma (makes sense for constrained tracks only)
o2::dataformats::GlobalTrackID gid
std::array< RecoLayer, o2::itsmft::ChipMappingITS::NLayers > layers
static constexpr float ladderWidthInv()
const ClusterNative * clusters[constants::MAXSECTOR][constants::MAXGLOBALPADROW]
unsigned int clusterOffset[constants::MAXSECTOR][constants::MAXGLOBALPADROW]
float refTimeOffset
additive time offset reference (\mus)
float refVDrift
reference vdrift for which factor was extracted
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
o2::InteractionRecord ir(0, 0)
char const *restrict const cmp