Belle II Software  release-05-02-19
CDCTrigger2DFinderModule.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/houghtracking/CDCTrigger2DFinderModule.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 /* defines */
21 #define CDC_SUPER_LAYERS 9
22 
23 using namespace std;
24 using namespace Belle2;
25 using namespace Belle2::CDC;
26 
27 //-----------------------------------------------------------------
28 // Register the Module
29 //-----------------------------------------------------------------
30 REG_MODULE(CDCTrigger2DFinder)
31 
32 //-----------------------------------------------------------------
33 // Implementation
34 //-----------------------------------------------------------------
35 
37 {
38  //Set module properties
39  setDescription("Hough tracking algorithm for CDC trigger.");
40  setPropertyFlags(c_ParallelProcessingCertified);
41 
42  // Define module parameters
43  addParam("hitCollectionName", m_hitCollectionName,
44  "Name of the input StoreArray of CDCTriggerSegmentHits.",
45  string(""));
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.",
51  string(""));
52  addParam("nCellsPhi", m_nCellsPhi,
53  "Number of Hough cells in phi (limits: [-180, 180]). Must be an even number.",
54  (unsigned)(160));
55  addParam("nCellsR", m_nCellsR,
56  "Number of Hough cells in 1/r. Must be an even number.",
57  (unsigned)(34));
58  addParam("minPt", m_minPt,
59  "Minimum Pt [GeV]. "
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);
65 
66  addParam("minHits", m_minHits,
67  "Minimum hits from different super layers required in a peak cell.",
68  (unsigned)(4));
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).",
72  (unsigned)(4));
73  addParam("minCells", m_minCells,
74  "Peaks with less than minCells connected cells are ignored.",
75  (unsigned)(2));
76  addParam("onlyLocalMax", m_onlyLocalMax,
77  "Switch to remove cells connected to a cell with higher super layer count.",
78  false);
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)",
83  (unsigned)(6));
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",
99  (unsigned)(3));
100  addParam("clusterSizeY", m_clusterSizeY,
101  "maximum number of 2 x 2 squares in cluster for pattern algorithm",
102  (unsigned)(3));
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);
107 
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(""));
111 
112  addParam("suppressClone", m_suppressClone,
113  "Switch to send only the first found track and suppress the "
114  "subsequent clones." , false);
115 }
116 
117 void
118 CDCTrigger2DFinderModule::initialize()
119 {
120  m_segmentHits.isRequired(m_hitCollectionName);
121  m_tracks.registerInDataStore(m_outputCollectionName);
122  m_clusters.registerInDataStore(m_clusterCollectionName);
123 
124  m_tracks.registerRelationTo(m_segmentHits);
125  m_tracks.registerRelationTo(m_clusters);
126 
127  if (m_storePlane > 0) m_houghPlane.registerInDataStore("HoughPlane");
128 
129  CDCGeometryPar& cdc = CDCGeometryPar::Instance();
130  int layerId = 3;
131  int nTS = 0;
132  for (int iSL = 0; iSL < 9; ++iSL) {
133  TSoffset[iSL] = nTS;
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);
138  }
139  layerId += (iSL > 0 ? 6 : 7);
140  }
141 
142  if (m_testFilename != "") {
143  testFile.open(m_testFilename);
144  }
145 
146  if (m_suppressClone) {
147  B2INFO("2D finder will exit with the first track candidate in time.");
148  }
149 }
150 
151 void
152 CDCTrigger2DFinderModule::event()
153 {
154  /* Clean hits */
155  hitMap.clear();
156  houghCand.clear();
157 
158  /* set default return value */
159  setReturnValue(true);
160 
161  if (m_testFilename != "") {
162  testFile << StoreObjPtr<EventMetaData>()->getEvent() << " "
163  << m_segmentHits.getEntries() << endl;
164  }
165 
166  if (m_segmentHits.getEntries() == 0) {
167  //B2WARNING("CDCTracking: tsHitsCollection is empty!");
168  return;
169  }
170 
171  /* get CDCHits coordinates in conformal space */
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;
177  }
178  if (iSL % 2) continue;
179  if (m_ignore2nd && m_segmentHits[iHit]->getPriorityPosition() < 3) continue;
180  double phi = m_segmentHits[iHit]->getSegmentID() - TSoffset[iSL];
181  if (m_usePriority) {
182  phi += 0.5 * (((m_segmentHits[iHit]->getPriorityPosition() >> 1) & 1)
183  - (m_segmentHits[iHit]->getPriorityPosition() & 1));
184  }
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)));
190  }
191 
192  /* Extent the Hough plane such that the cell number is a power of 2 (same for x and y).
193  * This is for the fast peak finder, which divides the Hough plane in half in each step.
194  * Peaks found outside of the actual limits are ignored. */
195  maxIterations = ceil(log2(max(m_nCellsPhi, m_nCellsR))) - 1;
196  nCells = pow(2, maxIterations + 1);
197  /* limits in phi: [-pi, pi] + extra cells */
198  double rectX = M_PI * nCells / m_nCellsPhi;
199  /* limits in R: [-R(minPt), R(minPt)] + extra cells + shift */
200  maxR = 0.5 * Const::speedOfLight * 1.5e-4 / m_minPt;
201  double rectY = maxR * nCells / m_nCellsR;
202  shiftR = 0;
203  if (m_shiftPt < 0) {
204  shiftR = -maxR / 2. / m_nCellsR;
205  } else if (m_shiftPt > 0) {
206  shiftR = maxR / 2. / m_nCellsR;
207  }
208 
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");
213 
214  /* prepare matrix for storing the Hough plane */
215  if (m_storePlane > 0) {
216  m_houghPlane.create();
217  m_houghPlane->ResizeTo(m_nCellsPhi, m_nCellsR);
218  }
219 
220  // hit map containing only the early hits
221  cdcMap fastHitMap;
222  if (m_suppressClone && !hitMap.empty()) {
223  // find the first track candidates in Hough plane
224  // only for z0 resolution study with single-track events
225  // This will surely fail with multi-track ones,
226  // in which case we really need tick-by-tick simulation for all hits.
227 
229  typedef pair<int, cdcPair> cdcHitPair;
230  // sequential hit map, ordered by TS found time
231  typedef vector<cdcHitPair> cdcSeqMap;
232  cdcSeqMap seqHitMap;
233  // copy hit map to sequential hit map and sort it by TS found time
234  for (auto hit : hitMap) {
235  seqHitMap.push_back(hit);
236  }
237  sort(seqHitMap.begin(), seqHitMap.end(), [this](cdcHitPair i, cdcHitPair j) {
238  return m_segmentHits[i.first]->foundTime() < m_segmentHits[j.first]->foundTime();
239  });
240  auto seqHitItr = seqHitMap.begin();
241  /* layer filter */
242  vector<bool> layerHit(CDC_SUPER_LAYERS, false);
243  // data clock cycle in unit of 2ns
244  short period = 16;
245  short firstTick = m_segmentHits[(*seqHitMap.begin()).first]->foundTime() / period + 1;
246  short lastTick = m_segmentHits[(*(seqHitMap.end() - 1)).first]->foundTime() / period + 1;
247  // add TS hits in every clock cycle until a track candidate is found
248  for (auto tick = firstTick * period; tick < lastTick * period; tick += period) {
249  int nHitInCycle = 0;
250  for (auto itr = seqHitItr; itr < seqHitMap.end(); ++itr) {
251  cdcHitPair currentHit = *itr;
252  // start from the first hit over SL threshold
253  if (count(layerHit.begin(), layerHit.end(), true) >= m_minHits &&
254  m_segmentHits[currentHit.first]->foundTime() > tick) {
255  break;
256  }
257  nHitInCycle++;
258  layerHit[m_segmentHits[currentHit.first]->getISuperLayer()] = true;
259  }
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");
266  break;
267  }
268  advance(seqHitItr, nHitInCycle);
269  }
270  } else {
271  /* find track candidates in Hough plane using all TS hits */
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");
275  }
276  }
277 
278  /* merge track candidates */
279  if (m_clusterPattern)
280  patternClustering(fastHitMap);
281  else
282  connectedRegions();
283 
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;
292  m_tracks[i]->getRelationsTo<CDCTriggerSegmentHit>(m_hitCollectionName);
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;
299  }
300  }
301  }
302 }
303 
304 void
305 CDCTrigger2DFinderModule::terminate()
306 {
307  if (m_testFilename != "") testFile.close();
308 }
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
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::CDCTrigger2DFinderModule
Definition: CDCTrigger2DFinderModule.h:93
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