diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertex.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertex.h index ce470361034f5..5343d26ec5ce5 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertex.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertex.h @@ -37,6 +37,11 @@ class PrimaryVertex : public Vertex> void setIR(const InteractionRecord& ir) { mIRMin = mIRMax = ir; } bool hasUniqueIR() const { return !mIRMin.isDummy() && (mIRMin == mIRMax); } + float getTMAD() const { return mTMAD; } + void setTMAD(float v) { mTMAD = v; } + float getZMAD() const { return mZMAD; } + void setZMAD(float v) { mZMAD = v; } + #ifndef GPUCA_ALIGPUCODE void print() const; std::string asString() const; @@ -45,8 +50,10 @@ class PrimaryVertex : public Vertex> protected: InteractionRecord mIRMin{}; ///< by default not assigned! InteractionRecord mIRMax{}; ///< by default not assigned! + float mTMAD = -1; ///< MAD estimator for Tsigma + float mZMAD = -1; ///< MAD estimator for Zsigma - ClassDefNV(PrimaryVertex, 1); + ClassDefNV(PrimaryVertex, 2); }; #ifndef GPUCA_ALIGPUCODE diff --git a/Detectors/GlobalTrackingWorkflow/src/PrimaryVertexingSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/PrimaryVertexingSpec.cxx index 36e3266c59514..e70122e7a01b5 100644 --- a/Detectors/GlobalTrackingWorkflow/src/PrimaryVertexingSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/PrimaryVertexingSpec.cxx @@ -81,6 +81,7 @@ void PrimaryVertexingSpec::init(InitContext& ic) throw std::runtime_error(fmt::format("directory {} for raw data dumps does not exist", dumpDir)); } mVertexer.setPoolDumpDirectory(dumpDir); + mVertexer.setTrackSources(mTrackSrc); } void PrimaryVertexingSpec::run(ProcessingContext& pc) @@ -102,14 +103,19 @@ void PrimaryVertexingSpec::run(ProcessingContext& pc) std::vector gids; auto maxTrackTimeError = PVertexerParams::Instance().maxTimeErrorMUS; auto trackMaxX = PVertexerParams::Instance().trackMaxX; + auto minIBHits = PVertexerParams::Instance().minIBHits; auto halfROFITS = 0.5 * mITSROFrameLengthMUS + mITSROFBiasMUS; auto hw2ErrITS = 2.f / std::sqrt(12.f) * mITSROFrameLengthMUS; // conversion from half-width to error for ITS - auto creator = [maxTrackTimeError, hw2ErrITS, halfROFITS, trackMaxX, &tracks, &gids](auto& _tr, GTrackID _origID, float t0, float terr) { + auto creator = [maxTrackTimeError, hw2ErrITS, halfROFITS, trackMaxX, minIBHits, &tracks, &gids, &recoData](auto& _tr, GTrackID _origID, float t0, float terr) { if constexpr (isBarrelTrack()) { if (!_origID.includesDet(DetID::ITS) || _tr.getX() > trackMaxX) { return true; // just in case this selection was not done on RecoContainer filling level } + auto itsID = recoData.getITSContributorGID(_origID); + if (!itsID.isSourceSet() || o2::math_utils::numberOfBitsSet(recoData.getITSTrack(itsID).getPattern() & 7) < minIBHits) { + return true; + } if constexpr (isITSTrack()) { t0 += halfROFITS; // ITS time is supplied in \mus as beginning of ROF terr *= hw2ErrITS; // error is supplied as a half-ROF duration, convert to \mus @@ -129,13 +135,17 @@ void PrimaryVertexingSpec::run(ProcessingContext& pc) recoData.fillTrackMCLabels(gids, tracksMCInfo); } mVertexer.setStartIR(recoData.startIR); - std::vector ft0Data; + static std::vector ft0Data; if (mValidateWithIR) { // select BCs for validation + ft0Data.clear(); const o2::ft0::InteractionTag& ft0Params = o2::ft0::InteractionTag::Instance(); auto ft0all = recoData.getFT0RecPoints(); for (const auto& ftRP : ft0all) { if (ft0Params.isSelected(ftRP)) { - ft0Data.push_back(ftRP.getInteractionRecord()); + ft0Data.emplace_back(InteractionCandidate{ftRP.getInteractionRecord(), + float(ftRP.getInteractionRecord().differenceInBC(recoData.startIR) * o2::constants::lhc::LHCBunchSpacingMUS), + float(ftRP.getTrigger().getAmplA() + ftRP.getTrigger().getAmplC()), + GTrackID::FT0}); } } } @@ -151,9 +161,9 @@ void PrimaryVertexingSpec::run(ProcessingContext& pc) } mTimer.Stop(); - LOGP(info, "Found {} PVs, Time CPU/Real:{:.3f}/{:.3f} (DBScan: {:.4f}, Finder:{:.4f}, Rej.Debris:{:.4f}, Reattach:{:.4f}) | {} trials for {} TZ-clusters, max.trials: {}, Slowest TZ-cluster: {} ms of mult {}", + LOGP(info, "Found {} PVs, Time CPU/Real:{:.3f}/{:.3f} (DBScan: {:.4f}, Finder:{:.4f}, MADSel:{:.4f}, Rej.Debris:{:.4f}, Reattach:{:.4f}) | {} trials for {} TZ-clusters, max.trials: {}, Slowest TZ-cluster: {} ms of mult {}", vertices.size(), mTimer.CpuTime() - timeCPU0, mTimer.RealTime() - timeReal0, - mVertexer.getTimeDBScan().CpuTime(), mVertexer.getTimeVertexing().CpuTime(), mVertexer.getTimeDebris().CpuTime(), mVertexer.getTimeReAttach().CpuTime(), + mVertexer.getTimeDBScan().CpuTime(), mVertexer.getTimeVertexing().CpuTime(), mVertexer.getTimeMADSel().CpuTime(), mVertexer.getTimeDebris().CpuTime(), mVertexer.getTimeReAttach().CpuTime(), mVertexer.getTotTrials(), mVertexer.getNTZClusters(), mVertexer.getMaxTrialsPerCluster(), mVertexer.getLongestClusterTimeMS(), mVertexer.getLongestClusterMult()); } diff --git a/Detectors/Vertexing/include/DetectorsVertexing/PVertexer.h b/Detectors/Vertexing/include/DetectorsVertexing/PVertexer.h index 2c71e5d5d9df1..e007068c0f56d 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/PVertexer.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/PVertexer.h @@ -31,6 +31,7 @@ #include "DetectorsVertexing/PVertexerParams.h" #include "ReconstructionDataFormats/GlobalTrackID.h" #include "DataFormatsCalibration/MeanVertexObject.h" +#include "ITSMFTBase/DPLAlpideParam.h" #include "gsl/span" #include #include @@ -59,7 +60,11 @@ class PVertexer void init(); void end(); template - int process(const TR& tracks, const gsl::span gids, const gsl::span bcData, + int process(const TR& tracks, const gsl::span gids, const gsl::span intCand, + std::vector& vertices, std::vector& vertexTrackIDs, std::vector& v2tRefs, + const gsl::span lblTracks, std::vector& lblVtx); + template + int process(const TR& tracks, const gsl::span gids, const gsl::span intRec, std::vector& vertices, std::vector& vertexTrackIDs, std::vector& v2tRefs, const gsl::span lblTracks, std::vector& lblVtx); @@ -82,6 +87,7 @@ class PVertexer void setBz(float bz) { mBz = bz; } void setValidateWithIR(bool v) { mValidateWithIR = v; } bool getValidateWithIR() const { return mValidateWithIR; } + void setTrackSources(GTrackID::mask_t s); auto& getTracksPool() const { return mTracksPool; } auto& getTimeZClusters() const { return mTimeZClusters; } @@ -117,6 +123,7 @@ class PVertexer TStopwatch& getTimeDBScan() { return mTimeDBScan; } TStopwatch& getTimeVertexing() { return mTimeVertexing; } TStopwatch& getTimeDebris() { return mTimeDebris; } + TStopwatch& getTimeMADSel() { return mTimeMADSel; } TStopwatch& getTimeReAttach() { return mTimeReAttach; } void setPoolDumpDirectory(const std::string& d) { mPoolDumpDirectory = d; } @@ -127,7 +134,7 @@ class PVertexer static constexpr int DBS_UNDEF = -2, DBS_NOISE = -1, DBS_INCHECK = -10; SeedHistoTZ buildHistoTZ(const VertexingInput& input); - int runVertexing(gsl::span gids, const gsl::span bcData, + int runVertexing(gsl::span gids, const gsl::span intCand, std::vector& vertices, std::vector& vertexTrackIDs, std::vector& v2tRefs, gsl::span lblTracks, std::vector& lblVtx); void createMCLabels(gsl::span lblTracks, const std::vector& trackIDs, const std::vector& v2tRefs, std::vector& lblVtx); @@ -144,6 +151,9 @@ class PVertexer bool upscaleSigma(VertexSeed& vtxSeed) const; bool relateTrackToMeanVertex(o2::track::TrackParCov& trc, float vtxErr2); bool relateTrackToVertex(o2::track::TrackParCov& trc, const o2d::VertexBase& vtxSeed) const; + void applyMADSelection(std::vector& vertices, std::vector& timeSort, const std::vector& v2tRefs, const std::vector& trackIDs); + void applITSOnlyFractionCut(std::vector& vertices, std::vector& timeSort, const std::vector& v2tRefs, const std::vector& trackIDs); + void applInteractionValidation(std::vector& vertices, std::vector& timeSort, const gsl::span intCand, int minContrib); template void createTracksPool(const TR& tracks, gsl::span gids); @@ -151,7 +161,7 @@ class PVertexer int findVertices(const VertexingInput& input, std::vector& vertices, std::vector& trackIDs, std::vector& v2tRefs); void reAttach(std::vector& vertices, std::vector& timeSort, std::vector& trackIDs, std::vector& v2tRefs); - std::pair getBestIR(const PVertex& vtx, const gsl::span bcData, int& currEntry) const; + std::pair getBestIR(const PVertex& vtx, const gsl::span intCand, int& currEntry) const; int dbscan_RangeQuery(int idxs, std::vector& cand, std::vector& status); void dbscan_clusterize(); @@ -173,7 +183,7 @@ class PVertexer float mBz = 0.; ///< mag.field at beam line float mDBScanDeltaT = 0.; ///< deltaT cut for DBScan check float mDBSMaxZ2InvCorePoint = 0; ///< inverse of max sigZ^2 of the track which can be core point in the DBScan - bool mValidateWithIR = false; ///< require vertex validation with InteractionRecords (if available) + bool mValidateWithIR = false; ///< require vertex validation with InteractionCandidates (if available) o2::InteractionRecord mStartIR{0, 0}; ///< IR corresponding to the start of the TF @@ -183,6 +193,8 @@ class PVertexer // ///========== Parameters to be set externally, e.g. from CCDB ==================== + GTrackID::mask_t mTrackSrc{}; + std::vector mSrcVec{}; const PVertexerParams* mPVParams = nullptr; float mTukey2I = 0; ///< 1./[Tukey parameter]^2 static constexpr float kDefTukey = 5.0f; ///< def.value for tukey constant @@ -192,12 +204,18 @@ class PVertexer size_t mNTZClustersIni = 0; size_t mTotTrials = 0; size_t mMaxTrialPerCluster = 0; + float mMaxTDiffDebris = 0; ///< when reducing debris, don't consider vertices separated by time > this value in \mus + float mMaxTDiffDebrisExtra = 0; ///< when reducing debris, don't consider vertices separated by time > this value in \mus (optional additional cut + float mMaxTDiffDebrisFiducial = 0; + float mMaxZDiffDebrisFiducial = 0; + float mMaxMultRatDebrisFiducial = 0; long mLongestClusterTimeMS = 0; int mLongestClusterMult = 0; bool mPoolDumpProduced = false; TStopwatch mTimeDBScan; TStopwatch mTimeVertexing; TStopwatch mTimeDebris; + TStopwatch mTimeMADSel; TStopwatch mTimeReAttach; std::string mPoolDumpDirectory{}; #ifdef _PV_DEBUG_TREE_ @@ -263,6 +281,7 @@ void PVertexer::createTracksPool(const TR& tracks, gsl::spanmeanVertexExtraErrSelection * mPVParams->meanVertexExtraErrSelection; + const auto& alpParams = o2::itsmft::DPLAlpideParam::Instance(); for (uint32_t i = 0; i < ntGlo; i++) { int id = sortedTrackID[i]; @@ -275,6 +294,18 @@ void PVertexer::createTracksPool(const TR& tracks, gsl::spanaddTimeSigma2, mPVParams->addZSigma2); if (tvf.wghHisto < 0) { mTracksPool.pop_back(); // discard bad track + continue; + } + if (gids[id].getSource() == GTrackID::ITSTPC) { // if the track was adjusted to ITS ROF boundary, flag it + float bcf = tvf.timeEst.getTimeStamp() / o2::constants::lhc::LHCBunchSpacingMUS + o2::constants::lhc::LHCMaxBunches; + int bcWrtROF = int(bcf - alpParams.roFrameBiasInBC) % alpParams.roFrameLengthInBC; + if (bcWrtROF == 0) { + float dbc = bcf - (int(bcf / alpParams.roFrameBiasInBC)) * alpParams.roFrameBiasInBC; + if (std::abs(dbc) < 1e-6) { + tvf.setITSTPCAdjusted(); + LOGP(debug, "Adjusted t={} -> bcf={} dbc = {}", tvf.timeEst.getTimeStamp(), bcf, dbc); + } + } } } @@ -288,12 +319,28 @@ void PVertexer::createTracksPool(const TR& tracks, gsl::span -int PVertexer::process(const TR& tracks, const gsl::span gids, const gsl::span bcData, +int PVertexer::process(const TR& tracks, const gsl::span gids, const gsl::span intCand, + std::vector& vertices, std::vector& vertexTrackIDs, std::vector& v2tRefs, + const gsl::span lblTracks, std::vector& lblVtx) +{ + createTracksPool(tracks, gids); + return runVertexing(gids, intCand, vertices, vertexTrackIDs, v2tRefs, lblTracks, lblVtx); +} + +//___________________________________________________________________ +template +int PVertexer::process(const TR& tracks, const gsl::span gids, const gsl::span intRec, std::vector& vertices, std::vector& vertexTrackIDs, std::vector& v2tRefs, const gsl::span lblTracks, std::vector& lblVtx) { createTracksPool(tracks, gids); - return runVertexing(gids, bcData, vertices, vertexTrackIDs, v2tRefs, lblTracks, lblVtx); + static std::vector intCand; + intCand.clear(); + intCand.reserve(intRec.size()); + for (const auto& ir : intRec) { + intCand.emplace_back(InteractionCandidate{ir, float(ir.differenceInBC(mStartIR) * o2::constants::lhc::LHCBunchSpacingMUS), 0, 0}); + } + return runVertexing(gids, intCand, vertices, vertexTrackIDs, v2tRefs, lblTracks, lblVtx); } //___________________________________________________________________ diff --git a/Detectors/Vertexing/include/DetectorsVertexing/PVertexerHelpers.h b/Detectors/Vertexing/include/DetectorsVertexing/PVertexerHelpers.h index 64274123e651e..ec11c6e289c71 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/PVertexerHelpers.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/PVertexerHelpers.h @@ -105,6 +105,8 @@ struct TrackVF { enum { kUsed, kNoVtx = -1, kDiscarded = kNoVtx - 1 }; + enum { kITSTPCAdjust = 0x1, + kDummyHBin = 0xffff }; float x; ///< reference X float y; ///< Y at X float z; ///< Z at X @@ -119,11 +121,14 @@ struct TrackVF { TimeEst timeEst; float wgh = 0.; ///< track weight wrt current vertex seed float wghHisto = 0.; // weight based on track errors, used for histogramming - int entry; ///< track entry in the input vector - int32_t bin = -1; // seeds histo bin - GTrackID gid{}; + int entry; ///< track entry in the input vector int vtxID = kNoVtx; ///< assigned vertex + GTrackID gid{}; + uint16_t bin = kDummyHBin; // seeds histo bin + uint16_t flags = 0; // + void setITSTPCAdjusted() { flags &= kITSTPCAdjust; } + bool isITSTPCAdjusted() const { return flags & kITSTPCAdjust; } bool canAssign() const { return wgh > 0. && vtxID == kNoVtx; } bool canUse() const { return vtxID == kNoVtx; } bool canUse(float zmin, float zmax) const @@ -215,7 +220,7 @@ struct TrackVF { struct SeedHistoTZ : public o2::dataformats::FlatHisto2D_f { using o2::dataformats::FlatHisto2D::FlatHisto2D; - int fillAndFlagBin(float x, float y, float w) + uint16_t fillAndFlagBin(float x, float y, float w) { uint32_t bin = getBin(x, y); if (isValidBin(bin)) { @@ -224,9 +229,9 @@ struct SeedHistoTZ : public o2::dataformats::FlatHisto2D_f { } fillBin(bin, w); nEntries++; - return bin; + return uint16_t(bin); } - return -1; + return 0xffff; } void clear() @@ -268,6 +273,13 @@ struct TrackVFDump { ClassDefNV(TrackVFDump, 1); }; +struct InteractionCandidate : public o2::InteractionRecord { + float time = 0; + float amplitude = 0; + uint32_t flag = 0; // origin, etc. + InteractionCandidate() = default; +}; + } // namespace vertexing } // namespace o2 diff --git a/Detectors/Vertexing/include/DetectorsVertexing/PVertexerParams.h b/Detectors/Vertexing/include/DetectorsVertexing/PVertexerParams.h index 4c926456e3ad0..5e7a25924f55d 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/PVertexerParams.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/PVertexerParams.h @@ -45,6 +45,7 @@ struct PVertexerParams : public o2::conf::ConfigurableParamHelper this value in \mus + + float maxTDiffDebris = 7.0; ///< when reducing debris, don't consider vertices separated by time > this value in \mus if >0, if <0: mult factor to ITS ROF float maxZDiffDebris = 1.0; ///< don't consider vertices separated by Z > this value in cm float maxMultRatDebris = 0.05; ///< don't consider vertices with multiplicity ratio above this float maxChi2TZDebris = 2000.; ///< don't consider vertices with mutual chi2 exceeding this (for pp should be ~10) float addTimeSigma2Debris = 0.05 * 0.05; ///< increment time error^2 by this amount when calculating vertex-to-vertex chi2 float addZSigma2Debris = 0.005 * 0.005; ///< increment z error^2 by this amount when calculating vertex-to-vertex chi2 + // extra debris reduction cut, ignored if maxTDiffDebrisExtra == 0. The maxTDiffDebrisExtra must not exceed maxTDiffDebris + float maxTDiffDebrisExtra = 0; ///< when reducing debris, don't consider vertices separated by time > this value in \mus if >0, if <0: mult factor to ITS ROF + float maxZDiffDebrisExtra = 1.0; ///< don't consider vertices separated by Z > this value in cm + float maxMultRatDebrisExtra = 0.05; ///< don't consider vertices with multiplicity ratio above this + float maxChi2TZDebrisExtra = 2000.; ///< don't consider vertices with mutual chi2 exceeding this (for pp should be ~10) + float addTimeSigma2DebrisExtra = 0.05 * 0.05; ///< increment time error^2 by this amount when calculating vertex-to-vertex chi2 + float addZSigma2DebrisExtra = 0.005 * 0.005; ///< increment z error^2 by this amount when calculating vertex-to-vertex chi2 + + float maxITSOnlyFraction = 1.0; ///< max ITS-only tracks fraction to accept, recommended value for PbPb = 0.85 + float minITSOnlyFraction = 0.0; ///< min ITS-only tracks fraction to accept + // + // MAD clean up (median abs. deviation in T and Z) + float maxTMAD = -1.; ///< max accepted tMAD, not tMAD cleanup if negative + float minTMAD = -1.; ///< min accepted tMAD + float maxZMAD = -1.; ///< max accepted zMAD, not zMAD cleanup if negative + float minZMAD = -1.; ///< min accepted zMAD + // validation with externally provided InteractionRecords (e.g. from FT0) + int minNContributorsForIRcutIni = -1; ///< min multiplicity to reject before cleanup if no matching interaction found (if >= 0) int minNContributorsForIRcut = 4; ///< do not apply IR cut to vertices below IR tagging efficiency threshold float maxTError = 0.2; ///< use min of vertex time error or this for nsigma evaluation float minTError = 0.003; ///< don't use error smaller than that (~BC/2/minNContributorsForFT0cut) diff --git a/Detectors/Vertexing/src/PVertexer.cxx b/Detectors/Vertexing/src/PVertexer.cxx index b933cd8867b10..16d80aaef91b7 100644 --- a/Detectors/Vertexing/src/PVertexer.cxx +++ b/Detectors/Vertexing/src/PVertexer.cxx @@ -17,6 +17,7 @@ #include "DetectorsBase/Propagator.h" #include "Math/SMatrix.h" #include "Math/SVector.h" +#include "MathUtils/fit.h" #include #include "CommonUtils/StringUtils.h" #include @@ -34,7 +35,7 @@ PVertexer::PVertexer() } //___________________________________________________________________ -int PVertexer::runVertexing(const gsl::span gids, const gsl::span bcData, +int PVertexer::runVertexing(const gsl::span gids, const gsl::span intCand, std::vector& vertices, std::vector& vertexTrackIDs, std::vector& v2tRefs, gsl::span lblTracks, std::vector& lblVtx) { @@ -82,16 +83,26 @@ int PVertexer::runVertexing(const gsl::span gids, const gsl: } #endif - mTimeDebris.Start(); + if (mValidateWithIR && mPVParams->minNContributorsForIRcutIni >= 0) { + applInteractionValidation(verticesLoc, vtTimeSortID, intCand, mPVParams->minNContributorsForIRcutIni); + } + + if (mPVParams->minITSOnlyFraction > 0.f || mPVParams->maxITSOnlyFraction < 1.0f) { + applITSOnlyFractionCut(verticesLoc, vtTimeSortID, v2tRefsLoc, trackIDs); + } + if (mPVParams->applyDebrisReduction) { + mTimeDebris.Start(); reduceDebris(verticesLoc, vtTimeSortID, lblVtxLoc); + mTimeDebris.Stop(); } - mTimeDebris.Stop(); - mTimeReAttach.Start(); + if (mPVParams->applyReattachment) { + mTimeReAttach.Start(); reAttach(verticesLoc, vtTimeSortID, trackIDs, v2tRefsLoc); + mTimeReAttach.Stop(); } - mTimeReAttach.Stop(); + if (lblTracks.size()) { // at this stage labels are needed just for the debug output createMCLabels(lblTracks, trackIDs, v2tRefsLoc, lblVtxLoc); } @@ -108,30 +119,31 @@ int PVertexer::runVertexing(const gsl::span gids, const gsl: vertexTrackIDs.reserve(trackIDs.size()); int trCopied = 0, count = 0, vtimeID = 0; - for (auto i : vtTimeSortID) { + for (auto& i : vtTimeSortID) { if (i < 0) { continue; // vertex was suppressed } auto& vtx = verticesLoc[i]; - - bool irSet = setCompatibleIR(vtx); - if (!irSet) { - continue; + if (!setCompatibleIR(vtx)) { + i = -1; } - // do we need to validate by Int. records ? - if (mValidateWithIR) { - auto bestMatch = getBestIR(vtx, bcData, vtimeID); - if (bestMatch.first >= 0) { - vtx.setFlags(PVertex::TimeValidated); - if (bestMatch.second == 1) { - vtx.setIR(bcData[bestMatch.first]); - } - LOG(debug) << "Validated with t0 " << bestMatch.first << " with " << bestMatch.second << " candidates"; - } else if (vtx.getNContributors() >= mPVParams->minNContributorsForIRcut) { - LOG(debug) << "Discarding " << vtx; - continue; // reject - } + } + // do we need to validate by Int. records ? + if (mValidateWithIR) { + applInteractionValidation(verticesLoc, vtTimeSortID, intCand, mPVParams->minNContributorsForIRcut); + } + // + if (mPVParams->maxTMAD > 0 || mPVParams->maxZMAD > 0) { + mTimeMADSel.Start(); + applyMADSelection(verticesLoc, vtTimeSortID, v2tRefsLoc, trackIDs); + mTimeMADSel.Stop(); + } + // + for (auto i : vtTimeSortID) { + if (i < 0) { + continue; // vertex was suppressed } + auto& vtx = verticesLoc[i]; vertices.push_back(vtx); if (lblVtxLoc.size()) { lblVtx.push_back(lblVtxLoc[i]); @@ -491,7 +503,7 @@ void PVertexer::reAttach(std::vector& vertices, std::vector& timeS if (dt > tRange) { // all following vertices will be have higher time than this track, stop checking break; } - bool useTime = trc.gid.getSource() != GTrackID::ITS && mPVParams->useTimeInChi2; // TODO Check if it is ok to not account time error for ITS tracks + bool useTime = (trc.gid.getSource() != GTrackID::ITS) && mPVParams->useTimeInChi2; // TODO Check if it is ok to not account time error for ITS tracks float wgh = 1.f - mTukey2I * trc.evalChi2ToVertex(vertices[vtvec[ivt].first], useTime); if (wgh > trc.wgh) { trc.wgh = wgh; @@ -501,7 +513,7 @@ void PVertexer::reAttach(std::vector& vertices, std::vector& timeS if (trc.vtxID != TrackVF::kNoVtx) { mTimeZClusters[trc.vtxID].trackIDs.push_back(itr); trc.vtxID = TrackVF::kNoVtx; // to enable for using in the fit - trc.bin = -1; + trc.bin = TrackVF::kDummyHBin; } } // refit vertices with reattached tracks @@ -533,6 +545,134 @@ void PVertexer::reAttach(std::vector& vertices, std::vector& timeS }); } +//___________________________________________________________________ +void PVertexer::applyMADSelection(std::vector& vertices, std::vector& timeSort, const std::vector& v2tRefs, const std::vector& trackIDs) +{ + int nv = timeSort.size(), nkill = 0; + std::vector dvec; + for (int ivt = 0; ivt < nv; ivt++) { + int iv = timeSort[ivt]; + if (iv < 0) { // already rejected? + continue; + } + auto& pv = vertices[iv]; + const auto trefs = v2tRefs[iv]; + int idMin = trefs.getFirstEntry(), idMax = idMin + trefs.getEntries(); + dvec.clear(); + dvec.reserve(pv.getNContributors()); + if (mPVParams->maxTMAD > 0) { + int nAdjusted = 0; + float tVtx = pv.getTimeStamp().getTimeStamp(); + for (int i = idMin; i < idMax; i++) { + const auto& trc = mTracksPool[trackIDs[i]]; + if (trc.gid.getSource() == GTrackID::ITS || // ITS should not be accounted in the mad MAD as it has no real errors + (trc.gid.getSource() == GTrackID::ITSTPC && trc.isITSTPCAdjusted() && (++nAdjusted) > 1)) { // adjusted tracks are counted only once since they are all the same + continue; + } else { + dvec.push_back(trc.timeEst.getTimeStamp() - tVtx); + } + } + float tmad = o2::math_utils::MAD2Sigma(dvec.size(), dvec.data()); + dvec.clear(); + if (tmad > mPVParams->maxTMAD || tmad < mPVParams->minTMAD) { + timeSort[ivt] = -1; // disable vertex + LOGP(debug, "Killing vertex {} with TMAD {}, {} of {} killed", iv, tmad, ++nkill, nv); + continue; + } + pv.setTMAD(tmad); + } + if (mPVParams->maxZMAD > 0) { + float zVtx = pv.getZ(); + for (int i = idMin; i < idMax; i++) { + const auto& trc = mTracksPool[trackIDs[i]]; + dvec.push_back(trc.z - zVtx); + } + float zmad = o2::math_utils::MAD2Sigma(dvec.size(), dvec.data()); + if (zmad > mPVParams->maxZMAD || zmad < mPVParams->minZMAD) { + timeSort[ivt] = -1; // disable vertex + LOGP(debug, "Killing vertex {} with ZMAD {}, {} of {} killed", iv, zmad, ++nkill, nv); + continue; + } + pv.setZMAD(zmad); + } + } +} + +//___________________________________________________________________ +void PVertexer::applITSOnlyFractionCut(std::vector& vertices, std::vector& timeSort, const std::vector& v2tRefs, const std::vector& trackIDs) +{ + int nv = timeSort.size(); + for (int ivt = 0; ivt < nv; ivt++) { + int iv = timeSort[ivt]; + if (iv < 0) { // already rejected? + continue; + } + const auto& pv = vertices[iv]; + float tVtx = pv.getTimeStamp().getTimeStamp(); + const auto trefs = v2tRefs[iv]; + int idMin = trefs.getFirstEntry(), idMax = idMin + trefs.getEntries(); + int nITS = 0; + for (int i = idMin; i < idMax; i++) { + const auto& trc = mTracksPool[trackIDs[i]]; + if (trc.gid.getSource() == GTrackID::Source::ITS) { // ITS should not be accounted in the mad MAD as it has no real errors + nITS++; + } + } + float frac = nITS / float(trefs.getEntries()); + if (frac > mPVParams->maxITSOnlyFraction || frac < mPVParams->minITSOnlyFraction) { + timeSort[ivt] = -1; // disable vertex + } + } +} + +//___________________________________________________________________ +void PVertexer::applInteractionValidation(std::vector& vertices, std::vector& timeSort, const gsl::span intCand, int minContrib) +{ + // get max error + int nv = timeSort.size(), nkill = 0; + float maxErr = std::max(mPVParams->minTError, mPVParams->nSigmaTimeCut * mPVParams->maxTError) + mPVParams->timeMarginVertexTime; + + int intCandMin = 0, nInt = intCand.size(); + for (int ivt = 0; ivt < nv; ivt++) { + int iv = timeSort[ivt]; + if (iv < 0) { // already rejected? + continue; + } + auto& pv = vertices[iv]; + while (intCandMin < nInt && intCand[intCandMin].time < (pv.getTimeStamp().getTimeStamp() + mPVParams->timeBiasMS - maxErr)) { + intCandMin++; // skip all times which have no chance to be matched + } + int i = intCandMin, nCompatible = 0, best = -1; + float bestDf = 1e12; + while (i < nInt) { + float tv = pv.getTimeStamp().getTimeStamp() + mPVParams->timeBiasMS; + float tvE = std::max(mPVParams->minTError, mPVParams->nSigmaTimeCut * std::min(mPVParams->maxTError, pv.getTimeStamp().getTimeStampError())) + mPVParams->timeMarginVertexTime; + if (intCand[i].time > (tv + tvE)) { + break; + } + if (intCand[i].time < (tv - tvE)) { + i++; + continue; + } + nCompatible++; + auto dfa = std::abs(tv - intCand[intCandMin].time); + if (dfa < bestDf) { + bestDf = dfa; + best = i; + } + i++; + } + if (best >= 0) { + pv.setFlags(PVertex::TimeValidated); + if (nCompatible == 1) { + pv.setIR(intCand[best]); + } + } else if (pv.getNContributors() >= minContrib) { + timeSort[ivt] = -1; // discard + } + } +} + //___________________________________________________________________ void PVertexer::reduceDebris(std::vector& vertices, std::vector& timeSort, const std::vector& lblVtx) { @@ -542,7 +682,7 @@ void PVertexer::reduceDebris(std::vector& vertices, std::vector& t std::vector multSort(nv); // sort time indices in multiplicity std::iota(multSort.begin(), multSort.end(), 0); std::sort(multSort.begin(), multSort.end(), [&timeSort, vertices](int i, int j) { - return vertices[timeSort[i]].getNContributors() > vertices[timeSort[j]].getNContributors(); + return timeSort[i] < 0 ? false : (timeSort[j] < 0 ? true : vertices[timeSort[i]].getNContributors() > vertices[timeSort[j]].getNContributors()); }); // suppress vertex pointed by j if needed, if time difference between i and j is too large, return true to stop checking @@ -550,37 +690,29 @@ void PVertexer::reduceDebris(std::vector& vertices, std::vector& t auto checkPair = [&vertices, &timeSort, &lblVtx, this](int i, int j) { auto &vtI = vertices[timeSort[i]], &vtJ = vertices[timeSort[j]]; auto tDiff = std::abs(vtI.getTimeStamp().getTimeStamp() - vtJ.getTimeStamp().getTimeStamp()); - if (tDiff > this->mPVParams->maxTDiffDebris) { + + if (tDiff > this->mMaxTDiffDebrisFiducial) { return true; // don't continue checking other neighbours in time } if (vtI.getNContributors() < vtJ.getNContributors()) { return false; // comparison goes from higher to lower mult vtx } - bool rej = true; + bool rej = false; float zDiff = std::abs(vtI.getZ() - vtJ.getZ()); - if (zDiff > this->mPVParams->maxZDiffDebris) { // cannot be reduced as too far in Z -#ifndef _PV_DEBUG_TREE_ - return false; -#endif - rej = false; - } - float multRat = float(vtJ.getNContributors()) / float(vtI.getNContributors()); - if (multRat > this->mPVParams->maxMultRatDebris) { -#ifndef _PV_DEBUG_TREE_ - return false; -#endif - rej = false; - } - float tiE = vtI.getTimeStamp().getTimeStampError(), tjE = vtJ.getTimeStamp().getTimeStampError(); - float chi2z = zDiff * zDiff / (vtI.getSigmaZ2() + vtJ.getSigmaZ2() + this->mPVParams->addZSigma2Debris); - float chi2t = tDiff * tDiff / (tiE * tiE + tjE * tjE + this->mPVParams->addTimeSigma2Debris); - if (chi2z + chi2t > this->mPVParams->maxChi2TZDebris) { -#ifndef _PV_DEBUG_TREE_ - return false; -#endif - rej = false; + if (zDiff < this->mMaxZDiffDebrisFiducial && float(vtJ.getNContributors() < float(vtI.getNContributors()) * this->mMaxMultRatDebrisFiducial)) { + if (this->mMaxTDiffDebrisExtra <= 0 || // if no extra cut requested, no need to recheck the differences + (zDiff < this->mPVParams->maxZDiffDebris && tDiff < this->mMaxTDiffDebris && float(vtJ.getNContributors()) < float(vtI.getNContributors()) * this->mPVParams->maxMultRatDebris)) { + float chi2z = zDiff * zDiff / (vtI.getSigmaZ2() + vtJ.getSigmaZ2() + this->mPVParams->addZSigma2Debris); + float chi2t = tDiff * tDiff / (vtI.getTimeStamp().getTimeStampError2() + vtJ.getTimeStamp().getTimeStampError2() + this->mPVParams->addTimeSigma2Debris); + rej = (chi2z + chi2t) < this->mPVParams->maxChi2TZDebris; + } + if (!rej && this->mMaxTDiffDebrisExtra > 0 && + (zDiff < this->mPVParams->maxZDiffDebrisExtra && tDiff < this->mMaxTDiffDebrisExtra && float(vtJ.getNContributors()) < float(vtI.getNContributors()) * this->mPVParams->maxMultRatDebrisExtra)) { + float chi2z = zDiff * zDiff / (vtI.getSigmaZ2() + vtJ.getSigmaZ2() + this->mPVParams->addZSigma2DebrisExtra); + float chi2t = tDiff * tDiff / (vtI.getTimeStamp().getTimeStampError2() + vtJ.getTimeStamp().getTimeStampError2() + this->mPVParams->addTimeSigma2DebrisExtra); + rej = (chi2z + chi2t) < this->mPVParams->maxChi2TZDebrisExtra; + } } - // all veto cuts passed, declare as fake! #ifdef _PV_DEBUG_TREE_ o2::MCEventLabel dummyLbl; this->mDebugDumpPVComp.emplace_back(PVtxCompDump{vtI, vtJ, chi2z, chi2t, rej}); @@ -677,6 +809,12 @@ void PVertexer::init() mDBScanDeltaT = mPVParams->dbscanDeltaT > 0.f ? mPVParams->dbscanDeltaT : mITSROFrameLengthMUS - mPVParams->dbscanDeltaT; mDBSMaxZ2InvCorePoint = mPVParams->dbscanMaxSigZCorPoint > 0 ? 1. / (mPVParams->dbscanMaxSigZCorPoint * mPVParams->dbscanMaxSigZCorPoint) : 1e6; + mMaxTDiffDebris = mPVParams->maxTDiffDebris < 0 ? mITSROFrameLengthMUS * (-mPVParams->maxTDiffDebris) : mPVParams->maxTDiffDebris; + mMaxTDiffDebrisExtra = mPVParams->maxTDiffDebrisExtra == 0 ? -1 : (mPVParams->maxTDiffDebrisExtra < 0 ? mITSROFrameLengthMUS * (-mPVParams->maxTDiffDebrisExtra) : mPVParams->maxTDiffDebrisExtra); + mMaxTDiffDebrisFiducial = mMaxTDiffDebrisExtra > 0 ? std::max(mMaxTDiffDebrisExtra, mMaxTDiffDebris) : mMaxTDiffDebris; + mMaxZDiffDebrisFiducial = mMaxTDiffDebrisExtra > 0 ? std::max(mPVParams->maxZDiffDebrisExtra, mPVParams->maxZDiffDebris) : mPVParams->maxZDiffDebris; + mMaxMultRatDebrisFiducial = mMaxTDiffDebrisExtra > 0 ? std::max(mPVParams->maxMultRatDebrisExtra, mPVParams->maxMultRatDebris) : mPVParams->maxMultRatDebris; + #ifdef _PV_DEBUG_TREE_ mDebugDumpFile = std::make_unique("pvtxDebug.root", "recreate"); @@ -737,7 +875,7 @@ void PVertexer::finalizeVertex(const VertexingInput& input, const PVertex& vtx, if (trc.canAssign()) { trackIDs.push_back(i); trc.vtxID = lastID; - if (trc.bin >= 0 && histo) { + if (trc.bin != TrackVF::kDummyHBin && histo) { histo->setBinContent(trc.bin, -1.f); // discard used bin } } @@ -981,22 +1119,22 @@ void PVertexer::dbscan_clusterize() } //___________________________________________________________________ -std::pair PVertexer::getBestIR(const PVertex& vtx, const gsl::span bcData, int& currEntry) const +std::pair PVertexer::getBestIR(const PVertex& vtx, const gsl::span intCand, int& currEntry) const { // select best matching interaction record - int best = -1, n = bcData.size(); - while (currEntry < n && bcData[currEntry] < vtx.getIRMin()) { + int best = -1, n = intCand.size(); + while (currEntry < n && intCand[currEntry] < vtx.getIRMin()) { currEntry++; // skip all times which have no chance to be matched } int i = currEntry, nCompatible = 0; float bestDf = 1e12; auto tVtxNS = (vtx.getTimeStamp().getTimeStamp() + mPVParams->timeBiasMS) * 1e3; // time in ns while (i < n) { - if (bcData[i] > vtx.getIRMax()) { + if (intCand[i] > vtx.getIRMax()) { break; } nCompatible++; - auto dfa = std::abs(bcData[i].differenceInBCNS(mStartIR) - tVtxNS); + auto dfa = std::abs(intCand[i].differenceInBCNS(mStartIR) - tVtxNS); if (dfa <= bestDf) { bestDf = dfa; best = i; @@ -1033,6 +1171,14 @@ SeedHistoTZ PVertexer::buildHistoTZ(const VertexingInput& input) float dz = hZMax - hZMin, dt = hTMax - hTMin; int nbz = 1 + int((dz) / mPVParams->histoBinZSize), nbt = 1 + int((dt) / mPVParams->histoBinTSize); + if (nbz * nbt > 0x7fff) { + float factor = float(nbz * nbt / 0x7fff); + if (nbt > nbz) { + nbt = (0x7fff - 1) / nbz; + } else { + nbz = (0x7fff - 1) / nbt; + } + } float dzh = 0.5f * (nbz * mPVParams->histoBinZSize - dz), dth = 0.5f * (nbt * mPVParams->histoBinTSize - dt); SeedHistoTZ seedHistoTZ(nbt, hTMin - dth, hTMax + dth, nbz, hZMin - dzh, hZMax + dzh); @@ -1202,7 +1348,7 @@ void PVertexer::dumpPool() int PVertexer::processFromExternalPool(const std::vector& pool, std::vector& vertices, std::vector& vertexTrackIDs, std::vector& v2tRefs) { // dummy inputs - std::vector bcData; + std::vector intCand; std::vector lblTracks; std::vector lblVtx; @@ -1213,5 +1359,17 @@ int PVertexer::processFromExternalPool(const std::vector& pool, std::ve mTracksPool.push_back(tr); gids.push_back(tr.gid); } - return runVertexing(gids, bcData, vertices, vertexTrackIDs, v2tRefs, lblTracks, lblVtx); + return runVertexing(gids, intCand, vertices, vertexTrackIDs, v2tRefs, lblTracks, lblVtx); +} + +//______________________________________________ +void PVertexer::setTrackSources(GTrackID::mask_t s) +{ + mTrackSrc = s; + // fill vector of sources to consider + for (int is = 0; is < GTrackID::NSources; is++) { + if (mTrackSrc[is]) { + mSrcVec.push_back(is); + } + } }