diff --git a/Detectors/TRD/qc/include/TRDQC/Tracking.h b/Detectors/TRD/qc/include/TRDQC/Tracking.h index f39c64286d0cc..c1c44c1a07dce 100644 --- a/Detectors/TRD/qc/include/TRDQC/Tracking.h +++ b/Detectors/TRD/qc/include/TRDQC/Tracking.h @@ -23,6 +23,7 @@ #include "DataFormatsTRD/Constants.h" #include "ReconstructionDataFormats/TrackTPCITS.h" #include "ReconstructionDataFormats/GlobalTrackID.h" +#include "DataFormatsTRD/TrackTriggerRecord.h" #include "DataFormatsTPC/TrackTPC.h" #include "DetectorsBase/Propagator.h" #include "GPUTRDRecoParam.h" @@ -103,6 +104,10 @@ class Tracking mLocalGain = localGain; } + // quantities necessary for pile-up correction + void setTriggeredBCFT0(std::vector t) { mTriggeredBCFT0 = t; } + void setFirstOrbit(uint32_t o) { mFirstOrbit = o; } + private: float mMaxSnp{o2::base::Propagator::MAX_SIN_PHI}; ///< max snp when propagating tracks float mMaxStep{o2::base::Propagator::MAX_STEP}; ///< maximum step for propagation @@ -115,12 +120,20 @@ class Tracking std::vector mTrackQC; // input from DPL - gsl::span mTracksITSTPC; ///< ITS-TPC seeding tracks - gsl::span mTracksTPC; ///< TPC seeding tracks - gsl::span mTracksITSTPCTRD; ///< TRD tracks reconstructed from TPC or ITS-TPC seeds - gsl::span mTracksTPCTRD; ///< TRD tracks reconstructed from TPC or TPC seeds - gsl::span mTrackletsRaw; ///< array of raw tracklets needed for TRD refit - gsl::span mTrackletsCalib; ///< array of calibrated tracklets needed for TRD refit + gsl::span mTracksITSTPC; ///< ITS-TPC seeding tracks + gsl::span mTracksTPC; ///< TPC seeding tracks + gsl::span mTracksITSTPCTRD; ///< TRD tracks reconstructed from TPC or ITS-TPC seeds + gsl::span mTracksTPCTRD; ///< TRD tracks reconstructed from TPC or TPC seeds + gsl::span mTrackTriggerRecordsITSTPCTRD; ///< TRD tracks reconstructed from TPC or ITS-TPC seeds + gsl::span mTrackTriggerRecordsTPCTRD; ///< TRD tracks reconstructed from TPC or TPC seeds + gsl::span mTrackletsRaw; ///< array of raw tracklets needed for TRD refit + gsl::span mTrackletsCalib; ///< array of calibrated tracklets needed for TRD refit + + // quantities necessary for pile-up correction + std::vector mTriggeredBCFT0; ///< array with the FT0 trigger times + int mCurrentTriggerRecord; + uint32_t mFirstOrbit; + int mCurrentTrackId; // corrections from ccdb, some need to be loaded only once hence an init flag o2::trd::LocalGainFactor mLocalGain; ///< local gain factors from krypton calibration diff --git a/Detectors/TRD/qc/src/Tracking.cxx b/Detectors/TRD/qc/src/Tracking.cxx index da2d05794e2d8..35f0734498a40 100644 --- a/Detectors/TRD/qc/src/Tracking.cxx +++ b/Detectors/TRD/qc/src/Tracking.cxx @@ -37,15 +37,23 @@ void Tracking::setInput(const o2::globaltracking::RecoContainer& input) mTracksTPCTRD = input.getTPCTRDTracks(); mTrackletsRaw = input.getTRDTracklets(); mTrackletsCalib = input.getTRDCalibratedTracklets(); + mTrackTriggerRecordsITSTPCTRD = input.getITSTPCTRDTriggers(); + mTrackTriggerRecordsTPCTRD = input.getTPCTRDTriggers(); } void Tracking::run() { + mCurrentTriggerRecord = 0; + mCurrentTrackId = 0; for (const auto& trkTrd : mTracksTPCTRD) { checkTrack(trkTrd, true); + mCurrentTrackId++; } + mCurrentTriggerRecord = 0; + mCurrentTrackId = 0; for (const auto& trkTrd : mTracksITSTPCTRD) { checkTrack(trkTrd, false); + mCurrentTrackId++; } } @@ -65,6 +73,59 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD) qcStruct.dEdxTotTPC = isTPCTRD ? mTracksTPC[id].getdEdx().dEdxTotTPC : mTracksTPC[mTracksITSTPC[id].getRefTPC()].getdEdx().dEdxTotTPC; } + // find corresponding track trigger record to get track timing + int triggeredBC = 0; + for (; mCurrentTriggerRecord < (isTPCTRD ? mTrackTriggerRecordsTPCTRD.size() : mTrackTriggerRecordsITSTPCTRD.size()); mCurrentTriggerRecord++) { + auto& tRecord = (isTPCTRD ? mTrackTriggerRecordsTPCTRD[mCurrentTriggerRecord] : mTrackTriggerRecordsITSTPCTRD[mCurrentTriggerRecord]); + if (mCurrentTrackId >= tRecord.getFirstTrack() && mCurrentTrackId < tRecord.getFirstTrack() + tRecord.getNumberOfTracks()) { + triggeredBC = tRecord.getBCData().differenceInBC({0, mFirstOrbit}); + break; + } + } + + // Find most probable BCs and RMS for pile-up correction and error. Same BC is assumed for all tracklets + float tCorrPileUp = 0.; + float tErrPileUp2 = 0; + float maxProb = 0.f; + // The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability + float sumCorr = 0.f; + float sumCorr2 = 0.f; + float sumProb = 0.f; + for (int iBC = 0; iBC < mTriggeredBCFT0.size(); iBC++) { + int deltaBC = roundf(mTriggeredBCFT0[iBC] - triggeredBC); + if (deltaBC <= mRecoParam.getPileUpRangeBefore()) { + continue; + } + if (deltaBC >= mRecoParam.getPileUpRangeAfter()) { + break; + } + // collect the charges + std::array q0; + std::array q1; + for (int iLy = 0; iLy < NLAYER; iLy++) { + int trkltId = trkTrd.getTrackletIndex(iLy); + if (trkltId < 0) { + q0[iLy] = -1; + q1[iLy] = -1; + } else { + q0[iLy] = mTrackletsRaw[trkltId].getQ0(); + q1[iLy] = mTrackletsRaw[trkltId].getQ1(); + } + } + // get pile-up probability + float probBC = mRecoParam.getPileUpProbTrack(deltaBC, q0, q1); + sumCorr += probBC * deltaBC; + sumCorr2 += probBC * deltaBC * deltaBC; + sumProb += probBC; + if (probBC > maxProb) { + maxProb = probBC; + tCorrPileUp = -deltaBC; + } + } + if (sumProb > 1e-6) { + tErrPileUp2 = sumCorr2 / sumProb - 2 * tCorrPileUp * sumCorr / sumProb + tCorrPileUp * tCorrPileUp; + } + for (int iLayer = 0; iLayer < NLAYER; ++iLayer) { int trkltId = trkTrd.getTrackletIndex(iLayer); if (trkltId < 0) { @@ -88,14 +149,24 @@ 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))) { 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; + + 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); qcStruct.trackProp[iLayer] = trk; diff --git a/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h b/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingQCSpec.h index b6a0cb9b78f2f..fc86d197f4df5 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,22 @@ 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)) { + triggeredBCFT0.push_back(f0rec.getInteractionRecord().differenceInBC({0, firstOrbit})); + } + } + } + mQC.setTriggeredBCFT0(triggeredBCFT0); mQC.run(); pc.outputs().snapshot(Output{"TRD", "TRACKINGQC", 0}, mQC.getTrackQC()); } @@ -94,6 +112,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 +152,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..90b60389799d7 100644 --- a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx @@ -415,6 +415,25 @@ void TRDGlobalTracking::run(ProcessingContext& pc) } LOGF(info, "%i tracks are loaded into the TRD tracker. Out of those %i ITS-TPC tracks and %i TPC tracks", nTracksLoadedITSTPC + nTracksLoadedTPC, nTracksLoadedITSTPC, nTracksLoadedTPC); + // Load the FT0 triggered BCs if this is requested + + if (mTrkMask[GTrackID::FT0]) { // pile-up tagging was requested + auto ft0recPoints = inputTracks.getFT0RecPoints(); + uint32_t firstOrbit = 0; + for (size_t ft0id = 0; ft0id < ft0recPoints.size(); ft0id++) { + const auto& f0rec = ft0recPoints[ft0id]; + if (ft0id == 0) { + firstOrbit = f0rec.getInteractionRecord().orbit; + } + if (o2::ft0::InteractionTag::Instance().isSelected(f0rec)) { + uint32_t currentOrbit = f0rec.getInteractionRecord().orbit; + mTriggeredBCFT0.push_back(f0rec.getInteractionRecord().bc + (currentOrbit - firstOrbit) * o2::constants::lhc::LHCMaxBunches); + } + } + } + + mTracker->SetFT0TriggeredBC(mTriggeredBCFT0.data(), mTriggeredBCFT0.size()); + // start the tracking // mTracker->DumpTracks(); mChainTracking->DoTRDGPUTracking(mTracker); @@ -797,6 +816,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 trkParam->resetCovariance(100); @@ -820,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()); @@ -827,9 +887,18 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, tiltCorrUp = 0.f; } - std::array trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp, zPosCorrUp}; + // conversion from slope in pad per time bin to slope in cm per BC = tracklets[trkltIdx].getSlopeFloat() * padWidth / BCperTimeBin + float slopeFactor = mTrackletsRaw[trkltId].getSlopeFloat() * pad->getWidthIPad() / 4.f; + float yCorrPileUp = tCorrPileUp * slopeFactor; + float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor; + + 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); if (!trkParam->update(trkltPosUp, trkltCovUp)) { diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx index f7adc2401df79..dfe524c38ac6f 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx @@ -24,14 +24,21 @@ 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; } #endif @@ -55,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 a0a8e71143d94..d561349a37a43 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h @@ -39,45 +39,184 @@ 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; /// Get tracklet z correction coefficient for track-eta based corraction GPUd() float getZCorrCoeffNRC() const { return mZCorrCoefNRC; } + /// Get BC intervals for pile-up + GPUd() int getPileUpRangeBefore() const { return mPileUpRangeBefore; } + GPUd() int getPileUpRangeAfter() const { return mPileUpRangeAfter; } + private: // tracklet error parameterization depends on the magnetic field float mLorentzAngle{0.f}; // 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 - ClassDefNV(GPUTRDRecoParam, 3); + // pile-up prob parametrization, depending on charges + // default parametrization, all tracklets + int mPileUpRangeBefore{-130}; ///< maximal number of BC for which pile-up from previous collision has an influence + int mPileUpMaxProb{0}; ///< number of BC with respect to triggered BC for the event with maximal probability + int mPileUpRangeAfter{70}; ///< maximal number of BC for which pile-up from next collision has an influence + // tracklets with Q0!=0 and Q1!=0 + int mPileUpRangeBefore11{-130}; ///< maximal number of BC for which pile-up from previous collision has an influence + int mPileUpMaxProb11{0}; ///< number of BC with respect to triggered BC for the event with maximal probability + int mPileUpRangeAfter11{30}; ///< maximal number of BC for which pile-up from next collision has an influence + // tracklets with Q0=0 and Q1!=0 + int mPileUpRangeBefore01{-80}; ///< maximal number of BC for which pile-up from previous collision has an influence + int mPileUpMaxProb01{30}; ///< number of BC with respect to triggered BC for the event with maximal probability + int mPileUpRangeAfter01{70}; ///< maximal number of BC for which pile-up from next collision has an influence + // tracklets with Q0!=0 and Q1=0 + int mPileUpRangeBefore10{-130}; ///< maximal number of BC for which pile-up from previous collision has an influence + int mPileUpMaxProb10{-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 + 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); }; +/// 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 + // 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/Definitions/GPUSettingsList.h b/GPU/GPUTracking/Definitions/GPUSettingsList.h index 606deb44d9528..78fb5f6d6971c 100644 --- a/GPU/GPUTracking/Definitions/GPUSettingsList.h +++ b/GPU/GPUTracking/Definitions/GPUSettingsList.h @@ -182,14 +182,16 @@ 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") 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") +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 80098ff151ebe..324acc9451f36 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -93,7 +93,7 @@ void* GPUTRDTracker_t::SetPointersTracks(void* base) } template -GPUTRDTracker_t::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mGenerateSpacePoints(false), mProcessPerTimeFrame(false), mNAngleHistogramBins(25), mAngleHistogramRange(50), mMemoryPermanent(-1), mMemoryTracklets(-1), mMemoryTracks(-1), mNMaxCollisions(0), mNMaxTracks(0), mNMaxSpacePoints(0), mTracks(nullptr), mTrackAttribs(nullptr), mNCandidates(1), mNTracks(0), mNEvents(0), mMaxBackendThreads(100), mTrackletIndexArray(nullptr), mHypothesis(nullptr), mCandidates(nullptr), mSpacePoints(nullptr), mGeo(nullptr), mRecoParam(nullptr), mDebugOutput(false), mMaxEta(0.84f), mRoadZ(18.f), mTPCVdrift(2.58f), mTPCTDriftOffset(0.f), mDebug(new GPUTRDTrackerDebug()) +GPUTRDTracker_t::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mGenerateSpacePoints(false), mProcessPerTimeFrame(false), mNAngleHistogramBins(25), mAngleHistogramRange(50), mMemoryPermanent(-1), mMemoryTracklets(-1), mMemoryTracks(-1), mNMaxCollisions(0), mNMaxTracks(0), mNMaxSpacePoints(0), mTracks(nullptr), mTrackAttribs(nullptr), mNCandidates(1), mNTracks(0), mNEvents(0), mMaxBackendThreads(100), mTrackletIndexArray(nullptr), mFT0TriggeredBC(nullptr), mNFT0BC(0), mHypothesis(nullptr), mCandidates(nullptr), mSpacePoints(nullptr), mGeo(nullptr), mRecoParam(nullptr), mDebugOutput(false), mMaxEta(0.84f), mRoadZ(18.f), mTPCVdrift(2.58f), mTPCTDriftOffset(0.f), mDebug(new GPUTRDTrackerDebug()) { //-------------------------------------------------------------------- // Default constructor @@ -351,6 +351,7 @@ GPUd() void GPUTRDTracker_t::DoTrackingThread(int32_t iTrk, int32_ } PROP prop(getPropagatorParam()); mTracks[iTrk].setChi2(Param().rec.trd.penaltyChi2); // TODO check if this should not be higher + auto trkStart = mTracks[iTrk]; for (int32_t iColl = 0; iColl < nCollisionIds; ++iColl) { // do track following for each collision candidate and keep best track @@ -436,6 +437,29 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK } mDebug->Reset(); t->setChi2(0.f); + + // Find compatible BC ids + int32_t nIdxBCMin = -1; + int32_t nIdxBCMax = -1; + + for (int32_t iBC = 0; iBC < mNFT0BC; iBC++) { + int32_t deltaBC = CAMath::Round(mFT0TriggeredBC[iBC] - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId] / o2::constants::lhc::LHCBunchSpacingMUS); + if (nIdxBCMin == -1 && deltaBC > mRecoParam->getPileUpRangeBefore()) { + nIdxBCMin = iBC; + } + if (deltaBC >= mRecoParam->getPileUpRangeAfter()) { + nIdxBCMax = iBC; + break; + } + if (iBC == mNFT0BC - 1) { + nIdxBCMax = iBC + 1; + if (nIdxBCMin == -1) { + // we did not find the correct BC, so we don't do any pile-up correction + nIdxBCMin = nIdxBCMax; + } + } + } + float zShiftTrk = 0.f; if (mProcessPerTimeFrame) { zShiftTrk = (mTrackAttribs[iTrk].mTime - GetConstantMem()->ioPtrs.trdTriggerTimes[collisionId]) * mTPCVdrift * mTrackAttribs[iTrk].mSide; @@ -585,22 +609,57 @@ 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, true, (tracklets[trkltIdx].GetQ0() != 0), (tracklets[trkltIdx].GetQ1() != 0)); + 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; + } + } + // 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(); - 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; + 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(); @@ -612,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); @@ -690,15 +749,52 @@ 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; } - 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, 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; + if (probBC > maxProb) { + maxProb = probBC; + yCorrPileUp = -slopeFactor * deltaBC; + } + } + if (sumProb > 1e-6f) { + 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 prop->setTrack(&trackNoUp); @@ -760,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 @@ -777,6 +873,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK if (iUpdate == 0 && mNCandidates > 1) { *t = mCandidates[2 * iUpdate + nextIdx]; } + } // end update loop if (!isOK) { @@ -939,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 @@ -947,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); @@ -955,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; } @@ -1018,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 f698e570d2158..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; } @@ -132,6 +133,11 @@ class GPUTRDTracker_t : public GPUProcessor GPUd() void SetRoadZ(float roadZ) { mRoadZ = roadZ; } GPUd() void SetTPCVdrift(float vDrift) { mTPCVdrift = vDrift; } GPUd() void SetTPCTDriftOffset(float t) { mTPCTDriftOffset = t; } + GPUd() void SetFT0TriggeredBC(int32_t* t, int32_t n) + { + mFT0TriggeredBC = t; + mNFT0BC = n; + } GPUd() bool GetIsDebugOutputOn() const { return mDebugOutput; } GPUd() float GetMaxEta() const { return mMaxEta; } @@ -170,6 +176,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 };