1 #include <framework/logging/Logger.h>
2 #include <trg/cdc/NeuroTrigger.h>
4 #include <cdc/geometry/CDCGeometryPar.h>
5 #include <framework/gearbox/Const.h>
6 #include <framework/gearbox/Unit.h>
7 #include <framework/datastore/StoreObjPtr.h>
8 #include <framework/datastore/StoreArray.h>
9 #include <trg/cdc/dbobjects/CDCTriggerNeuroConfig.h>
10 #include <trg/cdc/dataobjects/CDCTriggerTrack.h>
25 if (p.nHidden.size() != 1 && p.nHidden.size() != p.nMLP) {
26 B2ERROR(
"Number of nHidden lists should be 1 or " << p.nMLP);
29 if (p.outputScale.size() != 1 && p.outputScale.size() != p.nMLP) {
30 B2ERROR(
"Number of outputScale lists should be 1 or " << p.nMLP);
33 bool rangeProduct = (p.phiRange.size() * p.invptRange.size() * p.thetaRange.size() * p.SLpattern.size() == p.nMLP);
35 if (p.phiRange.size() != 1 && p.phiRange.size() != p.nMLP) {
36 B2ERROR(
"Number of phiRange lists should be 1 or " << p.nMLP);
39 if (p.invptRange.size() != 1 && p.invptRange.size() != p.nMLP) {
40 B2ERROR(
"Number of invptRange lists should be 1 or " << p.nMLP);
43 if (p.thetaRange.size() != 1 && p.thetaRange.size() != p.nMLP) {
44 B2ERROR(
"Number of thetaRange lists should be 1 or " << p.nMLP);
47 if (p.SLpattern.size() != 1 && p.SLpattern.size() != p.nMLP) {
48 B2ERROR(
"Number of SLpattern lists should be 1 or " << p.nMLP);
53 if (p.maxHitsPerSL.size() != 1 && p.maxHitsPerSL.size() != p.SLpattern.size()) {
54 B2ERROR(
"Number of maxHitsPerSL lists should be 1 or " << p.SLpattern.size());
57 if (p.SLpatternMask.size() != 1 && p.SLpatternMask.size() != p.SLpattern.size()) {
58 B2ERROR(
"Number of SLpatternMask lists should be 1 or " << p.SLpattern.size());
62 unsigned short nTarget = int(p.targetZ) + int(p.targetTheta);
64 B2ERROR(
"No outputs! Turn on either targetZ or targetTheta.");
68 for (
unsigned iPhi = 0; iPhi < p.phiRange.size(); ++iPhi) {
69 if (p.phiRange[iPhi].size() != 2) {
70 B2ERROR(
"phiRange should be exactly 2 values");
74 if (p.phiRange[iPhi][0] >= p.phiRange[iPhi][1]) {
75 B2ERROR(
"phiRange[0] should be smaller than phiRange[1]");
78 if (p.phiRange[iPhi][0] < -360. || p.phiRange[iPhi][1] > 360. ||
79 (p.phiRange[iPhi][1] - p.phiRange[iPhi][0]) > 360.) {
80 B2ERROR(
"phiRange should be in [-360, 360], with maximal width of 360");
84 for (
unsigned iPt = 0; iPt < p.invptRange.size(); ++iPt) {
85 if (p.invptRange[iPt].size() != 2) {
86 B2ERROR(
"invptRange should be exactly 2 values");
89 if (p.invptRange[iPt][0] >= p.invptRange[iPt][1]) {
90 B2ERROR(
"invptRange[0] should be smaller than invptRange[1]");
94 for (
unsigned iTheta = 0; iTheta < p.thetaRange.size(); ++iTheta) {
95 if (p.thetaRange[iTheta].size() != 2) {
96 B2ERROR(
"thetaRange should be exactly 2 values");
100 if (p.thetaRange[iTheta][0] >= p.thetaRange[iTheta][1]) {
101 B2ERROR(
"thetaRange[0] should be smaller than thetaRange[1]");
104 if (p.thetaRange[iTheta][0] < 0. || p.thetaRange[iTheta][1] > 180.) {
105 B2ERROR(
"thetaRange should be in [0, 180]");
109 for (
unsigned iScale = 0; iScale < p.outputScale.size(); ++iScale) {
110 if (p.outputScale[iScale].size() != 2 * nTarget) {
111 B2ERROR(
"outputScale should be exactly " << 2 * nTarget <<
" values");
116 if (p.phiRange.size() != p.phiRangeTrain.size()) {
117 B2ERROR(
"Number of phiRange lists and phiRangeTrain lists should be equal.");
120 for (
unsigned iPhi = 0; iPhi < p.phiRange.size(); ++iPhi) {
121 if (p.phiRangeTrain[iPhi].size() != 2) {
122 B2ERROR(
"phiRangeTrain should be exactly 2 values.");
124 }
else if (p.phiRangeTrain[iPhi][0] > p.phiRange[iPhi][0] ||
125 p.phiRangeTrain[iPhi][1] < p.phiRange[iPhi][1]) {
126 B2ERROR(
"phiRangeTrain should be wider than phiRange or equal.");
131 if (p.invptRange.size() != p.invptRangeTrain.size()) {
132 B2ERROR(
"Number of invptRange lists and invptRangeTrain lists should be equal.");
135 for (
unsigned iPt = 0; iPt < p.invptRange.size(); ++iPt) {
136 if (p.invptRangeTrain[iPt].size() != 2) {
137 B2ERROR(
"invptRangeTrain should be exactly 2 values.");
139 }
else if (p.invptRangeTrain[iPt][0] > p.invptRange[iPt][0] ||
140 p.invptRangeTrain[iPt][1] < p.invptRange[iPt][1]) {
141 B2ERROR(
"invptRangeTrain should be wider than invptRange or equal.");
146 if (p.thetaRange.size() != p.thetaRangeTrain.size()) {
147 B2ERROR(
"Number of thetaRange lists and thetaRangeTrain lists should be equal.");
150 for (
unsigned iTheta = 0; iTheta < p.thetaRange.size(); ++iTheta) {
151 if (p.thetaRangeTrain[iTheta].size() != 2) {
152 B2ERROR(
"thetaRangeTrain should be exactly 2 values.");
154 }
else if (p.thetaRangeTrain[iTheta][0] > p.thetaRange[iTheta][0] ||
155 p.thetaRangeTrain[iTheta][1] < p.thetaRange[iTheta][1]) {
156 B2ERROR(
"thetaRangeTrain should be wider than thetaRange or equal.");
166 for (
unsigned iMLP = 0; iMLP < p.nMLP; ++iMLP) {
168 vector<unsigned> indices = getRangeIndices(p, iMLP);
170 unsigned short maxHits = (p.maxHitsPerSL.size() == 1) ? p.maxHitsPerSL[0] : p.maxHitsPerSL[indices[3]];
171 unsigned long SLpattern = p.SLpattern[indices[3]];
172 unsigned long SLpatternMask = (p.SLpatternMask.size() == 1) ? p.SLpatternMask[0] : p.SLpatternMask[indices[3]];
173 unsigned short nInput = 27 * maxHits;
174 vector<float> nHidden = (p.nHidden.size() == 1) ? p.nHidden[0] : p.nHidden[iMLP];
175 vector<unsigned short> nNodes = {nInput};
176 for (
unsigned iHid = 0; iHid < nHidden.size(); ++iHid) {
177 if (p.multiplyHidden) {
178 nNodes.push_back(nHidden[iHid] * nNodes[0]);
180 nNodes.push_back(nHidden[iHid]);
183 nNodes.push_back(nTarget);
184 unsigned short targetVars = int(p.targetZ) + (int(p.targetTheta) << 1);
186 vector<float> phiRange = p.phiRangeTrain[indices[0]];
187 vector<float> invptRange = p.invptRangeTrain[indices[1]];
188 vector<float> thetaRange = p.thetaRangeTrain[indices[2]];
189 B2DEBUG(50,
"Ranges for sector " << iMLP
190 <<
": phiRange [" << phiRange[0] <<
", " << phiRange[1]
191 <<
"], invptRange [" << invptRange[0] <<
", " << invptRange[1]
192 <<
"], thetaRange [" << thetaRange[0] <<
", " << thetaRange[1]
193 <<
"], SLpattern " << SLpattern);
195 vector<float> outputScale = (p.outputScale.size() == 1) ? p.outputScale[0] : p.outputScale[iMLP];
202 outputScale[2 * int(p.targetZ)] *=
Unit::deg;
203 outputScale[2 * int(p.targetZ) + 1] *=
Unit::deg;
206 m_MLPs.push_back(
CDCTriggerMLP(nNodes, targetVars, outputScale,
207 phiRange, invptRange, thetaRange,
208 maxHits, SLpattern, SLpatternMask, p.tMax, p.T0fromHits,
218 std::vector<unsigned> indices = {0, 0, 0, 0};
219 if (p.phiRange.size() * p.invptRange.size() * p.thetaRange.size() * p.SLpattern.size() == p.nMLP) {
220 indices[0] = isector % p.phiRange.size();
221 indices[1] = (isector / p.phiRange.size()) % p.invptRange.size();
222 indices[2] = (isector / (p.phiRange.size() * p.invptRange.size())) % p.thetaRange.size();
223 indices[3] = isector / (p.phiRange.size() * p.invptRange.size() * p.thetaRange.size());
225 indices[0] = (p.phiRange.size() == 1) ? 0 : isector;
226 indices[1] = (p.invptRange.size() == 1) ? 0 : isector;
227 indices[2] = (p.thetaRange.size() == 1) ? 0 : isector;
228 indices[3] = (p.SLpattern.size() == 1) ? 0 : isector;
239 for (
int iSL = 0; iSL < 9; ++iSL) {
240 m_TSoffset[iSL] = nTS;
241 nTS += cdc.nWiresInLayer(layerId);
242 m_TSoffset[iSL + 1] = nTS;
243 for (
int priority = 0; priority < 2; ++ priority) {
244 m_radius[iSL][priority] = cdc.senseWireR(layerId + priority);
246 layerId += (iSL > 0 ? 6 : 7);
253 m_segmentHits.isRequired(hitCollectionName);
254 if (!((et_option ==
"fastestpriority") || (et_option ==
"zero") || (et_option ==
"fastest2d"))) {
255 m_eventTime.isRequired(eventTimeName);
257 m_hitCollectionName = hitCollectionName;
263 vector<int> indices = {};
265 if (m_MLPs.size() == 0) {
266 B2WARNING(
"Trying to select MLP before initializing MLPs.");
271 for (
unsigned isector = 0; isector < m_MLPs.size(); ++isector) {
272 if (m_MLPs[isector].inPhiRange(phi0) && m_MLPs[isector].inInvptRange(invpt)
273 && m_MLPs[isector].inThetaRange(theta)) {
274 indices.push_back(isector);
278 if (indices.size() == 0) {
279 B2DEBUG(150,
"Track does not match any sector.");
280 B2DEBUG(150,
"invpt=" << invpt <<
", phi=" << phi0 * 180. / M_PI <<
", theta=" << theta * 180. / M_PI);
290 if (MLPs.size() == 0) {
295 for (
unsigned i = 0; i < MLPs.size(); ++i) {
296 int isector = MLPs[i];
297 unsigned long sectorPattern = m_MLPs[isector].getSLpattern();
298 B2DEBUG(250,
"hitPattern " << pattern <<
" sectorPattern " << sectorPattern);
300 if (sectorPattern == 0) {
301 B2DEBUG(250,
"found match for general sector");
305 if (pattern == sectorPattern) {
306 B2DEBUG(250,
"found match for hit pattern " << pattern);
313 if (neurotrackinputmode) {
314 B2DEBUG(150,
"No sector found to match pattern, using sector 0" << pattern <<
".");
317 B2DEBUG(150,
"No sector found to match pattern " << pattern <<
".");
327 double omega = track.getOmega();
328 for (
int iSL = 0; iSL < 9; ++iSL) {
329 for (
int priority = 0; priority < 2; ++priority) {
330 m_alpha[iSL][priority] = (m_radius[iSL][priority] * abs(omega) < 2.) ?
331 asin(m_radius[iSL][priority] * omega / 2.) :
333 m_idRef[iSL][priority] = remainder(((track.getPhi0() - m_alpha[iSL][priority]) *
334 (m_TSoffset[iSL + 1] - m_TSoffset[iSL]) / 2. / M_PI),
335 (m_TSoffset[iSL + 1] - m_TSoffset[iSL]));
343 unsigned precisionPhi = m_precision[0];
344 unsigned precisionAlpha = m_precision[0];
345 unsigned precisionScale = m_precision[1];
346 unsigned precisionId = m_precision[2];
348 double omega = track.getOmega();
349 long phi = round(track.getPhi0() * (1 << precisionPhi));
351 for (
int iSL = 0; iSL < 9; ++iSL) {
352 for (
int priority = 0; priority < 2; ++priority) {
354 m_alpha[iSL][priority] =
355 (m_radius[iSL][priority] * abs(omega) < 2) ?
356 round(asin(m_radius[iSL][priority] * omega / 2) * (1 << precisionAlpha)) :
357 round(M_PI_2 * (1 << precisionAlpha));
358 long dphi = (precisionAlpha >= precisionPhi) ?
359 phi - (
long(m_alpha[iSL][priority]) >> (precisionAlpha - precisionPhi)) :
360 phi - (long(m_alpha[iSL][priority]) << (precisionPhi - precisionAlpha));
361 m_idRef[iSL][priority] =
362 long(dphi * round((m_TSoffset[iSL + 1] - m_TSoffset[iSL]) / 2. / M_PI * (1 << precisionScale)) /
363 (
long(1) << (precisionPhi + precisionScale - precisionId)));
365 m_alpha[iSL][priority] /= (1 << precisionAlpha);
366 m_idRef[iSL][priority] /= (1 << precisionId);
374 int iSL = hit.getISuperLayer();
375 int priority = hit.getPriorityPosition();
377 double relId = hit.getSegmentID() + 0.5 * (((priority >> 1) & 1) - (priority & 1))
378 - m_TSoffset[iSL] - m_idRef[iSL][int(priority < 3)];
379 relId = remainder(relId, (m_TSoffset[iSL + 1] - m_TSoffset[iSL]));
386 if (et_option != m_MLPs[isector].get_et_option()) {
387 B2WARNING(
"Used event time option is different to the one set in the MLP"
388 <<
LogVar(
"et_option", et_option) <<
LogVar(
"isector", isector)
389 <<
LogVar(
"et_option_mlp", m_MLPs[isector].get_et_option()));
391 if (et_option ==
"etf_or_fastestpriority") {
392 bool hasT0 = m_eventTime->hasBinnedEventT0(Const::CDC);
394 m_T0 = m_eventTime->getBinnedEventT0(Const::CDC);
401 for (
unsigned ihit = 0; ihit < axialHits.
size(); ++ihit) {
403 if (axialHits.
weight(ihit) < 0)
continue;
404 unsigned short iSL = axialHits[ihit]->getISuperLayer();
406 if ((iSL % 2 == 1) && !neuroinputmode)
continue;
408 double relId = getRelId(*axialHits[ihit]);
409 if (m_MLPs[isector].isRelevant(relId, iSL) &&
410 axialHits[ihit]->priorityTime() < m_T0) {
411 m_T0 = axialHits[ihit]->priorityTime();
414 if (!neuroinputmode) {
417 for (
int ihit = 0; ihit < hits.getEntries(); ++ihit) {
418 unsigned short iSL = hits[ihit]->getISuperLayer();
420 if (iSL % 2 == 0)
continue;
422 double relId = getRelId(*hits[ihit]);
423 if (m_MLPs[isector].isRelevant(relId, iSL) && hits[ihit]->priorityTime() < m_T0) {
424 m_T0 = hits[ihit]->priorityTime();
435 }
else if (et_option ==
"etf_or_fastest2d") {
436 bool hasT0 = m_eventTime->hasBinnedEventT0(Const::CDC);
438 m_T0 = m_eventTime->getBinnedEventT0(Const::CDC);
445 for (
unsigned ihit = 0; ihit < axialHits.
size(); ++ihit) {
447 if (axialHits.
weight(ihit) < 0)
continue;
448 unsigned short iSL = axialHits[ihit]->getISuperLayer();
450 if (iSL % 2 == 1)
continue;
452 double relId = getRelId(*axialHits[ihit]);
453 if (m_MLPs[isector].isRelevant(relId, iSL) &&
454 axialHits[ihit]->priorityTime() < m_T0) {
455 m_T0 = axialHits[ihit]->priorityTime();
465 }
else if (et_option ==
"fastestpriority") {
466 B2DEBUG(200,
"et_option is 'fastestpriority'");
471 B2DEBUG(200,
"looping over axials:");
472 for (
unsigned ihit = 0; ihit < axialHits.
size(); ++ihit) {
474 if (axialHits.
weight(ihit) < 0)
continue;
475 unsigned short iSL = axialHits[ihit]->getISuperLayer();
477 if ((iSL % 2 == 1) && !neuroinputmode)
continue;
479 B2DEBUG(200,
" check drifttime: SL" + std::to_string(iSL) +
",ID = " + std::to_string(axialHits[ihit]->getSegmentID()) +
", t = " +
480 std::to_string(axialHits[ihit]->priorityTime()));
481 double relId = getRelId(*axialHits[ihit]);
482 if (m_MLPs[isector].isRelevant(relId, iSL) &&
483 axialHits[ihit]->priorityTime() < m_T0) {
484 m_T0 = axialHits[ihit]->priorityTime();
485 B2DEBUG(200,
" new t0: " << std::to_string(m_T0));
488 if (!neuroinputmode) {
491 B2DEBUG(200,
"looping over stereos:");
492 for (
int ihit = 0; ihit < hits.getEntries(); ++ihit) {
493 unsigned short iSL = hits[ihit]->getISuperLayer();
495 if (iSL % 2 == 0)
continue;
497 B2DEBUG(200,
" check drifttime: SL" + std::to_string(iSL) +
",ID = " + std::to_string(hits[ihit]->getSegmentID()) +
", t = " +
498 std::to_string(hits[ihit]->priorityTime()));
499 double relId = getRelId(*hits[ihit]);
500 if (m_MLPs[isector].isRelevant(relId, iSL) && hits[ihit]->priorityTime() < m_T0) {
501 m_T0 = hits[ihit]->priorityTime();
502 B2DEBUG(200,
" new t0: " << std::to_string(m_T0));
513 }
else if (et_option ==
"fastest2d") {
518 for (
unsigned ihit = 0; ihit < axialHits.
size(); ++ihit) {
520 if (axialHits.
weight(ihit) < 0)
continue;
521 unsigned short iSL = axialHits[ihit]->getISuperLayer();
523 if (iSL % 2 == 1)
continue;
525 double relId = getRelId(*axialHits[ihit]);
526 if (m_MLPs[isector].isRelevant(relId, iSL) &&
527 axialHits[ihit]->priorityTime() < m_T0) {
528 m_T0 = axialHits[ihit]->priorityTime();
537 }
else if (et_option ==
"zero") {
540 }
else if (et_option ==
"etf_only") {
541 bool hasT0 = m_eventTime->hasBinnedEventT0(Const::CDC);
543 m_T0 = m_eventTime->getBinnedEventT0(Const::CDC);
546 B2ERROR(
"No ETF output, but forced to use ETF!");
550 }
else if (et_option ==
"etf_or_zero") {
551 bool hasT0 = m_eventTime->hasBinnedEventT0(Const::CDC);
553 m_T0 = m_eventTime->getBinnedEventT0(Const::CDC);
560 B2ERROR(
"No valid parameter for et_option (" << et_option <<
" )!");
567 B2WARNING(
"This function is deprecated, it used the flag 'T0fromHits'. "
568 <<
"Give additionally et_option as argument to use the new function. "
569 <<
"et_option is now set for T0fromHits=True to 'etf_fastestpriority' "
570 <<
"and for T0fromHits=false to 'etf_or_zero'.");
571 if (m_MLPs[isector].getT0fromHits()) {
572 getEventTime(isector, track,
"etf_or_fastestpriority");
574 getEventTime(isector, track,
"etf_or_zero");
584 unsigned long hitPattern = 0;
585 vector<unsigned> nHits;
590 for (
unsigned ihit = 0; ihit < axialHits.
size(); ++ ihit) {
592 if (axialHits.
weight(ihit) < 0) {
595 unsigned short iSL = axialHits[ihit]->getISuperLayer();
597 if ((!neurotrackinputmode) && (iSL % 2 == 1))
continue;
599 int t = (m_hasT0) ? axialHits[ihit]->priorityTime() - m_T0 : 0;
602 }
else if (t > expert.getTMax()) {
603 t = expert.getTMax();
605 double relId = getRelId(*axialHits[ihit]);
607 if (expert.isRelevant(relId, iSL)) {
608 if (nHits[iSL] < expert.getMaxHitsPerSL()) {
609 hitPattern |= 1 << (iSL + 9 * nHits[iSL]);
612 B2DEBUG(250,
"hit in SL " << iSL);
614 B2DEBUG(250,
"hit in SL " << iSL <<
" not relevant (relId = " << relId <<
")");
617 if (!neurotrackinputmode) {
619 for (
int ihit = 0; ihit < m_segmentHits.getEntries(); ++ ihit) {
620 unsigned short iSL = m_segmentHits[ihit]->getISuperLayer();
622 if (iSL % 2 == 0)
continue;
624 int t = (m_hasT0) ? m_segmentHits[ihit]->priorityTime() - m_T0 : 0;
627 }
else if (t > expert.getTMax()) {
628 t = expert.getTMax();
630 double relId = getRelId(*m_segmentHits[ihit]);
631 if (expert.isRelevant(relId, iSL)) {
632 if (nHits[iSL] < expert.getMaxHitsPerSL()) {
633 hitPattern |= 1 << (iSL + 9 * nHits[iSL]);
636 B2DEBUG(250,
"hit in SL " << iSL);
638 B2DEBUG(250,
"hit in SL " << iSL <<
" not relevant (relId = " << relId <<
")");
642 B2DEBUG(250,
"hitPattern " << hitPattern);
643 return hitPattern & expert.getSLpatternMask();
650 vector<unsigned> selectedHitIds = {};
653 tMin.assign(expert.nNodesLayer(0), expert.getTMax());
654 vector<bool> LRknown;
655 LRknown.assign(expert.nNodesLayer(0),
false);
657 hitIds.assign(expert.nNodesLayer(0), -1);
658 vector<unsigned> nHits;
664 B2DEBUG(250,
"start hit loop over all related hits");
665 for (
unsigned ihit = 0; ihit < allHits.
size(); ++ihit) {
667 if (allHits.
weight(ihit) < 0)
continue;
668 unsigned short iSL = allHits[ihit]->getISuperLayer();
669 if (expert.getSLpatternUnmasked() != 0 &&
670 !((expert.getSLpatternUnmasked() >> iSL) & 1)) {
671 B2DEBUG(250,
"skipping hit in SL " << iSL);
675 int t = (m_hasT0) ? allHits[ihit]->priorityTime() - m_T0 : 0;
678 }
else if (t > expert.getTMax()) {
679 t = expert.getTMax();
681 double relId = getRelId(*allHits[ihit]);
682 if (expert.isRelevant(relId, iSL)) {
684 unsigned short iRef = iSL;
685 if (expert.getMaxHitsPerSL() > 1) {
686 if (nHits[iSL] < expert.getMaxHitsPerSL() &&
687 (expert.getSLpatternUnmasked() >> (iSL + 9 * nHits[iSL])) & 1) {
688 iRef += 9 * nHits[iSL];
691 for (
unsigned compare = iSL; compare < iSL + 9 * nHits[iSL]; compare += 9) {
692 if ((LRknown[iRef] && !LRknown[compare]) ||
693 (LRknown[iRef] == LRknown[compare] && tMin[iRef] < tMin[compare]))
701 useHit = (allHits[ihit]->LRknown() && t <= tMin[iRef]);
703 useHit = (allHits[ihit]->LRknown() || t <= tMin[iRef]);
705 B2DEBUG(250,
"relevant wire SL " << iSL <<
" LR " << allHits[ihit]->getLeftRight()
706 <<
" t " << t <<
" iRef " << iRef <<
" useHit " << useHit);
709 LRknown[iRef] = allHits[ihit]->LRknown();
711 hitIds[iRef] = allHits[ihit]->getArrayIndex();
716 for (
unsigned iHit = 0; iHit < hitIds.size(); ++iHit) {
717 if (hitIds[iHit] >= 0) selectedHitIds.push_back(hitIds[iHit]);
719 return selectedHitIds;
723 bool returnAllRelevant)
726 vector<unsigned> selectedHitIds = {};
729 tMin.assign(expert.nNodesLayer(0), expert.getTMax());
730 vector<bool> LRknown;
731 LRknown.assign(expert.nNodesLayer(0),
false);
733 hitIds.assign(expert.nNodesLayer(0), -1);
734 vector<unsigned> nHits;
740 B2DEBUG(250,
"start axial hit loop");
741 for (
unsigned ihit = 0; ihit < axialHits.
size(); ++ihit) {
743 if (axialHits.
weight(ihit) < 0)
continue;
744 unsigned short iSL = axialHits[ihit]->getISuperLayer();
746 if (iSL % 2 == 1)
continue;
747 if (expert.getSLpatternUnmasked() != 0 &&
748 !((expert.getSLpatternUnmasked() >> iSL) & 1)) {
749 B2DEBUG(250,
"skipping hit in SL " << iSL);
754 int t = (m_hasT0) ? axialHits[ihit]->priorityTime() - m_T0 : 0;
757 }
else if (t > expert.getTMax()) {
758 t = expert.getTMax();
760 double relId = getRelId(*axialHits[ihit]);
761 if (expert.isRelevant(relId, iSL)) {
763 unsigned short iRef = iSL;
764 if (expert.getMaxHitsPerSL() > 1) {
765 if (nHits[iSL] < expert.getMaxHitsPerSL() &&
766 (expert.getSLpatternUnmasked() >> (iSL + 9 * nHits[iSL])) & 1) {
767 iRef += 9 * nHits[iSL];
770 for (
unsigned compare = iSL; compare < iSL + 9 * nHits[iSL]; compare += 9) {
771 if ((LRknown[iRef] && !LRknown[compare]) ||
772 (LRknown[iRef] == LRknown[compare] && tMin[iRef] < tMin[compare]))
780 useHit = (axialHits[ihit]->LRknown() && t <= tMin[iRef]);
782 useHit = (axialHits[ihit]->LRknown() || t <= tMin[iRef]);
784 B2DEBUG(250,
"relevant wire SL " << iSL <<
" LR " << axialHits[ihit]->getLeftRight()
785 <<
" t " << t <<
" iRef " << iRef <<
" useHit " << useHit);
788 LRknown[iRef] = axialHits[ihit]->LRknown();
790 hitIds[iRef] = axialHits[ihit]->getArrayIndex();
792 if (returnAllRelevant) selectedHitIds.push_back(axialHits[ihit]->getArrayIndex());
797 B2DEBUG(250,
"start stereo hit loop");
798 for (
int ihit = 0; ihit < m_segmentHits.getEntries(); ++ ihit) {
799 unsigned short iSL = m_segmentHits[ihit]->getISuperLayer();
801 if (iSL % 2 == 0)
continue;
802 if (expert.getSLpatternUnmasked() != 0 &&
803 !((expert.getSLpatternUnmasked() >> iSL) & 1)) {
804 B2DEBUG(250,
"skipping hit in SL " << iSL);
808 int t = (m_hasT0) ? m_segmentHits[ihit]->priorityTime() - m_T0 : 0;
811 }
else if (t > expert.getTMax()) {
812 t = expert.getTMax();
814 double relId = getRelId(*m_segmentHits[ihit]);
815 if (expert.isRelevant(relId, iSL)) {
817 unsigned short iRef = iSL;
818 if (expert.getMaxHitsPerSL() > 1) {
819 if (nHits[iSL] < expert.getMaxHitsPerSL() &&
820 (expert.getSLpatternUnmasked() >> (iSL + 9 * nHits[iSL])) & 1) {
821 iRef += 9 * nHits[iSL];
824 for (
unsigned compare = iSL; compare < iSL + 9 * nHits[iSL]; compare += 9) {
825 if ((LRknown[iRef] && !LRknown[compare]) ||
826 (LRknown[iRef] == LRknown[compare] && tMin[iRef] < tMin[compare]))
834 useHit = (m_segmentHits[ihit]->LRknown() && t <= tMin[iRef]);
836 useHit = (m_segmentHits[ihit]->LRknown() || t <= tMin[iRef]);
838 B2DEBUG(250,
"relevant wire SL " << iSL <<
" LR " << m_segmentHits[ihit]->getLeftRight()
839 <<
" t " << t <<
" iRef " << iRef <<
" useHit " << useHit);
842 LRknown[iRef] = m_segmentHits[ihit]->LRknown();
846 if (returnAllRelevant) selectedHitIds.push_back(ihit);
851 if (!returnAllRelevant) {
852 for (
unsigned iHit = 0; iHit < hitIds.size(); ++iHit) {
853 if (hitIds[iHit] >= 0) selectedHitIds.push_back(hitIds[iHit]);
856 return selectedHitIds;
864 vector<float> inputVector;
865 inputVector.assign(expert.nNodesLayer(0), 0.);
867 vector<unsigned> nHits;
869 for (
unsigned ii = 0; ii < hitIds.size(); ++ii) {
870 int ihit = hitIds[ii];
871 unsigned short iSL = m_segmentHits[ihit]->getISuperLayer();
872 unsigned short iRef = iSL + 9 * nHits[iSL];
874 int t = (m_hasT0) ? m_segmentHits[ihit]->priorityTime() - m_T0 : 0;
877 }
else if (t > expert.getTMax()) {
878 t = expert.getTMax();
880 int LR = m_segmentHits[ihit]->getLeftRight();
881 double relId = getRelId(*m_segmentHits[ihit]);
882 int priority = m_segmentHits[ihit]->getPriorityPosition();
884 inputVector[3 * iRef] = expert.scaleId(relId, iSL);
885 float scaleT = pow(2, floor(log2(1. / expert.getTMax())));
886 inputVector[3 * iRef + 1] = (((LR >> 1) & 1) - (LR & 1)) * t * scaleT;
887 inputVector[3 * iRef + 2] = m_alpha[iSL][int(priority < 3)] * 0.5;
896 vector<float> weights = expert.getWeights();
897 vector<float> layerinput = input;
898 vector<float> layeroutput = {};
900 for (
unsigned il = 1; il < expert.nLayers(); ++il) {
902 layerinput.push_back(1.);
905 layeroutput.assign(expert.nNodesLayer(il), 0.);
907 for (
unsigned io = 0; io < layeroutput.size(); ++io) {
909 for (
unsigned ii = 0; ii < layerinput.size(); ++ii) {
910 layeroutput[io] += layerinput[ii] * weights[iw++];
913 layeroutput[io] = tanh(layeroutput[io] / 2.);
916 layerinput = layeroutput;
918 return expert.unscaleTarget(layeroutput);
924 unsigned precisionInput = m_precision[3];
925 unsigned precisionWeights = m_precision[4];
926 unsigned precisionLUT = m_precision[5];
927 unsigned precisionTanh = m_precision[3];
928 unsigned dp = precisionInput + precisionWeights - precisionLUT;
932 vector<long> inputFix(input.size(), 0);
933 for (
unsigned ii = 0; ii < input.size(); ++ii) {
934 inputFix[ii] = long(input[ii] * (1 << precisionInput));
937 vector<float> weights = expert.getWeights();
938 vector<long> weightsFix(weights.size(), 0);
939 for (
unsigned iw = 0; iw < weights.size(); ++iw) {
940 weightsFix[iw] = long(round(weights[iw] * (1 << precisionWeights)));
943 unsigned xMax = unsigned(ceil(atanh(1. - 1. / (1 << (precisionTanh + 1))) *
944 (1 << (precisionLUT + 1))));
947 vector<long> layerinput = inputFix;
948 vector<long> layeroutput = {};
950 for (
unsigned il = 1; il < expert.nLayers(); ++il) {
952 layerinput.push_back(1 << precisionInput);
955 layeroutput.assign(expert.nNodesLayer(il), 0);
957 for (
unsigned io = 0; io < layeroutput.size(); ++io) {
959 for (
unsigned ii = 0; ii < layerinput.size(); ++ii) {
960 layeroutput[io] += layerinput[ii] * weightsFix[iw++];
963 unsigned long bin = abs(layeroutput[io]) >> dp;
965 float x = (bin + 0.5 - 1. / (1 << (dp + 1))) / (1 << precisionLUT);
966 long tanhLUT = (bin < xMax) ?
long(round(tanh(x / 2.) * (1 << precisionTanh))) : (1 << precisionTanh);
967 layeroutput[io] = (layeroutput[io] < 0) ? -tanhLUT : tanhLUT;
970 layerinput = layeroutput;
974 vector<float> output(layeroutput.size(), 0.);
975 for (
unsigned io = 0; io < output.size(); ++io) {
976 output[io] = layeroutput[io] / float(1 << precisionTanh);
978 return expert.unscaleTarget(output);
984 B2INFO(
"Saving networks to file " << filename <<
", array " << arrayname);
985 TFile datafile(filename.c_str(),
"UPDATE");
986 TObjArray* MLPs =
new TObjArray(m_MLPs.size());
987 for (
unsigned isector = 0; isector < m_MLPs.size(); ++isector) {
988 MLPs->Add(&m_MLPs[isector]);
990 MLPs->Write(arrayname.c_str(), TObject::kSingleKey | TObject::kOverwrite);
999 if (filename.size() < 1) {
1001 m_MLPs = m_cdctriggerneuroconfig->getMLPs();
1002 if (m_MLPs.size() == 0) {
1003 B2ERROR(
"Could not load Neurotrigger weights from database!");
1006 B2DEBUG(2,
"Loaded Neurotrigger MLP weights from database: " + m_cdctriggerneuroconfig->getNNName());
1007 B2DEBUG(100,
"loaded " << m_MLPs.size() <<
" networks from database");
1012 TFile datafile(filename.c_str(),
"READ");
1013 if (!datafile.IsOpen()) {
1014 B2WARNING(
"Could not open file " << filename);
1018 TObjArray* MLPs = (TObjArray*)datafile.Get(arrayname.c_str());
1021 B2WARNING(
"File " << filename <<
" does not contain key " << arrayname);
1025 for (
int isector = 0; isector < MLPs->GetEntriesFast(); ++isector) {
1027 if (expert) m_MLPs.push_back(*expert);
1028 else B2WARNING(
"Wrong type " << MLPs->At(isector)->ClassName() <<
", ignoring this entry.");
1033 B2DEBUG(100,
"loaded " << m_MLPs.size() <<
" networks");