forked from sPHENIX-Collaboration/coresoftware
/
PHTpcSeedFinder.cc
122 lines (105 loc) · 4.1 KB
/
PHTpcSeedFinder.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
/*!
* \file PHTpcSeedFinder.cc
* \brief
* \author Dmitry Arkhipkin <arkhipkin@gmail.com>
*/
#include "PHTpcSeedFinder.h"
#include "PHTpcTrackerUtil.h"
#include "externals/kdfinder.hpp" // for TrackCandidate, Circle, find_tr...
#include <phool/PHLog.h>
#include <log4cpp/CategoryStream.hh> // for CategoryStream
class TrkrClusterContainer;
class TrkrHitSetContainer;
float round(float var)
{
return (int(var * 100.0) / 100.0);
}
PHTpcSeedFinder::PHTpcSeedFinder()
: mMaxDistance1(3.0)
, mTripletAngle1(M_PI / 8)
, mMinHits1(10)
, mMaxDistance2(6.0)
, mTripletAngle2(M_PI / 8)
, mMinHits2(5)
, mNThreads(1)
, mRemoveLoopers(false)
, mMinLooperRadius(10.0)
, mMaxLooperRadius(70.0)
{
}
std::vector<kdfinder::TrackCandidate<double>*> PHTpcSeedFinder::findSeeds(TrkrClusterContainer* cluster_map, TrkrHitSetContainer* hitsets, double B)
{
LOG_WARN_IF("tracking.PHTpcSeedFnder.findSeeds", !cluster_map) << __FILE__ << "," << __LINE__ << " Can't find node TRKR_CLUSTER";
std::vector<kdfinder::TrackCandidate<double>*> result;
if (!cluster_map)
{
return result;
}
LOG_WARN_IF("tracking.PHTpcSeedFnder.findSeeds", !hitsets) << __FILE__ << "," << __LINE__ << " Can't find node TRKR_HITSET";
if (!hitsets)
{
return result;
}
//----- convert clusters to kdhits -----
std::vector<std::vector<double> > kdhits(PHTpcTrackerUtil::convert_clusters_to_hits(cluster_map, hitsets));
LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "imported " << kdhits.size() << " clusters";
//----- find unfitted kdtracks -----
std::vector<std::vector<double> > unused_hits;
std::vector<std::vector<std::vector<double> > > kdtracks;
kdtracks = kdfinder::find_tracks_iterative<double>(kdhits, unused_hits,
mMaxDistance1 /* max distance in cm*/, mTripletAngle1 /* triplet angle */, mMinHits1 /* min hits to keep track */, // first iteration
mMaxDistance2, mTripletAngle2, mMinHits2, // second iteration params
mNThreads /* nthreads */,
false); // print stats
LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "kdtracks: " << kdtracks.size();
//----- create fitted track candidates -----
result = kdfinder::get_track_candidates<double>(kdtracks, B);
LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "initial kdtrack candidates: " << result.size();
//----- filter out junk = bad tracks, <25 MeV Pt tracks etc.. -----
for (auto it = result.begin(); it != result.end();)
{
if (!(*it)->isFitted() || (*it)->Pt() < 0.025 || (*it)->Pt() > 200.0)
{
(*it)->deleteHits();
delete (*it);
it = result.erase(it);
}
else
{
++it;
}
}
LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "kdtrack candidates after filtering: " << result.size();
//----- filter out loopers, if enabled -----
if (mRemoveLoopers)
{
for (auto it = result.begin(); it != result.end();)
{
kdfinder::Circle<double>* c = (*it)->getCircleFit();
double x = c->a, y = c->b, r = c->r, xyr = std::sqrt(x * x + y * y);
if (((xyr - r) > mMinLooperRadius) && ((xyr + r) < mMaxLooperRadius))
{
(*it)->deleteHits();
delete (*it);
it = result.erase(it);
}
else
{
++it;
}
}
}
/*
//----- primitive yet very fast vertex seed finding -----
// FIXME: use it instead of Rave?
std::vector< std::pair< std::vector<double>, std::vector<size_t> > > vertices =
kdfinder::find_vertex_seeds<double>( result, 0, 0, 0.5, 2.0, 3, false );
LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "track candidate vertex seeds: " << vertices.size();
for( auto vertex: vertices ) {
std::vector<double> pos = vertex.first;
std::vector<size_t> ids = vertex.second;
LOG_DEBUG("tracking.PHTpcSeedFinder.findSeeds") << "vertex tracks: " << ids.size() << ", pos: " << pos[0] << ", " << pos[1] << ", " << pos[2];
}
*/
return result;
}