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"
26 bool diagonal,
int minhits,
int minhits_axial,
29 int mincells,
bool verbose,
30 string& axialFile,
string& stereoFile)
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;
46 B2DEBUG(55,
"initialized binnings");
48 initHitModAxSt(*m_parrayHitMod);
49 B2DEBUG(55,
"initialized HitMod, a map of tsid to (orient, relid, letter).");
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 ");
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)");
69 m_clusterer.setPlaneShape(m_planeShape);
95 m_nPhiOne = m_nPhiFull / m_params.phigeo;
99 m_nPhiUse = m_params.parcels * m_nPhiOne;
101 m_nPhiExp = m_params.parcelsExp * m_nPhiOne;
103 m_axbins.hitid = m_nAx;
104 m_stbins.hitid = m_nSt;
105 m_axbins.phi = m_nPhiUse;
106 m_stbins.phi = m_nPhiUse;
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;
115 m_expstbins.phi = m_nPhiExp;
117 m_fullbins.hitid = m_nTS;
118 m_fullbins.phi = m_nPhiFull;
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 }};
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);
138 m_nWires = {160, 160, 192, 224, 256, 288, 320, 352, 384};
140 m_planeShape = {m_nOmega, m_nPhiFull, m_nTheta};
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;
147 float ssPhi = (phiRange[1] - phiRange[0]) / m_nPhiOne;
148 float ssTheta = (thetaRange[1] - thetaRange[0]) / m_nTheta;
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);
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);
167 if (arrayFileGZ.is_open()) {
168 while (arrayStream >> uline) {
169 flatArray.push_back(uline);
173 B2ERROR(
"could not open array file: " << filename);
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];
194 unsigned short phigeo = m_params.phigeo;
195 std::vector<int> modSLs;
196 std::vector<int> toslid;
197 std::vector<int> sloffsets = {0};
198 std::vector<int> modoffsetsAxSt = {0, 0};
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);
204 sloffsets.push_back(sloffsets[isl] + m_nWires[isl]);
205 int last = modoffsetsAxSt[isl];
206 modoffsetsAxSt.push_back(last + m_nWires[isl] / phigeo);
208 for (
int its = 0; its < m_nTS; its++) {
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;
223 void NDFinder::squeezeOne(c5array& writeArray, c5array& readArray,
int outparcels,
int inparcels, c5index ihit, c5index iprio,
224 c5index itheta,
c5elem nomega)
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];
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];
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];
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);
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) {
271 for (c5index jtheta = 1; jtheta < zerobins.theta; jtheta++) {
272 expArray[ihit][iprio][iomega][phiCur][jtheta] = compArray[ihit][iprio][iomega][iphi + 2][itheta];
284 ushort minweightShow = 1)
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];
296 if (phiSlice == 0 and not started) {
298 }
else if (phiSlice != 0 and not started) {
300 }
else if (phiSlice != 0 and started) {
302 }
else if (phiSlice == 0 and started) {
308 B2DEBUG(55,
"printArray3D, phistart = " << phistart <<
", phiend = " << phiend);
309 auto d3shp = hitsToTracks.shape();
310 B2DEBUG(55,
"printArray shape: " << d3shp[0] <<
", " << d3shp[1] <<
", " << d3shp[2]);
312 if (phiend == phistart) {
316 for (c3index iomega = 0; iomega < bins.omega; iomega += om_inc) {
317 for (c3index iphi = phistart; iphi < phiend; iphi += phi_inc) {
319 ushort valRed = (ushort)((hitsToTracks[iomega][iphi][itheta]) / divide);
320 if (valRed <= minweightShow)
333 c2array& arrayHitMod = *m_parrayHitMod;
335 c2elem orient = arrayHitMod[m_hitIds[ihit]][0];
336 c2elem hitr = arrayHitMod[m_hitIds[ihit]][1];
337 c2elem letter = arrayHitMod[m_hitIds[ihit]][2];
339 unsigned short prio = m_prioPos[ ihit ];
340 short letterShort = (short) letter;
344 short DstartShort = (letterShort - 3) % m_params.phigeo * m_nPhiOne;
345 if (DstartShort < 0) {DstartShort = m_nPhiFull + DstartShort;}
348 short DstartComp = (letterShort - 5) % m_params.phigeo * m_nPhiOne;
349 if (DstartComp < 0) {DstartComp = m_nPhiFull + DstartComp;}
351 m_vecDstart.push_back(DstartShort);
352 m_hitOrients.push_back(orient);
355 addC3Comp(hitr, prio, *m_pcompAxial, DstartComp, m_compaxbins);
357 addC3Comp(hitr, prio, *m_pcompStereo, DstartComp, m_compstbins);
363 ushort prio,
const c5array& hitsToTracks,
367 c3array& houghPlane = *m_phoughPlane;
368 for (ushort itheta = 0; itheta < m_nTheta; itheta++) {
369 if (bins.theta > 1) {
372 for (ushort iomega = 0; iomega < m_nOmega; iomega++) {
373 ushort startfield = 0;
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];
391 for (
unsigned short ihit = 0; ihit < m_nHits; ihit++) {
395 B2DEBUG(55,
"m_houghPlane, 2");
396 printArray3D(*m_phoughPlane, m_fullbins);
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);
409 for (
unsigned long iclus = 0; iclus < clusters.size(); 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;
418 return hitsVsClusters;
421 vector<SimpleCluster>
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);
430 clusters[cid].add_hit(ihit, hitsVsClusters[cid][ihit], m_hitOrients[ihit]);
432 restHits.push_back(ihit);
437 vector<bool> processed(hitsVsClusters.size(),
false);
438 for (
long unsigned int iround = 0; iround <= (m_params.minhits - 2); iround++) {
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);
445 clusters[cid].add_hit(ihit, hitsVsClusters[cid][ihit], m_hitOrients[ihit]);
447 stillRest.push_back(ihit);
450 restHits = stillRest;
453 for (
long unsigned int iclus = 0; iclus < hitsVsClusters.size(); iclus++) {
454 if (processed[iclus])
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)) {
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;
467 processed[iclus] =
true;
468 for (
unsigned short ihirs : clu.
get_hits()) {
469 restHits.push_back(ihirs);
475 B2DEBUG(55,
"iround=" << iround <<
476 ", restHits: " << restHits.size() <<
", useclusters: " <<
477 useclusters.size() <<
", msg=" << msg.str());
487 c3array& houghPlane = *m_phoughPlane;
488 m_clusterer.setNewPlane(houghPlane);
489 std::vector<SimpleCluster> allClusters = m_clusterer.dbscan();
490 std::vector<SimpleCluster> clusters;
492 if (clu.getEntries().size() > m_params.mincells) {
493 clusters.push_back(clu);
497 std::vector<std::vector<unsigned short>> hitsVsClusters = getHitsVsClusters(clusters);
498 vector<SimpleCluster> useclusters = relateHitsToClusters(hitsVsClusters, clusters);
500 for (
unsigned long iclus = 0; iclus < useclusters.size(); 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});
511 vector<double> result = getWeightedMean(highWeight);
512 vector<double> estimate = getBinToVal(result);
513 vector<double> transformed = transform(estimate);
525 return - 1 / cdc_track_radius(1. / estVal);
527 }
else if (idx == 1) {
528 float phiMod = estVal;
532 return phiMod * TMath::DegToRad();
534 float thetRad = estVal * TMath::DegToRad();
535 return cos(thetRad) / sin(thetRad);
538 std::vector<double> NDFinder::transform(std::vector<double> estimate)
540 std::vector<double> transformed;
541 for (
int idx = 0; idx < 3; idx++) {
542 transformed.push_back(transformVar(estimate[idx], idx));
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);
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;
571 axomega /= weightSum;
573 axtheta /= weightSum;
574 vector<double> result = {axomega, axphi, axtheta};
581 std::vector<std::vector<ushort>> candWeights;
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);
587 if (candWeights.size() == 1) {
588 if (candWeights[0][1] > 0) {
591 }
else if (candWeights.size() > 1) {
593 bool operator()(vector<ushort> i, vector<ushort> j) {
return (i[1] > j[1]);}
595 sort(candWeights.begin(), candWeights.end(), mySortObject);
597 if (candWeights[0][1] > 0) {
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];
610 ushort cur_weight = 0;
611 cell_index cur_max_index = {0, 0, 0};
612 const c3array& houghPlane = *m_phoughPlane;
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;
620 return cur_max_index;
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) {
632 cur_elem.index = entry;
633 cur_elem.weight = cellWeight;
634 cellsAndWeight.push_back(cur_elem);
637 return cellsAndWeight;
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;
650 ushort iomega = peak[0];
651 ushort itheta = peak[2];
652 ushort orient = m_hitOrients[ihit];
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);
660 if (iphi < m_nPhiUse) {
662 contrib = (*m_parrayAxial)[hitr][prio][iomega][iphi][itheta];
664 contrib = (*m_parrayStereo)[hitr][prio][iomega][iphi][itheta];
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);
Store track parameters of found tracks.
std::vector< cellweight > getHighWeight(std::vector< cell_index > entries, float cutoff)
Candidate cells as seed for the clustering.
std::vector< SimpleCluster > relateHitsToClusters(std::vector< std::vector< unsigned short >> &hitsVsClusters, std::vector< SimpleCluster > &clusters)
Use confusion matrix to related hits to clusters.
void initBins()
Initialize the binnings and reserve the arrays.
void initHitModAxSt(c2array &hitMod)
Initialize hit modulo mappings.
void getCM()
NDFinder internal functions for track finding.
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.
void loadArray(const std::string &filename, ndbinning bins, c5array &hitsToTracks)
Load an NDFinder array of hit representations in track phase space.
void restoreZeros(ndbinning zerobins, ndbinning compbins, c5array &expArray, const c5array &compArray)
Restore non-zero suppressed hit curves.
ushort hitContrib(cell_index peak, ushort ihit)
Determine weight contribution of a single hit to a single cell.
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.
std::vector< double > getBinToVal(std::vector< double >)
Scale the weighted center to track parameter values.
void addC3Comp(ushort hitr, ushort prio, const c5array &hitsToTracks, short Dstart, ndbinning bins)
In place array addition to houghmap Comp: A = A + B.
void printParams()
Debug: print configured parameters.
void addLookup(unsigned short ihit)
Add a single axial or stereo hit to the houghmap.
int hitToCluster(std::vector< std::vector< unsigned short >> &hitsVsClusters, unsigned short hit)
Determine the best cluster for a single hit, given hitsVsClusters confusion matrix.
void findTracks()
main function for track finding
std::vector< std::vector< unsigned short > > getHitsVsClusters(std::vector< SimpleCluster > &clusters)
Create hits to clusters confusion matrix.
std::vector< double > getWeightedMean(std::vector< cellweight >)
Calculate the weighted center of a cluster.
float transformVar(float estVal, int idx)
Calculate physical units.
void printArray3D(c3array &hitsToTracks, ndbinning bins, ushort, ushort, ushort, ushort)
Debug Tool: Print part of the houghmap.
cell_index getMax(const std::vector< cell_index > &)
Peak cell in cluster.
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
unsigned long get_naxial()
Get number related axial hits.
std::vector< unsigned short > get_hits()
Get ids of related hits (indices of the TS StoreArray)
std::vector< cell_index > getEntries()
Get member cells in the cluster.
unsigned short c2elem
TS-Id to 1/32 phi-sector mapping is stored in a 2D array.
unsigned short c5elem
Store hit patterns in a 5D array (hitid, prio, omega, phi, theta)
Abstract base class for different kinds of events.
unsigned short minweight
minimum weight for a cluster cell
bool diagonal
Consider diagonal adjacent cells as neighbors.
unsigned short minpts
minimum number of neighbours for a cluster core cell
Default binning in a (7/32) phi-sector.