10 #include <trg/cdc/modules/houghtracking/CDCTrigger2DFinderModule.h>
12 #include <cdc/geometry/CDCGeometryPar.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/gearbox/Const.h>
21 #define CDC_SUPER_LAYERS 9
39 setDescription(
"Hough tracking algorithm for CDC trigger.");
40 setPropertyFlags(c_ParallelProcessingCertified);
43 addParam(
"hitCollectionName", m_hitCollectionName,
44 "Name of the input StoreArray of CDCTriggerSegmentHits.",
46 addParam(
"outputCollectionName", m_outputCollectionName,
47 "Name of the StoreArray holding the tracks found in the Hough tracking.",
48 string(
"TRGCDC2DFinderTracks"));
49 addParam(
"clusterCollectionName", m_clusterCollectionName,
50 "Name of the StoreArray holding the clusters formed in the Hough plane.",
52 addParam(
"nCellsPhi", m_nCellsPhi,
53 "Number of Hough cells in phi (limits: [-180, 180]). Must be an even number.",
55 addParam(
"nCellsR", m_nCellsR,
56 "Number of Hough cells in 1/r. Must be an even number.",
58 addParam(
"minPt", m_minPt,
60 "Hough plane limits in 1/r are [-1/r(minPt), 1/r(minPt)]", (
double)(0.3));
61 addParam(
"shiftPt", m_shiftPt,
62 "Shift the Hough plane by 1/4 cell size in 1/r to avoid "
63 "curvature 0 tracks (<0: shift in negative direction, "
64 "0: no shift, >0: shift in positive direction).", 0);
66 addParam(
"minHits", m_minHits,
67 "Minimum hits from different super layers required in a peak cell.",
69 addParam(
"minHitsShort", m_minHitsShort,
70 "Minimum hits required required in a peak cell for a short track"
71 " (must be in the first minHitsShort super layers).",
73 addParam(
"minCells", m_minCells,
74 "Peaks with less than minCells connected cells are ignored.",
76 addParam(
"onlyLocalMax", m_onlyLocalMax,
77 "Switch to remove cells connected to a cell with higher super layer count.",
79 addParam(
"connect", m_connect,
80 "Definition for connected cells. 4: direct (left/right/top/bottom), "
81 "6: direct + 2 diagonal (top right/bottom left), "
82 "8: direct + all diagonal (top right/top left/bottom right/bottom left)",
84 addParam(
"ignore2ndPriority", m_ignore2nd,
85 "Switch to skip second priority hits.",
false);
86 addParam(
"usePriorityPosition", m_usePriority,
87 "If true, use wire position of priority cell in track segment, "
88 "otherwise use wire position of center cell.",
true);
89 addParam(
"requireSL0", m_requireSL0,
90 "Switch to check separately for a hit in the innermost superlayer.",
false);
91 addParam(
"storeHoughPlane", m_storePlane,
92 "Switch for saving Hough plane as TMatrix in DataStore. "
93 "0: don't store anything, 1: store only peaks, 2: store full plane "
94 "(will increase runtime).", (
unsigned)(0));
95 addParam(
"clusterPattern", m_clusterPattern,
96 "use nested pattern algorithm to find clusters",
true);
97 addParam(
"clusterSizeX", m_clusterSizeX,
98 "maximum number of 2 x 2 squares in cluster for pattern algorithm",
100 addParam(
"clusterSizeY", m_clusterSizeY,
101 "maximum number of 2 x 2 squares in cluster for pattern algorithm",
103 addParam(
"hitRelationsFromCorners", m_hitRelationsFromCorners,
104 "Switch for creating relations to hits in the pattern algorithm. "
105 "If true, create relations from cluster corners, otherwise "
106 "from estimated cluster center (might not have relations).",
false);
108 addParam(
"testFilename", m_testFilename,
109 "If not empty, a file with input (hits) and output (tracks) "
110 "for each event is written (for firmware debugging).",
string(
""));
112 addParam(
"suppressClone", m_suppressClone,
113 "Switch to send only the first found track and suppress the "
114 "subsequent clones." ,
false);
118 CDCTrigger2DFinderModule::initialize()
120 m_segmentHits.isRequired(m_hitCollectionName);
121 m_tracks.registerInDataStore(m_outputCollectionName);
122 m_clusters.registerInDataStore(m_clusterCollectionName);
124 m_tracks.registerRelationTo(m_segmentHits);
125 m_tracks.registerRelationTo(m_clusters);
127 if (m_storePlane > 0) m_houghPlane.registerInDataStore(
"HoughPlane");
132 for (
int iSL = 0; iSL < 9; ++iSL) {
134 nTS += cdc.nWiresInLayer(layerId);
135 TSoffset[iSL + 1] = nTS;
136 for (
int priority = 0; priority < 2; ++ priority) {
137 radius[iSL][priority] = cdc.senseWireR(layerId + priority);
139 layerId += (iSL > 0 ? 6 : 7);
142 if (m_testFilename !=
"") {
143 testFile.open(m_testFilename);
146 if (m_suppressClone) {
147 B2INFO(
"2D finder will exit with the first track candidate in time.");
152 CDCTrigger2DFinderModule::event()
159 setReturnValue(
true);
161 if (m_testFilename !=
"") {
162 testFile << StoreObjPtr<EventMetaData>()->getEvent() <<
" "
163 << m_segmentHits.getEntries() << endl;
166 if (m_segmentHits.getEntries() == 0) {
172 for (
int iHit = 0; iHit < m_segmentHits.getEntries(); iHit++) {
173 unsigned short iSL = m_segmentHits[iHit]->getISuperLayer();
174 if (m_testFilename !=
"") {
175 testFile << iSL <<
" " << m_segmentHits[iHit]->getSegmentID() - TSoffset[iSL] <<
" "
176 << m_segmentHits[iHit]->getPriorityPosition() << endl;
178 if (iSL % 2)
continue;
179 if (m_ignore2nd && m_segmentHits[iHit]->getPriorityPosition() < 3)
continue;
180 double phi = m_segmentHits[iHit]->getSegmentID() - TSoffset[iSL];
182 phi += 0.5 * (((m_segmentHits[iHit]->getPriorityPosition() >> 1) & 1)
183 - (m_segmentHits[iHit]->getPriorityPosition() & 1));
185 phi = phi * 2. * M_PI / (TSoffset[iSL + 1] - TSoffset[iSL]);
186 double r = radius[iSL][int(m_usePriority &&
187 m_segmentHits[iHit]->getPriorityPosition() < 3)];
188 TVector2 pos(cos(phi) / r, sin(phi) / r);
189 hitMap.insert(std::make_pair(iHit, std::make_pair(iSL, pos)));
195 maxIterations = ceil(log2(max(m_nCellsPhi, m_nCellsR))) - 1;
196 nCells = pow(2, maxIterations + 1);
198 double rectX = M_PI * nCells / m_nCellsPhi;
200 maxR = 0.5 * Const::speedOfLight * 1.5e-4 / m_minPt;
201 double rectY = maxR * nCells / m_nCellsR;
204 shiftR = -maxR / 2. / m_nCellsR;
205 }
else if (m_shiftPt > 0) {
206 shiftR = maxR / 2. / m_nCellsR;
209 B2DEBUG(50,
"extending Hough plane to " << maxIterations <<
" iterations, "
210 << nCells <<
" cells: phi in ["
211 << -rectX * 180. / M_PI <<
", " << rectX * 180. / M_PI
212 <<
"] deg, 1/r in [" << -rectY + shiftR <<
", " << rectY + shiftR <<
"] /cm");
215 if (m_storePlane > 0) {
216 m_houghPlane.create();
217 m_houghPlane->ResizeTo(m_nCellsPhi, m_nCellsR);
222 if (m_suppressClone && !hitMap.empty()) {
229 typedef pair<int, cdcPair> cdcHitPair;
231 typedef vector<cdcHitPair> cdcSeqMap;
234 for (
auto hit : hitMap) {
235 seqHitMap.push_back(hit);
237 sort(seqHitMap.begin(), seqHitMap.end(), [
this](cdcHitPair i, cdcHitPair j) {
238 return m_segmentHits[i.first]->foundTime() < m_segmentHits[j.first]->foundTime();
240 auto seqHitItr = seqHitMap.begin();
242 vector<bool> layerHit(CDC_SUPER_LAYERS,
false);
245 short firstTick = m_segmentHits[(*seqHitMap.begin()).first]->foundTime() / period + 1;
246 short lastTick = m_segmentHits[(*(seqHitMap.end() - 1)).first]->foundTime() / period + 1;
248 for (
auto tick = firstTick * period; tick < lastTick * period; tick += period) {
250 for (
auto itr = seqHitItr; itr < seqHitMap.end(); ++itr) {
251 cdcHitPair currentHit = *itr;
253 if (count(layerHit.begin(), layerHit.end(),
true) >= m_minHits &&
254 m_segmentHits[currentHit.first]->foundTime() > tick) {
258 layerHit[m_segmentHits[currentHit.first]->getISuperLayer()] =
true;
260 copy_n(seqHitItr, nHitInCycle, inserter(fastHitMap, fastHitMap.end()));
261 fastInterceptFinder(fastHitMap, -rectX, rectX, -rectY + shiftR, rectY + shiftR, 0, 0, 0);
262 B2DEBUG(20,
"at tick " << tick <<
", number of candidates: " << houghCand.size());
263 if (!houghCand.empty()) {
264 B2DEBUG(10,
"found a track at clock " << tick <<
" with "
265 << fastHitMap.size() <<
"hits");
268 advance(seqHitItr, nHitInCycle);
272 fastInterceptFinder(hitMap, -rectX, rectX, -rectY + shiftR, rectY + shiftR, 0, 0, 0);
273 if (!houghCand.empty()) {
274 B2DEBUG(10,
"found a track with " << hitMap.size() <<
"hits");
279 if (m_clusterPattern)
280 patternClustering(fastHitMap);
284 if (m_testFilename !=
"") {
285 testFile << m_tracks.getEntries() << endl;
286 for (
int i = 0; i < m_tracks.getEntries(); ++i) {
287 float ix = (m_tracks[i]->getPhi0() - M_PI_4) * m_nCellsPhi / 2. / M_PI - 0.5;
288 float iy = (m_tracks[i]->getOmega() / 2. + maxR - shiftR) * m_nCellsR / 2. / maxR - 0.5;
289 testFile << round(2 * ix) / 2. <<
" " << round(2 * iy) / 2. <<
" "
290 << m_tracks[i]->getChargeSign() << endl;
293 testFile << hits.size() << endl;
294 for (
unsigned ihit = 0; ihit < hits.size(); ++ihit) {
295 unsigned short iSL = hits[ihit]->getISuperLayer();
296 testFile << iSL <<
" " << hits[ihit]->getSegmentID() - TSoffset[iSL] <<
" "
297 << hits[ihit]->getPriorityPosition() <<
" "
298 << hits.weight(ihit) << endl;
305 CDCTrigger2DFinderModule::terminate()
307 if (m_testFilename !=
"") testFile.close();