8#include <framework/logging/Logger.h>
9#include <trg/cdc/NeuroTrigger.h>
10#include <trg/cdc/NeuroTriggerParameters.h>
11#include <trg/cdc/dataobjects/CDCTriggerMLPData.h>
12#include <cdc/geometry/CDCGeometryPar.h>
13#include <framework/gearbox/Const.h>
14#include <framework/gearbox/Unit.h>
15#include <framework/datastore/StoreObjPtr.h>
16#include <framework/datastore/StoreArray.h>
17#include <trg/cdc/dbobjects/CDCTriggerNeuroConfig.h>
18#include <trg/cdc/dataobjects/CDCTriggerTrack.h>
22#include "boost/iostreams/filter/gzip.hpp"
23#include "boost/iostreams/filtering_streambuf.hpp"
24#include "boost/iostreams/filtering_stream.hpp"
25#include "boost/multi_array.hpp"
26#define BOOST_MULTI_ARRAY_NO_GENERATORS
40 if (p.nHidden.size() != 1 && p.nHidden.size() != p.nMLP) {
41 B2ERROR(
"Number of nHidden lists should be 1 or " << p.nMLP);
44 if (p.outputScale.size() != 1 && p.outputScale.size() != p.nMLP) {
45 B2ERROR(
"Number of outputScale lists should be 1 or " << p.nMLP);
48 bool rangeProduct = (p.phiRangeUse.size() * p.invptRangeUse.size() * p.thetaRangeUse.size() * p.SLpattern.size() == p.nMLP);
50 if (p.phiRangeUse.size() != 1 && p.phiRangeUse.size() != p.nMLP) {
51 B2ERROR(
"Number of phiRangeUse.lists should be 1 or " << p.nMLP);
54 if (p.invptRangeUse.size() != 1 && p.invptRangeUse.size() != p.nMLP) {
55 B2ERROR(
"Number of invptRangeUse.lists should be 1 or " << p.nMLP);
58 if (p.thetaRangeUse.size() != 1 && p.thetaRangeUse.size() != p.nMLP) {
59 B2ERROR(
"Number of thetaRangeUse.lists should be 1 or " << p.nMLP);
62 if (p.SLpattern.size() != 1 && p.SLpattern.size() != p.nMLP) {
63 B2ERROR(
"Number of SLpattern lists should be 1 or " << p.nMLP);
68 if (p.maxHitsPerSL.size() != 1 && p.maxHitsPerSL.size() != p.SLpattern.size()) {
69 B2ERROR(
"Number of maxHitsPerSL lists should be 1 or " << p.SLpattern.size());
72 if (p.SLpatternMask.size() != 1 && p.SLpatternMask.size() != p.SLpattern.size()) {
73 B2ERROR(
"Number of SLpatternMask lists should be 1 or " << p.SLpattern.size());
77 if (p.targetZ.isSet() || p.targetTheta.isSet()) {
78 unsigned short nTarget = int(p.targetZ) + int(p.targetTheta);
80 B2ERROR(
"No outputs! Turn on either targetZ or targetTheta.");
85 for (
unsigned iPhi = 0; iPhi < p.phiRangeUse.size(); ++iPhi) {
86 if (p.phiRangeUse[iPhi].size() != 2) {
87 B2ERROR(
"phiRangeUse should be exactly 2 values");
91 if (p.phiRangeUse[iPhi][0] >= p.phiRangeUse[iPhi][1]) {
92 B2ERROR(
"phiRangeUse should be smaller than phiRangeUse");
95 if (p.phiRangeUse[iPhi][0] < -360. || p.phiRangeUse[iPhi][1] > 360. ||
96 (p.phiRangeUse[iPhi][1] - p.phiRangeUse[iPhi][0]) > 360.) {
97 B2ERROR(
"phiRangeUse should be in [-360, 360], with maximal width of 360");
101 for (
unsigned iPt = 0; iPt < p.invptRangeUse.size(); ++iPt) {
102 if (p.invptRangeUse[iPt].size() != 2) {
103 B2ERROR(
"invptRangeUse should be exactly 2 values");
106 if (p.invptRangeUse[iPt][0] >= p.invptRangeUse[iPt][1]) {
107 B2ERROR(
"invptRangeUse should be smaller than invptRangeUse");
111 for (
unsigned iTheta = 0; iTheta < p.thetaRangeUse.size(); ++iTheta) {
112 if (p.thetaRangeUse[iTheta].size() != 2) {
113 B2ERROR(
"thetaRangeUse should be exactly 2 values");
117 if (p.thetaRangeUse[iTheta][0] >= p.thetaRangeUse[iTheta][1]) {
118 B2ERROR(
"thetaRangeUse should be smaller than thetaRangeUse");
121 if (p.thetaRangeUse[iTheta][0] < 0. || p.thetaRangeUse[iTheta][1] > 180.) {
122 B2ERROR(
"thetaRangeUse should be in [0, 180]");
126 int nTarget = (int) p.targetZ + (
int) p.targetTheta;
127 if (p.outputScale.size() == p.nMLP || p.outputScale.size() == 1) {
128 for (
unsigned iScale = 0; iScale < p.outputScale.size(); ++iScale) {
129 if (p.outputScale[iScale].size() != 2 *
static_cast<unsigned int>(nTarget)) {
130 B2ERROR(
"outputScale should be exactly " << 2 * nTarget <<
" values");
135 B2ERROR(
"the size of outputscale should be 1 or match the number of experts");
138 if (p.phiRangeUse.size() != p.phiRangeTrain.size()) {
139 B2ERROR(
"Number of phiRangeUse.lists and phiRangeTrain lists should be equal.");
142 for (
unsigned iPhi = 0; iPhi < p.phiRangeUse.size(); ++iPhi) {
143 if (p.phiRangeTrain[iPhi].size() != 2) {
144 B2ERROR(
"phiRangeTrain should be exactly 2 values.");
146 }
else if (p.phiRangeTrain[iPhi][0] > p.phiRangeUse[iPhi][0] ||
147 p.phiRangeTrain[iPhi][1] < p.phiRangeUse[iPhi][1]) {
148 B2ERROR(
"phiRangeTrain should be wider than phiRangeUse.or equal.");
153 if (p.invptRangeUse.size() != p.invptRangeTrain.size()) {
154 B2ERROR(
"Number of invptRangeUse.lists and invptRangeTrain lists should be equal.");
157 for (
unsigned iPt = 0; iPt < p.invptRangeUse.size(); ++iPt) {
158 if (p.invptRangeTrain[iPt].size() != 2) {
159 B2ERROR(
"invptRangeTrain should be exactly 2 values.");
161 }
else if (p.invptRangeTrain[iPt][0] > p.invptRangeUse[iPt][0] ||
162 p.invptRangeTrain[iPt][1] < p.invptRangeUse[iPt][1]) {
163 B2ERROR(
"invptRangeTrain should be wider than invptRangeUse.or equal.");
168 if (p.thetaRangeUse.size() != p.thetaRangeTrain.size()) {
169 B2ERROR(
"Number of thetaRangeUse.lists and thetaRangeTrain lists should be equal.");
172 for (
unsigned iTheta = 0; iTheta < p.thetaRangeUse.size(); ++iTheta) {
173 if (p.thetaRangeTrain[iTheta].size() != 2) {
174 B2ERROR(
"thetaRangeTrain should be exactly 2 values.");
176 }
else if (p.thetaRangeTrain[iTheta][0] > p.thetaRangeUse[iTheta][0] ||
177 p.thetaRangeTrain[iTheta][1] < p.thetaRangeUse[iTheta][1]) {
178 B2ERROR(
"thetaRangeTrain should be wider than thetaRangeUse.or equal.");
188 for (
unsigned iMLP = 0; iMLP < p.nMLP; ++iMLP) {
193 unsigned short maxHits = (p.maxHitsPerSL.size() == 1) ? p.maxHitsPerSL[0] : p.maxHitsPerSL[indices[3]];
194 unsigned long SLpattern = p.SLpattern[indices[3]];
195 unsigned long SLpatternMask = (p.SLpatternMask.size() == 1) ? p.SLpatternMask[0] : p.SLpatternMask[indices[3]];
196 unsigned short nInput = 27 * maxHits;
197 vector<
NNTParam<float>> nHidden = (p.nHidden.size() == 1) ? p.nHidden[0] : p.nHidden[iMLP];
198 vector<unsigned short> nNodes = {nInput};
199 for (
unsigned iHid = 0; iHid < nHidden.size(); ++iHid) {
200 if (p.multiplyHidden) {
201 nNodes.push_back(nHidden[iHid] * nNodes[0]);
203 nNodes.push_back(nHidden[iHid]);
206 nNodes.push_back(nTarget);
207 unsigned short targetVars = int(p.targetZ) + (int(p.targetTheta) << 1);
209 vector<float> phiRangeUse = p.tcastvector<
float>(p.phiRangeUse)[indices[0]];
210 vector<float> invptRangeUse = p.tcastvector<
float>(p.invptRangeUse)[indices[1]];
211 vector<float> thetaRangeUse = p.tcastvector<
float>(p.thetaRangeUse)[indices[2]];
212 vector<float> phiRangeTrain = p.tcastvector<
float>(p.phiRangeTrain)[indices[0]];
213 vector<float> invptRangeTrain = p.tcastvector<
float>(p.invptRangeTrain)[indices[1]];
214 vector<float> thetaRangeTrain = p.tcastvector<
float>(p.thetaRangeTrain)[indices[2]];
215 B2DEBUG(50,
"Ranges for sector " << iMLP
216 <<
": phiRange [" << phiRangeUse[0] <<
", " << phiRangeUse[1]
217 <<
"], invptRange [" << invptRangeUse[0] <<
", " << invptRangeUse[1]
218 <<
"], thetaRange [" << thetaRangeUse[0] <<
", " << thetaRangeUse[1]
219 <<
"], SLpattern " << SLpattern);
221 vector<float> outputScale = (p.outputScale.size() == 1) ? p.tcastvector<
float>(p.outputScale)[0] : p.tcastvector<
float>
222 (p.outputScale)[iMLP];
233 outputScale[2 * int(p.targetZ)] *=
Unit::deg;
234 outputScale[2 * int(p.targetZ) + 1] *=
Unit::deg;
238 phiRangeUse, invptRangeUse, thetaRangeUse,
239 phiRangeTrain, invptRangeTrain, thetaRangeTrain,
240 maxHits, SLpattern, SLpatternMask, p.tMax,
244 if (p.IDRanges.size() == p.nMLP) {
245 for (
auto exp : p.IDRanges) {
247 std::vector<float> irange = {exp.begin() + 1, exp.end()};
248 m_MLPs[
static_cast<int>(exp[0])].setRelID(irange);
250 }
else if (p.IDRanges.size() == 0) {
251 B2WARNING(
"idranges have not been initialized yet, did you forget it?");
253 B2ERROR(
"number of idranges should match the number of experts!");
263 std::vector<unsigned> indices = {0, 0, 0, 0};
264 if (p.phiRangeUse.size() * p.invptRangeUse.size() * p.thetaRangeUse.size() * p.SLpattern.size() == p.nMLP) {
265 indices[0] = isector % p.phiRangeUse.size();
266 indices[1] = (isector / p.phiRangeUse.size()) % p.invptRangeUse.size();
267 indices[2] = (isector / (p.phiRangeUse.size() * p.invptRangeUse.size())) % p.thetaRangeUse.size();
268 indices[3] = isector / (p.phiRangeUse.size() * p.invptRangeUse.size() * p.thetaRangeUse.size());
270 indices[0] = (p.phiRangeUse.size() == 1) ? 0 : isector;
271 indices[1] = (p.invptRangeUse.size() == 1) ? 0 : isector;
272 indices[2] = (p.thetaRangeUse.size() == 1) ? 0 : isector;
273 indices[3] = (p.SLpattern.size() == 1) ? 0 : isector;
285 for (
int iSL = 0; iSL < 9; ++iSL) {
287 nTS += cdc.nWiresInLayer(layerId);
289 for (
int priority = 0; priority < 2; ++ priority) {
290 m_radius[iSL][priority] = cdc.senseWireR(layerId + priority);
292 layerId += (iSL > 0 ? 6 : 7);
300 if (!((et_option ==
"fastestpriority") || (et_option ==
"etfhwin") || (et_option ==
"zero") || (et_option ==
"fastest2d"))) {
315 vector<int> indices = {};
318 B2WARNING(
"Trying to select MLP before initializing MLPs.");
323 for (
unsigned isector = 0; isector <
m_MLPs.size(); ++isector) {
324 if (
m_MLPs[isector].inPhiRangeTrain(phi0) &&
m_MLPs[isector].inInvptRangeTrain(invpt)
325 &&
m_MLPs[isector].inThetaRangeTrain(theta)) {
326 indices.push_back(isector);
330 if (indices.size() == 0) {
331 B2DEBUG(150,
"Track does not match any sector.");
332 B2DEBUG(150,
"invpt=" << invpt <<
", phi=" << phi0 * 180. / M_PI <<
", theta=" << theta * 180. / M_PI);
340 vector<int> indices = {};
343 B2WARNING(
"Trying to select MLP before initializing MLPs.");
348 for (
unsigned isector = 0; isector <
m_MLPs.size(); ++isector) {
349 if (
m_MLPs[isector].inPhiRangeUse(phi0) &&
m_MLPs[isector].inInvptRangeUse(invpt)
350 &&
m_MLPs[isector].inThetaRangeUse(theta)) {
351 indices.push_back(isector);
355 if (indices.size() == 0) {
356 B2DEBUG(150,
"Track does not match any sector.");
357 B2DEBUG(150,
"invpt=" << invpt <<
", phi=" << phi0 * 180. / M_PI <<
", theta=" << theta * 180. / M_PI);
367 if (MLPs.size() == 0) {
372 for (
unsigned i = 0; i < MLPs.size(); ++i) {
373 int isector = MLPs[i];
374 unsigned long sectorPattern =
m_MLPs[isector].getSLpattern();
375 B2DEBUG(250,
"hitPattern " << pattern <<
" sectorPattern " << sectorPattern);
377 if (sectorPattern == 0) {
378 B2DEBUG(250,
"found match for general sector");
382 if (pattern == sectorPattern) {
383 B2DEBUG(250,
"found match for hit pattern " << pattern);
390 if (neurotrackinputmode) {
391 B2DEBUG(150,
"No sector found to match pattern, using sector 0" << pattern <<
".");
394 B2DEBUG(150,
"No sector found to match pattern " << pattern <<
".");
404 double omega = track.getOmega();
405 for (
int iSL = 0; iSL < 9; ++iSL) {
406 for (
int priority = 0; priority < 2; ++priority) {
408 asin(
m_radius[iSL][priority] * omega / 2.) :
410 m_idRef[iSL][priority] = remainder(((track.getPhi0() -
m_alpha[iSL][priority]) *
426 double omega = track.getOmega();
427 long phi = round(track.getPhi0() * (1 << precisionPhiAlpha));
429 for (
int iSL = 0; iSL < 9; ++iSL) {
430 for (
int priority = 0; priority < 2; ++priority) {
433 (
m_radius[iSL][priority] * abs(omega) < 2) ?
434 round(asin(
m_radius[iSL][priority] * omega / 2) * (1 << precisionPhiAlpha)) :
435 round(M_PI_2 * (1 << precisionPhiAlpha));
436 long dphi = phi - (long(
m_alpha[iSL][priority]));
438 long(dphi * round((
m_TSoffset[iSL + 1] -
m_TSoffset[iSL]) / 2. / M_PI * (1 << precisionScale)) /
439 (
long(1) << (precisionPhiAlpha + precisionScale - precisionId)));
441 m_alpha[iSL][priority] /= (1 << precisionPhiAlpha);
442 m_idRef[iSL][priority] /= (1 << precisionId);
450 int iSL = hit.getISuperLayer();
451 int priority = hit.getPriorityPosition();
453 double relId = hit.getSegmentID() + 0.5 * (((priority >> 1) & 1) - (priority & 1))
464 B2DEBUG(200,
"looping over axials:");
465 for (
unsigned ihit = 0; ihit < Hits.
size(); ++ihit) {
467 if (Hits.
weight(ihit) < 0)
continue;
468 unsigned short iSL = Hits[ihit]->getISuperLayer();
469 if (iSL % 2 == 1 && onlyAxials) {
continue;}
471 B2DEBUG(200,
" check drifttime: SL" + std::to_string(iSL) +
",ID = " + std::to_string(Hits[ihit]->getSegmentID()) +
", t = " +
472 std::to_string(Hits[ihit]->priorityTime()));
473 double relId =
getRelId(*Hits[ihit]);
474 if (
m_MLPs[isector].isRelevant(relId, iSL) &&
475 Hits[ihit]->priorityTime() < tlow) {
476 tlow = Hits[ihit]->priorityTime();
477 B2DEBUG(200,
" new tlow: " << std::to_string(tlow));
489 B2DEBUG(20,
"Used event time option is different to the one set in the MLP"
490 <<
LogVar(
"et_option", et_option) <<
LogVar(
"isector", isector)
493 if (et_option ==
"fastestpriority") {
494 B2DEBUG(200,
"et_option is 'fastestpriority'");
507 }
else if (et_option ==
"fastest2d") {
519 }
else if (et_option ==
"zero") {
522 }
else if (et_option ==
"etf_only") {
528 B2ERROR(
"No ETF output, but forced to use ETF!");
532 }
else if (et_option ==
"etf_or_fastestpriority") {
538 getEventTime(isector, track,
"fastestpriority", neuroinputmode);
552 }
else if (et_option ==
"min_etf_fastestpriority") {
559 getEventTime(isector, track,
"fastestpriority", neuroinputmode);
563 }
else if (et_option ==
"etf_or_fastest2d") {
569 getEventTime(isector, track,
"fastest2d", neuroinputmode);
583 }
else if (et_option ==
"etf_or_zero") {
592 }
else if (et_option ==
"etfcc") {
593 if (track.getHasETFTime()) {
594 m_T0 = track.getETF_unpacked();
601 }
else if (et_option ==
"etfcc_or_zero") {
602 if (track.getHasETFTime()) {
603 m_T0 = track.getETF_unpacked();
610 }
else if (et_option ==
"etfcc_or_fastestpriority") {
611 if (track.getHasETFTime()) {
612 m_T0 = track.getETF_unpacked();
615 getEventTime(isector, track,
"fastestpriority", neuroinputmode);
618 }
else if (et_option ==
"etfhwin") {
619 m_T0 = track.getETF_recalced();
623 B2ERROR(
"No valid parameter for et_option (" << et_option <<
" )!");
632 unsigned long driftth = 0;
633 vector<unsigned> nHits;
638 for (
unsigned ihit = 0; ihit < axialHits.
size(); ++ ihit) {
640 if (axialHits.
weight(ihit) < 0) {
643 unsigned short iSL = axialHits[ihit]->getISuperLayer();
645 if ((!neurotrackinputmode) && (iSL % 2 == 1))
continue;
647 int t = (
m_hasT0) ? axialHits[ihit]->priorityTime() -
m_T0 : 0;
648 double relId =
getRelId(*axialHits[ihit]);
649 if (t < 0 || t > expert.getTMax()) {
650 if (expert.isRelevant(relId, iSL)) {
651 if (nHits[iSL] < expert.getMaxHitsPerSL()) {
652 driftth |= 1 << (8 - iSL + 9 * nHits[iSL]);
658 if (!neurotrackinputmode) {
663 if (iSL % 2 == 0)
continue;
667 if (t < 0 || t > expert.getTMax()) {
668 if (expert.isRelevant(relId, iSL)) {
669 if (nHits[iSL] < expert.getMaxHitsPerSL()) {
670 driftth |= 1 << (8 - iSL + 9 * nHits[iSL]);
682 unsigned long chitPattern = 0;
683 vector<unsigned> nHits;
688 for (
unsigned ihit = 0; ihit < axialHits.
size(); ++ ihit) {
690 if (axialHits.
weight(ihit) < 0) {
693 unsigned short iSL = axialHits[ihit]->getISuperLayer();
695 if ((!neurotrackinputmode) && (iSL % 2 == 1))
continue;
696 double relId =
getRelId(*axialHits[ihit]);
697 if (expert.isRelevant(relId, iSL)) {
698 if (nHits[iSL] < expert.getMaxHitsPerSL()) {
699 chitPattern |= 1 << (iSL + 9 * nHits[iSL]);
704 if (!neurotrackinputmode) {
709 if (iSL % 2 == 0)
continue;
712 if (expert.isRelevant(relId, iSL)) {
713 if (nHits[iSL] < expert.getMaxHitsPerSL()) {
714 chitPattern |= 1 << (iSL + 9 * nHits[iSL]);
727 unsigned long hitPattern = 0;
728 vector<unsigned> nHits;
733 for (
unsigned ihit = 0; ihit < axialHits.
size(); ++ ihit) {
735 if (axialHits.
weight(ihit) < 0) {
738 unsigned short iSL = axialHits[ihit]->getISuperLayer();
740 if ((!neurotrackinputmode) && (iSL % 2 == 1))
continue;
742 int t = (
m_hasT0) ? axialHits[ihit]->priorityTime() -
m_T0 : 0;
745 }
else if (t > expert.getTMax()) {
746 t = expert.getTMax();
748 double relId =
getRelId(*axialHits[ihit]);
750 if (expert.isRelevant(relId, iSL)) {
751 if (nHits[iSL] < expert.getMaxHitsPerSL()) {
752 hitPattern |= 1 << (iSL + 9 * nHits[iSL]);
755 B2DEBUG(250,
"hit in SL " << iSL);
757 B2DEBUG(250,
"hit in SL " << iSL <<
" not relevant (relId = " << relId <<
")");
760 if (!neurotrackinputmode) {
765 if (iSL % 2 == 0)
continue;
770 }
else if (t > expert.getTMax()) {
771 t = expert.getTMax();
774 if (expert.isRelevant(relId, iSL)) {
775 if (nHits[iSL] < expert.getMaxHitsPerSL()) {
776 hitPattern |= 1 << (iSL + 9 * nHits[iSL]);
779 B2DEBUG(250,
"hit in SL " << iSL);
781 B2DEBUG(250,
"hit in SL " << iSL <<
" not relevant (relId = " << relId <<
")");
785 B2DEBUG(250,
"hitPattern " << hitPattern);
786 return hitPattern & expert.getSLpatternMask();
793 vector<unsigned> selectedHitIds = {};
796 tMin.assign(expert.nNodesLayer(0), expert.getTMax());
797 vector<bool> LRknown;
798 LRknown.assign(expert.nNodesLayer(0),
false);
800 hitIds.assign(expert.nNodesLayer(0), -1);
801 vector<unsigned> nHits;
807 B2DEBUG(250,
"start hit loop over all related hits");
812 for (
unsigned ihit = 0; ihit < allHits.
size(); ++ihit) {
814 if (allHits.
weight(ihit) < 0)
continue;
815 unsigned short iSL = allHits[ihit]->getISuperLayer();
816 if (expert.getSLpatternUnmasked() != 0 &&
817 !((expert.getSLpatternUnmasked() >> iSL) & 1)) {
818 B2DEBUG(250,
"skipping hit in SL " << iSL);
822 int t = (
m_hasT0) ? allHits[ihit]->priorityTime() -
m_T0 : 0;
825 }
else if (t > expert.getTMax()) {
826 t = expert.getTMax();
828 double relId =
getRelId(*allHits[ihit]);
829 if (expert.isRelevant(relId, iSL)) {
831 unsigned short iRef = iSL;
832 if (expert.getMaxHitsPerSL() > 1) {
833 if (nHits[iSL] < expert.getMaxHitsPerSL() &&
834 (expert.getSLpatternUnmasked() >> (iSL + 9 * nHits[iSL])) & 1) {
835 iRef += 9 * nHits[iSL];
838 for (
unsigned compare = iSL; compare < iSL + 9 * nHits[iSL]; compare += 9) {
839 if ((LRknown[iRef] && !LRknown[compare]) ||
840 (LRknown[iRef] == LRknown[compare] && tMin[iRef] < tMin[compare]))
848 useHit = (allHits[ihit]->LRknown() && t <= tMin[iRef]);
850 useHit = (allHits[ihit]->LRknown() || t <= tMin[iRef]);
852 B2DEBUG(250,
"relevant wire SL " << iSL <<
" LR " << allHits[ihit]->getLeftRight()
853 <<
" t " << t <<
" iRef " << iRef <<
" useHit " << useHit);
856 LRknown[iRef] = allHits[ihit]->LRknown();
858 hitIds[iRef] = allHits[ihit]->getArrayIndex();
863 for (
unsigned iHit = 0; iHit < hitIds.size(); ++iHit) {
864 if (hitIds[iHit] >= 0) selectedHitIds.push_back(hitIds[iHit]);
866 return selectedHitIds;
870 bool returnAllRelevant)
873 vector<unsigned> selectedHitIds = {};
876 tMin.assign(expert.nNodesLayer(0), expert.getTMax());
877 vector<bool> LRknown;
878 LRknown.assign(expert.nNodesLayer(0),
false);
880 hitIds.assign(expert.nNodesLayer(0), -1);
881 vector<unsigned> nHits;
887 B2DEBUG(250,
"start axial hit loop");
888 for (
unsigned ihit = 0; ihit < axialHits.
size(); ++ihit) {
890 if (axialHits.
weight(ihit) < 0)
continue;
891 unsigned short iSL = axialHits[ihit]->getISuperLayer();
893 if (iSL % 2 == 1)
continue;
894 if (expert.getSLpatternUnmasked() != 0 &&
895 !((expert.getSLpatternUnmasked() >> iSL) & 1)) {
896 B2DEBUG(250,
"skipping hit in SL " << iSL);
901 int t = (
m_hasT0) ? axialHits[ihit]->priorityTime() -
m_T0 : 0;
904 }
else if (t > expert.getTMax()) {
905 t = expert.getTMax();
907 double relId =
getRelId(*axialHits[ihit]);
908 if (expert.isRelevant(relId, iSL)) {
910 unsigned short iRef = iSL;
911 if (expert.getMaxHitsPerSL() > 1) {
912 if (nHits[iSL] < expert.getMaxHitsPerSL() &&
913 (expert.getSLpatternUnmasked() >> (iSL + 9 * nHits[iSL])) & 1) {
914 iRef += 9 * nHits[iSL];
917 for (
unsigned compare = iSL; compare < iSL + 9 * nHits[iSL]; compare += 9) {
918 if ((LRknown[iRef] && !LRknown[compare]) ||
919 (LRknown[iRef] == LRknown[compare] && tMin[iRef] < tMin[compare]))
927 useHit = (axialHits[ihit]->LRknown() && t <= tMin[iRef]);
929 useHit = (axialHits[ihit]->LRknown() || t <= tMin[iRef]);
931 B2DEBUG(250,
"relevant wire SL " << iSL <<
" LR " << axialHits[ihit]->getLeftRight()
932 <<
" t " << t <<
" iRef " << iRef <<
" useHit " << useHit);
935 LRknown[iRef] = axialHits[ihit]->LRknown();
937 hitIds[iRef] = axialHits[ihit]->getArrayIndex();
939 if (returnAllRelevant) selectedHitIds.push_back(axialHits[ihit]->getArrayIndex());
944 B2DEBUG(250,
"start stereo hit loop");
948 if (iSL % 2 == 0)
continue;
949 if (expert.getSLpatternUnmasked() != 0 &&
950 !((expert.getSLpatternUnmasked() >> iSL) & 1)) {
951 B2DEBUG(250,
"skipping hit in SL " << iSL);
958 }
else if (t > expert.getTMax()) {
959 t = expert.getTMax();
962 if (expert.isRelevant(relId, iSL)) {
964 unsigned short iRef = iSL;
965 if (expert.getMaxHitsPerSL() > 1) {
966 if (nHits[iSL] < expert.getMaxHitsPerSL() &&
967 (expert.getSLpatternUnmasked() >> (iSL + 9 * nHits[iSL])) & 1) {
968 iRef += 9 * nHits[iSL];
971 for (
unsigned compare = iSL; compare < iSL + 9 * nHits[iSL]; compare += 9) {
972 if ((LRknown[iRef] && !LRknown[compare]) ||
973 (LRknown[iRef] == LRknown[compare] && tMin[iRef] < tMin[compare]))
981 useHit = (
m_segmentHits[ihit]->LRknown() && t <= tMin[iRef]);
983 useHit = (
m_segmentHits[ihit]->LRknown() || t <= tMin[iRef]);
985 B2DEBUG(250,
"relevant wire SL " << iSL <<
" LR " <<
m_segmentHits[ihit]->getLeftRight()
986 <<
" t " << t <<
" iRef " << iRef <<
" useHit " << useHit);
993 if (returnAllRelevant) selectedHitIds.push_back(ihit);
998 if (!returnAllRelevant) {
999 for (
unsigned iHit = 0; iHit < hitIds.size(); ++iHit) {
1000 if (hitIds[iHit] >= 0) selectedHitIds.push_back(hitIds[iHit]);
1003 return selectedHitIds;
1011 vector<float> inputVector;
1012 inputVector.assign(expert.nNodesLayer(0), 0.);
1014 vector<unsigned> nHits;
1016 for (
unsigned ii = 0; ii < hitIds.size(); ++ii) {
1017 int ihit = hitIds[ii];
1019 unsigned short iRef = iSL + 9 * nHits[iSL];
1025 }
else if (t > expert.getTMax()) {
1026 t = expert.getTMax();
1032 inputVector[3 * iRef] = expert.scaleId(relId, iSL);
1033 float scaleT = pow(2, floor(log2(1. / expert.getTMax())));
1034 inputVector[3 * iRef + 1] = (((LR >> 1) & 1) - (LR & 1)) * t * scaleT;
1035 inputVector[3 * iRef + 2] =
m_alpha[iSL][int(priority < 3)] * 0.5;
1044 vector<float> weights = expert.getWeights();
1045 vector<float> layerinput = input;
1046 vector<float> layeroutput = {};
1048 for (
unsigned il = 1; il < expert.nLayers(); ++il) {
1050 layerinput.push_back(1.);
1052 layeroutput.clear();
1053 layeroutput.assign(expert.nNodesLayer(il), 0.);
1055 for (
unsigned io = 0; io < layeroutput.size(); ++io) {
1057 for (
unsigned ii = 0; ii < layerinput.size(); ++ii) {
1058 layeroutput[io] += layerinput[ii] * weights[iw++];
1061 layeroutput[io] = tanh(layeroutput[io] / 2.);
1064 layerinput = layeroutput;
1066 return expert.unscaleTarget(layeroutput);
1076 unsigned dp = precisionInput + precisionWeights - precisionLUT;
1080 vector<long> inputFix(input.size(), 0);
1081 for (
unsigned ii = 0; ii < input.size(); ++ii) {
1082 inputFix[ii] = long(input[ii] * (1 << precisionInput));
1085 vector<float> weights = expert.getWeights();
1086 vector<long> weightsFix(weights.size(), 0);
1087 for (
unsigned iw = 0; iw < weights.size(); ++iw) {
1088 weightsFix[iw] = long(round(weights[iw] * (1 << precisionWeights)));
1091 unsigned xMax = unsigned(ceil(atanh(1. - 1. / (1 << (precisionTanh + 1))) *
1092 (1 << (precisionLUT + 1))));
1095 vector<long> layerinput = inputFix;
1096 vector<long> layeroutput = {};
1098 for (
unsigned il = 1; il < expert.nLayers(); ++il) {
1100 layerinput.push_back(1 << precisionInput);
1102 layeroutput.clear();
1103 layeroutput.assign(expert.nNodesLayer(il), 0);
1105 for (
unsigned io = 0; io < layeroutput.size(); ++io) {
1107 for (
unsigned ii = 0; ii < layerinput.size(); ++ii) {
1108 layeroutput[io] += layerinput[ii] * weightsFix[iw++];
1111 unsigned long bin = abs(layeroutput[io]) >> dp;
1113 float x = (bin + 0.5 - 1. / (1 << (dp + 1))) / (1 << precisionLUT);
1114 long tanhLUT = (bin < xMax) ?
long(round(tanh(x / 2.) * (1 << precisionTanh))) : (1 << precisionTanh);
1115 layeroutput[io] = (layeroutput[io] < 0) ? -tanhLUT : tanhLUT;
1118 layerinput = layeroutput;
1122 vector<float> output(layeroutput.size(), 0.);
1123 for (
unsigned io = 0; io < output.size(); ++io) {
1124 output[io] = layeroutput[io] / float(1 << precisionTanh);
1126 return expert.unscaleTarget(output);
1132 B2INFO(
"Saving networks to file " << filename <<
", array " << arrayname);
1133 TFile datafile(filename.c_str(),
"UPDATE");
1134 TObjArray* MLPs =
new TObjArray(
m_MLPs.size());
1135 for (
unsigned isector = 0; isector <
m_MLPs.size(); ++isector) {
1136 MLPs->Add(&
m_MLPs[isector]);
1138 MLPs->Write(arrayname.c_str(), TObject::kSingleKey | TObject::kOverwrite);
1145 std::ifstream gzipfile(filename, ios_base::in | ios_base::binary);
1146 boost::iostreams::filtering_istream arrayStream;
1147 arrayStream.push(boost::iostreams::gzip_decompressor());
1148 arrayStream.push(gzipfile);
1150 if (gzipfile.is_open()) {
1151 while (arrayStream >> hline) {
1152 for (
unsigned i = 0; i < 18; ++i) {
1153 m_MLPs[hline.exPert].relevantID[i] = hline.relID[i];
1156 }
else {
return false;}
1164 if (filename.size() < 1) {
1167 if (
m_MLPs.size() == 0) {
1168 B2ERROR(
"Could not load Neurotrigger weights from database!");
1172 B2DEBUG(100,
"loaded " <<
m_MLPs.size() <<
" networks from database");
1177 TFile datafile(filename.c_str(),
"READ");
1178 if (!datafile.IsOpen()) {
1179 B2WARNING(
"Could not open file " << filename);
1183 TObjArray* MLPs = (TObjArray*)datafile.Get(arrayname.c_str());
1186 B2WARNING(
"File " << filename <<
" does not contain key " << arrayname);
1190 for (
int isector = 0; isector < MLPs->GetEntriesFast(); ++isector) {
1192 if (expert)
m_MLPs.push_back(*expert);
1193 else B2WARNING(
"Wrong type " << MLPs->At(isector)->ClassName() <<
", ignoring this entry.");
1198 B2DEBUG(100,
"loaded " <<
m_MLPs.size() <<
" networks");
1200 B2INFO(
"Loaded Neurotrigger MLP weights from file: " + filename);
bool hasBinnedEventT0(const Const::EDetector detector) const
Check if one of the detectors in the given set has a binned t0 estimation.
int getBinnedEventT0(const Const::EDetector detector) const
Return the stored binned event t0 for the given detector or 0 if nothing stored.
Class to keep all parameters of an expert MLP for the neuro trigger.
Combination of several CDCHits to a track segment hit for the trigger.
Track created by the CDC trigger.
The Class for CDC Geometry Parameters.
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
Class to represent a complete set to describe a Neurotrigger.
unsigned long getPureDriftThreshold(unsigned isector, const CDCTriggerTrack &track, const bool neurotrackinputmode)
Get the drift threshold bits, where the time of the TS was outside of the accepted time window and th...
double m_alpha[9][2]
2D crossing angle of current track
std::vector< CDCTriggerMLP > m_MLPs
List of networks.
double getRelId(const CDCTriggerSegmentHit &hit)
Calculate phi position of a hit relative to 2D track (scaled to number of wires).
DBObjPtr< CDCTriggerNeuroConfig > m_cdctriggerneuroconfig
get NNT payload from database.
void initialize(const Parameters &p)
Set parameters and get some network independent parameters.
std::vector< unsigned > getRangeIndices(const NeuroTriggerParameters &p, unsigned isector)
Get indices for sector ranges in parameter lists.
std::vector< int > selectMLPs(float phi0, float invpt, float theta)
Select all matching expert MLPs based on the given track parameters.
void updateTrack(const CDCTriggerTrack &track)
Calculate 2D phi position and arclength for the given track and store them.
bool load(const std::string &filename, const std::string &arrayname="MLPs")
Load MLPs from file.
std::vector< unsigned > selectHits(unsigned isector, const CDCTriggerTrack &track, bool returnAllRelevant=false)
Select best hits for each super layer.
std::vector< float > runMLPFix(unsigned isector, std::vector< float > input)
Run an expert MLP with fixed point arithmetic.
void initializeCollections(std::string hitCollectionName, std::string eventTimeName, const std::string &et_option)
set the hit collection and event time to required and store the hit collection name
std::vector< float > getInputVector(unsigned isector, const std::vector< unsigned > &hitIds)
Calculate input values for MLP.
StoreArray< CDCTriggerSegmentHit > m_segmentHits
StoreArray containing the input track segment hits.
int selectMLPbyPattern(std::vector< int > &MLPs, unsigned long pattern, const bool neurotrackinputmode)
Select one MLP from a list of sector indices.
std::vector< unsigned > m_precision
Fixed point precision in bit after radix point.
std::string get_et_option()
Return value of m_et_option.
unsigned long getInputPattern(unsigned isector, const CDCTriggerTrack &track, const bool neurotrackinputmode)
Calculate input pattern for MLP.
int m_T0
Event time of current event / track.
void save(const std::string &filename, const std::string &arrayname="MLPs")
Save MLPs to file.
void updateTrackFix(const CDCTriggerTrack &track)
Calculate 2D phi position and arclength for the given track and store them.
std::vector< unsigned > selectHitsHWSim(unsigned isector, const CDCTriggerTrack &track)
Select hits for each super layer from the ones related to input track.
bool loadIDHist(const std::string &filename)
function to load idhist from file
std::vector< int > selectMLPsTrain(float phi0, float invpt, float theta)
Select all matching expert MLPs based on the given track parameters.
void setConstants()
Loads parameters from the geometry and precalculates some constants that will be needed.
unsigned m_TSoffset[10]
Number of track segments up to super layer.
double m_radius[9][2]
Radius of the CDC layers with priority wires (2 per super layer)
StoreObjPtr< BinnedEventT0 > m_eventTime
StoreObjPtr containing the event time.
unsigned long getCompleteHitPattern(unsigned isector, const CDCTriggerTrack &track, const bool neurotrackinputmode)
Get complete hit pattern of neurotrack.
void getEventTime(unsigned isector, const CDCTriggerTrack &track, std::string et_option, const bool)
Read out the event time and store it.
int getLowestTime(unsigned isector, RelationVector< CDCTriggerSegmentHit > Hits, bool onlyAxials)
helper function to get the fastest priority time of given ts array
double m_idRef[9][2]
2D phi position of current track scaled to number of wires
std::string m_hitCollectionName
Name of the StoreArray containing the input track segment hits.
std::vector< float > runMLP(unsigned isector, const std::vector< float > &input)
Run an expert MLP.
bool m_hasT0
Flag to show if stored event time is valid.
Class for type safe access to objects that are referred to in relations.
size_t size() const
Get number of relations.
float weight(int index) const
Get weight with index.
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
int getEntries() const
Get the number of objects in the array.
bool isValid() const
Check whether the object was created.
static const double deg
degree to radians
Class to store variables with their name which were sent to the logging service.
Abstract base class for different kinds of events.