Belle II Software  release-05-02-19
CDCTriggerHoughETFModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2012-2014 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Michael Schnell, Sara Neuhaus *
7  * *
8  **************************************************************************/
9 
10 #include <trg/cdc/modules/houghETF/CDCTriggerHoughETFModule.h>
11 
12 #include <cdc/geometry/CDCGeometryPar.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/gearbox/Const.h>
15 
16 #include <cmath>
17 #include <algorithm>
18 #include <iterator>
19 
20 #include <root/TMatrix.h>
21 
22 /* defines */
23 #define CDC_SUPER_LAYERS 9
24 
25 using namespace std;
26 using namespace Belle2;
27 using namespace Belle2::CDC;
28 
29 //-----------------------------------------------------------------
30 // Register the Module
31 //-----------------------------------------------------------------
32 REG_MODULE(CDCTriggerHoughETF)
33 
34 //-----------------------------------------------------------------
35 // Implementation
36 //-----------------------------------------------------------------
37 
39 {
40  //Set module properties
41  setDescription("Hough tracking algorithm for CDC trigger.");
42  setPropertyFlags(c_ParallelProcessingCertified);
43 
44  // Define module parameters
45  // 2D ETF specified parameters
46  // ----- PARAMETERS : GRL t0 finder with fastest priority timing (current at 2020.2.7)-----
47  // t0CalcMethod = 0
48  // usePriorityPosition = true
49  // arrivalOrder = 0
50  //
51  addParam("outputEventTimeName", m_EventTimeName,
52  "Name of the output StoreObjPtr.",
53  string(""));
54  addParam("storeTracks", m_storeTracks,
55  "store tracks",
56  false);
57  addParam("usePriorityTiming", m_usePriorityTiming,
58  "Use priority timing instead of fastest timing.",
59  true);
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.",
64  (unsigned)(0));
65  addParam("arrivalOrder", m_arrivalOrder,
66  "When t0CalcMethod == 0: Nth fastest ft is used as T0. (i.e. 0 is fastest)",
67  (unsigned)(0));
68  addParam("timeWindowStart", m_timeWindowStart,
69  "When t0CalcMethod == 2: start time of time window relative to median. (in ns)",
70  (short)(-64));
71  addParam("timeWindowEnd", m_timeWindowEnd,
72  "When t0CalcMethod == 2: end time of time window relative to median. (in ns)",
73  (short)(0));
74 
75  //common as CDCTrigger2DFinderModule
76  addParam("hitCollectionName", m_hitCollectionName,
77  "Name of the input StoreArray of CDCTriggerSegmentHits.",
78  string(""));
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.",
84  string(""));
85  addParam("nCellsPhi", m_nCellsPhi,
86  "Number of Hough cells in phi (limits: [-180, 180]). Must be an even number.",
87  (unsigned)(160));
88  addParam("nCellsR", m_nCellsR,
89  "Number of Hough cells in 1/r. Must be an even number.",
90  (unsigned)(34));
91  addParam("minPt", m_minPt,
92  "Minimum Pt [GeV]. "
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.",
100  (unsigned)(4));
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).",
104  (unsigned)(4));
105  addParam("minCells", m_minCells,
106  "Peaks with less than minCells connected cells are ignored.",
107  (unsigned)(2));
108  addParam("onlyLocalMax", m_onlyLocalMax,
109  "Switch to remove cells connected to a cell with higher super layer count.",
110  false);
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)",
115  (unsigned)(6));
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.",
122  false);
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",
133  (unsigned)(3));
134  addParam("clusterSizeY", m_clusterSizeY,
135  "maximum number of 2 x 2 squares in cluster for pattern algorithm",
136  (unsigned)(3));
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);
141 
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(""));
145 
146  addParam("suppressClone", m_suppressClone,
147  "Switch to send only the first found track and suppress the "
148  "subsequent clones." , false);
149 }
150 
151 void
152 CDCTriggerHoughETFModule::initialize()
153 {
154  // register DataStore elements
155  m_eventTime.registerInDataStore(m_EventTimeName);
156  m_segmentHits.isRequired(m_hitCollectionName);
157 
158  if (m_storeTracks) {
159  m_tracks.registerInDataStore(m_outputCollectionName);
160  m_clusters.registerInDataStore(m_clusterCollectionName);
161 
162  m_tracks.registerRelationTo(m_segmentHits);
163  m_tracks.registerRelationTo(m_clusters);
164  }
165 
166  if (m_storePlane > 0) m_houghPlane.registerInDataStore("HoughPlane");
167 
168  CDCGeometryPar& cdc = CDCGeometryPar::Instance();
169  int layerId = 3;
170  int nTS = 0;
171  for (int iSL = 0; iSL < 9; ++iSL) {
172  TSoffset[iSL] = nTS;
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);
177  }
178  layerId += (iSL > 0 ? 6 : 7);
179  }
180 
181  if (m_testFilename != "") {
182  testFile.open(m_testFilename);
183  }
184 
185  if (m_suppressClone) {
186  B2INFO("2D finder will exit with the first track candidate in time.");
187  }
188 
189  if (m_t0CalcMethod != 0 && m_t0CalcMethod != 1) B2WARNING("t0CalcMethod is invalid. calculate as default.");
190 }
191 
192 void
193 CDCTriggerHoughETFModule::event()
194 {
195  /* Clean hits */
196  hitMap.clear();
197  houghCand.clear();
198  associatedTSHitsList.clear();
199 
200  if (!m_eventTime.isValid()) m_eventTime.create();
201  /* set default return value */
202  setReturnValue(true);
203 
204  if (m_testFilename != "") {
205  testFile << StoreObjPtr<EventMetaData>()->getEvent() << " "
206  << m_segmentHits.getEntries() << endl;
207  }
208 
209  if (m_segmentHits.getEntries() == 0) {
210  //B2WARNING("CDCTracking: tsHitsCollection is empty!");
211  return;
212  }
213 
214  /* get CDCHits coordinates in conformal space */
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;
220  }
221  if (iSL % 2) continue; //ysue
222  if (m_ignore2nd && m_segmentHits[iHit]->getPriorityPosition() < 3) continue;
223  double phi = m_segmentHits[iHit]->getSegmentID() - TSoffset[iSL];
224  if (m_usePriority) {
225  phi += 0.5 * (((m_segmentHits[iHit]->getPriorityPosition() >> 1) & 1)
226  - (m_segmentHits[iHit]->getPriorityPosition() & 1));
227  }
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)));
233  }
234 
235  /* Extent the Hough plane such that the cell number is a power of 2 (same for x and y).
236  * This is for the fast peak finder, which divides the Hough plane in half in each step.
237  * Peaks found outside of the actual limits are ignored. */
238  maxIterations = ceil(log2(max(m_nCellsPhi, m_nCellsR))) - 1;
239  nCells = pow(2, maxIterations + 1);
240  /* limits in phi: [-pi, pi] + extra cells */
241  double rectX = M_PI * nCells / m_nCellsPhi;
242  /* limits in R: [-R(minPt), R(minPt)] + extra cells + shift */
243  maxR = 0.5 * Const::speedOfLight * 1.5e-4 / m_minPt;
244  double rectY = maxR * nCells / m_nCellsR;
245  shiftR = 0;
246  if (m_shiftPt < 0) {
247  shiftR = -maxR / 2. / m_nCellsR;
248  } else if (m_shiftPt > 0) {
249  shiftR = maxR / 2. / m_nCellsR;
250  }
251 
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");
256 
257  /* prepare matrix for storing the Hough plane */
258  if (m_storePlane > 0) {
259  m_houghPlane.create();
260  m_houghPlane->ResizeTo(m_nCellsPhi, m_nCellsR);
261  }
262 
263  // hit map containing only the early hits
264  cdcMap fastHitMap;
265  if (m_suppressClone && !hitMap.empty()) {
266  // find the first track candidates in Hough plane
267  // only for z0 resolution study with single-track events
268  // This will surely fail with multi-track ones,
269  // in which case we really need tick-by-tick simulation for all hits.
270 
272  typedef pair<int, cdcPair> cdcHitPair;
273  // sequential hit map, ordered by TS found time
274  typedef vector<cdcHitPair> cdcSeqMap;
275  cdcSeqMap seqHitMap;
276  // copy hit map to sequential hit map and sort it by TS found time
277  for (auto hit : hitMap) {
278  seqHitMap.push_back(hit);
279  }
280  sort(seqHitMap.begin(), seqHitMap.end(), [this](cdcHitPair i, cdcHitPair j) {
281  return m_segmentHits[i.first]->foundTime() < m_segmentHits[j.first]->foundTime();
282  });
283  auto seqHitItr = seqHitMap.begin();
284  /* layer filter */
285  vector<bool> layerHit(CDC_SUPER_LAYERS, false);
286  // data clock cycle in unit of 2ns
287  short period = 16;
288  short firstTick = m_segmentHits[(*seqHitMap.begin()).first]->foundTime() / period + 1;
289  short lastTick = m_segmentHits[(*(seqHitMap.end() - 1)).first]->foundTime() / period + 1;
290  // add TS hits in every clock cycle until a track candidate is found
291  for (auto tick = firstTick * period; tick < lastTick * period; tick += period) {
292  int nHitInCycle = 0;
293  for (auto itr = seqHitItr; itr < seqHitMap.end(); ++itr) {
294  cdcHitPair currentHit = *itr;
295  // start from the first hit over SL threshold
296  if (count(layerHit.begin(), layerHit.end(), true) >= m_minHits &&
297  m_segmentHits[currentHit.first]->foundTime() > tick) {
298  break;
299  }
300  nHitInCycle++;
301  layerHit[m_segmentHits[currentHit.first]->getISuperLayer()] = true;
302  }
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");
309  break;
310  }
311  advance(seqHitItr, nHitInCycle);
312  }
313  } else {
314  /* find track candidates in Hough plane using all TS hits */
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");
318  }
319  }
320 
321  /* merge track candidates */
322  if (m_clusterPattern) {
323  if (patternClustering(fastHitMap))
324  m_eventTime->addBinnedEventT0(calcEventTiming(), Const::CDC);
325  } else {
326  connectedRegions();
327  }
328 
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;
337  m_tracks[i]->getRelationsTo<CDCTriggerSegmentHit>(m_hitCollectionName);
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;
344  }
345  }
346  }
347 }
348 
349 void
350 CDCTriggerHoughETFModule::terminate()
351 {
352  if (m_testFilename != "") testFile.close();
353 }
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::CDCTriggerHoughETFModule
Definition: CDCTriggerHoughETFModule.h:94
Belle2::cdcMap
std::map< int, cdcPair > cdcMap
Map of <counter, cdcPair>, for hits with indices.
Definition: CDCTriggerHoughETFModule.h:44
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::CDC::CDCGeometryPar
The Class for CDC Geometry Parameters.
Definition: CDCGeometryPar.h:75
Belle2::RelationVector
Class for type safe access to objects that are referred to in relations.
Definition: DataStore.h:38
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDCTriggerSegmentHit
Combination of several CDCHits to a track segment hit for the trigger.
Definition: CDCTriggerSegmentHit.h:16
Belle2::CDC
Definition: CrudeT0CalibrationAlgorithm.h:28