Project
Loading...
Searching...
No Matches
PVertexer.h
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
15
16#ifndef O2_PVERTEXER_H
17#define O2_PVERTEXER_H
18
19#include <array>
20#include <utility>
26#include "MathUtils/Utils.h"
35#include "gsl/span"
36#include <numeric>
37#include <TTree.h>
38#include <TFile.h>
39#include <TStopwatch.h>
40
41//#define _PV_DEBUG_TREE_ // if enabled, produce dbscan and vertex comparison dump
42
43namespace o2
44{
45namespace vertexing
46{
47
48namespace o2d = o2::dataformats;
49
51{
52
53 public:
54 enum class FitStatus : int { Failure,
58 OK };
59 PVertexer();
60 void init();
61 void end();
62 template <typename TR>
63 int process(const TR& tracks, const gsl::span<o2d::GlobalTrackID> gids, const gsl::span<InteractionCandidate> intCand,
64 std::vector<PVertex>& vertices, std::vector<o2d::VtxTrackIndex>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs,
65 const gsl::span<const o2::MCCompLabel> lblTracks, std::vector<o2::MCEventLabel>& lblVtx);
66 template <typename TR>
67 int process(const TR& tracks, const gsl::span<o2d::GlobalTrackID> gids, const gsl::span<o2::InteractionRecord> intRec,
68 std::vector<PVertex>& vertices, std::vector<o2d::VtxTrackIndex>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs,
69 const gsl::span<const o2::MCCompLabel> lblTracks, std::vector<o2::MCEventLabel>& lblVtx);
70
71 int processFromExternalPool(const std::vector<TrackVF>& pool, std::vector<PVertex>& vertices, std::vector<o2d::VtxTrackIndex>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs);
72
73 bool findVertex(const VertexingInput& input, PVertex& vtx);
74
75 void setStartIR(const o2::InteractionRecord& ir) { mStartIR = ir; }
76
77 void setTukey(float t)
78 {
79 mTukey2I = t > 0.f ? 1.f / (t * t) : 1.f / (PVertexerParams::kDefTukey * PVertexerParams::kDefTukey);
80 }
81 float getTukey() const;
82
83 bool setCompatibleIR(PVertex& vtx);
84
85 void setBunchFilling(const o2::BunchFilling& bf);
86
87 void setBz(float bz) { mBz = bz; }
88 void setValidateWithIR(bool v) { mValidateWithIR = v; }
89 bool getValidateWithIR() const { return mValidateWithIR; }
91
92 auto& getTracksPool() const { return mTracksPool; }
93 auto& getTimeZClusters() const { return mTimeZClusters; }
94
95 auto& getMeanVertex() const { return mMeanVertex; }
97 {
98 if (!v) {
99 return;
100 }
101 mMeanVertex = *v;
102 mMeanVertexSeed = *v;
103 initMeanVertexConstraint();
104 }
105
107 {
108 mITSROFrameLengthMUS = v;
109 }
110
111 // special methods used to refit already found vertex skipping certain number of tracks.
112 template <typename TR>
113 bool prepareVertexRefit(const TR& tracks, const o2d::VertexBase& vtxSeed);
114
115 PVertex refitVertex(const std::vector<bool> useTrack, const o2d::VertexBase& vtxSeed);
116
117 auto getNTZClusters() const { return mNTZClustersIni; }
118 auto getTotTrials() const { return mTotTrials; }
119 auto getMaxTrialsPerCluster() const { return mMaxTrialPerCluster; }
120 auto getLongestClusterMult() const { return mLongestClusterMult; }
121 auto getLongestClusterTimeMS() const { return mLongestClusterTimeMS; }
122 auto getNKilledBCValid() const { return mNKilledBCValid; }
123 auto getNKilledIntCand() const { return mNKilledIntCand; }
124 auto getNKilledDebris() const { return mNKilledDebris; }
125 auto getNKilledQuality() const { return mNKilledQuality; }
126 auto getNKilledITSOnly() const { return mNKilledITSOnly; }
127 auto getNIniFound() const { return mNIniFound; }
128
129 TStopwatch& getTimeDBScan() { return mTimeDBScan; }
130 TStopwatch& getTimeVertexing() { return mTimeVertexing; }
131 TStopwatch& getTimeDebris() { return mTimeDebris; }
132 TStopwatch& getTimeMADSel() { return mTimeMADSel; }
133 TStopwatch& getTimeReAttach() { return mTimeReAttach; }
134
135 void setPoolDumpDirectory(const std::string& d) { mPoolDumpDirectory = d; }
136
137 void printInpuTracksStatus(const VertexingInput& input) const;
138
139 private:
140 static constexpr int DBS_UNDEF = -2, DBS_NOISE = -1, DBS_INCHECK = -10;
141
142 SeedHistoTZ buildHistoTZ(const VertexingInput& input);
143 int runVertexing(gsl::span<o2d::GlobalTrackID> gids, const gsl::span<InteractionCandidate> intCand,
144 std::vector<PVertex>& vertices, std::vector<o2d::VtxTrackIndex>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs,
145 gsl::span<const o2::MCCompLabel> lblTracks, std::vector<o2::MCEventLabel>& lblVtx);
146 void createMCLabels(gsl::span<const o2::MCCompLabel> lblTracks, const std::vector<uint32_t>& trackIDs, const std::vector<V2TRef>& v2tRefs, std::vector<o2::MCEventLabel>& lblVtx);
147 void reduceDebris(std::vector<PVertex>& vertices, std::vector<int>& timeSort, const std::vector<o2::MCEventLabel>& lblVtx);
148 FitStatus fitIteration(const VertexingInput& input, VertexSeed& vtxSeed);
149 void finalizeVertex(const VertexingInput& input, const PVertex& vtx, std::vector<PVertex>& vertices, std::vector<V2TRef>& v2tRefs, std::vector<uint32_t>& trackIDs, SeedHistoTZ* histo = nullptr);
150 void accountTrack(TrackVF& trc, VertexSeed& vtxSeed) const;
151 bool solveVertex(VertexSeed& vtxSeed) const;
152 FitStatus evalIterations(VertexSeed& vtxSeed, PVertex& vtx) const;
153 TimeEst timeEstimate(const VertexingInput& input) const;
154 float findZSeedHistoPeak() const;
155 void initMeanVertexConstraint();
156 void applyConstraint(VertexSeed& vtxSeed) const;
157 bool upscaleSigma(VertexSeed& vtxSeed) const;
158 bool relateTrackToMeanVertex(o2::track::TrackParCov& trc, float vtxErr2);
159 bool relateTrackToVertex(o2::track::TrackParCov& trc, const o2d::VertexBase& vtxSeed) const;
160 void applyMADSelection(std::vector<PVertex>& vertices, std::vector<int>& timeSort, const std::vector<V2TRef>& v2tRefs, const std::vector<uint32_t>& trackIDs);
161 void applITSOnlyFractionCut(std::vector<PVertex>& vertices, std::vector<int>& timeSort, const std::vector<V2TRef>& v2tRefs, const std::vector<uint32_t>& trackIDs);
162 void applInteractionValidation(std::vector<PVertex>& vertices, std::vector<int>& timeSort, const gsl::span<InteractionCandidate> intCand, int minContrib);
163
164 template <typename TR>
165 void createTracksPool(const TR& tracks, gsl::span<const o2d::GlobalTrackID> gids);
166
167 int findVertices(const VertexingInput& input, std::vector<PVertex>& vertices, std::vector<uint32_t>& trackIDs, std::vector<V2TRef>& v2tRefs);
168 void reAttach(std::vector<PVertex>& vertices, std::vector<int>& timeSort, std::vector<uint32_t>& trackIDs, std::vector<V2TRef>& v2tRefs);
169
170 std::pair<int, int> getBestIR(const PVertex& vtx, const gsl::span<InteractionCandidate> intCand, int& currEntry) const;
171
172 int dbscan_RangeQuery(int idxs, std::vector<int>& cand, std::vector<int>& status);
173 void dbscan_clusterize();
174 void doDBScanDump(const VertexingInput& input, gsl::span<const o2::MCCompLabel> lblTracks);
175 void doVtxDump(std::vector<PVertex>& vertices, std::vector<uint32_t> trackIDsLoc, std::vector<V2TRef>& v2tRefsLoc, gsl::span<const o2::MCCompLabel> lblTracks);
176 void doDBGPoolDump(gsl::span<const o2::MCCompLabel> lblTracks);
177 void dumpPool();
178
179 o2::BunchFilling mBunchFilling;
180 std::array<int16_t, o2::constants::lhc::LHCMaxBunches> mClosestBunchAbove{-1}; // closest filled bunch from above, 1st element -1 to disable usage by default
181 std::array<int16_t, o2::constants::lhc::LHCMaxBunches> mClosestBunchBelow{-1}; // closest filled bunch from below, 1st element -1 to disable usage by default
182 o2d::MeanVertexObject mMeanVertex{}; // calibrated mean vertex object
183 o2d::VertexBase mMeanVertexSeed{}; // mean vertex at particular Z (accounting for slopes
184 std::array<float, 3> mXYConstraintInvErr = {1.0f, 0.f, 1.0f};
185 //
186 std::vector<TrackVF> mTracksPool;
187 std::vector<TimeZCluster> mTimeZClusters;
188 float mITSROFrameLengthMUS = 0;
189 float mBz = 0.;
190 float mDBScanDeltaT = 0.;
191 float mDBSMaxZ2InvCorePoint = 0;
192 bool mValidateWithIR = false;
193 o2::InteractionRecord mStartIR{0, 0};
194 // structure for the vertex refit
195 o2d::VertexBase mVtxRefitOrig{};
196 std::vector<int> mRefitTrackIDs{};
197 //
198
200 GTrackID::mask_t mTrackSrc{};
201 std::vector<int> mSrcVec{};
202 const PVertexerParams* mPVParams = nullptr;
203 float mTukey2I = 0;
204 static constexpr float kDefTukey = 5.0f;
205 static constexpr float kHugeF = 1.e12;
206 static constexpr float kAlmost0F = 1e-12;
207 static constexpr double kAlmost0D = 1e-16;
208 int mNIniFound = 0;
209 int mNKilledBCValid = 0;
210 int mNKilledIntCand = 0;
211 int mNKilledDebris = 0;
212 int mNKilledQuality = 0;
213 int mNKilledITSOnly = 0;
214 size_t mNTZClustersIni = 0;
215 size_t mTotTrials = 0;
216 size_t mMaxTrialPerCluster = 0;
217 float mMaxTDiffDebris = 0;
218 float mMaxTDiffDebrisExtra = 0;
219 float mMaxTDiffDebrisFiducial = 0;
220 float mMaxZDiffDebrisFiducial = 0;
221 float mMaxMultRatDebrisFiducial = 0;
222 long mLongestClusterTimeMS = 0;
223 int mLongestClusterMult = 0;
224 bool mPoolDumpProduced = false;
225 bool mITSOnly = false;
226 TStopwatch mTimeDBScan;
227 TStopwatch mTimeVertexing;
228 TStopwatch mTimeDebris;
229 TStopwatch mTimeMADSel;
230 TStopwatch mTimeReAttach;
231 std::string mPoolDumpDirectory{};
232#ifdef _PV_DEBUG_TREE_
233 std::unique_ptr<TFile> mDebugDumpFile;
234 std::unique_ptr<TTree> mDebugDBScanTree;
235 std::unique_ptr<TTree> mDebugPoolTree;
236 std::unique_ptr<TTree> mDebugVtxTree;
237 std::unique_ptr<TTree> mDebugVtxCompTree;
238
239 std::vector<TrackVFDump> mDebugDumpDBSTrc;
240 std::vector<GTrackID> mDebugDumpDBSGID;
241 std::vector<o2::MCCompLabel> mDebugDumpDBSTrcMC;
242
243 PVertex mDebugDumpVtx;
244 std::vector<TrackVFDump> mDebugDumpVtxTrc;
245 std::vector<GTrackID> mDebugDumpVtxGID;
246 std::vector<o2::MCCompLabel> mDebugDumpVtxTrcMC;
247
248 std::vector<PVtxCompDump> mDebugDumpPVComp;
249 std::vector<o2::MCEventLabel> mDebugDumpPVCompLbl0; // for some reason the added as a class member
250 std::vector<o2::MCEventLabel> mDebugDumpPVCompLbl1; // gets stored as simple uint
251#endif
252};
253
254//___________________________________________________________________
255inline void PVertexer::applyConstraint(VertexSeed& vtxSeed) const
256{
257 // impose meanVertex constraint, i.e. account terms
258 // (V_i-Constrain_i)^2/sig2constr_i for i=X,Y in the fit chi2 definition
259 vtxSeed.cxx += mXYConstraintInvErr[0];
260 vtxSeed.cxy += mXYConstraintInvErr[1];
261 vtxSeed.cyy += mXYConstraintInvErr[2];
262 float xv = mMeanVertex.getXAtZ(vtxSeed.getZ()), yv = mMeanVertex.getYAtZ(vtxSeed.getZ());
263 vtxSeed.cx0 += mXYConstraintInvErr[0] * xv + mXYConstraintInvErr[1] * yv;
264 vtxSeed.cy0 += mXYConstraintInvErr[1] * xv + mXYConstraintInvErr[2] * yv;
265}
266
267//___________________________________________________________________
268inline bool PVertexer::upscaleSigma(VertexSeed& vtxSeed) const
269{
270 // scale upward the scaleSigma2 if needes
271 if (vtxSeed.scaleSigma2 < mPVParams->maxScale2) {
272 auto s = vtxSeed.scaleSigma2 * mPVParams->upscaleFactor;
273 vtxSeed.setScale(s > mPVParams->maxScale2 ? mPVParams->maxScale2 : s, mTukey2I);
274 return true;
275 }
276 return false;
277}
278
279//___________________________________________________________________
280template <typename TR>
281void PVertexer::createTracksPool(const TR& tracks, gsl::span<const o2d::GlobalTrackID> gids)
282{
283 // create pull of all candidate tracks in a global array ordered in time
284 mTracksPool.clear();
285 auto ntGlo = tracks.size();
286 std::vector<int> sortedTrackID(ntGlo);
287 mTracksPool.reserve(ntGlo);
288 std::iota(sortedTrackID.begin(), sortedTrackID.end(), 0);
289 std::sort(sortedTrackID.begin(), sortedTrackID.end(), [&tracks](int i, int j) {
290 return tracks[i].timeEst.getTimeStamp() < tracks[j].timeEst.getTimeStamp();
291 });
292
293 // check all containers
294 float vtxErr2 = 0.5 * (mMeanVertex.getSigmaX2() + mMeanVertex.getSigmaY2()) + mPVParams->meanVertexExtraErrSelection * mPVParams->meanVertexExtraErrSelection;
296
297 for (uint32_t i = 0; i < ntGlo; i++) {
298 int id = sortedTrackID[i];
300 trc.updateCov(mPVParams->sysErrY2, o2::track::kSigY2);
301 trc.updateCov(mPVParams->sysErrZ2, o2::track::kSigZ2);
302 if (!relateTrackToMeanVertex(trc, vtxErr2)) {
303 continue;
304 }
305 auto& tvf = mTracksPool.emplace_back(trc, tracks[id].getTimeMUS(), id, gids[id], mPVParams->addTimeSigma2, mPVParams->addZSigma2);
306 if (tvf.wghHisto < 0) {
307 mTracksPool.pop_back(); // discard bad track
308 continue;
309 }
310 if (gids[id].getSource() == GTrackID::ITSTPC) { // if the track was adjusted to ITS ROF boundary, flag it
311 float bcf = tvf.timeEst.getTimeStamp() / o2::constants::lhc::LHCBunchSpacingMUS + o2::constants::lhc::LHCMaxBunches;
312 int bcWrtROF = int(bcf - alpParams.roFrameBiasInBC) % alpParams.roFrameLengthInBC;
313 if (bcWrtROF == 0) {
314 float dbc = bcf - (int(bcf / alpParams.roFrameBiasInBC)) * alpParams.roFrameBiasInBC;
315 if (std::abs(dbc) < 1e-6) {
316 tvf.setITSTPCAdjusted();
317 LOGP(debug, "Adjusted t={} -> bcf={} dbc = {}", tvf.timeEst.getTimeStamp(), bcf, dbc);
318 }
319 }
320 }
321 }
322
323 if (mTracksPool.empty()) {
324 return;
325 }
326 //
327 auto tMin = mTracksPool.front().timeEst.getTimeStamp();
328 auto tMax = mTracksPool.back().timeEst.getTimeStamp();
329}
330
331//___________________________________________________________________
332template <typename TR>
333int PVertexer::process(const TR& tracks, const gsl::span<o2d::GlobalTrackID> gids, const gsl::span<InteractionCandidate> intCand,
334 std::vector<PVertex>& vertices, std::vector<o2d::VtxTrackIndex>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs,
335 const gsl::span<const o2::MCCompLabel> lblTracks, std::vector<o2::MCEventLabel>& lblVtx)
336{
337 createTracksPool(tracks, gids);
338 return runVertexing(gids, intCand, vertices, vertexTrackIDs, v2tRefs, lblTracks, lblVtx);
339}
340
341//___________________________________________________________________
342template <typename TR>
343int PVertexer::process(const TR& tracks, const gsl::span<o2d::GlobalTrackID> gids, const gsl::span<o2::InteractionRecord> intRec,
344 std::vector<PVertex>& vertices, std::vector<o2d::VtxTrackIndex>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs,
345 const gsl::span<const o2::MCCompLabel> lblTracks, std::vector<o2::MCEventLabel>& lblVtx)
346{
347 createTracksPool(tracks, gids);
348 static std::vector<InteractionCandidate> intCand;
349 intCand.clear();
350 intCand.reserve(intRec.size());
351 for (const auto& ir : intRec) {
352 intCand.emplace_back(InteractionCandidate{ir, float(ir.differenceInBC(mStartIR) * o2::constants::lhc::LHCBunchSpacingMUS), 0, 0});
353 }
354 return runVertexing(gids, intCand, vertices, vertexTrackIDs, v2tRefs, lblTracks, lblVtx);
355}
356
357//___________________________________________________________________
358template <typename TR>
359bool PVertexer::prepareVertexRefit(const TR& tracks, const o2d::VertexBase& vtxSeed)
360{
361 // Create pull of all tracks for the refitting starting from the tracks of already found vertex.
362 // This method should be called once before calling refitVertex method for given vtxSeed (which can be
363 // called multiple time with different useTrack request)
364 mTracksPool.clear();
365 mTracksPool.reserve(tracks.size());
366 for (uint32_t i = 0; i < tracks.size(); i++) {
367 o2::track::TrackParCov trc = tracks[i];
368 if (!relateTrackToVertex(trc, vtxSeed)) {
369 continue;
370 }
371 auto& tvf = mTracksPool.emplace_back(trc, TimeEst{0.f, 1.f}, i, GTrackID{}, 1., mPVParams->addZSigma2);
372 tvf.entry = i;
373 }
374 if (int(mTracksPool.size()) < mPVParams->minTracksPerVtx) {
375 return false;
376 }
377 mRefitTrackIDs.resize(tracks.size());
378 std::iota(mRefitTrackIDs.begin(), mRefitTrackIDs.end(), 0);
379 mVtxRefitOrig = vtxSeed;
380 return true;
381}
382
383} // namespace vertexing
384} // namespace o2
385#endif
General auxilliary methods.
Base track model for the Barrel, params only, w/o covariance.
int32_t i
Global index for barrel track: provides provenance (detectors combination), index in respective array...
Header to collect LHC related constants.
Primary vertex finder helper classes.
uint32_t j
Definition RawData.h:0
std::ostringstream debug
void setTrackSources(GTrackID::mask_t s)
void setBz(float bz)
Definition PVertexer.h:87
auto getNIniFound() const
Definition PVertexer.h:127
int process(const TR &tracks, const gsl::span< o2d::GlobalTrackID > gids, const gsl::span< InteractionCandidate > intCand, std::vector< PVertex > &vertices, std::vector< o2d::VtxTrackIndex > &vertexTrackIDs, std::vector< V2TRef > &v2tRefs, const gsl::span< const o2::MCCompLabel > lblTracks, std::vector< o2::MCEventLabel > &lblVtx)
Definition PVertexer.h:333
auto getTotTrials() const
Definition PVertexer.h:118
void setStartIR(const o2::InteractionRecord &ir)
set InteractionRecods for the beginning of the TF
Definition PVertexer.h:75
auto getLongestClusterMult() const
Definition PVertexer.h:120
bool getValidateWithIR() const
Definition PVertexer.h:89
auto & getTracksPool() const
Definition PVertexer.h:92
auto getLongestClusterTimeMS() const
Definition PVertexer.h:121
TStopwatch & getTimeDebris()
Definition PVertexer.h:131
bool prepareVertexRefit(const TR &tracks, const o2d::VertexBase &vtxSeed)
Definition PVertexer.h:359
auto & getTimeZClusters() const
Definition PVertexer.h:93
auto getMaxTrialsPerCluster() const
Definition PVertexer.h:119
auto getNKilledITSOnly() const
Definition PVertexer.h:126
TStopwatch & getTimeVertexing()
Definition PVertexer.h:130
auto & getMeanVertex() const
Definition PVertexer.h:95
TStopwatch & getTimeMADSel()
Definition PVertexer.h:132
void printInpuTracksStatus(const VertexingInput &input) const
auto getNKilledQuality() const
Definition PVertexer.h:125
bool findVertex(const VertexingInput &input, PVertex &vtx)
PVertex refitVertex(const std::vector< bool > useTrack, const o2d::VertexBase &vtxSeed)
void setPoolDumpDirectory(const std::string &d)
Definition PVertexer.h:135
auto getNTZClusters() const
Definition PVertexer.h:117
auto getNKilledDebris() const
Definition PVertexer.h:124
int processFromExternalPool(const std::vector< TrackVF > &pool, std::vector< PVertex > &vertices, std::vector< o2d::VtxTrackIndex > &vertexTrackIDs, std::vector< V2TRef > &v2tRefs)
void setBunchFilling(const o2::BunchFilling &bf)
bool setCompatibleIR(PVertex &vtx)
void setValidateWithIR(bool v)
Definition PVertexer.h:88
TStopwatch & getTimeReAttach()
Definition PVertexer.h:133
TStopwatch & getTimeDBScan()
Definition PVertexer.h:129
auto getNKilledBCValid() const
Definition PVertexer.h:122
void setTukey(float t)
Definition PVertexer.h:77
void setMeanVertex(const o2d::MeanVertexObject *v)
Definition PVertexer.h:96
void setITSROFrameLength(float v)
Definition PVertexer.h:106
auto getNKilledIntCand() const
Definition PVertexer.h:123
const GLdouble * v
Definition glcorearb.h:832
GLuint id
Definition glcorearb.h:650
constexpr double LHCBunchSpacingMUS
constexpr int LHCMaxBunches
Definition of a container to keep/associate and arbitrary number of labels associated to an index wit...
TrackParCovF TrackParCov
Definition Track.h:33
struct o2::upgrades_utils::@453 tracks
structure to keep trigger-related info
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
int64_t differenceInBC(const InteractionRecord &other) const
float meanVertexExtraErrSelection
extra error to meanvertex sigma used when selecting tracks
int minTracksPerVtx
min N tracks per vertex
float sysErrY2
systematic error on track Y error
float addTimeSigma2
increment time error^2 by this amount when calculating histo weight
float maxScale2
max slaling factor^2
float addZSigma2
increment z error^2 by this amount when calculating histo weight
float sysErrZ2
systematic error on track Z error
static constexpr float kDefTukey
def.value for tukey constant
float upscaleFactor
factor for upscaling if not candidate is found
weights and scaling params for current vertex
o2::InteractionRecord ir(0, 0)