From 8eca6cb72900c3dbd38e8477fc9bcbfdc201edf5 Mon Sep 17 00:00:00 2001 From: Gauthier Legras Date: Wed, 15 Apr 2026 16:30:09 +0200 Subject: [PATCH 1/9] TRD: correction for pile-up --- Detectors/TRD/qc/include/TRDQC/Tracking.h | 25 +++-- Detectors/TRD/qc/src/Tracking.cxx | 71 +++++++++++++- .../TRDWorkflow/TRDGlobalTrackingQCSpec.h | 24 ++++- .../TRDWorkflow/TRDGlobalTrackingSpec.h | 1 + .../workflow/src/TRDGlobalTrackingSpec.cxx | 65 ++++++++++++- GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx | 68 ++++++++++++++ GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h | 30 +++++- GPU/GPUTracking/DataTypes/GPUTRDTrack.h | 2 +- GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx | 92 ++++++++++++++++++- GPU/GPUTracking/TRDTracking/GPUTRDTracker.h | 3 + .../TRDTracking/GPUTRDTrackletWord.h | 4 + 11 files changed, 370 insertions(+), 15 deletions(-) diff --git a/Detectors/TRD/qc/include/TRDQC/Tracking.h b/Detectors/TRD/qc/include/TRDQC/Tracking.h index f39c64286d0cc..d5a402949efb1 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" @@ -102,6 +103,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 @@ -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..faa4aecad778c 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++; } } @@ -64,6 +72,60 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) if (mPID) { 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); @@ -93,9 +155,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..b84fcb8a5e4f4 100644 --- a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx @@ -414,6 +414,23 @@ void TRDGlobalTracking::run(ProcessingContext& pc) LOGF(debug, "Loaded TPC track %i with time %f. Window from %f to %f", nTracksLoadedTPC, trkAttribs.mTime, trkAttribs.mTime - trkAttribs.mTimeSubMax, trkAttribs.mTime + trkAttribs.mTimeAddMax); } 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(); @@ -796,6 +813,46 @@ 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 @@ -826,10 +883,16 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, if (!((trkParam->getSigmaZ2() < (padLength * padLength / 12.f)) && (std::fabs(mTrackletsCalib[trkltId].getZ() - trkParam->getZ()) < padLength))) { tiltCorrUp = 0.f; } + + // 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, zPosCorrUp}; + 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..088a6f18467ff 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..c17012a28aa83 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h @@ -56,9 +56,15 @@ 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 @@ -74,8 +80,30 @@ class GPUTRDRecoParam float mCorrYDyC{0.f}; float mZCorrCoefNRC{1.4f}; ///< tracklet z-position depends linearly on track dip angle + + // 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, 3); + ClassDefNV(GPUTRDRecoParam, 4); }; } // namespace gpu diff --git a/GPU/GPUTracking/DataTypes/GPUTRDTrack.h b/GPU/GPUTracking/DataTypes/GPUTRDTrack.h index b358e8b82d480..cbb84e8fff306 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDTrack.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDTrack.h @@ -91,7 +91,7 @@ class GPUTRDTrack_t : public T GPUd() bool getIsFindable(int32_t iLayer) const { return (mFlags >> iLayer) & 0x1; } GPUd() int32_t getNmissingConsecLayers(int32_t iLayer) const; GPUd() int32_t getIsPenaltyAdded(int32_t iLayer) const { return getIsFindable(iLayer) && getTrackletIndex(iLayer) < 0; } - + // setters GPUd() void setRefGlobalTrackIdRaw(uint32_t id) { mRefGlobalTrackId = id; } // This method is only defined in TrackTRD.h and is intended to be used only with that TRD track type diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index 80098ff151ebe..94395bea59882 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,30 @@ 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; @@ -447,6 +472,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK const GPUTRDpadPlane* pad = nullptr; const GPUTRDTrackletWord* tracklets = GetConstantMem()->ioPtrs.trdTracklets; const GPUTRDSpacePoint* spacePoints = GetConstantMem()->ioPtrs.trdSpacePoints; + #ifdef ENABLE_GPUTRDDEBUG TRDTRK trackNoUp(*t); @@ -585,9 +611,37 @@ 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 +650,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 +751,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 +862,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..3023d90c5c65a 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h @@ -132,6 +132,7 @@ 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 +171,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 }; From df845378b68d4b7eb46681bb5ed02f14e0e817f2 Mon Sep 17 00:00:00 2001 From: Gauthier Legras Date: Wed, 15 Apr 2026 17:56:00 +0200 Subject: [PATCH 2/9] clang format --- Detectors/TRD/qc/include/TRDQC/Tracking.h | 20 +++++----- Detectors/TRD/qc/src/Tracking.cxx | 27 +++++++------ .../workflow/src/TRDGlobalTrackingSpec.cxx | 26 ++++++------- GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx | 24 ++++++------ GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h | 26 ++++++------- GPU/GPUTracking/DataTypes/GPUTRDTrack.h | 2 +- GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx | 39 +++++++++---------- GPU/GPUTracking/TRDTracking/GPUTRDTracker.h | 6 ++- 8 files changed, 86 insertions(+), 84 deletions(-) diff --git a/Detectors/TRD/qc/include/TRDQC/Tracking.h b/Detectors/TRD/qc/include/TRDQC/Tracking.h index d5a402949efb1..c1c44c1a07dce 100644 --- a/Detectors/TRD/qc/include/TRDQC/Tracking.h +++ b/Detectors/TRD/qc/include/TRDQC/Tracking.h @@ -103,7 +103,7 @@ class Tracking { mLocalGain = localGain; } - + // quantities necessary for pile-up correction void setTriggeredBCFT0(std::vector t) { mTriggeredBCFT0 = t; } void setFirstOrbit(uint32_t o) { mFirstOrbit = o; } @@ -120,17 +120,17 @@ 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 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 + 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 + std::vector mTriggeredBCFT0; ///< array with the FT0 trigger times int mCurrentTriggerRecord; uint32_t mFirstOrbit; int mCurrentTrackId; diff --git a/Detectors/TRD/qc/src/Tracking.cxx b/Detectors/TRD/qc/src/Tracking.cxx index faa4aecad778c..514488a3c2d0a 100644 --- a/Detectors/TRD/qc/src/Tracking.cxx +++ b/Detectors/TRD/qc/src/Tracking.cxx @@ -72,7 +72,7 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) if (mPID) { 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++) { @@ -80,10 +80,10 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) 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; + 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; @@ -108,8 +108,7 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) if (trkltId < 0) { q0[iLy] = -1; q1[iLy] = -1; - } - else { + } else { q0[iLy] = mTrackletsRaw[trkltId].getQ0(); q1[iLy] = mTrackletsRaw[trkltId].getQ1(); } @@ -121,11 +120,11 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) sumProb += probBC; if (probBC > maxProb) { maxProb = probBC; - tCorrPileUp = - deltaBC; + tCorrPileUp = -deltaBC; } - } - if (sumProb > 1e-6) tErrPileUp2 = sumCorr2 / sumProb - 2 * tCorrPileUp * sumCorr / sumProb + tCorrPileUp * tCorrPileUp; - + } + 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); @@ -155,12 +154,12 @@ 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; } - - // 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; + + // 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); diff --git a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx index b84fcb8a5e4f4..7404b71c59ad5 100644 --- a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx @@ -414,7 +414,7 @@ void TRDGlobalTracking::run(ProcessingContext& pc) LOGF(debug, "Loaded TPC track %i with time %f. Window from %f to %f", nTracksLoadedTPC, trkAttribs.mTime, trkAttribs.mTime - trkAttribs.mTimeSubMax, trkAttribs.mTime + trkAttribs.mTimeAddMax); } 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 @@ -422,14 +422,15 @@ void TRDGlobalTracking::run(ProcessingContext& pc) 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 (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 @@ -813,7 +814,7 @@ 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; @@ -835,8 +836,7 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, if (trkltId < 0) { q0[iLy] = -1; q1[iLy] = -1; - } - else { + } else { q0[iLy] = mTrackletsRaw[trkltId].getQ0(); q1[iLy] = mTrackletsRaw[trkltId].getQ1(); } @@ -848,11 +848,11 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, sumProb += probBC; if (probBC > maxProb) { maxProb = probBC; - tCorrPileUp = - deltaBC; + tCorrPileUp = -deltaBC; } - } - if (sumProb > 1e-6) tErrPileUp2 = sumCorr2 / sumProb - 2 * tCorrPileUp * sumCorr / sumProb + tCorrPileUp * tCorrPileUp; - + } + if (sumProb > 1e-6) + tErrPileUp2 = sumCorr2 / sumProb - 2 * tCorrPileUp * sumCorr / sumProb + tCorrPileUp * tCorrPileUp; if (inwards) { // reset covariance to something big for inwards refit @@ -883,9 +883,9 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, if (!((trkParam->getSigmaZ2() < (padLength * padLength / 12.f)) && (std::fabs(mTrackletsCalib[trkltId].getZ() - trkParam->getZ()) < padLength))) { tiltCorrUp = 0.f; } - - // 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; + + // 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; diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx index 088a6f18467ff..11af412f7782b 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx @@ -78,7 +78,6 @@ void GPUTRDRecoParam::recalcTrkltCov(const float tilt, const float snp, const fl 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 @@ -106,14 +105,14 @@ float GPUTRDRecoParam::getPileUpProbTracklet(int nBC, int Q0, int Q1) const 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) + if (nBC <= minBC || nBC >= maxBC) return 0.; - float maxProb = 2./(maxBC - minBC); - if (nBC > minBC && nBC <= maxProbBC) + float maxProb = 2. / (maxBC - minBC); + if (nBC > minBC && nBC <= maxProbBC) return maxProb / (maxProbBC - minBC) * (nBC - minBC); - else + else return maxProb / (maxProbBC - maxBC) * (nBC - maxBC); } @@ -126,19 +125,20 @@ float GPUTRDRecoParam::getPileUpProbTrack(int nBC, std::array Q0, std::a // 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 - + // 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.; - + 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; + if (Q0[i] < 0 || Q1[i] < 0) + continue; float probTracklet = GPUTRDRecoParam::getPileUpProbTracklet(nBC, Q0[i], Q1[i]); probTrack *= probTracklet / probNoInfo; } diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h index c17012a28aa83..851159cfcf8ab 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h @@ -61,7 +61,7 @@ class GPUTRDRecoParam /// 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; } @@ -80,28 +80,28 @@ class GPUTRDRecoParam float mCorrYDyC{0.f}; float mZCorrCoefNRC{1.4f}; ///< tracklet z-position depends linearly on track dip angle - + // 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 + 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 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 + 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 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 + 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); }; diff --git a/GPU/GPUTracking/DataTypes/GPUTRDTrack.h b/GPU/GPUTracking/DataTypes/GPUTRDTrack.h index cbb84e8fff306..b358e8b82d480 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDTrack.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDTrack.h @@ -91,7 +91,7 @@ class GPUTRDTrack_t : public T GPUd() bool getIsFindable(int32_t iLayer) const { return (mFlags >> iLayer) & 0x1; } GPUd() int32_t getNmissingConsecLayers(int32_t iLayer) const; GPUd() int32_t getIsPenaltyAdded(int32_t iLayer) const { return getIsFindable(iLayer) && getTrackletIndex(iLayer) < 0; } - + // setters GPUd() void setRefGlobalTrackIdRaw(uint32_t id) { mRefGlobalTrackId = id; } // This method is only defined in TrackTRD.h and is intended to be used only with that TRD track type diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index 94395bea59882..2f338800bafb7 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -351,7 +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 @@ -437,11 +437,11 @@ 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()) { @@ -460,7 +460,6 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK } } - float zShiftTrk = 0.f; if (mProcessPerTimeFrame) { zShiftTrk = (mTrackAttribs[iTrk].mTime - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId]) * mTPCVdrift * mTrackAttribs[iTrk].mSide; @@ -472,7 +471,6 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK const GPUTRDpadPlane* pad = nullptr; const GPUTRDTrackletWord* tracklets = GetConstantMem()->ioPtrs.trdTracklets; const GPUTRDSpacePoint* spacePoints = GetConstantMem()->ioPtrs.trdSpacePoints; - #ifdef ENABLE_GPUTRDDEBUG TRDTRK trackNoUp(*t); @@ -611,8 +609,8 @@ 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, + + // 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; @@ -623,8 +621,8 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK 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; + // 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()); @@ -633,10 +631,11 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK sumProb += probBC; if (probBC > maxProb) { maxProb = probBC; - yCorrPileUp = - slopeFactor * deltaBC; + yCorrPileUp = -slopeFactor * deltaBC; } } - if (sumProb > 1e-6f) yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp; + 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) @@ -751,8 +750,8 @@ 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; } - - // 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, + + // 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; @@ -763,8 +762,8 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK 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; + // 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()); @@ -773,13 +772,13 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK sumProb += probBC; if (probBC > maxProb) { maxProb = probBC; - yCorrPileUp = - slopeFactor * deltaBC; + yCorrPileUp = -slopeFactor * deltaBC; } } - if (sumProb > 1e-6f) yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp; + 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); @@ -862,7 +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 3023d90c5c65a..d491539b8fcf8 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h @@ -132,7 +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() void SetFT0TriggeredBC(int32_t* t, int32_t n) + { + mFT0TriggeredBC = t; + mNFT0BC = n; + } GPUd() bool GetIsDebugOutputOn() const { return mDebugOutput; } GPUd() float GetMaxEta() const { return mMaxEta; } From 65e927c10c4ea5051308d9e5df6ca3bcd9bd9f25 Mon Sep 17 00:00:00 2001 From: Gauthier Legras Date: Fri, 17 Apr 2026 14:40:50 +0200 Subject: [PATCH 3/9] small update --- Detectors/TRD/qc/src/Tracking.cxx | 3 +- .../TRDWorkflow/TRDGlobalTrackingQCSpec.h | 3 +- GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx | 36 +++++++++++++++---- GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h | 2 +- GPU/GPUTracking/Definitions/GPUSettingsList.h | 4 +-- 5 files changed, 35 insertions(+), 13 deletions(-) diff --git a/Detectors/TRD/qc/src/Tracking.cxx b/Detectors/TRD/qc/src/Tracking.cxx index 514488a3c2d0a..3360ff2c97c3f 100644 --- a/Detectors/TRD/qc/src/Tracking.cxx +++ b/Detectors/TRD/qc/src/Tracking.cxx @@ -78,8 +78,7 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) 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; + triggeredBC = tRecord.getBCData().differenceInBC({0, mFirstOrbit});; break; } } diff --git a/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h b/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h index e52b3bcc23842..25afc48291db1 100644 --- a/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h +++ b/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h @@ -80,8 +80,7 @@ class TRDGlobalTrackingQC : public Task 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); + triggeredBCFT0.push_back(f0rec.getInteractionRecord().differenceInBC({0,firstOrbit})); } } } diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx index 11af412f7782b..3908d80dc7b9f 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx @@ -27,11 +27,15 @@ void GPUTRDRecoParam::init(float bz, const GPUSettingsRec* rec) float resRPhiIdeal2 = 1.6e-3f; if (rec) { resRPhiIdeal2 = rec->trd.trkltResRPhiIdeal * rec->trd.trkltResRPhiIdeal; + mPileUpRangeBefore = -rec->trd.pileupBwdNBC; + mPileUpRangeAfter = rec->trd.pileupFwdNBC; } #ifndef GPUCA_STANDALONE else { const auto& rtrd = GPU_GET_CONFIG(GPUSettingsRecTRD); resRPhiIdeal2 = rtrd.trkltResRPhiIdeal * rtrd.trkltResRPhiIdeal; + mPileUpRangeBefore = -rtrd.pileupBwdNBC; + mPileUpRangeAfter = rtrd.pileupFwdNBC; } #endif @@ -78,13 +82,19 @@ void GPUTRDRecoParam::recalcTrkltCov(const float tilt, const float snp, const fl cov[2] = c2 * (t2 * sy2 + sz2); } -float GPUTRDRecoParam::getPileUpProbTracklet(int nBC, int Q0, int Q1) const +float getPileUpProbTracklet(int nBC, int Q0, int Q1) { // 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 + + float prob = 0.; + int maxBC = mPileUpRangeAfter; int minBC = mPileUpRangeBefore; int maxProbBC = mPileUpMaxProb; + if (nBC <= mPileUpRangeBefore || nBC >= mPileUpRangeAfter) + return prob; + if (Q0 > 0 && Q1 > 0) { maxBC = mPileUpRangeAfter11; minBC = mPileUpRangeBefore11; @@ -99,6 +109,16 @@ float GPUTRDRecoParam::getPileUpProbTracklet(int nBC, int Q0, int Q1) const maxBC = mPileUpRangeAfter10; minBC = mPileUpRangeBefore10; maxProbBC = mPileUpMaxProb10; + + // if Q1 = 0, there is a second maximum at nBC=0, probably due to tracklets with low energy loss in the drift/TR regions + // so we enlarge the probability around there + if (nBC > maxProbBC && nBC <= 0) { + prob += 2. / (maxBC - minBC) / (0 - maxProbBC) * (nBC - maxProbBC); + } + if (nBC > 0 && nBC < maxBC) { + prob += 2. / (maxBC - minBC) / (0 - maxBC) * (nBC - maxBC); + } + } if (Q0 == 0 && Q1 == 0) { maxBC = mPileUpRangeAfter00; @@ -107,13 +127,17 @@ float GPUTRDRecoParam::getPileUpProbTracklet(int nBC, int Q0, int Q1) const } // 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) + 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); + if (nBC > minBC && nBC <= maxProbBC) { + prob += maxProb / (maxProbBC - minBC) * (nBC - minBC); + } + else { + prob += maxProb / (maxProbBC - maxBC) * (nBC - maxBC); + } + return prob; } float GPUTRDRecoParam::getPileUpProbTrack(int nBC, std::array Q0, std::array Q1) const diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h index 851159cfcf8ab..b5c9d3c77c16e 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h @@ -96,7 +96,7 @@ class GPUTRDRecoParam 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 mPileUpMaxProb10{-60}; ///< 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 diff --git a/GPU/GPUTracking/Definitions/GPUSettingsList.h b/GPU/GPUTracking/Definitions/GPUSettingsList.h index 606deb44d9528..df611cc156951 100644 --- a/GPU/GPUTracking/Definitions/GPUSettingsList.h +++ b/GPU/GPUTracking/Definitions/GPUSettingsList.h @@ -188,8 +188,8 @@ AddOptionRTC(addDeflectionInChi2, uint8_t, 0, "", 0, "Set to 1 to add the deflec AddOptionRTC(stopTrkAfterNMissLy, uint8_t, 6, "", 0, "Abandon track following after N layers without a TRD match") AddOptionRTC(nTrackletsMin, uint8_t, 3, "", 0, "Tracks with less attached tracklets are discarded after the tracking") AddOptionRTC(matCorrType, uint8_t, 2, "", 0, "Material correction to use: 0 - none, 1 - TGeo, 2 - matLUT") -AddOptionRTC(pileupFwdNBC, uint8_t, 80, "", 0, "Post-trigger Pile-up integration time in BCs") -AddOptionRTC(pileupBwdNBC, uint8_t, 80, "", 0, "Pre-trigger Pile-up integration time in BCs") +AddOptionRTC(pileupFwdNBC, uint8_t, 70, "", 0, "Post-trigger Pile-up integration time in BCs") +AddOptionRTC(pileupBwdNBC, uint8_t, 130, "", 0, "Pre-trigger Pile-up integration time in BCs") AddHelp("help", 'h') EndConfig() From 54be31a6cd396f2dfb1ae71739b718ba4535f093 Mon Sep 17 00:00:00 2001 From: Gauthier Legras Date: Fri, 17 Apr 2026 14:42:25 +0200 Subject: [PATCH 4/9] clang format --- Detectors/TRD/qc/src/Tracking.cxx | 3 ++- .../include/TRDWorkflow/TRDGlobalTrackingQCSpec.h | 2 +- GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx | 14 ++++++-------- 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/Detectors/TRD/qc/src/Tracking.cxx b/Detectors/TRD/qc/src/Tracking.cxx index 3360ff2c97c3f..3acb47bf36caa 100644 --- a/Detectors/TRD/qc/src/Tracking.cxx +++ b/Detectors/TRD/qc/src/Tracking.cxx @@ -78,7 +78,8 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) for (; mCurrentTriggerRecord < (isTPCTRD ? mTrackTriggerRecordsTPCTRD.size() : mTrackTriggerRecordsITSTPCTRD.size()); mCurrentTriggerRecord++) { auto& tRecord = (isTPCTRD ? mTrackTriggerRecordsTPCTRD[mCurrentTriggerRecord] : mTrackTriggerRecordsITSTPCTRD[mCurrentTriggerRecord]); if (mCurrentTrackId >= tRecord.getFirstTrack() && mCurrentTrackId < tRecord.getFirstTrack() + tRecord.getNumberOfTracks()) { - triggeredBC = tRecord.getBCData().differenceInBC({0, mFirstOrbit});; + triggeredBC = tRecord.getBCData().differenceInBC({0, mFirstOrbit}); + ; break; } } diff --git a/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h b/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h index 25afc48291db1..fc86d197f4df5 100644 --- a/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h +++ b/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h @@ -80,7 +80,7 @@ class TRDGlobalTrackingQC : public Task mQC.setFirstOrbit(firstOrbit); } if (o2::ft0::InteractionTag::Instance().isSelected(f0rec)) { - triggeredBCFT0.push_back(f0rec.getInteractionRecord().differenceInBC({0,firstOrbit})); + triggeredBCFT0.push_back(f0rec.getInteractionRecord().differenceInBC({0, firstOrbit})); } } } diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx index 3908d80dc7b9f..bc3b3b0877570 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx @@ -86,15 +86,15 @@ float getPileUpProbTracklet(int nBC, int Q0, int Q1) { // 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 - + float prob = 0.; - + int maxBC = mPileUpRangeAfter; int minBC = mPileUpRangeBefore; int maxProbBC = mPileUpMaxProb; - if (nBC <= mPileUpRangeBefore || nBC >= mPileUpRangeAfter) + if (nBC <= mPileUpRangeBefore || nBC >= mPileUpRangeAfter) return prob; - + if (Q0 > 0 && Q1 > 0) { maxBC = mPileUpRangeAfter11; minBC = mPileUpRangeBefore11; @@ -109,7 +109,7 @@ float getPileUpProbTracklet(int nBC, int Q0, int Q1) maxBC = mPileUpRangeAfter10; minBC = mPileUpRangeBefore10; maxProbBC = mPileUpMaxProb10; - + // if Q1 = 0, there is a second maximum at nBC=0, probably due to tracklets with low energy loss in the drift/TR regions // so we enlarge the probability around there if (nBC > maxProbBC && nBC <= 0) { @@ -118,7 +118,6 @@ float getPileUpProbTracklet(int nBC, int Q0, int Q1) if (nBC > 0 && nBC < maxBC) { prob += 2. / (maxBC - minBC) / (0 - maxBC) * (nBC - maxBC); } - } if (Q0 == 0 && Q1 == 0) { maxBC = mPileUpRangeAfter00; @@ -133,8 +132,7 @@ float getPileUpProbTracklet(int nBC, int Q0, int Q1) float maxProb = 2. / (maxBC - minBC); if (nBC > minBC && nBC <= maxProbBC) { prob += maxProb / (maxProbBC - minBC) * (nBC - minBC); - } - else { + } else { prob += maxProb / (maxProbBC - maxBC) * (nBC - maxBC); } return prob; From 83911f50c21c95230cbd1788faea92a50c1c64c9 Mon Sep 17 00:00:00 2001 From: Gauthier Legras Date: Fri, 17 Apr 2026 15:58:03 +0200 Subject: [PATCH 5/9] small fix --- Detectors/TRD/qc/src/Tracking.cxx | 1 - GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/Detectors/TRD/qc/src/Tracking.cxx b/Detectors/TRD/qc/src/Tracking.cxx index 3acb47bf36caa..8ef0b9dd61a24 100644 --- a/Detectors/TRD/qc/src/Tracking.cxx +++ b/Detectors/TRD/qc/src/Tracking.cxx @@ -79,7 +79,6 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) auto& tRecord = (isTPCTRD ? mTrackTriggerRecordsTPCTRD[mCurrentTriggerRecord] : mTrackTriggerRecordsITSTPCTRD[mCurrentTriggerRecord]); if (mCurrentTrackId >= tRecord.getFirstTrack() && mCurrentTrackId < tRecord.getFirstTrack() + tRecord.getNumberOfTracks()) { triggeredBC = tRecord.getBCData().differenceInBC({0, mFirstOrbit}); - ; break; } } diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx index bc3b3b0877570..5f64b227ceb89 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx @@ -82,7 +82,7 @@ void GPUTRDRecoParam::recalcTrkltCov(const float tilt, const float snp, const fl cov[2] = c2 * (t2 * sy2 + sz2); } -float getPileUpProbTracklet(int nBC, int Q0, int Q1) +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 From 85c0920785cefdecd62ead4fbc3e9cd785fde81b Mon Sep 17 00:00:00 2001 From: Gauthier Legras Date: Wed, 22 Apr 2026 12:06:06 +0200 Subject: [PATCH 6/9] move getPileUpProbTrack to the header --- GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx | 86 ------------------ GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h | 91 ++++++++++++++++++- GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx | 4 +- 3 files changed, 92 insertions(+), 89 deletions(-) diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx index 5f64b227ceb89..d82039ee22b08 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx @@ -81,89 +81,3 @@ 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 - - float prob = 0.; - - int maxBC = mPileUpRangeAfter; - int minBC = mPileUpRangeBefore; - int maxProbBC = mPileUpMaxProb; - if (nBC <= mPileUpRangeBefore || nBC >= mPileUpRangeAfter) - return prob; - - 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 Q1 = 0, there is a second maximum at nBC=0, probably due to tracklets with low energy loss in the drift/TR regions - // so we enlarge the probability around there - if (nBC > maxProbBC && nBC <= 0) { - prob += 2. / (maxBC - minBC) / (0 - maxProbBC) * (nBC - maxProbBC); - } - if (nBC > 0 && nBC < maxBC) { - prob += 2. / (maxBC - minBC) / (0 - maxBC) * (nBC - maxBC); - } - } - 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) { - prob += maxProb / (maxProbBC - minBC) * (nBC - minBC); - } else { - prob += maxProb / (maxProbBC - maxBC) * (nBC - maxBC); - } - return prob; -} - -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 b5c9d3c77c16e..68c26e20adc30 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h @@ -56,7 +56,7 @@ 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 getPileUpProbTracklet(int nBC, bool withChargeInfo, bool Q0 = true, bool Q1 = true) const; GPUd() float getPileUpProbTrack(int nBC, std::array Q0, std::array Q1) const; /// Get tracklet z correction coefficient for track-eta based corraction @@ -106,6 +106,95 @@ class GPUTRDRecoParam ClassDefNV(GPUTRDRecoParam, 4); }; +GPUdi() float GPUTRDRecoParam::getPileUpProbTracklet(int nBC, bool withChargeInfo, bool Q0, bool 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 (bool is false) or not (bool is true) + + float prob = 0.; + + int maxBC = mPileUpRangeAfter; + int minBC = mPileUpRangeBefore; + int maxProbBC = mPileUpMaxProb; + if (nBC <= mPileUpRangeBefore || nBC >= mPileUpRangeAfter) { + return prob; + } + + if (withChargeInfo) { + if (Q0 && Q1) { + maxBC = mPileUpRangeAfter11; + minBC = mPileUpRangeBefore11; + maxProbBC = mPileUpMaxProb11; + } + if (!Q0 && Q1) { + maxBC = mPileUpRangeAfter01; + minBC = mPileUpRangeBefore01; + maxProbBC = mPileUpMaxProb01; + } + if (Q0 && !Q1) { + maxBC = mPileUpRangeAfter10; + minBC = mPileUpRangeBefore10; + maxProbBC = mPileUpMaxProb10; + + // if Q1 = 0, there is a second maximum at nBC=0, probably due to tracklets with low energy loss in the drift/TR regions + // so we enlarge the probability around there + if (nBC > maxProbBC && nBC <= 0) { + prob += 2. / (maxBC - minBC) / (0 - maxProbBC) * (nBC - maxProbBC); + } + if (nBC > 0 && nBC < maxBC) { + prob += 2. / (maxBC - minBC) / (0 - maxBC) * (nBC - maxBC); + } + } + if (!Q0 && !Q1) { + 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) { + prob += maxProb / (maxProbBC - minBC) * (nBC - minBC); + } else { + prob += maxProb / (maxProbBC - maxBC) * (nBC - maxBC); + } + return prob; +} + +GPUdi() 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, false); + + 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, true, (Q0[i] != 0), (Q1[i] != 0)); + probTrack *= probTracklet / probNoInfo; + } + + return probTrack; +} + } // namespace gpu } // namespace o2 diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index 2f338800bafb7..b8307760cb993 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -625,7 +625,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK 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()); + float probBC = mRecoParam->getPileUpProbTracklet(deltaBC, true, (tracklets[trkltIdx].GetQ0() != 0), (tracklets[trkltIdx].GetQ1() != 0)); sumCorr += probBC * slopeFactor * deltaBC; sumCorr2 += probBC * slopeFactor * deltaBC * slopeFactor * deltaBC; sumProb += probBC; @@ -766,7 +766,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK 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()); + float probBC = mRecoParam->getPileUpProbTracklet(deltaBC, true, (tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetQ0() != 0), (tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetQ1() != 0)); sumCorr += probBC * slopeFactor * deltaBC; sumCorr2 += probBC * slopeFactor * deltaBC * slopeFactor * deltaBC; sumProb += probBC; From db38fa1093c4b32e45778b8f467c1f1d7ea59d73 Mon Sep 17 00:00:00 2001 From: Gauthier Legras Date: Thu, 23 Apr 2026 09:47:25 +0200 Subject: [PATCH 7/9] fix checkcode errors --- Detectors/TRD/qc/src/Tracking.cxx | 3 ++- Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx | 6 ++++-- GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx | 6 ++++-- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/Detectors/TRD/qc/src/Tracking.cxx b/Detectors/TRD/qc/src/Tracking.cxx index 8ef0b9dd61a24..a296eb9038434 100644 --- a/Detectors/TRD/qc/src/Tracking.cxx +++ b/Detectors/TRD/qc/src/Tracking.cxx @@ -122,8 +122,9 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) tCorrPileUp = -deltaBC; } } - if (sumProb > 1e-6) + 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); diff --git a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx index 7404b71c59ad5..871c450813d03 100644 --- a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx @@ -422,8 +422,9 @@ void TRDGlobalTracking::run(ProcessingContext& pc) uint32_t firstOrbit = 0; for (size_t ft0id = 0; ft0id < ft0recPoints.size(); ft0id++) { const auto& f0rec = ft0recPoints[ft0id]; - if (ft0id == 0) + 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); @@ -851,8 +852,9 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, tCorrPileUp = -deltaBC; } } - if (sumProb > 1e-6) + if (sumProb > 1e-6) { tErrPileUp2 = sumCorr2 / sumProb - 2 * tCorrPileUp * sumCorr / sumProb + tCorrPileUp * tCorrPileUp; + } if (inwards) { // reset covariance to something big for inwards refit diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index b8307760cb993..3e8dd7c3464e8 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -634,8 +634,9 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK yCorrPileUp = -slopeFactor * deltaBC; } } - if (sumProb > 1e-6f) + 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) @@ -775,8 +776,9 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK yCorrPileUp = -slopeFactor * deltaBC; } } - if (sumProb > 1e-6f) + 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}; From 501f33d466f070146599fb6297e4f2815aeb65da Mon Sep 17 00:00:00 2001 From: Gauthier Legras Date: Tue, 16 Jun 2026 10:49:45 +0200 Subject: [PATCH 8/9] update error parametrization --- Detectors/TRD/qc/src/Tracking.cxx | 5 +- .../workflow/src/TRDGlobalTrackingSpec.cxx | 6 +- GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx | 22 +++---- GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h | 59 +++++++++++++------ GPU/GPUTracking/Definitions/GPUSettingsList.h | 2 + GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx | 49 ++++++++++----- GPU/GPUTracking/TRDTracking/GPUTRDTracker.h | 7 ++- 7 files changed, 101 insertions(+), 49 deletions(-) diff --git a/Detectors/TRD/qc/src/Tracking.cxx b/Detectors/TRD/qc/src/Tracking.cxx index a296eb9038434..db171d4d7621e 100644 --- a/Detectors/TRD/qc/src/Tracking.cxx +++ b/Detectors/TRD/qc/src/Tracking.cxx @@ -149,6 +149,7 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) const PadPlane* pad = Geometry::instance()->getPadPlane(trkltDet); float tilt = tan(TMath::DegToRad() * pad->getTiltingAngle()); // tilt is signed! and returned in degrees float tiltCorrUp = tilt * (mTrackletsCalib[trkltId].getZ() - trk.getZ()); + float dyTiltCorr = tilt * trk.getTgl() * Geometry::instance()->cdrHght(); float zPosCorrUp = mTrackletsCalib[trkltId].getZ() + mRecoParam.getZCorrCoeffNRC() * trk.getTgl(); float padLength = pad->getRowSize(tracklet.getPadRow()); if (!((trk.getSigmaZ2() < (padLength * padLength / 12.f)) && (std::fabs(mTrackletsCalib[trkltId].getZ() - trk.getZ()) < padLength))) { @@ -159,10 +160,12 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) float slopeFactor = mTrackletsRaw[trkltId].getSlopeFloat() * pad->getWidthIPad() / 4.f; float yCorrPileUp = tCorrPileUp * slopeFactor; float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor; + + float angularPull = (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trk.getSnp())) / std::sqrt(mRecoParam.getDyRes(trk.getSnp(), 0)); std::array trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp, zPosCorrUp}; std::array trkltCovUp; - mRecoParam.recalcTrkltCov(tilt, trk.getSnp(), pad->getRowSize(tracklet.getPadRow()), trkltCovUp); + mRecoParam.recalcTrkltCov(tilt, trk.getSnp(), pad->getRowSize(tracklet.getPadRow()), trkltCovUp, angularPull, 0); trkltCovUp[0] += yAddErrPileUp2; auto chi2trklt = trk.getPredictedChi2(trkltPosUp, trkltCovUp); diff --git a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx index 871c450813d03..b6249fe1c6029 100644 --- a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx @@ -879,6 +879,7 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, } const PadPlane* pad = Geometry::instance()->getPadPlane(trkltDet); float tilt = tan(TMath::DegToRad() * pad->getTiltingAngle()); // tilt is signed! and returned in degrees + float dyTiltCorr = tilt * trkParam->getTgl() * Geometry::instance()->cdrHght(); float tiltCorrUp = tilt * (mTrackletsCalib[trkltId].getZ() - trkParam->getZ()); float zPosCorrUp = mTrackletsCalib[trkltId].getZ() + mRecoParam.getZCorrCoeffNRC() * trkParam->getTgl(); float padLength = pad->getRowSize(mTrackletsRaw[trkltId].getPadRow()); @@ -890,10 +891,13 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, float slopeFactor = mTrackletsRaw[trkltId].getSlopeFloat() * pad->getWidthIPad() / 4.f; float yCorrPileUp = tCorrPileUp * slopeFactor; float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor; + + int nTrackletsChamber = mTracker->GetNtrackletsChamber(trk.getCollisionId(), trkltDet); + float angularPull = (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trkParam->getSnp())) / std::sqrt(mRecoParam.getDyRes(trkParam->getSnp(), nTrackletsChamber)); std::array trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp, zPosCorrUp}; std::array trkltCovUp; - mRecoParam.recalcTrkltCov(tilt, trkParam->getSnp(), pad->getRowSize(mTrackletsRaw[trkltId].getPadRow()), trkltCovUp); + mRecoParam.recalcTrkltCov(tilt, trkParam->getSnp(), pad->getRowSize(mTrackletsRaw[trkltId].getPadRow()), trkltCovUp, (mRec->GetParam().rec.trd.useAngularPull != 0 ? angularPull : 0.), nTrackletsChamber); trkltCovUp[0] += yAddErrPileUp2; chi2 += trkParam->getPredictedChi2(trkltPosUp, trkltCovUp); diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx index d82039ee22b08..dfe524c38ac6f 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx @@ -24,16 +24,19 @@ using namespace o2::gpu; // error parameterizations taken from http://cds.cern.ch/record/2724259 Appendix A void GPUTRDRecoParam::init(float bz, const GPUSettingsRec* rec) { - float resRPhiIdeal2 = 1.6e-3f; + float resRPhiIdeal = 0.04f; + float resVsTanPhiMisalign = 0.f; if (rec) { - resRPhiIdeal2 = rec->trd.trkltResRPhiIdeal * rec->trd.trkltResRPhiIdeal; + resRPhiIdeal = rec->trd.trkltResRPhiIdeal; + resVsTanPhiMisalign = rec->trd.trkltResVsTanPhiMisalign; mPileUpRangeBefore = -rec->trd.pileupBwdNBC; mPileUpRangeAfter = rec->trd.pileupFwdNBC; } #ifndef GPUCA_STANDALONE else { const auto& rtrd = GPU_GET_CONFIG(GPUSettingsRecTRD); - resRPhiIdeal2 = rtrd.trkltResRPhiIdeal * rtrd.trkltResRPhiIdeal; + resRPhiIdeal = rtrd.trkltResRPhiIdeal; + resVsTanPhiMisalign = rtrd.trkltResVsTanPhiMisalign; mPileUpRangeBefore = -rtrd.pileupBwdNBC; mPileUpRangeAfter = rtrd.pileupFwdNBC; } @@ -59,23 +62,22 @@ void GPUTRDRecoParam::init(float bz, const GPUSettingsRec* rec) LOGP(warning, "No error parameterization available for Bz= {}. Keeping default value (sigma_y = const. = 1cm)", bz); } - mRPhiA2 = resRPhiIdeal2; + mRPhiA = resRPhiIdeal; + mRPhiATgp = resVsTanPhiMisalign; mLorentzAngle = -0.02f + 0.13f * bz / 5.f; mDyA2 = 6e-3f; mDyC2 = 0.3f; - mCorrYDyA = 0.27f; - mCorrYDyC = -0.44f; - LOGP(info, "Loaded parameterizations for Bz={}: PhiRes:[{},{},{}] DyRes:[{},{},{}] CorrYDy:[{},{},{}]", - bz, mRPhiA2, mLorentzAngle, mRPhiC2, mDyA2, mLorentzAngle, mDyC2, mCorrYDyA, mLorentzAngle, mCorrYDyC); + LOGP(info, "Loaded parameterizations for Bz={}: PhiRes:[{},{},{},{}] DyRes:[{},{},{}]", + bz, mRPhiA, mRPhiATgp, mLorentzAngle, mRPhiC2, mDyA2, mLorentzAngle, mDyC2); } -void GPUTRDRecoParam::recalcTrkltCov(const float tilt, const float snp, const float rowSize, float* cov) const +void GPUTRDRecoParam::recalcTrkltCov(const float tilt, const float snp, const float rowSize, float* cov, const float pull, const int occupancy) const { float t2 = tilt * tilt; // tan^2 (tilt) float c2 = 1.f / (1.f + t2); // cos^2 (tilt) - float sy2 = getRPhiRes(snp); + float sy2 = getRPhiRes(snp, CAMath::Abs(pull), occupancy); float sz2 = rowSize * rowSize / 12.f; cov[0] = c2 * (sy2 + t2 * sz2); cov[1] = c2 * tilt * (sz2 - sy2); diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h index 68c26e20adc30..3f7b91a7ba9bd 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h @@ -39,23 +39,17 @@ class GPUTRDRecoParam #if !defined(GPUCA_GPUCODE_DEVICE) /// Recalculate tracklet covariance based on phi angle of related track - GPUd() void recalcTrkltCov(const float tilt, const float snp, const float rowSize, std::array& cov) const + GPUd() void recalcTrkltCov(const float tilt, const float snp, const float rowSize, std::array& cov, const float pull = 0., const int occupancy = 0) const { - recalcTrkltCov(tilt, snp, rowSize, cov.data()); + recalcTrkltCov(tilt, snp, rowSize, cov.data(), pull, occupancy); } #endif - GPUd() void recalcTrkltCov(const float tilt, const float snp, const float rowSize, float* cov) const; - - /// Get tracklet r-phi resolution for given phi angle - /// Resolution depends on the track angle sin(phi) = snp and is approximated by the formula - /// sigma_y(snp) = sqrt(a^2 + c^2 * (snp - b)^2) - /// more details are given in http://cds.cern.ch/record/2724259 in section 5.3.3 - /// \param phi angle of related track - /// \return sigma_y^2 of tracklet - GPUd() float getRPhiRes(float snp) const { return (mRPhiA2 + mRPhiC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle)); } - 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() void recalcTrkltCov(const float tilt, const float snp, const float rowSize, float* cov, const float pull = 0., const int occupancy = 0) const; + + GPUd() float getRPhiRes(float snp, float pull = 0.f, int occupancy = 0) const; + GPUd() float getDyRes(float snp, int occupancy = 0) const { return mDyA2 + mDyC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle) + mOccDyA * occupancy; } // 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() const { return mCorrYDy; } GPUd() float getPileUpProbTracklet(int nBC, bool withChargeInfo, bool Q0 = true, bool Q1 = true) const; GPUd() float getPileUpProbTrack(int nBC, std::array Q0, std::array Q1) const; @@ -70,14 +64,21 @@ class GPUTRDRecoParam // tracklet error parameterization depends on the magnetic field float mLorentzAngle{0.f}; // rphi - float mRPhiA2{1.f}; ///< parameterization for tracklet position resolution - float mRPhiC2{0.f}; ///< parameterization for tracklet position resolution + float mRPhiA{1.f}; ///< parameterization for tracklet position resolution + float mRPhiATgp{1.f}; ///< parameterization for tracklet position resolution + float mRPhiC2{0.f}; ///< parameterization for tracklet position resolution // angle float mDyA2{1.225e-3f}; ///< parameterization for tracklet angular resolution float mDyC2{0.f}; ///< parameterization for tracklet angular resolution - // correlation coefficient between y residual and dy residual - float mCorrYDyA{0.f}; - float mCorrYDyC{0.f}; + // variation in y when dy variates by one sigma (= cov / sigma_dy = corr * sigma_y) (valid within 2sigma of dy) + float mCorrYDy{0.13f}; + // error parametrization vs angular pull (pol2) + float mPullA{6.8e-3f}; + float mPullB{0.049f}; + // error parametrization of y position vs occupancy defined as ntracklets within chamber (prop to sqrt(occupancy)) + float mOccA{3.3e-4f}; + // error parametrization for dy vs occupancy defined as ntracklets within chamber (prop to sqrt(occupancy)) + float mOccDyA{2.5e-4f}; float mZCorrCoefNRC{1.4f}; ///< tracklet z-position depends linearly on track dip angle @@ -106,6 +107,26 @@ class GPUTRDRecoParam ClassDefNV(GPUTRDRecoParam, 4); }; +/// Get tracklet r-phi resolution for given phi angle +/// Resolution depends on the track angle sin(phi) = snp and is approximated by the formula +/// sigma_y(snp) = sqrt(a^2 + c^2 * (snp - b)^2) +/// more details are given in http://cds.cern.ch/record/2724259 in section 5.3.3 +/// \param phi angle of related track +/// \return sigma_y^2 of tracklet +/// also depend on absolute pull and on chamber occupancy +GPUdi() float GPUTRDRecoParam::getRPhiRes(float snp, float pull, int occupancy) const { + // flat uncertainty + radial-alignment uncertainty depending on tan(phi) + float tgp = (CAMath::Abs(snp) < 0.99999f) ? CAMath::Abs(snp) / CAMath::Sqrt(1 - snp * snp) : 1e6; + float resIdeal = mRPhiA + mRPhiATgp * tgp; + if (pull > 10) { + // parametrization does not really work well for such large pull values + pull = 10.f; + } + float resPull = mPullA * pull * pull + mPullB * pull; // parametrization as pol2 summed in quadrature + float resOccupancy = mOccA * occupancy; // parametrization as sqrt() summed in quadrature + return (resIdeal * resIdeal + mRPhiC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle) + resPull * resPull + resOccupancy); +} + GPUdi() float GPUTRDRecoParam::getPileUpProbTracklet(int nBC, bool withChargeInfo, bool Q0, bool 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 diff --git a/GPU/GPUTracking/Definitions/GPUSettingsList.h b/GPU/GPUTracking/Definitions/GPUSettingsList.h index df611cc156951..78fb5f6d6971c 100644 --- a/GPU/GPUTracking/Definitions/GPUSettingsList.h +++ b/GPU/GPUTracking/Definitions/GPUSettingsList.h @@ -182,6 +182,7 @@ AddOptionRTC(addTimeRoadITSTPC, float, 2.5f, "", 0, "Increase time search road b AddOptionRTC(extraRoadY, float, 5.f, "", 0, "Addition to search road around track prolongation along Y in cm") AddOptionRTC(extraRoadZ, float, 10.f, "", 0, "Addition to search road around track prolongation along Z in cm") AddOptionRTC(trkltResRPhiIdeal, float, 1.f, "", 0, "Optimal tracklet rphi resolution in cm (in case phi of track = lorentz angle)") +AddOptionRTC(trkltResVsTanPhiMisalign, float, 0.f, "", 0, "tan(phi) dependence of tracklet error due to radial misalignment (centered at 0 angle)") AddOptionRTC(maxChi2Red, float, 99.f, "", 0, "maximum chi2 per attached tracklet for TRD tracks TODO: currently effectively disabled, requires tuning") AddOptionRTC(applyDeflectionCut, uint8_t, 0, "", 0, "Set to 1 to enable tracklet selection based on deflection") AddOptionRTC(addDeflectionInChi2, uint8_t, 0, "", 0, "Set to 1 to add the deflection in the chi2 calculation for matching") @@ -190,6 +191,7 @@ AddOptionRTC(nTrackletsMin, uint8_t, 3, "", 0, "Tracks with less attached trackl AddOptionRTC(matCorrType, uint8_t, 2, "", 0, "Material correction to use: 0 - none, 1 - TGeo, 2 - matLUT") AddOptionRTC(pileupFwdNBC, uint8_t, 70, "", 0, "Post-trigger Pile-up integration time in BCs") AddOptionRTC(pileupBwdNBC, uint8_t, 130, "", 0, "Pre-trigger Pile-up integration time in BCs") +AddOptionRTC(useAngularPull, uint8_t, 1, "", 0, "0 = don't use angular pull; 1 = additional error based on angular pull for refit only; 2 = add error also for chi2") AddHelp("help", 'h') EndConfig() diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index 3e8dd7c3464e8..643b0daeb10bc 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -638,6 +638,9 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp; } } + // number of tracklets within the chamber is the current TRD occupancy estimator + int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + currDet + 1] - mTrackletIndexArray[trkltIdxOffset + currDet]; + float angularPull = GetAngularPull(spacePoints[trkltIdx].getDy() + dyTiltCorr, trkWork->getSnp(), nTrackletsChamber); // 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(); @@ -645,17 +648,18 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK 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; + float trkltPosTmpYZ[2] = {yPosCorr, zPosCorr}; float trkltCovTmp[3] = {0.f}; 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); + // tracklet is in window: get predicted chi2 for update and store tracklet index if best guess + RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[trkltIdx].GetZbin()), (Param().rec.trd.useAngularPull == 2 ? angularPull : 0.f), nTrackletsChamber, 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 float trkltCovTmpWithDy[6] = {trkltCovTmp[0], trkltCovTmp[1], trkltCovTmp[2], 0.f, 0.f, 0.f}; - RecalcTrkltCovDy(tilt, trkWork->getSnp(), trkltCovTmpWithDy); + RecalcTrkltCovDy(tilt, trkWork->getSnp(), (Param().rec.trd.useAngularPull == 2 ? angularPull : 0.f), nTrackletsChamber, trkltCovTmpWithDy); trkltCovTmpWithDy[0] += trkWork->getSigmaY2(); trkltCovTmpWithDy[1] += trkWork->getSigmaZY(); trkltCovTmpWithDy[2] += trkWork->getSigmaZ2(); @@ -667,7 +671,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK } } // TODO cut on angular pull should be made stricter when proper v-drift calibration for the TRD tracklets is implemented - if ((chi2 > Param().rec.trd.maxChi2) || (Param().rec.trd.applyDeflectionCut && CAMath::Abs(GetAngularPull(spacePoints[trkltIdx].getDy() + dyTiltCorr, trkWork->getSnp())) > 4)) { + if ((chi2 > Param().rec.trd.maxChi2) || (Param().rec.trd.applyDeflectionCut && CAMath::Abs(angularPull) > 4)) { continue; } Hypothesis hypo(trkWork->getNlayersFindable(), iCandidate, trkltIdx, trkWork->getChi2() + chi2); @@ -745,11 +749,13 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK pad = mGeo->GetPadPlane(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector()); float tiltCorrUp = tilt * (spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() - trkWork->getZ()); + float dyTiltCorr = tilt * trkWork->getTgl() * mGeo->GetCdrHght(); float zPosCorrUp = spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() + mRecoParam->getZCorrCoeffNRC() * trkWork->getTgl(); zPosCorrUp -= zShiftTrk; float padLength = pad->GetRowSize(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin()); if (!((trkWork->getSigmaZ2() < (padLength * padLength / 12.f)) && (CAMath::Abs(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() - trkWork->getZ()) < padLength))) { tiltCorrUp = 0.f; + 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, @@ -780,10 +786,14 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp; } } + + const auto currDet = tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector(); + int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + currDet + 1] - mTrackletIndexArray[trkltIdxOffset + currDet]; 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); + float angularPull = GetAngularPull(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getDy() + dyTiltCorr, trkWork->getSnp(), nTrackletsChamber); + RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin()), ((Param().rec.trd.useAngularPull != 0) ? angularPull : 0.f), nTrackletsChamber, trkltCovUp); trkltCovUp[0] += yAddErrPileUp2; #ifdef ENABLE_GPUTRDDEBUG @@ -846,7 +856,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK trkWork->setIsCrossingNeighbor(iLayer); trkWork->setHasPadrowCrossing(); } - const auto currDet = tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector(); + // Mark tracklets as Padrow crossing if they have a neighboring tracklet. for (int32_t trkltIdx = glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet]; trkltIdx < glbTrkltIdxOffset + mTrackletIndexArray[trkltIdxOffset + currDet + 1]; ++trkltIdx) { // skip orig tracklet @@ -1026,7 +1036,7 @@ GPUd() float GPUTRDTracker_t::GetAlphaOfSector(const int32_t sec) } template -GPUd() void GPUTRDTracker_t::RecalcTrkltCov(const float tilt, const float snp, const float rowSize, float (&cov)[3]) +GPUd() void GPUTRDTracker_t::RecalcTrkltCov(const float tilt, const float snp, const float rowSize, const float pull, const int occupancy, float (&cov)[3]) { //-------------------------------------------------------------------- // recalculate tracklet covariance taking track phi angle into account @@ -1034,7 +1044,7 @@ GPUd() void GPUTRDTracker_t::RecalcTrkltCov(const float tilt, cons //-------------------------------------------------------------------- float t2 = tilt * tilt; // tan^2 (tilt) float c2 = 1.f / (1.f + t2); // cos^2 (tilt) - float sy2 = mRecoParam->getRPhiRes(snp); + float sy2 = mRecoParam->getRPhiRes(snp, CAMath::Abs(pull), occupancy); float sz2 = rowSize * rowSize / 12.f; cov[0] = c2 * (sy2 + t2 * sz2); cov[1] = c2 * tilt * (sz2 - sy2); @@ -1042,14 +1052,14 @@ GPUd() void GPUTRDTracker_t::RecalcTrkltCov(const float tilt, cons } template -GPUd() void GPUTRDTracker_t::RecalcTrkltCovDy(const float tilt, const float snp, float (&cov)[6]) +GPUd() void GPUTRDTracker_t::RecalcTrkltCovDy(const float tilt, const float snp, const float pull, const int occupancy, float (&cov)[6]) { float t2 = tilt * tilt; // tan^2 (tilt) float c2 = 1.f / (1.f + t2); // cos^2 (tilt) - float sy2 = mRecoParam->getRPhiRes(snp); - float sdy2 = mRecoParam->getDyRes(snp); - cov[3] = mRecoParam->getCorrYDy(snp) * CAMath::Sqrt(sdy2 * c2 * sy2); - cov[4] = -tilt * mRecoParam->getCorrYDy(snp) * CAMath::Sqrt(sdy2 * c2 * sy2); + float sy2 = mRecoParam->getRPhiRes(snp, CAMath::Abs(pull), occupancy); + float sdy2 = mRecoParam->getDyRes(snp, occupancy); + cov[3] = mRecoParam->getCorrYDy() * CAMath::Sqrt(sdy2 * c2); + cov[4] = -tilt * mRecoParam->getCorrYDy() * CAMath::Sqrt(sdy2 * c2); cov[5] = sdy2; } @@ -1105,16 +1115,25 @@ GPUd() bool GPUTRDTracker_t::InvertCov(float (&cov)[6]) } template -GPUd() float GPUTRDTracker_t::GetAngularPull(float dYtracklet, float snp) const +GPUd() float GPUTRDTracker_t::GetAngularPull(float dYtracklet, float snp, int occupancy) const { float dYtrack = mRecoParam->convertAngleToDy(snp); - float dYresolution = mRecoParam->getDyRes(snp); + float dYresolution = mRecoParam->getDyRes(snp, occupancy); if (dYresolution < 1e-6f) { return 999.f; } return (dYtracklet - dYtrack) / CAMath::Sqrt(dYresolution); } +template +GPUd() int GPUTRDTracker_t::GetNtrackletsChamber(int collisionId, int detector) const +{ + // get the number of tracklets for a given chamber and trigger + int32_t trkltIdxOffset = collisionId * (kNChambers + 1); + int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + detector + 1] - mTrackletIndexArray[trkltIdxOffset + detector]; + return nTrackletsChamber; +} + template GPUd() void GPUTRDTracker_t::FindChambersInRoad(const TRDTRK* t, const float roadY, const float roadZ, const int32_t iLayer, int32_t* det, const float zMax, const float alpha, const float zShiftTrk) const { diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h index d491539b8fcf8..0f94732f7d536 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h @@ -115,13 +115,14 @@ class GPUTRDTracker_t : public GPUProcessor GPUd() bool AdjustSector(PROP* prop, TRDTRK* t) const; GPUd() int32_t GetSector(float alpha) const; GPUd() float GetAlphaOfSector(const int32_t sec) const; - GPUd() float GetAngularPull(float dYtracklet, float snp) const; - GPUd() void RecalcTrkltCov(const float tilt, const float snp, const float rowSize, float (&cov)[3]); - GPUd() void RecalcTrkltCovDy(const float tilt, const float snp, float (&cov)[6]); + GPUd() float GetAngularPull(float dYtracklet, float snp, int occupancy) const; + GPUd() void RecalcTrkltCov(const float tilt, const float snp, const float rowSize, const float pull, const int occupancy, float (&cov)[3]); + GPUd() void RecalcTrkltCovDy(const float tilt, const float snp, const float pull, const int occupancy, float (&cov)[6]); GPUd() bool InvertCov(float (&cov)[6]); GPUd() void FindChambersInRoad(const TRDTRK* t, const float roadY, const float roadZ, const int32_t iLayer, int32_t* det, const float zMax, const float alpha, const float zShiftTrk) const; GPUd() bool IsGeoFindable(const TRDTRK* t, const int32_t layer, const float alpha, const float zShiftTrk) const; GPUd() void InsertHypothesis(Hypothesis hypo, int32_t& nCurrHypothesis, int32_t idxOffset); + GPUd() int GetNtrackletsChamber(int32_t collisionId, int32_t detector) const; // settings GPUd() void SetGenerateSpacePoints(bool flag) { mGenerateSpacePoints = flag; } From b42cec9ef94b99b98dde65f32f69653fde80eb58 Mon Sep 17 00:00:00 2001 From: Gauthier Legras Date: Tue, 16 Jun 2026 10:52:02 +0200 Subject: [PATCH 9/9] clang format --- Detectors/TRD/qc/src/Tracking.cxx | 2 +- .../TRD/workflow/src/TRDGlobalTrackingSpec.cxx | 2 +- GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h | 15 ++++++++------- GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx | 6 +++--- 4 files changed, 13 insertions(+), 12 deletions(-) diff --git a/Detectors/TRD/qc/src/Tracking.cxx b/Detectors/TRD/qc/src/Tracking.cxx index db171d4d7621e..35f0734498a40 100644 --- a/Detectors/TRD/qc/src/Tracking.cxx +++ b/Detectors/TRD/qc/src/Tracking.cxx @@ -160,7 +160,7 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) float slopeFactor = mTrackletsRaw[trkltId].getSlopeFloat() * pad->getWidthIPad() / 4.f; float yCorrPileUp = tCorrPileUp * slopeFactor; float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor; - + float angularPull = (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trk.getSnp())) / std::sqrt(mRecoParam.getDyRes(trk.getSnp(), 0)); std::array trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp, zPosCorrUp}; diff --git a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx index b6249fe1c6029..90b60389799d7 100644 --- a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx @@ -891,7 +891,7 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, float slopeFactor = mTrackletsRaw[trkltId].getSlopeFloat() * pad->getWidthIPad() / 4.f; float yCorrPileUp = tCorrPileUp * slopeFactor; float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor; - + int nTrackletsChamber = mTracker->GetNtrackletsChamber(trk.getCollisionId(), trkltDet); float angularPull = (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trkParam->getSnp())) / std::sqrt(mRecoParam.getDyRes(trkParam->getSnp(), nTrackletsChamber)); diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h index 3f7b91a7ba9bd..d561349a37a43 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h @@ -48,7 +48,7 @@ class GPUTRDRecoParam GPUd() float getRPhiRes(float snp, float pull = 0.f, int occupancy = 0) const; GPUd() float getDyRes(float snp, int occupancy = 0) const { return mDyA2 + mDyC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle) + mOccDyA * occupancy; } // 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 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() const { return mCorrYDy; } GPUd() float getPileUpProbTracklet(int nBC, bool withChargeInfo, bool Q0 = true, bool Q1 = true) const; GPUd() float getPileUpProbTrack(int nBC, std::array Q0, std::array Q1) const; @@ -64,9 +64,9 @@ class GPUTRDRecoParam // tracklet error parameterization depends on the magnetic field float mLorentzAngle{0.f}; // rphi - float mRPhiA{1.f}; ///< parameterization for tracklet position resolution - float mRPhiATgp{1.f}; ///< parameterization for tracklet position resolution - float mRPhiC2{0.f}; ///< parameterization for tracklet position resolution + float mRPhiA{1.f}; ///< parameterization for tracklet position resolution + float mRPhiATgp{1.f}; ///< parameterization for tracklet position resolution + float mRPhiC2{0.f}; ///< parameterization for tracklet position resolution // angle float mDyA2{1.225e-3f}; ///< parameterization for tracklet angular resolution float mDyC2{0.f}; ///< parameterization for tracklet angular resolution @@ -114,7 +114,8 @@ class GPUTRDRecoParam /// \param phi angle of related track /// \return sigma_y^2 of tracklet /// also depend on absolute pull and on chamber occupancy -GPUdi() float GPUTRDRecoParam::getRPhiRes(float snp, float pull, int occupancy) const { +GPUdi() float GPUTRDRecoParam::getRPhiRes(float snp, float pull, int occupancy) const +{ // flat uncertainty + radial-alignment uncertainty depending on tan(phi) float tgp = (CAMath::Abs(snp) < 0.99999f) ? CAMath::Abs(snp) / CAMath::Sqrt(1 - snp * snp) : 1e6; float resIdeal = mRPhiA + mRPhiATgp * tgp; @@ -123,8 +124,8 @@ GPUdi() float GPUTRDRecoParam::getRPhiRes(float snp, float pull, int occupancy) pull = 10.f; } float resPull = mPullA * pull * pull + mPullB * pull; // parametrization as pol2 summed in quadrature - float resOccupancy = mOccA * occupancy; // parametrization as sqrt() summed in quadrature - return (resIdeal * resIdeal + mRPhiC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle) + resPull * resPull + resOccupancy); + float resOccupancy = mOccA * occupancy; // parametrization as sqrt() summed in quadrature + return (resIdeal * resIdeal + mRPhiC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle) + resPull * resPull + resOccupancy); } GPUdi() float GPUTRDRecoParam::getPileUpProbTracklet(int nBC, bool withChargeInfo, bool Q0, bool Q1) const diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index 643b0daeb10bc..324acc9451f36 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -786,9 +786,9 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp; } } - + const auto currDet = tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector(); - int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + currDet + 1] - mTrackletIndexArray[trkltIdxOffset + currDet]; + int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + currDet + 1] - mTrackletIndexArray[trkltIdxOffset + currDet]; float trkltPosUp[2] = {spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getY() - tiltCorrUp + yCorrPileUp, zPosCorrUp}; float trkltCovUp[3] = {0.f}; @@ -1130,7 +1130,7 @@ GPUd() int GPUTRDTracker_t::GetNtrackletsChamber(int collisionId, { // get the number of tracklets for a given chamber and trigger int32_t trkltIdxOffset = collisionId * (kNChambers + 1); - int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + detector + 1] - mTrackletIndexArray[trkltIdxOffset + detector]; + int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + detector + 1] - mTrackletIndexArray[trkltIdxOffset + detector]; return nTrackletsChamber; }