10 #include <trg/cdc/modules/houghETF/CDCTriggerHoughETFModule.h>
12 #include <cdc/geometry/CDCGeometryPar.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/gearbox/Const.h>
20 #include <root/TMatrix.h>
23 #define CDC_SUPER_LAYERS 9
41 setDescription(
"Hough tracking algorithm for CDC trigger.");
42 setPropertyFlags(c_ParallelProcessingCertified);
51 addParam(
"outputEventTimeName", m_EventTimeName,
52 "Name of the output StoreObjPtr.",
54 addParam(
"storeTracks", m_storeTracks,
57 addParam(
"usePriorityTiming", m_usePriorityTiming,
58 "Use priority timing instead of fastest timing.",
60 addParam(
"t0CalcMethod", m_t0CalcMethod,
61 "0: Nth fastest fastest time."
62 "1: median of all timings."
63 "2: median of timings in timing window.",
65 addParam(
"arrivalOrder", m_arrivalOrder,
66 "When t0CalcMethod == 0: Nth fastest ft is used as T0. (i.e. 0 is fastest)",
68 addParam(
"timeWindowStart", m_timeWindowStart,
69 "When t0CalcMethod == 2: start time of time window relative to median. (in ns)",
71 addParam(
"timeWindowEnd", m_timeWindowEnd,
72 "When t0CalcMethod == 2: end time of time window relative to median. (in ns)",
76 addParam(
"hitCollectionName", m_hitCollectionName,
77 "Name of the input StoreArray of CDCTriggerSegmentHits.",
79 addParam(
"outputCollectionName", m_outputCollectionName,
80 "Name of the StoreArray holding the tracks found in the Hough tracking.",
81 string(
"TRGCDCETFTracks"));
82 addParam(
"clusterCollectionName", m_clusterCollectionName,
83 "Name of the StoreArray holding the clusters formed in the Hough plane.",
85 addParam(
"nCellsPhi", m_nCellsPhi,
86 "Number of Hough cells in phi (limits: [-180, 180]). Must be an even number.",
88 addParam(
"nCellsR", m_nCellsR,
89 "Number of Hough cells in 1/r. Must be an even number.",
91 addParam(
"minPt", m_minPt,
93 "Hough plane limits in 1/r are [-1/r(minPt), 1/r(minPt)]", (
double)(0.3));
94 addParam(
"shiftPt", m_shiftPt,
95 "Shift the Hough plane by 1/4 cell size in 1/r to avoid "
96 "curvature 0 tracks (<0: shift in negative direction, "
97 "0: no shift, >0: shift in positive direction).", 0);
98 addParam(
"minHits", m_minHits,
99 "Minimum hits from different super layers required in a peak cell.",
101 addParam(
"minHitsShort", m_minHitsShort,
102 "Minimum hits required required in a peak cell for a short track"
103 " (must be in the first minHitsShort super layers).",
105 addParam(
"minCells", m_minCells,
106 "Peaks with less than minCells connected cells are ignored.",
108 addParam(
"onlyLocalMax", m_onlyLocalMax,
109 "Switch to remove cells connected to a cell with higher super layer count.",
111 addParam(
"connect", m_connect,
112 "Definition for connected cells. 4: direct (left/right/top/bottom), "
113 "6: direct + 2 diagonal (top right/bottom left), "
114 "8: direct + all diagonal (top right/top left/bottom right/bottom left)",
116 addParam(
"ignore2ndPriority", m_ignore2nd,
117 "Switch to skip second priority hits.",
false);
118 addParam(
"usePriorityPosition", m_usePriority,
119 "If true, use wire position of priority cell in track segment, "
120 "otherwise use wire position of center cell."
121 "Note that the real houghETF firmware will not receive the information of priority position.",
123 addParam(
"requireSL0", m_requireSL0,
124 "Switch to check separately for a hit in the innermost superlayer.",
false);
125 addParam(
"storeHoughPlane", m_storePlane,
126 "Switch for saving Hough plane as TMatrix in DataStore. "
127 "0: don't store anything, 1: store only peaks, 2: store full plane "
128 "(will increase runtime).", (
unsigned)(0));
129 addParam(
"clusterPattern", m_clusterPattern,
130 "use nested pattern algorithm to find clusters",
true);
131 addParam(
"clusterSizeX", m_clusterSizeX,
132 "maximum number of 2 x 2 squares in cluster for pattern algorithm",
134 addParam(
"clusterSizeY", m_clusterSizeY,
135 "maximum number of 2 x 2 squares in cluster for pattern algorithm",
137 addParam(
"hitRelationsFromCorners", m_hitRelationsFromCorners,
138 "Switch for creating relations to hits in the pattern algorithm. "
139 "If true, create relations from cluster corners, otherwise "
140 "from estimated cluster center (might not have relations).",
false);
142 addParam(
"testFilename", m_testFilename,
143 "If not empty, a file with input (hits) and output (tracks) "
144 "for each event is written (for firmware debugging).",
string(
""));
146 addParam(
"suppressClone", m_suppressClone,
147 "Switch to send only the first found track and suppress the "
148 "subsequent clones." ,
false);
152 CDCTriggerHoughETFModule::initialize()
155 m_eventTime.registerInDataStore(m_EventTimeName);
156 m_segmentHits.isRequired(m_hitCollectionName);
159 m_tracks.registerInDataStore(m_outputCollectionName);
160 m_clusters.registerInDataStore(m_clusterCollectionName);
162 m_tracks.registerRelationTo(m_segmentHits);
163 m_tracks.registerRelationTo(m_clusters);
166 if (m_storePlane > 0) m_houghPlane.registerInDataStore(
"HoughPlane");
171 for (
int iSL = 0; iSL < 9; ++iSL) {
173 nTS += cdc.nWiresInLayer(layerId);
174 TSoffset[iSL + 1] = nTS;
175 for (
int priority = 0; priority < 2; ++ priority) {
176 radius[iSL][priority] = cdc.senseWireR(layerId + priority);
178 layerId += (iSL > 0 ? 6 : 7);
181 if (m_testFilename !=
"") {
182 testFile.open(m_testFilename);
185 if (m_suppressClone) {
186 B2INFO(
"2D finder will exit with the first track candidate in time.");
189 if (m_t0CalcMethod != 0 && m_t0CalcMethod != 1) B2WARNING(
"t0CalcMethod is invalid. calculate as default.");
193 CDCTriggerHoughETFModule::event()
198 associatedTSHitsList.clear();
200 if (!m_eventTime.isValid()) m_eventTime.create();
202 setReturnValue(
true);
204 if (m_testFilename !=
"") {
205 testFile << StoreObjPtr<EventMetaData>()->getEvent() <<
" "
206 << m_segmentHits.getEntries() << endl;
209 if (m_segmentHits.getEntries() == 0) {
215 for (
int iHit = 0; iHit < m_segmentHits.getEntries(); iHit++) {
216 unsigned short iSL = m_segmentHits[iHit]->getISuperLayer();
217 if (m_testFilename !=
"") {
218 testFile << iSL <<
" " << m_segmentHits[iHit]->getSegmentID() - TSoffset[iSL] <<
" "
219 << m_segmentHits[iHit]->getPriorityPosition() << endl;
221 if (iSL % 2)
continue;
222 if (m_ignore2nd && m_segmentHits[iHit]->getPriorityPosition() < 3)
continue;
223 double phi = m_segmentHits[iHit]->getSegmentID() - TSoffset[iSL];
225 phi += 0.5 * (((m_segmentHits[iHit]->getPriorityPosition() >> 1) & 1)
226 - (m_segmentHits[iHit]->getPriorityPosition() & 1));
228 phi = phi * 2. * M_PI / (TSoffset[iSL + 1] - TSoffset[iSL]);
229 double r = radius[iSL][int(m_usePriority &&
230 m_segmentHits[iHit]->getPriorityPosition() < 3)];
231 TVector2 pos(cos(phi) / r, sin(phi) / r);
232 hitMap.insert(std::make_pair(iHit, std::make_pair(iSL, pos)));
238 maxIterations = ceil(log2(max(m_nCellsPhi, m_nCellsR))) - 1;
239 nCells = pow(2, maxIterations + 1);
241 double rectX = M_PI * nCells / m_nCellsPhi;
243 maxR = 0.5 * Const::speedOfLight * 1.5e-4 / m_minPt;
244 double rectY = maxR * nCells / m_nCellsR;
247 shiftR = -maxR / 2. / m_nCellsR;
248 }
else if (m_shiftPt > 0) {
249 shiftR = maxR / 2. / m_nCellsR;
252 B2DEBUG(50,
"extending Hough plane to " << maxIterations <<
" iterations, "
253 << nCells <<
" cells: phi in ["
254 << -rectX * 180. / M_PI <<
", " << rectX * 180. / M_PI
255 <<
"] deg, 1/r in [" << -rectY + shiftR <<
", " << rectY + shiftR <<
"] /cm");
258 if (m_storePlane > 0) {
259 m_houghPlane.create();
260 m_houghPlane->ResizeTo(m_nCellsPhi, m_nCellsR);
265 if (m_suppressClone && !hitMap.empty()) {
272 typedef pair<int, cdcPair> cdcHitPair;
274 typedef vector<cdcHitPair> cdcSeqMap;
277 for (
auto hit : hitMap) {
278 seqHitMap.push_back(hit);
280 sort(seqHitMap.begin(), seqHitMap.end(), [
this](cdcHitPair i, cdcHitPair j) {
281 return m_segmentHits[i.first]->foundTime() < m_segmentHits[j.first]->foundTime();
283 auto seqHitItr = seqHitMap.begin();
285 vector<bool> layerHit(CDC_SUPER_LAYERS,
false);
288 short firstTick = m_segmentHits[(*seqHitMap.begin()).first]->foundTime() / period + 1;
289 short lastTick = m_segmentHits[(*(seqHitMap.end() - 1)).first]->foundTime() / period + 1;
291 for (
auto tick = firstTick * period; tick < lastTick * period; tick += period) {
293 for (
auto itr = seqHitItr; itr < seqHitMap.end(); ++itr) {
294 cdcHitPair currentHit = *itr;
296 if (count(layerHit.begin(), layerHit.end(),
true) >= m_minHits &&
297 m_segmentHits[currentHit.first]->foundTime() > tick) {
301 layerHit[m_segmentHits[currentHit.first]->getISuperLayer()] =
true;
303 copy_n(seqHitItr, nHitInCycle, inserter(fastHitMap, fastHitMap.end()));
304 fastInterceptFinder(fastHitMap, -rectX, rectX, -rectY + shiftR, rectY + shiftR, 0, 0, 0);
305 B2DEBUG(20,
"at tick " << tick <<
", number of candidates: " << houghCand.size());
306 if (!houghCand.empty()) {
307 B2DEBUG(10,
"found a track at clock " << tick <<
" with "
308 << fastHitMap.size() <<
"hits");
311 advance(seqHitItr, nHitInCycle);
315 fastInterceptFinder(hitMap, -rectX, rectX, -rectY + shiftR, rectY + shiftR, 0, 0, 0);
316 if (!houghCand.empty()) {
317 B2DEBUG(10,
"found a track with " << hitMap.size() <<
"hits");
322 if (m_clusterPattern) {
323 if (patternClustering(fastHitMap))
324 m_eventTime->addBinnedEventT0(calcEventTiming(), Const::CDC);
329 if (m_testFilename !=
"") {
330 testFile << m_tracks.getEntries() << endl;
331 for (
int i = 0; i < m_tracks.getEntries(); ++i) {
332 float ix = (m_tracks[i]->getPhi0() - M_PI_4) * m_nCellsPhi / 2. / M_PI - 0.5;
333 float iy = (m_tracks[i]->getOmega() / 2. + maxR - shiftR) * m_nCellsR / 2. / maxR - 0.5;
334 testFile << round(2 * ix) / 2. <<
" " << round(2 * iy) / 2. <<
" "
335 << m_tracks[i]->getChargeSign() << endl;
338 testFile << hits.size() << endl;
339 for (
unsigned ihit = 0; ihit < hits.size(); ++ihit) {
340 unsigned short iSL = hits[ihit]->getISuperLayer();
341 testFile << iSL <<
" " << hits[ihit]->getSegmentID() - TSoffset[iSL] <<
" "
342 << hits[ihit]->getPriorityPosition() <<
" "
343 << hits.weight(ihit) << endl;
350 CDCTriggerHoughETFModule::terminate()
352 if (m_testFilename !=
"") testFile.close();