5 namespace l1tVertexFinder {
7 void VertexFinder::computeAndSetVertexParameters(
RecoVertex<>& vertex,
8 const std::vector<float>& bin_centers,
9 const std::vector<unsigned int>&
counts) {
14 double highestPt = 0.;
15 unsigned int numHighPtTracks = 0;
21 std::vector<double> bin_pt(bin_centers.size(), 0.0);
22 unsigned int ibin = 0;
23 unsigned int itrack = 0;
27 trackPt =
track->pt();
28 if (trackPt > settings_->vx_TrackMaxPt()) {
31 highestPt = (trackPt > highestPt) ? trackPt : highestPt;
32 if (settings_->vx_TrackMaxPtBehavior() == 0)
34 else if (settings_->vx_TrackMaxPtBehavior() == 1)
35 trackPt = settings_->vx_TrackMaxPt();
38 pt +=
std::pow(trackPt, settings_->vx_weightedmean());
39 if (bin_centers.empty() && counts.empty()) {
40 SumZ +=
track->z0() *
std::pow(trackPt, settings_->vx_weightedmean());
43 bin_pt[ibin] +=
std::pow(trackPt, settings_->vx_weightedmean());
44 if (itrack == counts[ibin]) {
45 SumZ += bin_centers[ibin] * bin_pt[ibin];
46 z0square += bin_centers[ibin] * bin_centers[ibin];
53 z0 = SumZ / ((settings_->vx_weightedmean() > 0) ? pt : vertex.
numTracks());
57 vertex.
setParameters(pt, z0, z0width, highPt, numHighPtTracks, highestPt);
60 void VertexFinder::GapClustering() {
64 for (
unsigned int i = 0;
i < fitTracks_.size(); ++
i) {
67 if ((i + 1 < fitTracks_.size() and fitTracks_[i + 1].z0() - fitTracks_[
i].z0() > settings_->vx_distance())
or
68 i == fitTracks_.size() - 1) {
69 if (Vertex.
numTracks() >= settings_->vx_minTracks()) {
70 computeAndSetVertexParameters(Vertex, {}, {});
71 vertices_.push_back(Vertex);
82 if (
std::abs(track0->z0() - track1->z0()) > distance) {
83 distance =
std::abs(track0->z0() - track1->z0());
95 if (
std::abs(track0->z0() - track1->z0()) < distance) {
96 distance =
std::abs(track0->z0() - track1->z0());
105 float distanceSum = 0;
109 distanceSum +=
std::abs(track0->z0() - track1->z0());
118 computeAndSetVertexParameters(cluster0, {}, {});
119 computeAndSetVertexParameters(cluster1, {}, {});
124 void VertexFinder::agglomerativeHierarchicalClustering() {
129 std::vector<RecoVertex<>> vClusters;
130 vClusters.resize(fitTracks_.size());
132 for (
unsigned int i = 0;
i < fitTracks_.size(); ++
i) {
133 vClusters[
i].insert(&fitTracks_[
i]);
138 float MinimumScore = 9999;
140 unsigned int clusterId0 = 0;
141 unsigned int clusterId1 = 0;
142 for (
unsigned int iClust = 0; iClust < vClusters.size() - 1; iClust++) {
146 if (settings_->vx_distanceType() == 0)
147 M = maxDistance(vClusters[iClust], vClusters[iClust + 1]);
148 else if (settings_->vx_distanceType() == 1)
149 M = minDistance(vClusters[iClust], vClusters[iClust + 1]);
150 else if (settings_->vx_distanceType() == 2)
151 M = meanDistance(vClusters[iClust], vClusters[iClust + 1]);
153 M = centralDistance(vClusters[iClust], vClusters[iClust + 1]);
155 if (M < MinimumScore) {
158 clusterId1 = iClust + 1;
161 if (MinimumScore > settings_->vx_distance()
or vClusters[clusterId1].tracks().empty())
164 vClusters[clusterId1].insert(
track);
166 vClusters.erase(vClusters.begin() + clusterId0);
170 if (clust.numTracks() >= settings_->vx_minTracks()) {
171 computeAndSetVertexParameters(clust, {}, {});
172 vertices_.push_back(clust);
177 void VertexFinder::DBSCAN() {
179 std::vector<unsigned int>
visited;
180 std::vector<unsigned int> saved;
185 for (
unsigned int i = 0;
i < fitTracks_.size(); ++
i) {
186 if (
find(visited.begin(), visited.end(),
i) != visited.end())
190 visited.push_back(
i);
191 std::set<unsigned int> neighbourTrackIds;
192 unsigned int numDensityTracks = 0;
193 if (fitTracks_[
i].
pt() > settings_->vx_dbscan_pt())
197 for (
unsigned int k = 0;
k < fitTracks_.size(); ++
k) {
199 if (
k !=
i and
std::abs(fitTracks_[
k].z0() - fitTracks_[
i].z0()) < settings_->vx_distance()) {
200 neighbourTrackIds.insert(
k);
201 if (fitTracks_[
k].
pt() > settings_->vx_dbscan_pt()) {
207 if (numDensityTracks < settings_->vx_dbscan_mintracks()) {
213 for (
unsigned int id : neighbourTrackIds) {
214 if (
find(visited.begin(), visited.end(),
id) == visited.end()) {
215 visited.push_back(
id);
216 std::vector<unsigned int> neighbourTrackIds2;
217 for (
unsigned int k = 0;
k < fitTracks_.size(); ++
k) {
219 if (
std::abs(fitTracks_[
k].z0() - fitTracks_[
id].z0()) < settings_->vx_distance()) {
220 neighbourTrackIds2.push_back(
k);
225 for (
unsigned int id2 : neighbourTrackIds2) {
226 neighbourTrackIds.insert(id2);
230 if (
find(saved.begin(), saved.end(),
id) == saved.end())
231 vertex.
insert(&fitTracks_[
id]);
233 computeAndSetVertexParameters(vertex, {}, {});
234 if (vertex.
numTracks() >= settings_->vx_minTracks())
235 vertices_.push_back(vertex);
241 void VertexFinder::PVR() {
246 acceptedTracks.push_back(
track);
249 while (discardedTracks.size() >= settings_->vx_minTracks()
or start ==
true) {
251 bool removing =
true;
252 discardedTracks.clear();
254 float oldDistance = 0.;
256 if (settings_->debug() > 2)
257 edm::LogInfo(
"VertexFinder") <<
"PVR::AcceptedTracks " << acceptedTracks.size();
261 z0start +=
track.z0();
265 z0start /= acceptedTracks.size();
266 if (settings_->debug() > 2)
267 edm::LogInfo(
"VertexFinder") <<
"PVR::z0 vertex " << z0start;
268 FitTrackCollection::iterator badTrackIt = acceptedTracks.end();
271 for (FitTrackCollection::iterator it = acceptedTracks.begin(); it < acceptedTracks.end(); ++it) {
274 if (
std::abs(track->
z0() - z0start) > settings_->vx_distance() and
275 std::abs(track->
z0() - z0start) > oldDistance) {
277 oldDistance =
std::abs(track->
z0() - z0start);
283 const L1Track badTrack = *badTrackIt;
284 if (settings_->debug() > 2)
285 edm::LogInfo(
"VertexFinder") <<
"PVR::Removing track " << badTrack.
z0() <<
" at distance " << oldDistance;
286 discardedTracks.push_back(badTrack);
287 acceptedTracks.erase(badTrackIt);
291 if (acceptedTracks.size() >= settings_->vx_minTracks()) {
296 computeAndSetVertexParameters(vertex, {}, {});
297 vertices_.push_back(vertex);
299 if (settings_->debug() > 2)
300 edm::LogInfo(
"VertexFinder") <<
"PVR::DiscardedTracks size " << discardedTracks.size();
301 acceptedTracks.clear();
302 acceptedTracks = discardedTracks;
306 void VertexFinder::adaptiveVertexReconstruction() {
312 discardedTracks.push_back(
track);
315 while (discardedTracks.size() >= settings_->vx_minTracks()
or start ==
true) {
317 discardedTracks2.clear();
318 FitTrackCollection::iterator it = discardedTracks.begin();
320 acceptedTracks.push_back(track);
321 float z0sum = track.
z0();
323 for (FitTrackCollection::iterator it2 = discardedTracks.begin(); it2 < discardedTracks.end(); ++it2) {
325 const L1Track secondTrack = *it2;
327 z0sum += secondTrack.
z0();
328 float z0vertex = z0sum / (acceptedTracks.size() + 1);
332 for (
const L1Track& accTrack : acceptedTracks) {
334 float Residual = accTrack.z0() - z0vertex;
344 chi2 += Residual * Residual;
345 dof = (acceptedTracks.size() + 1) * 2 - 1;
347 if (chi2 / dof < settings_->vx_chi2cut()) {
348 acceptedTracks.push_back(secondTrack);
350 discardedTracks2.push_back(secondTrack);
351 z0sum -= secondTrack.
z0();
356 if (acceptedTracks.size() >= settings_->vx_minTracks()) {
358 for (
const L1Track& track : acceptedTracks) {
361 computeAndSetVertexParameters(vertex, {}, {});
362 vertices_.push_back(vertex);
365 acceptedTracks.clear();
366 discardedTracks.clear();
367 discardedTracks = discardedTracks2;
371 void VertexFinder::HPV() {
379 if (
track.pt() < 50.) {
383 vertex.insert(&
track);
386 vertex.insert(&
track);
391 computeAndSetVertexParameters(vertex, {}, {});
393 vertices_.push_back(vertex);
396 void VertexFinder::Kmeans() {
397 unsigned int NumberOfClusters = settings_->vx_kmeans_nclusters();
399 vertices_.resize(NumberOfClusters);
400 float ClusterSeparation = 30. / NumberOfClusters;
402 for (
unsigned int i = 0;
i < NumberOfClusters; ++
i) {
403 float ClusterCentre = -15. + ClusterSeparation * (
i + 0.5);
404 vertices_[
i].setZ0(ClusterCentre);
406 unsigned int iterations = 0;
408 while (iterations < settings_->vx_kmeans_iterations()) {
409 for (
unsigned int i = 0;
i < NumberOfClusters; ++
i) {
410 vertices_[
i].clear();
415 if (iterations == settings_->vx_kmeans_iterations() - 3)
416 distance = settings_->vx_distance() * 2;
417 if (iterations > settings_->vx_kmeans_iterations() - 3)
418 distance = settings_->vx_distance();
419 unsigned int ClusterId;
421 for (
unsigned int id = 0;
id < NumberOfClusters; ++
id) {
429 vertices_[ClusterId].insert(&
track);
432 for (
unsigned int i = 0;
i < NumberOfClusters; ++
i) {
433 if (vertices_[
i].numTracks() >= settings_->vx_minTracks())
434 computeAndSetVertexParameters(vertices_[
i], {}, {});
440 void VertexFinder::findPrimaryVertex() {
444 for (
unsigned int i = 0;
i < vertices_.size(); ++
i) {
445 if (vertices_[
i].
pt() > vertexPt) {
446 vertexPt = vertices_[
i].pt();
452 void VertexFinder::associatePrimaryVertex(
double trueZ0) {
454 for (
unsigned int id = 0;
id < vertices_.size(); ++
id) {
456 distance =
std::abs(trueZ0 - vertices_[
id].z0());
462 void VertexFinder::fastHistoLooseAssociation() {
466 for (
float z = settings_->vx_histogram_min(); z < settings_->vx_histogram_max();
467 z += settings_->vx_histogram_binwidth()) {
474 computeAndSetVertexParameters(vertex, {}, {});
476 if (vertex.
pt() > vxPt) {
477 leading_vertex = vertex;
482 vertices_.emplace_back(leading_vertex);
489 std::ceil((settings_->vx_histogram_max() - settings_->vx_histogram_min()) / settings_->vx_histogram_binwidth());
490 std::vector<RecoVertex<>>
hist(nbins);
491 std::vector<RecoVertex<>> sums(nbins - settings_->vx_windowSize());
492 std::vector<float> bounds(nbins + 1);
495 settings_->vx_histogram_min(),
496 settings_->vx_histogram_binwidth());
500 if ((
track.z0() < settings_->vx_histogram_min()) || (
track.z0() > settings_->vx_histogram_max()))
502 if (
track.getTTTrackPtr()->chi2() > settings_->vx_TrackMaxChi2())
504 if (
track.pt() < settings_->vx_TrackMinPt())
508 float nPS = 0., nstubs = 0;
511 const auto& theStubs =
track.getTTTrackPtr()->getStubRefs();
512 if (theStubs.empty()) {
513 edm::LogWarning(
"VertexFinder") <<
"fastHisto::Could not retrieve the vector of stubs.";
518 for (
const auto& stub : theStubs) {
521 DetId detId(stub->getDetId());
531 if (nstubs < settings_->vx_NStubMin())
533 if (nPS < settings_->vx_NStubPSMin())
537 int trk_nstub = (int)
track.getTTTrackPtr()->getStubRefs().size();
538 float chi2dof =
track.getTTTrackPtr()->chi2() / (2 * trk_nstub - 4);
540 if (settings_->vx_DoPtComp()) {
541 float trk_consistency =
track.getTTTrackPtr()->stubPtConsistency();
542 if (trk_nstub == 4) {
549 if (settings_->vx_DoTightChi2()) {
550 if (
track.pt() > 10.0 && chi2dof > 5.0)
557 hist.at(index).insert(&
track);
562 std::vector<float> bin_centers(settings_->vx_windowSize(), 0.0);
563 std::vector<unsigned int>
counts(settings_->vx_windowSize(), 0);
564 for (
unsigned int i = 0;
i < sums.size();
i++) {
565 for (
unsigned int j = 0;
j < settings_->vx_windowSize();
j++) {
566 bin_centers[
j] = settings_->vx_histogram_min() + ((
i +
j) * settings_->vx_histogram_binwidth()) +
567 (0.5 * settings_->vx_histogram_binwidth());
569 sums.at(
i) += hist.at(
i +
j);
571 computeAndSetVertexParameters(sums.at(
i), bin_centers,
counts);
575 float sigma_max = -999;
577 std::vector<int>
found;
578 found.reserve(settings_->vx_nvtx());
579 for (
unsigned int ivtx = 0; ivtx < settings_->vx_nvtx(); ivtx++) {
582 for (
unsigned int i = 0;
i < sums.size();
i++) {
584 if (
find(found.begin(), found.end(),
i) != found.end())
586 if (sums.at(
i).pt() > sigma_max) {
587 sigma_max = sums.at(
i).pt();
592 found.push_back(imax);
593 vertices_.emplace_back(sums.at(imax));
static double constexpr NA
Avogadro's number.
constexpr int32_t ceil(float num)
unsigned int tidRing(const DetId &id) const
uint16_t *__restrict__ id
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventIDconst &, edm::Timestampconst & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
void setZ0(double z)
Set z0 position [cm].
std::vector< L1Track > FitTrackCollection
__host__ __device__ constexpr RandomIt upper_bound(RandomIt first, RandomIt last, const T &value, Compare comp={})
auto const & tracks
cannot be loose
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
unsigned int numTracks() const
Number of tracks originating from this vertex.
Simple wrapper class for TTTrack.
Abs< T >::type abs(const T &t)
static constexpr auto TOB
double z0() const
Vertex z0 position [cm].
const std::vector< const T * > & tracks() const
Tracks in the vertex.
double pt() const
Sum of fitted tracks transverse momentum [GeV].
void clear()
Clear track vector.
__host__ __device__ constexpr RandomIt lower_bound(RandomIt first, RandomIt last, const T &value, Compare comp={})
Log< level::Warning, false > LogWarning
static constexpr auto TID
Power< A, B >::type pow(const A &a, const B &b)
void setParameters(double pt, double z0, double width=-999., bool highPt=false, unsigned int nHighPt=-999, double highestPt=-999., bool pv=false)
Set the vertex parameters.
unsigned int tobLayer(const DetId &id) const
void insert(const T *fitTrack)
Assign fitted track to this vertex.