Belle II Software  release-08-01-10
NDFinder.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <string>
10 #include <cmath>
11 #include <iostream>
12 #include <fstream>
13 #include "framework/logging/Logger.h"
14 #include "boost/iostreams/filter/gzip.hpp"
15 #include "boost/iostreams/filtering_streambuf.hpp"
16 #include "boost/iostreams/filtering_stream.hpp"
17 #include "trg/cdc/NDFinder.h"
18 
19 #include <TMath.h>
20 
21 using namespace Belle2;
22 using namespace std;
23 
24 
25 void NDFinder::init(int minweight, int minpts,
26  bool diagonal, int minhits, int minhits_axial,
27  double thresh,
28  double minassign,
29  int mincells, bool verbose,
30  string& axialFile, string& stereoFile)
31 {
32  m_params.minhits = (long unsigned int) minhits;
33  m_params.minhits_axial = (long unsigned int) minhits_axial;
34  m_params.thresh = (float) thresh;
35  m_params.minassign = (float) minassign;
36  m_params.mincells = (long unsigned int) mincells;
37  m_params.axialFile = axialFile;
38  m_params.stereoFile = stereoFile;
39  m_clusterer_params.minweight = (unsigned short) minweight;
40  m_clusterer_params.minpts = (unsigned short) minpts;
41  m_clusterer_params.diagonal = diagonal;
42  m_verbose = verbose;
43 
44 
45  initBins();
46  B2DEBUG(55, "initialized binnings");
47 
48  initHitModAxSt(*m_parrayHitMod);
49  B2DEBUG(55, "initialized HitMod, a map of tsid to (orient, relid, letter).");
50 
52  loadArray(m_params.axialFile, m_compaxbins, *m_pcompAxial);
53  B2DEBUG(55, "loaded zero suppressed axial array ");
54  loadArray(m_params.stereoFile, m_compstbins, *m_pcompStereo);
55  B2DEBUG(55, "loaded zero suppressed stereo array ");
56 
58  restoreZeros(m_expaxbins, m_compaxbins, *m_parrayAxialExp, *m_pcompAxial);
59  B2DEBUG(55, "restored expanded axial array (11/32 phi) ");
60  restoreZeros(m_expstbins, m_compstbins, *m_parrayStereoExp, *m_pcompStereo);
61  B2DEBUG(55, "restored expanded stereo array (11/32 phi) ");
62  squeezeAll(m_axbins, *m_parrayAxial, *m_parrayAxialExp, m_params.parcels, m_params.parcelsExp);
63  B2DEBUG(55, "squeezed axial array (11/32 phi) --> (7/32 phi): ");
64  squeezeAll(m_stbins, *m_parrayStereo, *m_parrayStereoExp, m_params.parcels, m_params.parcelsExp);
65  B2DEBUG(55, "squeezed stereo array (11/32 phi) --> (7/32 phi)");
66 
67  reset();
68  m_clusterer = Clusterizend(m_clusterer_params);
69  m_clusterer.setPlaneShape(m_planeShape);
70 }
71 
73 {
84  m_nTS = 2336;
85  m_nSL = 9;
86  m_params.phigeo = 32;
87 
88  m_nAx = 41;
89  m_nSt = 32;
90  m_nPrio = 3;
91 
92  m_nOmega = 40;
93  m_nPhiFull = 384;
94  m_nTheta = 9;
95  m_nPhiOne = m_nPhiFull / m_params.phigeo; // 384/32 = 12
97  m_nPhiComp = 15;
99  m_nPhiUse = m_params.parcels * m_nPhiOne;
101  m_nPhiExp = m_params.parcelsExp * m_nPhiOne;
102 
103  m_axbins.hitid = m_nAx;
104  m_stbins.hitid = m_nSt;
105  m_axbins.phi = m_nPhiUse; //84
106  m_stbins.phi = m_nPhiUse; //84
107  m_compaxbins.hitid = m_nAx;
108  m_compstbins.hitid = m_nSt;
109  m_compaxbins.phi = m_nPhiComp;
110  m_compstbins.phi = m_nPhiComp;
111  m_compaxbins.theta = 1;
112  m_expaxbins.hitid = m_nAx;
113  m_expstbins.hitid = m_nSt;
114  m_expaxbins.phi = m_nPhiExp; //132
115  m_expstbins.phi = m_nPhiExp; //132
116 
117  m_fullbins.hitid = m_nTS;
118  m_fullbins.phi = m_nPhiFull;
119 
120  m_pc5shapeax = {{ m_nAx, m_nPrio, m_nOmega, m_nPhiUse, m_nTheta }};
121  m_pc5shapest = {{ m_nSt, m_nPrio, m_nOmega, m_nPhiUse, m_nTheta }};
122  m_pc3shape = {{ m_nOmega, m_nPhiFull, m_nTheta }};
123  m_pc2shapeHitMod = {{ m_nTS, m_nPrio}};
124  m_pc5shapeCompAx = {{ m_nAx, m_nPrio, m_nOmega, m_nPhiComp, 1 }};
125  m_pc5shapeCompSt = {{ m_nSt, m_nPrio, m_nOmega, m_nPhiComp, m_nTheta }};
126  m_pc5shapeExpAx = {{ m_nAx, m_nPrio, m_nOmega, m_nPhiExp, m_nTheta }};
127  m_pc5shapeExpSt = {{ m_nSt, m_nPrio, m_nOmega, m_nPhiExp, m_nTheta }};
128 
129  m_parrayAxial = new c5array(m_pc5shapeax);
130  m_parrayStereo = new c5array(m_pc5shapest);
131  m_phoughPlane = new c3array(m_pc3shape);
132  m_parrayHitMod = new c2array(m_pc2shapeHitMod);
133  m_pcompAxial = new c5array(m_pc5shapeCompAx);
134  m_pcompStereo = new c5array(m_pc5shapeCompSt);
135  m_parrayAxialExp = new c5array(m_pc5shapeExpAx);
136  m_parrayStereoExp = new c5array(m_pc5shapeExpSt);
137 
138  m_nWires = {160, 160, 192, 224, 256, 288, 320, 352, 384};
139 
140  m_planeShape = {m_nOmega, m_nPhiFull, m_nTheta}; //{40, 384, 9};
141 
143  vector<float> omegaRange = { -5., 5.};
144  vector<float> phiRange = {0., 11.25};
145  vector<float> thetaRange = {19., 140.};
146  float ssOmega = (omegaRange[1] - omegaRange[0]) / m_nOmega; //40;
147  float ssPhi = (phiRange[1] - phiRange[0]) / m_nPhiOne; //12;
148  float ssTheta = (thetaRange[1] - thetaRange[0]) / m_nTheta; //9;
149  m_acceptRanges.push_back(omegaRange);
150  m_acceptRanges.push_back(phiRange);
151  m_acceptRanges.push_back(thetaRange);
152  m_slotsizes.push_back(ssOmega);
153  m_slotsizes.push_back(ssPhi);
154  m_slotsizes.push_back(ssTheta);
155 }
156 
157 
158 
159 void NDFinder::loadArray(const std::string& filename, ndbinning bins, c5array& hitsToTracks)
160 {
161  vector<c5elem> flatArray;
162  ifstream arrayFileGZ(filename, ios_base::in | ios_base::binary);
163  boost::iostreams::filtering_istream arrayStream;
164  arrayStream.push(boost::iostreams::gzip_decompressor());
165  arrayStream.push(arrayFileGZ);
166  c5elem uline;
167  if (arrayFileGZ.is_open()) {
168  while (arrayStream >> uline) {
169  flatArray.push_back(uline);
170  }
171  arrayFileGZ.close();
172  } else {
173  B2ERROR("could not open array file: " << filename);
174  }
175  B2DEBUG(55, "loaded array from file " << filename);
176  unsigned long icount = 0;
177  for (c5index ihit = 0; ihit < bins.hitid; ihit++) {
178  for (c5index iprio = 0; iprio < bins.prio; iprio++) {
179  for (c5index iomega = 0; iomega < bins.omega; iomega++) {
180  for (c5index iphi = 0; iphi < bins.phi; iphi++) {
181  for (c5index itheta = 0; itheta < bins.theta; itheta++) {
182  hitsToTracks[ihit][iprio][iomega][iphi][itheta] = flatArray[icount];
183  icount++;
184  }
185  }
186  }
187  }
188  }
189 }
190 
191 
192 void NDFinder::initHitModAxSt(c2array& hitMod)
193 {
194  unsigned short phigeo = m_params.phigeo; // 32
195  std::vector<int> modSLs; //nWires per SL in 1/32 sector (9)
196  std::vector<int> toslid; //tsid to sl map (2336)
197  std::vector<int> sloffsets = {0}; // integrated number of wires in inner SLs (9);
198  std::vector<int> modoffsetsAxSt = {0, 0}; // integrated number of wires in inner SLs in 1/32 sector, separate for axial and stereo SLs (9);
199  for (int isl = 0; isl < m_nSL; isl++) {
200  modSLs.push_back(m_nWires[isl] / phigeo);
201  for (int itsrel = 0; itsrel < m_nWires[isl]; itsrel++) {
202  toslid.push_back(isl);
203  }
204  sloffsets.push_back(sloffsets[isl] + m_nWires[isl]);
205  int last = modoffsetsAxSt[isl];
206  modoffsetsAxSt.push_back(last + m_nWires[isl] / phigeo);
207  }
208  for (int its = 0; its < m_nTS; its++) { // 2336
209  int sl = toslid[its];
210  bool isAxial = (sl % 2 == 0);
211  int idInSL = its - sloffsets[sl];
212  int modSL = modSLs[sl];
213  int modId = idInSL % modSL;
214  int letter = (int) floor(idInSL / modSL);
215  int relId = (int)(modoffsetsAxSt[sl] + modId);
216  hitMod[its][0] = (int)(isAxial);
217  hitMod[its][1] = relId;
218  hitMod[its][2] = letter;
219  }
220 }
221 
222 
223 void NDFinder::squeezeOne(c5array& writeArray, c5array& readArray, int outparcels, int inparcels, c5index ihit, c5index iprio,
224  c5index itheta, c5elem nomega)
225 {
226  int outnphi = (int)(12 * outparcels);
227  c5index trafstart = (c5index)((inparcels - outparcels) / 2 * 12);
228  c5index trafend = (c5index)(trafstart + outnphi);
229  for (c5index iomega = 0; iomega < nomega; iomega++) {
230  for (c5index iphi = 0; iphi < outnphi; iphi++) {
231  c5index readPhi = trafstart + iphi;
232  writeArray[ihit][iprio][iomega][iphi][itheta] = readArray[ihit][iprio][iomega][readPhi][itheta];
233  }
234  for (c5index iphi = 0; iphi < trafstart; iphi++) {
235  c5index writePhi = (c5index)(outnphi - trafstart + iphi);
236  writeArray[ihit][iprio][iomega][writePhi][itheta] += readArray[ihit][iprio][iomega][iphi][itheta];
237  }
238  for (c5index iphi = 0; iphi < trafstart; iphi++) {
239  c5index readPhi = trafend + iphi;
240  writeArray[ihit][iprio][iomega][iphi][itheta] += readArray[ihit][iprio][iomega][readPhi][itheta];
241  }
242  }
243 }
244 
245 
246 void NDFinder::squeezeAll(ndbinning writebins, c5array& writeArray, c5array& readArray, int outparcels, int inparcels)
247 {
248  for (c5index ihit = 0; ihit < writebins.hitid; ihit++) {
249  for (c5index iprio = 0; iprio < writebins.prio; iprio++) {
250  for (c5index itheta = 0; itheta < writebins.theta; itheta++) {
251  squeezeOne(writeArray, readArray, outparcels, inparcels, ihit, iprio, itheta, writebins.omega);
252  }
253  }
254  }
255 }
256 
257 
258 void NDFinder::restoreZeros(ndbinning zerobins, ndbinning compbins, c5array& expArray, const c5array& compArray)
259 {
260  B2DEBUG(55, "restoreZeros: zerobins.theta " << zerobins.theta << ", combins.theta " << compbins.theta);
261  for (c5index ihit = 0; ihit < compbins.hitid; ihit++) {
262  for (c5index iprio = 0; iprio < compbins.prio; iprio++) {
263  for (c5index iomega = 0; iomega < compbins.omega; iomega++) {
264  for (c5index itheta = 0; itheta < compbins.theta; itheta++) {
265  c5elem phiStart = compArray[ihit][iprio][iomega][0][itheta];
266  c5elem phiWidth = compArray[ihit][iprio][iomega][1][itheta];
267  for (c5index iphi = 0; iphi < phiWidth; iphi++) {
268  c5elem phiCur = phiStart + iphi;
269  expArray[ihit][iprio][iomega][phiCur][itheta] = compArray[ihit][iprio][iomega][iphi + 2][itheta];
270  if (compbins.theta == 1) { // case axial, expand in theta
271  for (c5index jtheta = 1; jtheta < zerobins.theta; jtheta++) {
272  expArray[ihit][iprio][iomega][phiCur][jtheta] = compArray[ihit][iprio][iomega][iphi + 2][itheta];
273  }
274  }
275  }
276  }
277  }
278  }
279  }
280 }
281 
282 
283 void NDFinder::printArray3D(c3array& hitsToTracks, ndbinning bins, ushort phi_inc = 1, ushort om_inc = 1, ushort divide = 4,
284  ushort minweightShow = 1)
285 {
286  ushort phistart = 0;
287  ushort phiend = bins.phi;
288  bool started = false;
289  unsigned long phiSlice = 0;
290  for (c3index iphi = 0; iphi < bins.phi; iphi++) {
291  for (c3index itheta = 0; itheta < bins.theta; itheta++) {
292  for (c3index iomega = 0; iomega < bins.omega; iomega++) {
293  phiSlice += hitsToTracks[iomega][iphi][itheta];
294  }
295  }
296  if (phiSlice == 0 and not started) {
297  phistart = iphi + 1;
298  } else if (phiSlice != 0 and not started) {
299  started = true;
300  } else if (phiSlice != 0 and started) {
301  phiend = iphi;
302  } else if (phiSlice == 0 and started) {
303  phiend = iphi;
304  break;
305  }
306  }
307 
308  B2DEBUG(55, "printArray3D, phistart = " << phistart << ", phiend = " << phiend);
309  auto d3shp = hitsToTracks.shape();
310  B2DEBUG(55, "printArray shape: " << d3shp[0] << ", " << d3shp[1] << ", " << d3shp[2]);
311 
312  if (phiend == phistart) {
313  return;
314  }
315  c3index itheta = 4; // only print itheta = 4
316  for (c3index iomega = 0; iomega < bins.omega; iomega += om_inc) {
317  for (c3index iphi = phistart; iphi < phiend; iphi += phi_inc) {
318  // reduce printed weight
319  ushort valRed = (ushort)((hitsToTracks[iomega][iphi][itheta]) / divide);
320  if (valRed <= minweightShow) // skip small weights
321  cout << " ";
322  else
323  cout << valRed;
324  }
325  cout << "|" << endl;
326  }
327 }
328 
329 
331 void NDFinder::addLookup(unsigned short ihit)
332 {
333  c2array& arrayHitMod = *m_parrayHitMod;
334 
335  c2elem orient = arrayHitMod[m_hitIds[ihit]][0];
336  c2elem hitr = arrayHitMod[m_hitIds[ihit]][1];
337  c2elem letter = arrayHitMod[m_hitIds[ihit]][2];
338 
339  unsigned short prio = m_prioPos[ ihit ];
340  short letterShort = (short) letter;
341 
342  // Get hit contribution to cluster:
343  // 7 of 32 phi parcels: center = 3
344  short DstartShort = (letterShort - 3) % m_params.phigeo * m_nPhiOne;
345  if (DstartShort < 0) {DstartShort = m_nPhiFull + DstartShort;}
346  // Add hit to hough plane
347  // 11 of 32 phi parcels: center = 5
348  short DstartComp = (letterShort - 5) % m_params.phigeo * m_nPhiOne;
349  if (DstartComp < 0) {DstartComp = m_nPhiFull + DstartComp;}
350 
351  m_vecDstart.push_back(DstartShort);
352  m_hitOrients.push_back(orient);
353 
354  if (orient == 1) {
355  addC3Comp(hitr, prio, *m_pcompAxial, DstartComp, m_compaxbins);
356  } else {
357  addC3Comp(hitr, prio, *m_pcompStereo, DstartComp, m_compstbins);
358  }
359 }
360 
361 
362 void NDFinder::addC3Comp(ushort hitr,
363  ushort prio, const c5array& hitsToTracks,
364  short Dstart, ndbinning bins)
365 {
366  ushort ntheta = 0;
367  c3array& houghPlane = *m_phoughPlane;
368  for (ushort itheta = 0; itheta < m_nTheta; itheta++) { //9
369  if (bins.theta > 1) { //stereo
370  ntheta = itheta;
371  }
372  for (ushort iomega = 0; iomega < m_nOmega; iomega++) { //40
373  ushort startfield = 0;
374  ushort lenfield = 1;
375  ushort xstart = hitsToTracks[hitr][prio][iomega][startfield][ntheta];
376  ushort xlen = hitsToTracks[hitr][prio][iomega][lenfield][ntheta];
377  for (ushort iphiz = 0; iphiz < xlen; iphiz++) {
378  ushort iphix = iphiz + 2;
379  ushort iphi = iphiz + xstart;
380  ushort iHoughPhi = (iphi + Dstart) % m_nPhiFull;
381  houghPlane[iomega][iHoughPhi][itheta] += hitsToTracks[hitr][prio][iomega][iphix][ntheta];
382  }
383  }
384  }
385 }
386 
387 
389 {
391  for (unsigned short ihit = 0; ihit < m_nHits; ihit++) {
392  addLookup(ihit);
393  }
394  if (m_verbose) {
395  B2DEBUG(55, "m_houghPlane, 2");
396  printArray3D(*m_phoughPlane, m_fullbins);
397  }
398  getCM();
399 }
400 
401 
402 std::vector<std::vector<unsigned short>> NDFinder::getHitsVsClusters(std::vector<SimpleCluster>& clusters)
403 {
404  unsigned short nClusters = clusters.size();
405  unsigned short default_value = 0;
406  std::vector<unsigned short> hitElem(m_nHits, default_value);
407  std::vector<std::vector<unsigned short>> hitsVsClusters(nClusters, hitElem);
408 
409  for (unsigned long iclus = 0; iclus < clusters.size(); iclus++) {
410  SimpleCluster cli = clusters[iclus];
411  vector<cell_index> entries = cli.getEntries();
412  cell_index maxid = getMax(entries);
413  for (unsigned long ihit = 0; ihit < m_hitIds.size(); ihit++) {
414  ushort contrib = hitContrib(maxid, ihit);
415  hitsVsClusters[iclus][ihit] = contrib;
416  }
417  }
418  return hitsVsClusters;
419 }
420 
421 vector<SimpleCluster>
422 NDFinder::relateHitsToClusters(std::vector<std::vector<unsigned short>>& hitsVsClusters, std::vector<SimpleCluster>& clusters)
423 {
424  vector<SimpleCluster> useclusters;
425  vector<unsigned long> restHits;
426  if (hitsVsClusters.size() > 0) {
427  for (unsigned long ihit = 0; ihit < m_hitIds.size(); ihit++) {
428  int cid = hitToCluster(hitsVsClusters, ihit);
429  if (cid != -1) { // add hit to cluster
430  clusters[cid].add_hit(ihit, hitsVsClusters[cid][ihit], m_hitOrients[ihit]);
431  } else {
432  restHits.push_back(ihit);
433  }
434  }
435 
436  // select clusters, try to reassign rest hits
437  vector<bool> processed(hitsVsClusters.size(), false);
438  for (long unsigned int iround = 0; iround <= (m_params.minhits - 2); iround++) {
439  if (iround > 0) {
440  vector<unsigned long> stillRest;
441  for (unsigned long irhit = 0; irhit < restHits.size(); irhit++) {
442  unsigned long ihit = restHits[irhit];
443  int cid = hitToCluster(hitsVsClusters, ihit);
444  if (cid != -1) { // add hit to cluster
445  clusters[cid].add_hit(ihit, hitsVsClusters[cid][ihit], m_hitOrients[ihit]);
446  } else {
447  stillRest.push_back(ihit);
448  }
449  }
450  restHits = stillRest;
451  }
452  stringstream msg;
453  for (long unsigned int iclus = 0; iclus < hitsVsClusters.size(); iclus++) {
454  if (processed[iclus])
455  continue;
456  SimpleCluster& clu = clusters[iclus];
457  //if (clu.get_hits().size() >= m_params.minhits) {
458  if (clu.get_hits().size() >= m_params.minhits && clu.get_naxial() >= m_params.minhits_axial) {
459  useclusters.push_back(clusters[iclus]);
460  processed[iclus] = true;
461  msg << ", add iclu=" << iclus << "(nHits=" << clu.get_hits().size() << ",round=" << iround << ")";
462  } else if (clu.get_hits().size() < ((m_params.minhits - 2) + iround)) { //remove small clusters
463  msg << ", skip iclu=" << iclus << "(nHits=" << clu.get_hits().size() << ",round=" << iround << ")";
464  for (unsigned long int ihics = 0; ihics < hitsVsClusters[0].size(); ihics++) {
465  hitsVsClusters[iclus][ihics] = 0;
466  }
467  processed[iclus] = true;
468  for (unsigned short ihirs : clu.get_hits()) {
469  restHits.push_back(ihirs);
470  }
471  }
472  }
473  if (m_verbose) {
475  B2DEBUG(55, "iround=" << iround <<
476  ", restHits: " << restHits.size() << ", useclusters: " <<
477  useclusters.size() << ", msg=" << msg.str());
478  }
479  }
480  }
481  return useclusters;
482 }
483 
484 
486 {
487  c3array& houghPlane = *m_phoughPlane;
488  m_clusterer.setNewPlane(houghPlane);
489  std::vector<SimpleCluster> allClusters = m_clusterer.dbscan();
490  std::vector<SimpleCluster> clusters;
491  for (SimpleCluster& clu : allClusters) {
492  if (clu.getEntries().size() > m_params.mincells) {
493  clusters.push_back(clu);
494  }
495  }
496 
497  std::vector<std::vector<unsigned short>> hitsVsClusters = getHitsVsClusters(clusters);
498  vector<SimpleCluster> useclusters = relateHitsToClusters(hitsVsClusters, clusters);
499 
500  for (unsigned long iclus = 0; iclus < useclusters.size(); iclus++) {
501  SimpleCluster cli = useclusters[iclus];
502  vector<cell_index> entries = cli.getEntries();
503  cell_index maxid = getMax(entries);
504  ushort maxval = houghPlane[maxid[0]][maxid[1]][maxid[2]];
505  float cutoff = m_params.thresh * maxval;
506  vector<cellweight> highWeight = getHighWeight(entries, cutoff);
507  vector<vector<long int>> flatHW;
508  for (auto cx : highWeight) {
509  flatHW.push_back({cx.index[0], cx.index[1], cx.index[2], (unsigned short int)cx.weight});
510  }
511  vector<double> result = getWeightedMean(highWeight);
512  vector<double> estimate = getBinToVal(result);
513  vector<double> transformed = transform(estimate);
515  m_NDFinderTracks.push_back(NDFinderTrack(transformed, cli));
516  }
517 }
518 
519 float NDFinder::transformVar(float estVal, int idx)
520 {
521  if (idx == 0) { //omega
522  if (estVal == 0.) {
523  return estVal;
524  } else {
525  return - 1 / cdc_track_radius(1. / estVal);
526  }
527  } else if (idx == 1) { // phi
528  float phiMod = estVal;
529  if (estVal > 180) {
530  phiMod -= 360.;
531  }
532  return phiMod * TMath::DegToRad();
533  } else { // theta
534  float thetRad = estVal * TMath::DegToRad();
535  return cos(thetRad) / sin(thetRad);
536  }
537 }
538 std::vector<double> NDFinder::transform(std::vector<double> estimate)
539 {
540  std::vector<double> transformed;
541  for (int idx = 0; idx < 3; idx++) {
542  transformed.push_back(transformVar(estimate[idx], idx));
543  }
544  return transformed;
545 }
546 
547 
548 std::vector<double> NDFinder::getBinToVal(std::vector<double> thisAv)
549 {
550  vector<double> estimate;
551  for (ushort idim = 0; idim < 3; idim++) {
552  double trafd = m_acceptRanges[idim][0] + (thisAv[idim] + 0.5) * m_slotsizes[idim];
553  estimate.push_back(trafd);
554  }
555  return estimate;
556 }
557 
558 std::vector<double> NDFinder::getWeightedMean(std::vector<cellweight> highWeight)
559 {
560  double axomega = 0.;
561  double axphi = 0.;
562  double axtheta = 0.;
563  long weightSum = 0;
564 
565  for (cellweight& elem : highWeight) {
566  axomega += elem.index[0] * elem.weight;
567  axphi += elem.index[1] * elem.weight;
568  axtheta += elem.index[2] * elem.weight;
569  weightSum += elem.weight;
570  }
571  axomega /= weightSum;
572  axphi /= weightSum;
573  axtheta /= weightSum;
574  vector<double> result = {axomega, axphi, axtheta};
575  return result;
576 }
577 
578 int NDFinder::hitToCluster(std::vector<std::vector<unsigned short>>& hitsVsClusters, unsigned short ihit)
579 {
580  int cur_clus = -1;
581  std::vector<std::vector<ushort>> candWeights;
582 
583  for (unsigned long iclus = 0; iclus < hitsVsClusters.size(); iclus++) {
584  std::vector<ushort> cwEntry = {(ushort) iclus, hitsVsClusters[iclus][ihit]};
585  candWeights.push_back(cwEntry);
586  }
587  if (candWeights.size() == 1) { // single cluster case
588  if (candWeights[0][1] > 0) { // hit contributes non-zero weight
589  cur_clus = 0;
590  }
591  } else if (candWeights.size() > 1) {
592  struct mySortClass {
593  bool operator()(vector<ushort> i, vector<ushort> j) {return (i[1] > j[1]);}
594  } mySortObject;
595  sort(candWeights.begin(), candWeights.end(), mySortObject);
596 
597  if (candWeights[0][1] > 0) { // the hit contributes non-zero weight
598  float crit = fabs((float)candWeights[0][1] - (float)candWeights[1][1]) / ((float)candWeights[0][1]);
599  if (crit > m_params.minassign) {
600  cur_clus = candWeights[0][0];
601  }
602  }
603  }
604  return cur_clus;
605 }
606 
607 
608 cell_index NDFinder::getMax(const std::vector<cell_index>& entries)
609 {
610  ushort cur_weight = 0;
611  cell_index cur_max_index = {0, 0, 0};
612  const c3array& houghPlane = *m_phoughPlane;
613 
614  for (const cell_index& entry : entries) {
615  if (houghPlane[entry[0]][entry[1]][entry[2]] > cur_weight) {
616  cur_weight = houghPlane[entry[0]][entry[1]][entry[2]];
617  cur_max_index = entry;
618  }
619  }
620  return cur_max_index;
621 }
622 
623 
624 vector<cellweight> NDFinder::getHighWeight(std::vector<cell_index> entries, float cutoff)
625 {
626  vector<cellweight> cellsAndWeight;
627  const c3array& houghPlane = *m_phoughPlane;
628  for (cell_index& entry : entries) {
629  ushort cellWeight = houghPlane[entry[0]][entry[1]][entry[2]];
630  if (cellWeight > cutoff) {
631  cellweight cur_elem;
632  cur_elem.index = entry;
633  cur_elem.weight = cellWeight;
634  cellsAndWeight.push_back(cur_elem);
635  }
636  }
637  return cellsAndWeight;
638 }
639 
640 
641 ushort NDFinder::hitContrib(cell_index peak, ushort ihit)
642 {
643  ushort contrib = 0;
644  ushort iHoughPhi = peak[1];
645  short Dstart = m_vecDstart[ihit];
646  ushort iphi = iHoughPhi - Dstart;
647  if (Dstart > iHoughPhi && Dstart > 300) {
648  iphi = m_nPhiFull - Dstart + iHoughPhi;
649  }
650  ushort iomega = peak[0];
651  ushort itheta = peak[2];
652  ushort orient = m_hitOrients[ihit];
653 
654  const c2array& arrayHitMod = *m_parrayHitMod;
655  c2elem hitr = arrayHitMod[m_hitIds[ihit]][1];
656  unsigned short prio = m_prioPos[ ihit ];
657  if (Dstart > m_nPhiFull) {
658  B2ERROR("phi overflow: iHoughPhi = " << iHoughPhi << ", Dstart = " << Dstart << ", iphi=" << iphi);
659  }
660  if (iphi < m_nPhiUse) { // hit covers current phi area, get contribution
661  if (orient == 1) { // axial
662  contrib = (*m_parrayAxial)[hitr][prio][iomega][iphi][itheta];
663  } else { // stereo
664  contrib = (*m_parrayStereo)[hitr][prio][iomega][iphi][itheta];
665  }
666  }
667  return contrib;
668 }
669 
670 
671 void
673 {
674  clusterer_params cpa = m_clusterer.getParams();
675  B2DEBUG(55, "clusterer_params minweight=" << cpa.minweight << ", minpts=" << cpa.minpts << ", diagonal=" << cpa.diagonal);
676  B2DEBUG(55, "ndFinderParams minhits=" << m_params.minhits <<
677  ", minhits_axial=" << m_params.minhits_axial << ", thresh=" << m_params.thresh << ", minassign=" <<
678  m_params.minassign <<
679  ", mincells=" << m_params.mincells);
680 }
Clustering module.
Definition: Clusterizend.h:114
Store track parameters of found tracks.
Definition: NDFinder.h:28
std::vector< cellweight > getHighWeight(std::vector< cell_index > entries, float cutoff)
Candidate cells as seed for the clustering.
Definition: NDFinder.cc:624
std::vector< SimpleCluster > relateHitsToClusters(std::vector< std::vector< unsigned short >> &hitsVsClusters, std::vector< SimpleCluster > &clusters)
Use confusion matrix to related hits to clusters.
Definition: NDFinder.cc:422
void initBins()
Initialize the binnings and reserve the arrays.
Definition: NDFinder.cc:72
void initHitModAxSt(c2array &hitMod)
Initialize hit modulo mappings.
Definition: NDFinder.cc:192
void getCM()
NDFinder internal functions for track finding.
Definition: NDFinder.cc:485
void squeezeAll(ndbinning writebins, c5array &writeArray, c5array &readArray, int outparcels, int inparcels)
Loop over all hits and theta bins and squeeze all 2D (omega,phi) planes.
Definition: NDFinder.cc:246
void loadArray(const std::string &filename, ndbinning bins, c5array &hitsToTracks)
Load an NDFinder array of hit representations in track phase space.
Definition: NDFinder.cc:159
void restoreZeros(ndbinning zerobins, ndbinning compbins, c5array &expArray, const c5array &compArray)
Restore non-zero suppressed hit curves.
Definition: NDFinder.cc:258
ushort hitContrib(cell_index peak, ushort ihit)
Determine weight contribution of a single hit to a single cell.
Definition: NDFinder.cc:641
void squeezeOne(c5array &writeArray, c5array &readArray, int outparcels, int inparcels, c5index ihit, c5index iprio, c5index itheta, c5elem nomega)
Squeeze phi-axis in a 2D (omega,phi) plane.
Definition: NDFinder.cc:223
std::vector< double > getBinToVal(std::vector< double >)
Scale the weighted center to track parameter values.
Definition: NDFinder.cc:548
void addC3Comp(ushort hitr, ushort prio, const c5array &hitsToTracks, short Dstart, ndbinning bins)
In place array addition to houghmap Comp: A = A + B.
Definition: NDFinder.cc:362
void printParams()
Debug: print configured parameters.
Definition: NDFinder.cc:672
void addLookup(unsigned short ihit)
Add a single axial or stereo hit to the houghmap.
Definition: NDFinder.cc:331
int hitToCluster(std::vector< std::vector< unsigned short >> &hitsVsClusters, unsigned short hit)
Determine the best cluster for a single hit, given hitsVsClusters confusion matrix.
Definition: NDFinder.cc:578
void findTracks()
main function for track finding
Definition: NDFinder.cc:388
std::vector< std::vector< unsigned short > > getHitsVsClusters(std::vector< SimpleCluster > &clusters)
Create hits to clusters confusion matrix.
Definition: NDFinder.cc:402
std::vector< double > getWeightedMean(std::vector< cellweight >)
Calculate the weighted center of a cluster.
Definition: NDFinder.cc:558
float transformVar(float estVal, int idx)
Calculate physical units.
Definition: NDFinder.cc:519
void printArray3D(c3array &hitsToTracks, ndbinning bins, ushort, ushort, ushort, ushort)
Debug Tool: Print part of the houghmap.
Definition: NDFinder.cc:283
cell_index getMax(const std::vector< cell_index > &)
Peak cell in cluster.
Definition: NDFinder.cc:608
void init(int minweight, int minpts, bool diagonal, int minhits, int minhits_axial, double thresh, double minassign, int mincells, bool verbose, std::string &axialFile, std::string &stereoFile)
initialization
Definition: NDFinder.cc:25
Type for found clusters.
Definition: Clusterizend.h:34
unsigned long get_naxial()
Get number related axial hits.
Definition: Clusterizend.h:79
std::vector< unsigned short > get_hits()
Get ids of related hits (indices of the TS StoreArray)
Definition: Clusterizend.h:89
std::vector< cell_index > getEntries()
Get member cells in the cluster.
Definition: Clusterizend.h:62
unsigned short c2elem
TS-Id to 1/32 phi-sector mapping is stored in a 2D array.
Definition: NDFinderDefs.h:34
unsigned short c5elem
Store hit patterns in a 5D array (hitid, prio, omega, phi, theta)
Definition: NDFinderDefs.h:23
Abstract base class for different kinds of events.
unsigned short minweight
minimum weight for a cluster cell
Definition: Clusterizend.h:22
bool diagonal
Consider diagonal adjacent cells as neighbors.
Definition: Clusterizend.h:26
unsigned short minpts
minimum number of neighbours for a cluster core cell
Definition: Clusterizend.h:24
Default binning in a (7/32) phi-sector.
Definition: NDFinderDefs.h:41