diff --git a/Detectors/TRD/qc/include/TRDQC/Tracking.h b/Detectors/TRD/qc/include/TRDQC/Tracking.h index f39c64286d0cc..c1c44c1a07dce 100644 --- a/Detectors/TRD/qc/include/TRDQC/Tracking.h +++ b/Detectors/TRD/qc/include/TRDQC/Tracking.h @@ -23,6 +23,7 @@ #include "DataFormatsTRD/Constants.h" #include "ReconstructionDataFormats/TrackTPCITS.h" #include "ReconstructionDataFormats/GlobalTrackID.h" +#include "DataFormatsTRD/TrackTriggerRecord.h" #include "DataFormatsTPC/TrackTPC.h" #include "DetectorsBase/Propagator.h" #include "GPUTRDRecoParam.h" @@ -103,6 +104,10 @@ class Tracking mLocalGain = localGain; } + // quantities necessary for pile-up correction + void setTriggeredBCFT0(std::vector t) { mTriggeredBCFT0 = t; } + void setFirstOrbit(uint32_t o) { mFirstOrbit = o; } + private: float mMaxSnp{o2::base::Propagator::MAX_SIN_PHI}; ///< max snp when propagating tracks float mMaxStep{o2::base::Propagator::MAX_STEP}; ///< maximum step for propagation @@ -115,12 +120,20 @@ class Tracking std::vector mTrackQC; // input from DPL - gsl::span mTracksITSTPC; ///< ITS-TPC seeding tracks - gsl::span mTracksTPC; ///< TPC seeding tracks - gsl::span mTracksITSTPCTRD; ///< TRD tracks reconstructed from TPC or ITS-TPC seeds - gsl::span mTracksTPCTRD; ///< TRD tracks reconstructed from TPC or TPC seeds - gsl::span mTrackletsRaw; ///< array of raw tracklets needed for TRD refit - gsl::span mTrackletsCalib; ///< array of calibrated tracklets needed for TRD refit + gsl::span mTracksITSTPC; ///< ITS-TPC seeding tracks + gsl::span mTracksTPC; ///< TPC seeding tracks + gsl::span mTracksITSTPCTRD; ///< TRD tracks reconstructed from TPC or ITS-TPC seeds + gsl::span mTracksTPCTRD; ///< TRD tracks reconstructed from TPC or TPC seeds + gsl::span mTrackTriggerRecordsITSTPCTRD; ///< TRD tracks reconstructed from TPC or ITS-TPC seeds + gsl::span mTrackTriggerRecordsTPCTRD; ///< TRD tracks reconstructed from TPC or TPC seeds + gsl::span mTrackletsRaw; ///< array of raw tracklets needed for TRD refit + gsl::span mTrackletsCalib; ///< array of calibrated tracklets needed for TRD refit + + // quantities necessary for pile-up correction + std::vector mTriggeredBCFT0; ///< array with the FT0 trigger times + int mCurrentTriggerRecord; + uint32_t mFirstOrbit; + int mCurrentTrackId; // corrections from ccdb, some need to be loaded only once hence an init flag o2::trd::LocalGainFactor mLocalGain; ///< local gain factors from krypton calibration diff --git a/Detectors/TRD/qc/src/Tracking.cxx b/Detectors/TRD/qc/src/Tracking.cxx index da2d05794e2d8..514488a3c2d0a 100644 --- a/Detectors/TRD/qc/src/Tracking.cxx +++ b/Detectors/TRD/qc/src/Tracking.cxx @@ -37,15 +37,23 @@ void Tracking::setInput(const o2::globaltracking::RecoContainer& input) mTracksTPCTRD = input.getTPCTRDTracks(); mTrackletsRaw = input.getTRDTracklets(); mTrackletsCalib = input.getTRDCalibratedTracklets(); + mTrackTriggerRecordsITSTPCTRD = input.getITSTPCTRDTriggers(); + mTrackTriggerRecordsTPCTRD = input.getTPCTRDTriggers(); } void Tracking::run() { + mCurrentTriggerRecord = 0; + mCurrentTrackId = 0; for (const auto& trkTrd : mTracksTPCTRD) { checkTrack(trkTrd, true); + mCurrentTrackId++; } + mCurrentTriggerRecord = 0; + mCurrentTrackId = 0; for (const auto& trkTrd : mTracksITSTPCTRD) { checkTrack(trkTrd, false); + mCurrentTrackId++; } } @@ -65,6 +73,59 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) qcStruct.dEdxTotTPC = isTPCTRD ? mTracksTPC[id].getdEdx().dEdxTotTPC : mTracksTPC[mTracksITSTPC[id].getRefTPC()].getdEdx().dEdxTotTPC; } + // find corresponding track trigger record to get track timing + int triggeredBC = 0; + for (; mCurrentTriggerRecord < (isTPCTRD ? mTrackTriggerRecordsTPCTRD.size() : mTrackTriggerRecordsITSTPCTRD.size()); mCurrentTriggerRecord++) { + auto& tRecord = (isTPCTRD ? mTrackTriggerRecordsTPCTRD[mCurrentTriggerRecord] : mTrackTriggerRecordsITSTPCTRD[mCurrentTriggerRecord]); + if (mCurrentTrackId >= tRecord.getFirstTrack() && mCurrentTrackId < tRecord.getFirstTrack() + tRecord.getNumberOfTracks()) { + uint32_t currentOrbit = tRecord.getBCData().orbit; + triggeredBC = tRecord.getBCData().bc + (currentOrbit - mFirstOrbit) * o2::constants::lhc::LHCMaxBunches; + break; + } + } + + // Find most probable BCs and RMS for pile-up correction and error. Same BC is assumed for all tracklets + float tCorrPileUp = 0.; + float tErrPileUp2 = 0; + float maxProb = 0.f; + // The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability + float sumCorr = 0.f; + float sumCorr2 = 0.f; + float sumProb = 0.f; + for (int iBC = 0; iBC < mTriggeredBCFT0.size(); iBC++) { + int deltaBC = roundf(mTriggeredBCFT0[iBC] - triggeredBC); + if (deltaBC <= mRecoParam.getPileUpRangeBefore()) { + continue; + } + if (deltaBC >= mRecoParam.getPileUpRangeAfter()) { + break; + } + // collect the charges + std::array q0; + std::array q1; + for (int iLy = 0; iLy < NLAYER; iLy++) { + int trkltId = trkTrd.getTrackletIndex(iLy); + if (trkltId < 0) { + q0[iLy] = -1; + q1[iLy] = -1; + } else { + q0[iLy] = mTrackletsRaw[trkltId].getQ0(); + q1[iLy] = mTrackletsRaw[trkltId].getQ1(); + } + } + // get pile-up probability + float probBC = mRecoParam.getPileUpProbTrack(deltaBC, q0, q1); + sumCorr += probBC * deltaBC; + sumCorr2 += probBC * deltaBC * deltaBC; + sumProb += probBC; + if (probBC > maxProb) { + maxProb = probBC; + tCorrPileUp = -deltaBC; + } + } + if (sumProb > 1e-6) + tErrPileUp2 = sumCorr2 / sumProb - 2 * tCorrPileUp * sumCorr / sumProb + tCorrPileUp * tCorrPileUp; + for (int iLayer = 0; iLayer < NLAYER; ++iLayer) { int trkltId = trkTrd.getTrackletIndex(iLayer); if (trkltId < 0) { @@ -93,9 +154,16 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) if (!((trk.getSigmaZ2() < (padLength * padLength / 12.f)) && (std::fabs(mTrackletsCalib[trkltId].getZ() - trk.getZ()) < padLength))) { tiltCorrUp = 0.f; } - std::array trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp, zPosCorrUp}; + + // conversion from slope in pad per time bin to slope in cm per BC = tracklets[trkltIdx].getSlopeFloat() * padWidth / BCperTimeBin + float slopeFactor = mTrackletsRaw[trkltId].getSlopeFloat() * pad->getWidthIPad() / 4.f; + float yCorrPileUp = tCorrPileUp * slopeFactor; + float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor; + + std::array trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp, zPosCorrUp}; std::array trkltCovUp; mRecoParam.recalcTrkltCov(tilt, trk.getSnp(), pad->getRowSize(tracklet.getPadRow()), trkltCovUp); + trkltCovUp[0] += yAddErrPileUp2; auto chi2trklt = trk.getPredictedChi2(trkltPosUp, trkltCovUp); qcStruct.trackProp[iLayer] = trk; diff --git a/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h b/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h index b6a0cb9b78f2f..e52b3bcc23842 100644 --- a/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h +++ b/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h @@ -30,6 +30,8 @@ #include "DataFormatsParameters/GRPObject.h" #include "ReconstructionDataFormats/GlobalTrackID.h" #include "DataFormatsGlobalTracking/RecoContainer.h" +#include "DataFormatsFT0/RecPoints.h" +#include "FT0Reconstruction/InteractionTag.h" #include "TRDQC/Tracking.h" #include @@ -47,7 +49,7 @@ namespace trd class TRDGlobalTrackingQC : public Task { public: - TRDGlobalTrackingQC(std::shared_ptr dr, std::shared_ptr gr, bool tpcAvailable) : mDataRequest(dr), mGGCCDBRequest(gr), mTPCavailable(tpcAvailable) {} + TRDGlobalTrackingQC(std::shared_ptr dr, std::shared_ptr gr, bool tpcAvailable, o2::dataformats::GlobalTrackID::mask_t src) : mDataRequest(dr), mGGCCDBRequest(gr), mTPCavailable(tpcAvailable), mTrkMask(src) {} ~TRDGlobalTrackingQC() override = default; void init(InitContext& ic) final { @@ -67,6 +69,23 @@ class TRDGlobalTrackingQC : public Task updateTimeDependentParams(pc); // Make sure this is called after recoData.collectData, which may load some conditions mQC.reset(); mQC.setInput(recoData); + std::vector triggeredBCFT0; + if (mTrkMask[GTrackID::FT0]) { // pile-up tagging was requested + auto ft0recPoints = recoData.getFT0RecPoints(); + uint32_t firstOrbit = 0; + for (size_t ft0id = 0; ft0id < ft0recPoints.size(); ft0id++) { + const auto& f0rec = ft0recPoints[ft0id]; + if (ft0id == 0) { + firstOrbit = f0rec.getInteractionRecord().orbit; + mQC.setFirstOrbit(firstOrbit); + } + if (o2::ft0::InteractionTag::Instance().isSelected(f0rec)) { + uint32_t currentOrbit = f0rec.getInteractionRecord().orbit; + triggeredBCFT0.push_back(f0rec.getInteractionRecord().bc + (currentOrbit - firstOrbit) * o2::constants::lhc::LHCMaxBunches); + } + } + } + mQC.setTriggeredBCFT0(triggeredBCFT0); mQC.run(); pc.outputs().snapshot(Output{"TRD", "TRACKINGQC", 0}, mQC.getTrackQC()); } @@ -94,6 +113,7 @@ class TRDGlobalTrackingQC : public Task } } + o2::dataformats::GlobalTrackID::mask_t mTrkMask; ///< seeding track sources (TPC, ITS-TPC) std::shared_ptr mDataRequest; std::shared_ptr mGGCCDBRequest; bool mTPCavailable{false}; @@ -133,7 +153,7 @@ DataProcessorSpec getTRDGlobalTrackingQCSpec(o2::dataformats::GlobalTrackID::mas "trd-tracking-qc", dataRequest->inputs, outputs, - AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, isTPCavailable)}, + AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, isTPCavailable, src)}, Options{}}; } diff --git a/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingSpec.h b/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingSpec.h index 93f07dd58445e..bfe82a05cc0cb 100644 --- a/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingSpec.h +++ b/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingSpec.h @@ -109,6 +109,7 @@ class TRDGlobalTracking : public o2::framework::Task #endif std::array mCovDiagInner{}; ///< total cov.matrix extra diagonal error from TrackTuneParams std::array mCovDiagOuter{}; ///< total cov.matrix extra diagonal error from TrackTuneParams + std::vector mTriggeredBCFT0; ///< array with the FT0 trigger times // PID PIDPolicy mPolicy{PIDPolicy::DEFAULT}; ///< Model to load an evaluate std::unique_ptr mBase; ///< PID engine diff --git a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx index 0f578efd3aa5b..7404b71c59ad5 100644 --- a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx @@ -415,6 +415,24 @@ void TRDGlobalTracking::run(ProcessingContext& pc) } LOGF(info, "%i tracks are loaded into the TRD tracker. Out of those %i ITS-TPC tracks and %i TPC tracks", nTracksLoadedITSTPC + nTracksLoadedTPC, nTracksLoadedITSTPC, nTracksLoadedTPC); + // Load the FT0 triggered BCs if this is requested + + if (mTrkMask[GTrackID::FT0]) { // pile-up tagging was requested + auto ft0recPoints = inputTracks.getFT0RecPoints(); + uint32_t firstOrbit = 0; + for (size_t ft0id = 0; ft0id < ft0recPoints.size(); ft0id++) { + const auto& f0rec = ft0recPoints[ft0id]; + if (ft0id == 0) + firstOrbit = f0rec.getInteractionRecord().orbit; + if (o2::ft0::InteractionTag::Instance().isSelected(f0rec)) { + uint32_t currentOrbit = f0rec.getInteractionRecord().orbit; + mTriggeredBCFT0.push_back(f0rec.getInteractionRecord().bc + (currentOrbit - firstOrbit) * o2::constants::lhc::LHCMaxBunches); + } + } + } + + mTracker->SetFT0TriggeredBC(mTriggeredBCFT0.data(), mTriggeredBCFT0.size()); + // start the tracking // mTracker->DumpTracks(); mChainTracking->DoTRDGPUTracking(mTracker); @@ -797,6 +815,45 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, } } + // Find most probable BCs and RMS for pile-up correction and error. Same BC is assumed for all tracklets + float tCorrPileUp = 0.; + float tErrPileUp2 = 0; + float maxProb = 0.f; + // The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability + float sumCorr = 0.f; + float sumCorr2 = 0.f; + float sumProb = 0.f; + for (int iBC = 0; iBC < mTriggeredBCFT0.size(); iBC++) { + int deltaBC = roundf(mTriggeredBCFT0[iBC] - mChainTracking->mIOPtrs.trdTriggerTimes[trk.getCollisionId()] / o2::constants::lhc::LHCBunchSpacingMUS); + if (deltaBC <= mRecoParam.getPileUpRangeBefore() || deltaBC >= mRecoParam.getPileUpRangeAfter()) { + continue; + } + // collect the charges + std::array q0; + std::array q1; + for (int iLy = 0; iLy < NLAYER; iLy++) { + int trkltId = trk.getTrackletIndex(iLy); + if (trkltId < 0) { + q0[iLy] = -1; + q1[iLy] = -1; + } else { + q0[iLy] = mTrackletsRaw[trkltId].getQ0(); + q1[iLy] = mTrackletsRaw[trkltId].getQ1(); + } + } + // get pile-up probability + float probBC = mRecoParam.getPileUpProbTrack(deltaBC, q0, q1); + sumCorr += probBC * deltaBC; + sumCorr2 += probBC * deltaBC * deltaBC; + sumProb += probBC; + if (probBC > maxProb) { + maxProb = probBC; + tCorrPileUp = -deltaBC; + } + } + if (sumProb > 1e-6) + tErrPileUp2 = sumCorr2 / sumProb - 2 * tCorrPileUp * sumCorr / sumProb + tCorrPileUp * tCorrPileUp; + if (inwards) { // reset covariance to something big for inwards refit trkParam->resetCovariance(100); @@ -827,9 +884,15 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, tiltCorrUp = 0.f; } - std::array trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp, zPosCorrUp}; + // conversion from slope in pad per time bin to slope in cm per BC = tracklets[trkltIdx].getSlopeFloat() * padWidth / BCperTimeBin + float slopeFactor = mTrackletsRaw[trkltId].getSlopeFloat() * pad->getWidthIPad() / 4.f; + float yCorrPileUp = tCorrPileUp * slopeFactor; + float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor; + + std::array trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp, zPosCorrUp}; std::array trkltCovUp; mRecoParam.recalcTrkltCov(tilt, trkParam->getSnp(), pad->getRowSize(mTrackletsRaw[trkltId].getPadRow()), trkltCovUp); + trkltCovUp[0] += yAddErrPileUp2; chi2 += trkParam->getPredictedChi2(trkltPosUp, trkltCovUp); if (!trkParam->update(trkltPosUp, trkltCovUp)) { diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx index f7adc2401df79..11af412f7782b 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx @@ -77,3 +77,71 @@ void GPUTRDRecoParam::recalcTrkltCov(const float tilt, const float snp, const fl cov[1] = c2 * tilt * (sz2 - sy2); cov[2] = c2 * (t2 * sy2 + sz2); } + +float GPUTRDRecoParam::getPileUpProbTracklet(int nBC, int Q0, int Q1) const +{ + // get the probability that the tracklet with charges Q0 and Q1 belongs to a given BC, with a (signed) distance nBC from the TRD-triggered BC + // parametrization depends on whether charges are 0 or not + int maxBC = mPileUpRangeAfter; + int minBC = mPileUpRangeBefore; + int maxProbBC = mPileUpMaxProb; + if (Q0 > 0 && Q1 > 0) { + maxBC = mPileUpRangeAfter11; + minBC = mPileUpRangeBefore11; + maxProbBC = mPileUpMaxProb11; + } + if (Q0 == 0 && Q1 > 0) { + maxBC = mPileUpRangeAfter01; + minBC = mPileUpRangeBefore01; + maxProbBC = mPileUpMaxProb01; + } + if (Q0 > 0 && Q1 == 0) { + maxBC = mPileUpRangeAfter10; + minBC = mPileUpRangeBefore10; + maxProbBC = mPileUpMaxProb10; + } + if (Q0 == 0 && Q1 == 0) { + maxBC = mPileUpRangeAfter00; + minBC = mPileUpRangeBefore00; + maxProbBC = mPileUpMaxProb00; + } + + // prob is 0 if the BC is too far, maximal for a given nBC, and with two linear functions in between. The maximum is chosen so that the integral is 1. + if (nBC <= minBC || nBC >= maxBC) + return 0.; + float maxProb = 2. / (maxBC - minBC); + if (nBC > minBC && nBC <= maxProbBC) + return maxProb / (maxProbBC - minBC) * (nBC - minBC); + else + return maxProb / (maxProbBC - maxBC) * (nBC - maxBC); +} + +float GPUTRDRecoParam::getPileUpProbTrack(int nBC, std::array Q0, std::array Q1) const +{ + // get the probability that the track belongs to a given BC, with a (signed) distance nBC from the TRD-triggered BC + // it depends on the individual probabilities for every of its tracklets. + // + // If P(BC|L0,L1,...) is the probability that the track belongs to a given BC, given the information on the tracklet charges in L0,L1, ... + // P(BC|L0,L1,...) proportional to P(BC)*P(L0,L1,...|BC), prop to P(BC)*P(L0|BC)*P(L1|BC)*... since for a given track and BC, charge in different layers are independent + // prop to P(BC) * P(BC|L0)/P(BC) * P(BC|L1)/P(BC) * ... + // + // P(BC) is the probability with no charge information: we start from this probability, and each tracklet adds new information on pileup probability + + // basic probability, if we had no info on the charges + float probNoInfo = GPUTRDRecoParam::getPileUpProbTracklet(nBC, -1, -1); + + float probTrack = probNoInfo; + if (probNoInfo < 1e-6f) + return 0.; + + // For each tracklet, we add the info on its charge + for (int i = 0; i < 6; i++) { + // negative charge values if the tracklet is not present + if (Q0[i] < 0 || Q1[i] < 0) + continue; + float probTracklet = GPUTRDRecoParam::getPileUpProbTracklet(nBC, Q0[i], Q1[i]); + probTrack *= probTracklet / probNoInfo; + } + + return probTrack; +} diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h index a0a8e71143d94..851159cfcf8ab 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h @@ -56,10 +56,16 @@ class GPUTRDRecoParam GPUd() float getDyRes(float snp) const { return mDyA2 + mDyC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle); } // a^2 + c^2 * (snp - b)^2 GPUd() float convertAngleToDy(float snp) const { return 3.f * snp / CAMath::Sqrt(1 - snp * snp); } // when calibrated, sin(phi) = (dy / xDrift) / sqrt(1+(dy/xDrift)^2) works well GPUd() float getCorrYDy(float snp) const { return mCorrYDyA + mCorrYDyC * (snp - mLorentzAngle) * (snp - mLorentzAngle); } // a + c * (snp - b)^2 + GPUd() float getPileUpProbTracklet(int nBC, int Q0, int Q1) const; + GPUd() float getPileUpProbTrack(int nBC, std::array Q0, std::array Q1) const; /// Get tracklet z correction coefficient for track-eta based corraction GPUd() float getZCorrCoeffNRC() const { return mZCorrCoefNRC; } + /// Get BC intervals for pile-up + GPUd() int getPileUpRangeBefore() const { return mPileUpRangeBefore; } + GPUd() int getPileUpRangeAfter() const { return mPileUpRangeAfter; } + private: // tracklet error parameterization depends on the magnetic field float mLorentzAngle{0.f}; @@ -75,7 +81,29 @@ class GPUTRDRecoParam float mZCorrCoefNRC{1.4f}; ///< tracklet z-position depends linearly on track dip angle - ClassDefNV(GPUTRDRecoParam, 3); + // pile-up prob parametrization, depending on charges + // default parametrization, all tracklets + int mPileUpRangeBefore{-130}; ///< maximal number of BC for which pile-up from previous collision has an influence + int mPileUpMaxProb{0}; ///< number of BC with respect to triggered BC for the event with maximal probability + int mPileUpRangeAfter{70}; ///< maximal number of BC for which pile-up from next collision has an influence + // tracklets with Q0!=0 and Q1!=0 + int mPileUpRangeBefore11{-130}; ///< maximal number of BC for which pile-up from previous collision has an influence + int mPileUpMaxProb11{0}; ///< number of BC with respect to triggered BC for the event with maximal probability + int mPileUpRangeAfter11{30}; ///< maximal number of BC for which pile-up from next collision has an influence + // tracklets with Q0=0 and Q1!=0 + int mPileUpRangeBefore01{-80}; ///< maximal number of BC for which pile-up from previous collision has an influence + int mPileUpMaxProb01{30}; ///< number of BC with respect to triggered BC for the event with maximal probability + int mPileUpRangeAfter01{70}; ///< maximal number of BC for which pile-up from next collision has an influence + // tracklets with Q0!=0 and Q1=0 + int mPileUpRangeBefore10{-130}; ///< maximal number of BC for which pile-up from previous collision has an influence + int mPileUpMaxProb10{-30}; ///< number of BC with respect to triggered BC for the event with maximal probability + int mPileUpRangeAfter10{30}; ///< maximal number of BC for which pile-up from next collision has an influence + // tracklets with Q0=0 and Q1=0 + int mPileUpRangeBefore00{-10}; ///< maximal number of BC for which pile-up from previous collision has an influence + int mPileUpMaxProb00{22}; ///< number of BC with respect to triggered BC for the event with maximal probability + int mPileUpRangeAfter00{40}; ///< maximal number of BC for which pile-up from next collision has an influence + + ClassDefNV(GPUTRDRecoParam, 4); }; } // namespace gpu diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index 80098ff151ebe..2f338800bafb7 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -93,7 +93,7 @@ void* GPUTRDTracker_t::SetPointersTracks(void* base) } template -GPUTRDTracker_t::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mGenerateSpacePoints(false), mProcessPerTimeFrame(false), mNAngleHistogramBins(25), mAngleHistogramRange(50), mMemoryPermanent(-1), mMemoryTracklets(-1), mMemoryTracks(-1), mNMaxCollisions(0), mNMaxTracks(0), mNMaxSpacePoints(0), mTracks(nullptr), mTrackAttribs(nullptr), mNCandidates(1), mNTracks(0), mNEvents(0), mMaxBackendThreads(100), mTrackletIndexArray(nullptr), mHypothesis(nullptr), mCandidates(nullptr), mSpacePoints(nullptr), mGeo(nullptr), mRecoParam(nullptr), mDebugOutput(false), mMaxEta(0.84f), mRoadZ(18.f), mTPCVdrift(2.58f), mTPCTDriftOffset(0.f), mDebug(new GPUTRDTrackerDebug()) +GPUTRDTracker_t::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mGenerateSpacePoints(false), mProcessPerTimeFrame(false), mNAngleHistogramBins(25), mAngleHistogramRange(50), mMemoryPermanent(-1), mMemoryTracklets(-1), mMemoryTracks(-1), mNMaxCollisions(0), mNMaxTracks(0), mNMaxSpacePoints(0), mTracks(nullptr), mTrackAttribs(nullptr), mNCandidates(1), mNTracks(0), mNEvents(0), mMaxBackendThreads(100), mTrackletIndexArray(nullptr), mFT0TriggeredBC(nullptr), mNFT0BC(0), mHypothesis(nullptr), mCandidates(nullptr), mSpacePoints(nullptr), mGeo(nullptr), mRecoParam(nullptr), mDebugOutput(false), mMaxEta(0.84f), mRoadZ(18.f), mTPCVdrift(2.58f), mTPCTDriftOffset(0.f), mDebug(new GPUTRDTrackerDebug()) { //-------------------------------------------------------------------- // Default constructor @@ -351,6 +351,7 @@ GPUd() void GPUTRDTracker_t::DoTrackingThread(int32_t iTrk, int32_ } PROP prop(getPropagatorParam()); mTracks[iTrk].setChi2(Param().rec.trd.penaltyChi2); // TODO check if this should not be higher + auto trkStart = mTracks[iTrk]; for (int32_t iColl = 0; iColl < nCollisionIds; ++iColl) { // do track following for each collision candidate and keep best track @@ -436,6 +437,29 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK } mDebug->Reset(); t->setChi2(0.f); + + // Find compatible BC ids + int32_t nIdxBCMin = -1; + int32_t nIdxBCMax = -1; + + for (int32_t iBC = 0; iBC < mNFT0BC; iBC++) { + int32_t deltaBC = CAMath::Round(mFT0TriggeredBC[iBC] - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId] / o2::constants::lhc::LHCBunchSpacingMUS); + if (nIdxBCMin == -1 && deltaBC > mRecoParam->getPileUpRangeBefore()) { + nIdxBCMin = iBC; + } + if (deltaBC >= mRecoParam->getPileUpRangeAfter()) { + nIdxBCMax = iBC; + break; + } + if (iBC == mNFT0BC - 1) { + nIdxBCMax = iBC + 1; + if (nIdxBCMin == -1) { + // we did not find the correct BC, so we don't do any pile-up correction + nIdxBCMin = nIdxBCMax; + } + } + } + float zShiftTrk = 0.f; if (mProcessPerTimeFrame) { zShiftTrk = (mTrackAttribs[iTrk].mTime - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId]) * mTPCVdrift * mTrackAttribs[iTrk].mSide; @@ -585,9 +609,38 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK tiltCorr = 0.f; // will be zero also for TPC tracks which are shifted in z dyTiltCorr = 0.f; } + + // Correction for pile-up: if the track comes from a pile-up event, it should not be extrapolated to the anode plane at t=0, + // contrary to the assumption when tracklet is reconstructed, so we cancel this extrapolation. + // The correction is extracted from the most probable trigger. There is also an additional error depending on all the compatible triggers + float yCorrPileUp = 0.f; + float yAddErrPileUp2 = 0.f; + if (nIdxBCMax - nIdxBCMin >= 2) { + float maxProb = 0.f; + // The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability + float sumCorr = 0.f; + float sumCorr2 = 0.f; + float sumProb = 0.f; + // conversion from slope in pad per time bin to slope in cm per BC = tracklets[trkltIdx].getSlopeFloat() * padWidth / BCperTimeBin + float slopeFactor = tracklets[trkltIdx].GetSlopeFloat() * mGeo->GetPadPlaneWidthIPad(tracklets[trkltIdx].GetDetector()) / 4.f; + for (int32_t iBC = nIdxBCMin; iBC < nIdxBCMax; iBC++) { + int32_t deltaBC = CAMath::Round(mFT0TriggeredBC[iBC] - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId] / o2::constants::lhc::LHCBunchSpacingMUS); + float probBC = mRecoParam->getPileUpProbTracklet(deltaBC, tracklets[trkltIdx].GetQ0(), tracklets[trkltIdx].GetQ1()); + sumCorr += probBC * slopeFactor * deltaBC; + sumCorr2 += probBC * slopeFactor * deltaBC * slopeFactor * deltaBC; + sumProb += probBC; + if (probBC > maxProb) { + maxProb = probBC; + yCorrPileUp = -slopeFactor * deltaBC; + } + } + if (sumProb > 1e-6f) + yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp; + } + // correction for mean z position of tracklet (is not the center of the pad if track eta != 0) float zPosCorr = spacePoints[trkltIdx].getZ() + mRecoParam->getZCorrCoeffNRC() * trkWork->getTgl(); - float yPosCorr = spacePoints[trkltIdx].getY() - tiltCorr; + float yPosCorr = spacePoints[trkltIdx].getY() - tiltCorr + yCorrPileUp; zPosCorr -= zShiftTrk; // shift tracklet instead of track in order to avoid having to do a re-fit for each collision float deltaY = yPosCorr - projY; float deltaZ = zPosCorr - projZ; @@ -596,6 +649,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK if ((CAMath::Abs(deltaY) < roadY) && (CAMath::Abs(deltaZ) < roadZ)) { // TODO: check if this is still necessary after the cut before propagation of track // tracklet is in windwow: get predicted chi2 for update and store tracklet index if best guess RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[trkltIdx].GetZbin()), trkltCovTmp); + trkltCovTmp[0] += yAddErrPileUp2; float chi2 = prop->getPredictedChi2(trkltPosTmpYZ, trkltCovTmp); if (Param().rec.trd.addDeflectionInChi2 && (trkWork->getSnp() < 1.f - 1e-6f) && (trkWork->getSnp() > -1.f + 1e-6f)) { // we add the slope in the chi2 calculation @@ -696,9 +750,39 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK if (!((trkWork->getSigmaZ2() < (padLength * padLength / 12.f)) && (CAMath::Abs(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() - trkWork->getZ()) < padLength))) { tiltCorrUp = 0.f; } - float trkltPosUp[2] = {spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getY() - tiltCorrUp, zPosCorrUp}; + + // Correction for pile-up: if the track comes from a pile-up event, it should not be extrapolated to the anode plane at t=0, + // contrary to the assumption when tracklet is reconstructed, so we cancel this extrapolation. + // The correction is extracted from the most probable trigger. There is also an additional error depending on all the compatible triggers + float yCorrPileUp = 0.f; + float yAddErrPileUp2 = 0.f; + if (nIdxBCMax - nIdxBCMin >= 2) { + float maxProb = 0.f; + // The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability + float sumCorr = 0.f; + float sumCorr2 = 0.f; + float sumProb = 0.f; + // conversion from slope in pad per time bin to slope in cm per BC = tracklets[trkltIdx].getSlopeFloat() * padWidth / BCperTimeBin + float slopeFactor = tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetSlopeFloat() * mGeo->GetPadPlaneWidthIPad(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector()) / 4.f; + for (int32_t iBC = nIdxBCMin; iBC < nIdxBCMax; iBC++) { + int32_t deltaBC = CAMath::Round(mFT0TriggeredBC[iBC] - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId] / o2::constants::lhc::LHCBunchSpacingMUS); + float probBC = mRecoParam->getPileUpProbTracklet(deltaBC, tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetQ0(), tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetQ1()); + sumCorr += probBC * slopeFactor * deltaBC; + sumCorr2 += probBC * slopeFactor * deltaBC * slopeFactor * deltaBC; + sumProb += probBC; + if (probBC > maxProb) { + maxProb = probBC; + yCorrPileUp = -slopeFactor * deltaBC; + } + } + if (sumProb > 1e-6f) + yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp; + } + + float trkltPosUp[2] = {spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getY() - tiltCorrUp + yCorrPileUp, zPosCorrUp}; float trkltCovUp[3] = {0.f}; RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin()), trkltCovUp); + trkltCovUp[0] += yAddErrPileUp2; #ifdef ENABLE_GPUTRDDEBUG prop->setTrack(&trackNoUp); @@ -777,6 +861,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK if (iUpdate == 0 && mNCandidates > 1) { *t = mCandidates[2 * iUpdate + nextIdx]; } + } // end update loop if (!isOK) { diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h index f698e570d2158..d491539b8fcf8 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h @@ -132,6 +132,11 @@ class GPUTRDTracker_t : public GPUProcessor GPUd() void SetRoadZ(float roadZ) { mRoadZ = roadZ; } GPUd() void SetTPCVdrift(float vDrift) { mTPCVdrift = vDrift; } GPUd() void SetTPCTDriftOffset(float t) { mTPCTDriftOffset = t; } + GPUd() void SetFT0TriggeredBC(int32_t* t, int32_t n) + { + mFT0TriggeredBC = t; + mNFT0BC = n; + } GPUd() bool GetIsDebugOutputOn() const { return mDebugOutput; } GPUd() float GetMaxEta() const { return mMaxEta; } @@ -170,6 +175,8 @@ class GPUTRDTracker_t : public GPUProcessor // the array has (kNChambers + 1) * numberOfCollisions entries // note, that for collision iColl one has to add an offset corresponding to the index of the first tracklet of iColl to the index stored in mTrackletIndexArray int32_t* mTrackletIndexArray; + int32_t* mFT0TriggeredBC; // arrays with the FT0 triggered BCs, in number of BCs since the beginning of the TF + int32_t mNFT0BC; // number of FT0 BCs Hypothesis* mHypothesis; // array with multiple track hypothesis TRDTRK* mCandidates; // array of tracks for multiple hypothesis tracking GPUTRDSpacePoint* mSpacePoints; // array with tracklet coordinates in global tracking frame diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTrackletWord.h b/GPU/GPUTracking/TRDTracking/GPUTRDTrackletWord.h index cd7dfb9432b93..c69621c1684c1 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTrackletWord.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTrackletWord.h @@ -97,6 +97,10 @@ class GPUTRDTrackletWord : private o2::trd::Tracklet64 GPUd() float GetdY() const { return getUncalibratedDy(); } GPUd() int32_t GetDetector() const { return getDetector(); } GPUd() int32_t GetHCId() const { return getHCID(); } + GPUd() float GetSlopeFloat() const { return getSlopeFloat(); } + GPUd() int GetQ0() const { return getQ0(); } + GPUd() int GetQ1() const { return getQ1(); } + GPUd() int GetQ2() const { return getQ2(); } // IMPORTANT: Do not add members, this class must keep the same memory layout as o2::trd::Tracklet64 };