From 6c08a49cb7a0bd8cff03f3e838f91f0536261c99 Mon Sep 17 00:00:00 2001 From: rodwyer100 Date: Fri, 20 Sep 2024 06:37:57 -0700 Subject: [PATCH] Implementation of TS clustering and tracking as envisioned in LUT firmware. (#1461) * First commit of firmware changes * First commit of firmware changes * Update TrigScint/include/TrigScint/TrigScintFirmwareTracker.h Co-authored-by: Tamas Vami * Update TrigScint/src/TrigScint/TrigScintFirmwareTracker.cxx Co-authored-by: Tamas Vami * Update TrigScint/exampleConfigs/firmwareEx.py Co-authored-by: Tamas Vami * Update TrigScint/exampleConfigs/firmwareEx.py Co-authored-by: Tamas Vami * Commiting removal of print out statements * Commiting removal of superfluous firmware ap int libraries * Commiting removal of hard coded channel numbers * Forgot to add some changes that we unstagged w.r.t. the comments and CMakeFix. I haven't included all of the class name changes (capitalizations) I may eventually include * Committed smart pointer (for real this time). * Cleaning firmware objects to align with capitalization scheme * A few more things * Trying to do clang fixes * Reverting to old ldmx-env.sh * Changing to array structure, haven't validated this yes please don't push just now. * Apply clang-format * yay! all done * Apply clang-format Attempts were made to use std::array instead of passing around bare pointers. There are some C-style fixed length arrays (T arr[N]) mixed with STL fixed-length arrays (std::array) but we have avoided spurious `new` floating around. --------- Co-authored-by: rodwyer@stanford.edu Co-authored-by: Tamas Vami Co-authored-by: github-actions[bot] --- TrigScint/CMakeLists.txt | 4 +- TrigScint/exampleConfigs/firmwareEx.py | 168 ++++++++++ TrigScint/exampleConfigs/meganEx.py | 167 ++++++++++ .../TrigScint/Firmware/clusterproducer.h | 12 + TrigScint/include/TrigScint/Firmware/objdef.h | 134 ++++++++ .../include/TrigScint/Firmware/testutils.h | 21 ++ .../TrigScint/Firmware/trackproducer.h | 15 + .../TrigScint/TrigScintFirmwareTracker.h | 117 +++++++ TrigScint/python/trigScint.py | 36 ++ .../TrigScint/Firmware/clusterproducer_sw.cxx | 77 +++++ .../TrigScint/Firmware/trackproducer_hw.cxx | 159 +++++++++ .../TrigScint/TrigScintFirmwareTracker.cxx | 308 ++++++++++++++++++ 12 files changed, 1217 insertions(+), 1 deletion(-) create mode 100644 TrigScint/exampleConfigs/firmwareEx.py create mode 100644 TrigScint/exampleConfigs/meganEx.py create mode 100644 TrigScint/include/TrigScint/Firmware/clusterproducer.h create mode 100755 TrigScint/include/TrigScint/Firmware/objdef.h create mode 100755 TrigScint/include/TrigScint/Firmware/testutils.h create mode 100755 TrigScint/include/TrigScint/Firmware/trackproducer.h create mode 100644 TrigScint/include/TrigScint/TrigScintFirmwareTracker.h create mode 100644 TrigScint/src/TrigScint/Firmware/clusterproducer_sw.cxx create mode 100755 TrigScint/src/TrigScint/Firmware/trackproducer_hw.cxx create mode 100644 TrigScint/src/TrigScint/TrigScintFirmwareTracker.cxx diff --git a/TrigScint/CMakeLists.txt b/TrigScint/CMakeLists.txt index 2d1a84735..6344dccc2 100644 --- a/TrigScint/CMakeLists.txt +++ b/TrigScint/CMakeLists.txt @@ -40,10 +40,12 @@ if(BUILD_EVENT_ONLY) endif() +setup_library(module TrigScint name Firmware) +target_include_directories(TrigScint_Firmware PUBLIC ../Trigger/HLS_arbitrary_Precision_Types/include) setup_library(module TrigScint dependencies Framework::Framework Recon::Event DetDescr::DetDescr Tools::Tools SimCore::Event + Tools::Tools SimCore::Event TrigScint::Firmware ) - setup_python(package_name LDMX/TrigScint) diff --git a/TrigScint/exampleConfigs/firmwareEx.py b/TrigScint/exampleConfigs/firmwareEx.py new file mode 100644 index 000000000..fc9d4a93d --- /dev/null +++ b/TrigScint/exampleConfigs/firmwareEx.py @@ -0,0 +1,168 @@ +#!/bin/python + +import sys +import os +import json + +# we need the ldmx configuration package to construct the object + +from LDMX.Framework import ldmxcfg + +# set a 'pass name' +passName="sim" +p=ldmxcfg.Process(passName) + +#import all processors +from LDMX.SimCore import generators +from LDMX.SimCore import simulator +from LDMX.Biasing import filters + +from LDMX.Detectors.makePath import * +from LDMX.SimCore import simcfg + +#pull in command line options +nEle=4 # simulated beam electrons +runNum=10 +version="ldmx-det-v14" +outputNameString= "ldmxdetv14gap10mm_firmware.root" #sample identifier +outDir= "" #sample identifier + +# +# Instantiate the simulator. +# +sim = simulator.simulator("test") + +# +# Set the path to the detector to use (pulled from job config) +# +sim.setDetector( version, True ) +sim.scoringPlanes = makeScoringPlanesPath(version) + +outname=outputNameString #+".root" +print("NAME = " + outname) + +# +# Set run parameters. These are all pulled from the job config +# +p.run = runNum +p.maxEvents = 100 +nElectrons = nEle +beamEnergy = 4.0; #in GeV + +sim.description = "Inclusive "+str(beamEnergy)+" GeV electron events, "+str(nElectrons)+"e" +#sim.randomSeeds = [ SEED1 , SEED2 ] +sim.beamSpotSmear = [20., 80., 0] + + +mpgGen = generators.multi( "mgpGen" ) # this is the line that actually creates the generator +mpgGen.vertex = [ -44., 0., -880. ] # mm +mpgGen.nParticles = nElectrons +mpgGen.pdgID = 11 +mpgGen.enablePoisson = False #True + +import math +theta = math.radians(5.45) +beamEnergyMeV=1000*beamEnergy +px = beamEnergyMeV*math.sin(theta) +py = 0.; +pz= beamEnergyMeV*math.cos(theta) +mpgGen.momentum = [ px, py, pz ] + +# +# Set the multiparticle gun as generator +# +sim.generators = [ mpgGen ] + +#reconstruction and vetoes + +#Ecal and Hcal hardwired/geometry stuff +#import LDMX.Ecal.EcalGeometry +import LDMX.Ecal.ecal_hardcoded_conditions +from LDMX.Ecal import EcalGeometry +#egeom = EcalGeometry.EcalGeometryProvider.getInstance() +#Hcal hardwired/geometry stuff +from LDMX.Hcal import HcalGeometry +import LDMX.Hcal.hcal_hardcoded_conditions +#hgeom = HcalGeometry.HcalGeometryProvider.getInstance() + + +from LDMX.Ecal import digi as eDigi +from LDMX.Ecal import vetos +from LDMX.Hcal import digi as hDigi +from LDMX.Hcal import hcal + +from LDMX.Recon.simpleTrigger import TriggerProcessor + +from LDMX.TrigScint.trigScint import TrigScintDigiProducer +from LDMX.TrigScint.trigScint import TrigScintClusterProducer +from LDMX.TrigScint.trigScint import trigScintTrack +from LDMX.TrigScint.trigScint import TrigScintFirmwareTracker + +tsSimColls=[ "TriggerPad2SimHits", "TriggerPad3SimHits", "TriggerPad1SimHits" ] + +# ecal digi chain +# ecalDigi =eDigi.EcalDigiProducer('EcalDigis') +# ecalReco =eDigi.EcalRecProducer('ecalRecon') +# ecalVeto =vetos.EcalVetoProcessor('ecalVetoBDT') + +# #hcal digi chain +# hcalDigi =hDigi.HcalDigiProducer('hcalDigis') +# hcalReco =hDigi.HcalRecProducer('hcalRecon') +# hcalVeto =hcal.HcalVetoProcessor('hcalVeto') +# #hcalDigi.inputCollName="HcalSimHits" +#hcalDigi.inputPassName=passName + +# TS digi + clustering + track chain +tsDigisTag =TrigScintDigiProducer.pad2() +tsDigisTag.input_collection = tsSimColls[0]# +"_"+passName +tsDigisTag.input_pass_name = "sim" +tsDigisUp =TrigScintDigiProducer.pad3() +tsDigisUp.input_collection = tsSimColls[1]# +"_"+passName +tsDigisUp.input_pass_name = "sim" +tsDigisDown=TrigScintDigiProducer.pad1() +tsDigisDown.input_collection = tsSimColls[2]# +"_"+passName +tsDigisDown.input_pass_name = "sim" + +tsClustersTag =TrigScintClusterProducer.pad2() +tsClustersUp =TrigScintClusterProducer.pad1() +tsClustersDown =TrigScintClusterProducer.pad3() + + +tsDigisUp.verbosity=0 +tsClustersUp.verbosity=1 +trigScintTrack.verbosity=1 + +trigScintTrack.delta_max = 0.75 + +trigFirm = TrigScintFirmwareTracker( "trigFirm" ) +trigFirm.input_pass_name = "sim" +trigFirm.digis1_collection = "trigScintDigisPad1" +trigFirm.digis2_collection = "trigScintDigisPad2" +trigFirm.digis3_collection = "trigScintDigisPad3" +trigFirm.output_collection = "TriggerPadTracksFirmware" + +from LDMX.Recon.electronCounter import ElectronCounter +eCount = ElectronCounter( nElectrons, "ElectronCounter") # first argument is number of electrons in simulation +eCount.use_simulated_electron_number = False +eCount.input_collection="TriggerPadTracks" +eCount.input_pass_name=passName + +# # p.sequence=[ sim, ecalDigi, ecalReco, ecalVeto, hcalDigi, hcalReco, hcalVeto, tsDigisTag, tsDigisUp, tsDigisDown, tsClustersTag, tsClustersUp, tsClustersDown, trigScintTrack, eCount ] +# #hcal digi keeps crashing in config step +p.sequence=[ sim, tsDigisTag, tsDigisUp, tsDigisDown, tsClustersTag, tsClustersUp, tsClustersDown, trigScintTrack, trigFirm, eCount] +# p.sequence=[sim] + +p.outputFiles=[outname] + +p.termLogLevel = 0 # default is 2 (WARNING); but then logFrequency is ignored. level 1 = INFO. + +#print this many events to stdout (independent on number of events, edge case: round-off effects when not divisible. so can go up by a factor 2 or so) +logEvents=20 +if p.maxEvents < logEvents : + logEvents = p.maxEvents +p.logFrequency = int( p.maxEvents/logEvents ) + +#json.dumps(p.parameterDump(), indent=2) + +with open('parameterDump.json', 'w') as outfile: + json.dump(p.parameterDump(), outfile, indent=4) diff --git a/TrigScint/exampleConfigs/meganEx.py b/TrigScint/exampleConfigs/meganEx.py new file mode 100644 index 000000000..e7e4bc203 --- /dev/null +++ b/TrigScint/exampleConfigs/meganEx.py @@ -0,0 +1,167 @@ +#!/bin/python + +import sys +import os +import json + +# we need the ldmx configuration package to construct the object + +from LDMX.Framework import ldmxcfg + +# set a 'pass name' +passName="sim" +p=ldmxcfg.Process(passName) + +#import all processors +from LDMX.SimCore import generators +from LDMX.SimCore import simulator +from LDMX.Biasing import filters + +from LDMX.Detectors.makePath import * +from LDMX.SimCore import simcfg + +#pull in command line options +nEle=4 # simulated beam electrons +runNum=10 +version="ldmx-det-v14" +outputNameString= "ldmxdetv14gap10mm.root" #sample identifier +outDir= "" #sample identifier + +# +# Instantiate the simulator. +# +sim = simulator.simulator("test") + +# +# Set the path to the detector to use (pulled from job config) +# +sim.setDetector( version, True ) +sim.scoringPlanes = makeScoringPlanesPath(version) + +outname=outputNameString #+".root" +print("NAME = " + outname) + +# +# Set run parameters. These are all pulled from the job config +# +p.run = runNum +p.maxEvents = 100 +nElectrons = nEle +beamEnergy = 4.0; #in GeV + +sim.description = "Inclusive "+str(beamEnergy)+" GeV electron events, "+str(nElectrons)+"e" +#sim.randomSeeds = [ SEED1 , SEED2 ] +sim.beamSpotSmear = [20., 80., 0] + + +mpgGen = generators.multi( "mgpGen" ) # this is the line that actually creates the generator +mpgGen.vertex = [ -44., 0., -880. ] # mm +mpgGen.nParticles = nElectrons +mpgGen.pdgID = 11 +mpgGen.enablePoisson = False #True + +import math +theta = math.radians(5.45) +beamEnergyMeV=1000*beamEnergy +px = beamEnergyMeV*math.sin(theta) +py = 0.; +pz= beamEnergyMeV*math.cos(theta) +mpgGen.momentum = [ px, py, pz ] + +# +# Set the multiparticle gun as generator +# +sim.generators = [ mpgGen ] + +#reconstruction and vetoes + +#Ecal and Hcal hardwired/geometry stuff +#import LDMX.Ecal.EcalGeometry +import LDMX.Ecal.ecal_hardcoded_conditions +from LDMX.Ecal import EcalGeometry +#egeom = EcalGeometry.EcalGeometryProvider.getInstance() +#Hcal hardwired/geometry stuff +from LDMX.Hcal import HcalGeometry +import LDMX.Hcal.hcal_hardcoded_conditions +#hgeom = HcalGeometry.HcalGeometryProvider.getInstance() + + +from LDMX.Ecal import digi as eDigi +from LDMX.Ecal import vetos +from LDMX.Hcal import digi as hDigi +from LDMX.Hcal import hcal + +from LDMX.Recon.simpleTrigger import TriggerProcessor + +from LDMX.TrigScint.trigScint import TrigScintDigiProducer +from LDMX.TrigScint.trigScint import TrigScintClusterProducer +from LDMX.TrigScint.trigScint import trigScintTrack + +if "v12" in version : + tsSimColls=[ "TriggerPadTagSimHits", "TriggerPadUpSimHits", "TriggerPadDnSimHits" ] +else : + tsSimColls=[ "TriggerPad2SimHits", "TriggerPad3SimHits", "TriggerPad1SimHits" ] + +# ecal digi chain +# ecalDigi =eDigi.EcalDigiProducer('EcalDigis') +# ecalReco =eDigi.EcalRecProducer('ecalRecon') +# ecalVeto =vetos.EcalVetoProcessor('ecalVetoBDT') + +# #hcal digi chain +# hcalDigi =hDigi.HcalDigiProducer('hcalDigis') +# hcalReco =hDigi.HcalRecProducer('hcalRecon') +# hcalVeto =hcal.HcalVetoProcessor('hcalVeto') +# #hcalDigi.inputCollName="HcalSimHits" +#hcalDigi.inputPassName=passName + +# TS digi + clustering + track chain +tsDigisTag =TrigScintDigiProducer.pad2() +tsDigisTag.input_collection = tsSimColls[0]# +"_"+passName +tsDigisTag.input_pass_name = "sim" +tsDigisUp =TrigScintDigiProducer.pad3() +tsDigisUp.input_collection = tsSimColls[1]# +"_"+passName +tsDigisUp.input_pass_name = "sim" +tsDigisDown=TrigScintDigiProducer.pad1() +tsDigisDown.input_collection = tsSimColls[2]# +"_"+passName +tsDigisDown.input_pass_name = "sim" + +tsClustersTag =TrigScintClusterProducer.pad2() +tsClustersUp =TrigScintClusterProducer.pad1() +tsClustersDown =TrigScintClusterProducer.pad3() + +if "v12" in version : + tsClustersTag.pad_time = -2. + tsClustersUp.pad_time = 0. + tsClustersDown.pad_time = 0. + +tsDigisUp.verbosity=0 +tsClustersUp.verbosity=1 +trigScintTrack.verbosity=1 + +trigScintTrack.delta_max = 0.75 + +from LDMX.Recon.electronCounter import ElectronCounter +eCount = ElectronCounter( nElectrons, "ElectronCounter") # first argument is number of electrons in simulation +eCount.use_simulated_electron_number = False +eCount.input_collection="TriggerPadTracks" +eCount.input_pass_name=passName + +# # p.sequence=[ sim, ecalDigi, ecalReco, ecalVeto, hcalDigi, hcalReco, hcalVeto, tsDigisTag, tsDigisUp, tsDigisDown, tsClustersTag, tsClustersUp, tsClustersDown, trigScintTrack, eCount ] +# #hcal digi keeps crashing in config step +p.sequence=[ sim, tsDigisTag, tsDigisUp, tsDigisDown, tsClustersTag, tsClustersUp, tsClustersDown, trigScintTrack, eCount] +# p.sequence=[sim] + +p.outputFiles=[outname] + +p.termLogLevel = 0 # default is 2 (WARNING); but then logFrequency is ignored. level 1 = INFO. + +#print this many events to stdout (independent on number of events, edge case: round-off effects when not divisible. so can go up by a factor 2 or so) +logEvents=20 +if p.maxEvents < logEvents : + logEvents = p.maxEvents +p.logFrequency = int( p.maxEvents/logEvents ) + +json.dumps(p.parameterDump(), indent=2) + +with open('parameterDump.json', 'w') as outfile: + json.dump(p.parameterDump(), outfile, indent=4) diff --git a/TrigScint/include/TrigScint/Firmware/clusterproducer.h b/TrigScint/include/TrigScint/Firmware/clusterproducer.h new file mode 100644 index 000000000..94fd1ccd6 --- /dev/null +++ b/TrigScint/include/TrigScint/Firmware/clusterproducer.h @@ -0,0 +1,12 @@ +#ifndef CLUSTERPRODUCER_H +#define CLUSTERPRODUCER_H + +#include "objdef.h" + +void copyHit1(Hit One, Hit Two); +void copyHit2(Hit One, Hit Two); +void clusterproducer_ref(Hit inHit[NHITS], Cluster outClus[NCLUS]); +std::array clusterproducer_sw(Hit inHit[NHITS]); +void clusterproducer_hw(Hit inHit[NHITS], Cluster outClus[NCLUS]); + +#endif diff --git a/TrigScint/include/TrigScint/Firmware/objdef.h b/TrigScint/include/TrigScint/Firmware/objdef.h new file mode 100755 index 000000000..39c80cda4 --- /dev/null +++ b/TrigScint/include/TrigScint/Firmware/objdef.h @@ -0,0 +1,134 @@ +#ifndef OBJDEF_H +#define OBJDEF_H + +#include "ap_int.h" +#define NTIMES 6 +#define NHITS 25 +#define NCLUS 25 +#define NCHAN 50 +#define NTRK 10 +#define W = 10 + +#define NCENT 99 + +#define NDIGIS 14 +#define COMBO 9 + +// 2*NCHAN*NTIMES are the number of bytes per event plus 4+4+4+3+1 bytes for the +// header + +#define NSAMPLES 6 + +// NSAMPLES/8 is the number of 64 bit words + +#define NWORDS 72 + +struct Digi { + int mID, bID; + int adc0, adc1, adc2, adc3, adc4, adc5; + int tdc0, tdc1, tdc2, tdc3, tdc4, tdc5; +}; +inline void clearDigi(Digi& c) { + c.mID = 0; + c.bID = 0; + c.adc0 = 0; + c.adc1 = 0; + c.adc2 = 0; + c.adc3 = 0; + c.adc4 = 0; + c.adc5 = 0; + c.tdc0 = 0; + c.tdc1 = 0; + c.tdc2 = 0; + c.tdc3 = 0; + c.tdc4 = 0; + c.tdc5 = 0; +} +struct Hit { + ap_int<12> mID, bID; + ap_int<12> Amp, Time; // TrigTime; +}; +inline void clearHit(Hit& c) { + c.mID = 0; + c.bID = -1; + c.Amp = 0; + c.Time = 0; // c.TrigTime=0.0; +} +inline void cpyHit(Hit& c1, Hit& c2) { + c1.mID = c2.mID; + c1.bID = c2.bID; + c1.Amp = c2.Amp; + c1.Time = c2.Time; +} + +struct Cluster { + Hit Seed; + Hit Sec; + ap_int<11> Cent; + // int nhits, mID, SeedID; + // float CentX, CentY, CentZ, Amp, Time, TrigTime; +}; +inline void clearClus(Cluster& c) { + clearHit(c.Seed); + clearHit(c.Sec); + c.Cent = (ap_int<11>)(0); // clearHit(c.For); +} +inline void calcCent(Cluster& c) { + if (c.Seed.Amp > 0) { + c.Cent = (ap_int<12>)(10. * + ((float)(c.Seed.Amp * c.Seed.bID + + c.Sec.Amp * c.Sec.bID)) / + ((float)(c.Seed.Amp + c.Sec.Amp))); + } else { + c.Cent = (ap_int<12>)(0); + } +} +inline void cpyCluster(Cluster& c1, Cluster& c2) { + cpyHit(c1.Seed, c2.Seed); + cpyHit(c1.Sec, c2.Sec); +} + +struct Track { + Cluster Pad1; + Cluster Pad2; + Cluster Pad3; + ap_int<12> resid; +}; +inline void clearTrack(Track& c) { + clearClus(c.Pad1); + clearClus(c.Pad2); + clearClus(c.Pad3); + c.resid = 5000; +} +inline ap_int<12> calcTCent(Track& c) { + calcCent(c.Pad1); + calcCent(c.Pad2); + calcCent(c.Pad3); + float one = (float)c.Pad1.Cent; + float two = (float)c.Pad2.Cent; + float three = (float)c.Pad3.Cent; + float mean = (one + two + three) / 3.0; + ap_int<12> Cent = (ap_int<10>)((int)(mean)); + return Cent; +} +inline void calcResid(Track& c) { + calcCent(c.Pad1); + calcCent(c.Pad2); + calcCent(c.Pad3); + float one = (float)c.Pad1.Cent; + float two = (float)c.Pad2.Cent; + float three = (float)c.Pad3.Cent; + float mean = (one + two + three) / 3.0; + c.resid = (ap_int<12>)((int)(((one - mean) * (one - mean) + + (two - mean) * (two - mean) + + (three - mean) * (three - mean)) / + 3.0)); +} +inline void cpyTrack(Track& c1, Track& c2) { + cpyCluster(c1.Pad1, c2.Pad1); + cpyCluster(c1.Pad2, c2.Pad2); + cpyCluster(c1.Pad3, c2.Pad3); + c1.resid = c2.resid; +} + +#endif diff --git a/TrigScint/include/TrigScint/Firmware/testutils.h b/TrigScint/include/TrigScint/Firmware/testutils.h new file mode 100755 index 000000000..b6b1e2240 --- /dev/null +++ b/TrigScint/include/TrigScint/Firmware/testutils.h @@ -0,0 +1,21 @@ +#ifndef TESTUTILS_H +#define TESTUTILS_H +#include "objdef.h" + +bool compareHit(Hit Hit1, Hit Hit2) { + return ((Hit1.mID == Hit2.mID) and (Hit1.bID == Hit2.bID) and + (Hit1.Amp == Hit2.Amp) and + (Hit1.Time == Hit2.Time)); // and(Hit1.TrigTime==Hit2.TrigTime)); +} + +bool compareClus(Cluster clus1[NHITS], Cluster clus2[NHITS]) { + for (int i = 0; i < NHITS; ++i) { + if (not((compareHit(clus1[i].Seed, clus2[i].Seed)) and + (compareHit(clus1[i].Sec, clus2[i].Sec)))) { + return false; + } + } + return true; +} + +#endif diff --git a/TrigScint/include/TrigScint/Firmware/trackproducer.h b/TrigScint/include/TrigScint/Firmware/trackproducer.h new file mode 100755 index 000000000..c3d2f3cbc --- /dev/null +++ b/TrigScint/include/TrigScint/Firmware/trackproducer.h @@ -0,0 +1,15 @@ +#ifndef TRACKPRODUCER_H +#define TRACKPRODUCER_H + +#include "objdef.h" + +void copyCluster1(Cluster One, Cluster Two); +void copyCluster2(Cluster One, Cluster Two); +void trackproducer_ref(Cluster Pad1[NTRK], Cluster Pad2[NCLUS], + Cluster Pad3[NCLUS], Track outTrk[NTRK], + ap_int<12> lookup[NCENT][COMBO][2]); +void trackproducer_hw(Cluster Pad1[NTRK], Cluster Pad2[NCLUS], + Cluster Pad3[NCLUS], Track outTrk[NTRK], + ap_int<12> lookup[NCENT][COMBO][2]); + +#endif diff --git a/TrigScint/include/TrigScint/TrigScintFirmwareTracker.h b/TrigScint/include/TrigScint/TrigScintFirmwareTracker.h new file mode 100644 index 000000000..974ccc8ad --- /dev/null +++ b/TrigScint/include/TrigScint/TrigScintFirmwareTracker.h @@ -0,0 +1,117 @@ +/** + * @file TrigScintFirmwareTracker.h + * @brief Tracker made to emulate and stage real firmware, emulates existing + * ldmx software but has LUT structure. + * @author Rory O'Dwyer, Stanford University + */ + +#ifndef TRIGSCINT_TRIGSCINTFIRMWARETRACKER_H +#define TRIGSCINT_TRIGSCINTFIRMWARETRACKER_H + +// LDMX Framework +#include "Framework/Configure/Parameters.h" // Needed to import parameters from configuration file +#include "Framework/Event.h" +#include "Framework/EventProcessor.h" //Needed to declare processor +#include "Recon/Event/EventConstants.h" +#include "TrigScint/Event/TrigScintHit.h" +#include "TrigScint/Event/TrigScintTrack.h" +#include "TrigScint/Firmware/objdef.h" +#include "TrigScint/TrigScintFirmwareTracker.h" + +namespace trigscint { + +/** + * @class TrigScintFirmwareTracker + * @brief + */ +class TrigScintFirmwareTracker : public framework::Producer { + public: + TrigScintFirmwareTracker(const std::string& name, framework::Process& process) + : Producer(name, process) {} + + void configure(framework::config::Parameters& ps) override; + + void produce(framework::Event& event) override; + + ldmx::TrigScintTrack makeTrack(Track outTrk); + + /** + * add a hit at index idx to a cluster + */ + + private: + // collection of clusters produced + std::vector digis1_; + + // collection of clusters produced + std::vector digis2_; + + // collection of clusters produced + std::vector digis3_; + + // min threshold for adding a hit to a cluster + double minThr_{0.}; + + // max number of neighboring hits to combine when forming a cluster + int maxWidth_{2}; + + // specific verbosity of this producer + int verbose_{0}; + + // expected arrival time of hits in the pad [ns] + double padTime_{0.}; + + // maximum allowed delay for hits to be considered for clustering + double timeTolerance_{0.}; + + // output collection (clusters) + std::string output_collection_; + + // input collection (hits) + std::string digis1_collection_; + std::string digis2_collection_; + std::string digis3_collection_; + + std::vector tracks_; + + // specific pass name to use for track making + std::string passName_{""}; + + // vertical bar start index + int vertBarStartIdx_{52}; + + // cluster channel nb centroid (will be content weighted) + float centroid_{0.}; + + // cluster channel nb horizontal centroid (will be content weighted) + float centroidX_{-1}; + + // cluster channel nb vertical centroid (will be content weighted) + float centroidY_{-1}; + + // energy (edep), PE, or sth + float val_{0.}; + + // edep content, only; leave val_ for PE + float valE_{0.}; + + // book keep which channels have already been added to the cluster at hand + std::vector v_addedIndices_; + + // book keep which channels have already been added to any cluster + std::vector v_usedIndices_; + + // fraction of cluster energy deposition associated with beam electron sim + // hits + float beamE_{0.}; + + // cluster time (energy weighted based on hit time) + float time_{0.}; + + // empty map container + std::map hitChannelMap_; +}; + +} // namespace trigscint + +#endif /* TRIGSCINT_TRIGSCINTFIRMWARETRACKER_H */ diff --git a/TrigScint/python/trigScint.py b/TrigScint/python/trigScint.py index 7d8ccb312..f8fcbec77 100644 --- a/TrigScint/python/trigScint.py +++ b/TrigScint/python/trigScint.py @@ -285,6 +285,42 @@ def __init__(self,name) : trigScintTrack = TrigScintTrackProducer( "trigScintTrack" ) +class TrigScintTrackProducer(ldmxcfg.Producer) : + """Configuration for track producer for Trigger Scintillators""" + + def __init__(self,name) : + super().__init__(name,'trigscint::','TrigScint') + + self.delta_max = 0.75 + self.tracking_threshold = 0. #to add in neighboring channels + self.seeding_collection = "TriggerPad1Clusters" + self.further_input_collections = ["TriggerPad2Clusters","TriggerPad3Clusters"] + self.allow_skip_last_collection = False + self.vertical_bar_start_index = 52 + self.number_horizontal_bars = 24 #16 for x,y segmented geometry only + self.number_vertical_bars = 0 #8 for x,y segmented geometry only + self.horizontal_bar_width = 3. + self.horizontal_bar_gap = 0.3 + self.vertical_bar_width = 3. + self.vertical_bar_gap = 0.3 + self.input_pass_name="" #take any pass + self.output_collection="TriggerPadTracks" + self.verbosity = 0 + +class TrigScintFirmwareTracker(ldmxcfg.Producer) : + """Configuration for the track producer from the Firmware Tracker""" + def __init__(self,name) : + super().__init__(name,'trigscint::TrigScintFirmwareTracker','TrigScint') + self.clustering_threshold=40.0 + self.digis1_collection='trigScintDigisPad1' + self.digis2_collection='trigScintDigisPad2' + self.digis3_collection='trigScintDigisPad3' + self.input_pass_name="" + self.output_collection="TriggerPadTracks" + self.verbosity = 0 + self.time_tolerance = 50.0 + self.pad_time = -1.5 + class QIEAnalyzer(ldmxcfg.Analyzer) : """Configuration for linearized QIE analyzer for Trigger Scintillators""" diff --git a/TrigScint/src/TrigScint/Firmware/clusterproducer_sw.cxx b/TrigScint/src/TrigScint/Firmware/clusterproducer_sw.cxx new file mode 100644 index 000000000..140011d93 --- /dev/null +++ b/TrigScint/src/TrigScint/Firmware/clusterproducer_sw.cxx @@ -0,0 +1,77 @@ +#include + +#include + +#include "TrigScint/Firmware/clusterproducer.h" +#include "TrigScint/Firmware/objdef.h" + +std::array clusterproducer_sw(Hit inHit[NHITS]) { + ap_int<12> SEEDTHR = 30; + ap_int<12> CLUSTHR = 30; + + ap_int<12> mapL1[NCHAN]; + + std::array outClus; + + for (int i = 0; i < NCLUS; ++i) { + clearClus(outClus[i]); + } + + // CLEAR THE MAP + for (int i = 0; i < NCHAN; ++i) { + mapL1[i] = -1; + } + // MAP TO CHANNELS + for (int j = 0; j < NHITS; ++j) { + if (inHit[j].bID > -1) { + mapL1[inHit[j].bID] = j; + } + } + // NOW WE JUST LOOK FOR HITS EXCEEDING SEED, IF THEY DO WE PAIR 'EM. + for (int k = 0; k < NCLUS; ++k) { + bool doNextCluster = true; + if ((mapL1[2 * k] > -1)) { + if (inHit[mapL1[2 * k]].Amp > SEEDTHR) { + clearClus(outClus[k]); + outClus[k].Seed.mID = inHit[mapL1[2 * k]].mID; + outClus[k].Seed.bID = inHit[mapL1[2 * k]].bID; + outClus[k].Seed.Amp = inHit[mapL1[2 * k]].Amp; + outClus[k].Seed.Time = inHit[mapL1[2 * k]].Time; + if (mapL1[2 * k + 1] > -1) { + if (inHit[mapL1[2 * k + 1]].Amp > CLUSTHR) { + outClus[k].Sec.mID = inHit[mapL1[2 * k + 1]].mID; + outClus[k].Sec.bID = inHit[mapL1[2 * k + 1]].bID; + outClus[k].Sec.Amp = inHit[mapL1[2 * k + 1]].Amp; + outClus[k].Sec.Time = inHit[mapL1[2 * k + 1]].Time; + doNextCluster = false; + // You can comment this line to turn it into Serialized + clearHit(inHit[mapL1[2 * k + 1]]); + } + } + } + } + if ((mapL1[2 * k + 1] > -1) and (doNextCluster)) { + if (inHit[mapL1[2 * k + 1]].Amp > SEEDTHR) { + clearClus(outClus[k]); + outClus[k].Seed.mID = inHit[mapL1[2 * k + 1]].mID; + outClus[k].Seed.bID = inHit[mapL1[2 * k + 1]].bID; + outClus[k].Seed.Amp = inHit[mapL1[2 * k + 1]].Amp; + outClus[k].Seed.Time = inHit[mapL1[2 * k + 1]].Time; + if (k < NCLUS - 1) { + if (mapL1[2 * k + 2] > -1) { + if (inHit[mapL1[2 * k + 2]].Amp > CLUSTHR) { + outClus[k].Sec.mID = inHit[mapL1[2 * k + 2]].mID; + outClus[k].Sec.bID = inHit[mapL1[2 * k + 2]].bID; + outClus[k].Sec.Amp = inHit[mapL1[2 * k + 2]].Amp; + outClus[k].Sec.Time = inHit[mapL1[2 * k + 2]].Time; + // You can comment this line to turn it into Serialized + clearHit(inHit[mapL1[2 * k + 2]]); + } + } + } + } + } + } + + return outClus; +} diff --git a/TrigScint/src/TrigScint/Firmware/trackproducer_hw.cxx b/TrigScint/src/TrigScint/Firmware/trackproducer_hw.cxx new file mode 100755 index 000000000..b894a27cf --- /dev/null +++ b/TrigScint/src/TrigScint/Firmware/trackproducer_hw.cxx @@ -0,0 +1,159 @@ +#include + +#include + +#include "TrigScint/Firmware/objdef.h" +#include "TrigScint/Firmware/trackproducer.h" + +void trackproducer_hw(Cluster Pad1[NTRK], Cluster Pad2[NCLUS], + Cluster Pad3[NCLUS], Track outTrk[NTRK], + ap_int<12> lookup[NCENT][COMBO][2]) { +#pragma HLS ARRAY_PARTITION variable = Pad1 dim = 0 complete +#pragma HLS ARRAY_PARTITION variable = Pad2 dim = 0 complete +#pragma HLS ARRAY_PARTITION variable = Pad3 dim = 0 complete +#pragma HLS ARRAY_PARTITION variable = outTrk dim = 0 complete +#pragma HLS ARRAY_PARTITION variable = lookup dim = 0 complete +#pragma HLS PIPELINE II = 10 + Track test; +#pragma HLS ARRAY_PARTITION variable = test complete + + // This firmware module loops over first the Pad1 seeds (NTRK) and then the + // patterns (COMBO) For each seed it check 9 combinations of tracks. These + // combinations, which depend on alignment essentially consist of the clusters + // that have channels immediattely above or below the Pad1 cluster in the + // first layer, which you may observe from the LUT if you printed it. I would + // only need to check the pattern without all these continue statements, but + // the continue statements further reduce the pattern collection size by only + // applying certain patterns iff a secondary hit is there Thats why this looks + // complicated at all: the continues just include logic on whether a pattern + // should have a secondary hit. It also checks the track residual, only + // keeping one pattern for each pad1 cluster. + for (int i = 0; i < NTRK; i++) { + if (2 * Pad1[i].Seed.bID > 2 * NCHAN) { + continue; + } + for (int I = 0; I < COMBO; I++) { + clearTrack(test); + if (not(Pad1[i].Seed.Amp > 0)) { + continue; + } // Continue if Seed not Satisfied + ap_int<12> centroid = 2 * Pad1[i].Seed.bID; + if (Pad1[i].Sec.Amp > 0) { + centroid += 1; + } + cpyCluster(test.Pad1, Pad1[i]); + if ((lookup[centroid][I][0] == -1) or (lookup[centroid][I][1] == -1)) { + continue; + } // Pattern Empty + if (not(Pad2[lookup[centroid][I][0] / 4].Seed.Amp > 0)) { + continue; + } // Continue if Seed not Satisfied + if ((lookup[centroid][I][0] % 4 == 0) and + ((Pad2[lookup[centroid][I][0] / 4].Sec.bID >= 0) or + (Pad2[lookup[centroid][I][0] / 4].Seed.bID % 2 == 1))) { + continue; + } // Continue if Sec is not Expected, and not Empty + if ((lookup[centroid][I][0] % 4 == 1) and + ((Pad2[lookup[centroid][I][0] / 4].Sec.bID < 0) or + (Pad2[lookup[centroid][I][0] / 4].Seed.bID % 2 == 1))) { + continue; + } // Continue if Sec is Expected, and Empty + if ((lookup[centroid][I][0] % 4 == 2) and + ((Pad2[lookup[centroid][I][0] / 4].Sec.bID >= 0) or + (Pad2[lookup[centroid][I][0] / 4].Seed.bID % 2 == 0))) { + continue; + } // Continue if Sec is not Expected, and not Empty + if ((lookup[centroid][I][0] % 4 == 3) and + ((Pad2[lookup[centroid][I][0] / 4].Sec.bID < 0) or + (Pad2[lookup[centroid][I][0] / 4].Seed.bID % 2 == 0))) { + continue; + } // Continue if Sec is Expected, and Empty + if (not(Pad3[lookup[centroid][I][1] / 4].Seed.Amp > 0)) { + continue; + } // Continue if Seed not Satisfied + if ((lookup[centroid][I][1] % 4 == 0) and + ((Pad3[lookup[centroid][I][1] / 4].Sec.bID >= 0) or + (Pad3[lookup[centroid][I][1] / 4].Seed.bID % 2 == 1))) { + continue; + } // Continue if Sec is not Expected, and not Empty + if ((lookup[centroid][I][1] % 4 == 1) and + ((Pad3[lookup[centroid][I][1] / 4].Sec.bID < 0) or + (Pad3[lookup[centroid][I][1] / 4].Seed.bID % 2 == 1))) { + continue; + } // Continue if Sec is Expected, and Empty + if ((lookup[centroid][I][1] % 4 == 2) and + ((Pad3[lookup[centroid][I][1] / 4].Sec.bID >= 0) or + (Pad3[lookup[centroid][I][1] / 4].Seed.bID % 2 == 0))) { + continue; + } // Continue if Sec is not Expected, and not Empty + if ((lookup[centroid][I][1] % 4 == 3) and + ((Pad3[lookup[centroid][I][1] / 4].Sec.bID < 0) or + (Pad3[lookup[centroid][I][1] / 4].Seed.bID % 2 == 0))) { + continue; + } // Continue if Sec is Expected, and Empty + cpyCluster(test.Pad2, Pad2[lookup[centroid][I][0] / 4]); + cpyCluster(test.Pad3, Pad3[lookup[centroid][I][1] / 4]); + calcResid(test); + if (test.resid < outTrk[i].resid) { + cpyTrack(outTrk[i], test); + } + } + } + // While we ultimately envision having the firmware do duplicate track removal + // in the other two layers in a separate firmware module, they are done here + // so as to not have track over counting and to validate the processor. Thats + // what occurs here below. + for (int i = 1; i < NTRK - 1; i++) { + if ((outTrk[i - 1].Pad2.Seed.bID == outTrk[i].Pad2.Seed.bID) and + (outTrk[i].Pad2.Seed.bID >= 0)) { + if (outTrk[i - 1].resid <= outTrk[i].resid) { + clearTrack(outTrk[i]); + } else { + clearTrack(outTrk[i - 1]); + } + } + if ((outTrk[i].Pad2.Seed.bID == outTrk[i + 1].Pad2.Seed.bID) and + (outTrk[i + 1].Pad2.Seed.bID >= 0)) { + if (outTrk[i + 1].resid <= outTrk[i].resid) { + clearTrack(outTrk[i]); + } else { + clearTrack(outTrk[i + 1]); + } + } + if ((outTrk[i - 1].Pad2.Seed.bID == outTrk[i + 1].Pad2.Seed.bID) and + (outTrk[i + 1].Pad2.Seed.bID >= 0)) { + if (outTrk[i - 1].resid <= outTrk[i + 1].resid) { + clearTrack(outTrk[i + 1]); + } else { + clearTrack(outTrk[i - 1]); + } + } + } + for (int i = 1; i < NTRK - 1; i++) { + if ((outTrk[i - 1].Pad3.Seed.bID == outTrk[i].Pad3.Seed.bID) and + (outTrk[i].Pad3.Seed.bID >= 0)) { + if (outTrk[i - 1].resid <= outTrk[i].resid) { + clearTrack(outTrk[i]); + } else { + clearTrack(outTrk[i - 1]); + } + } + if ((outTrk[i].Pad3.Seed.bID == outTrk[i + 1].Pad3.Seed.bID) and + (outTrk[i + 1].Pad3.Seed.bID >= 0)) { + if (outTrk[i + 1].resid <= outTrk[i].resid) { + clearTrack(outTrk[i]); + } else { + clearTrack(outTrk[i + 1]); + } + } + if ((outTrk[i - 1].Pad3.Seed.bID == outTrk[i + 1].Pad3.Seed.bID) and + (outTrk[i + 1].Pad3.Seed.bID >= 0)) { + if (outTrk[i - 1].resid <= outTrk[i + 1].resid) { + clearTrack(outTrk[i + 1]); + } else { + clearTrack(outTrk[i - 1]); + } + } + } + return; +} diff --git a/TrigScint/src/TrigScint/TrigScintFirmwareTracker.cxx b/TrigScint/src/TrigScint/TrigScintFirmwareTracker.cxx new file mode 100644 index 000000000..6af83a407 --- /dev/null +++ b/TrigScint/src/TrigScint/TrigScintFirmwareTracker.cxx @@ -0,0 +1,308 @@ + +#include "TrigScint/TrigScintFirmwareTracker.h" + +#include +#include + +#include "TrigScint/Firmware/clusterproducer.h" +#include "TrigScint/Firmware/objdef.h" +#include "TrigScint/Firmware/trackproducer.h" + +namespace trigscint { + +void TrigScintFirmwareTracker::configure(framework::config::Parameters &ps) { + minThr_ = ps.getParameter("clustering_threshold"); + digis1_collection_ = ps.getParameter("digis1_collection"); + digis2_collection_ = ps.getParameter("digis2_collection"); + digis3_collection_ = ps.getParameter("digis3_collection"); + passName_ = ps.getParameter("input_pass_name"); + output_collection_ = ps.getParameter("output_collection"); + verbose_ = ps.getParameter("verbosity"); + timeTolerance_ = ps.getParameter("time_tolerance"); + padTime_ = ps.getParameter("pad_time"); + if (verbose_) { + ldmx_log(info) << "In TrigScintFirmwareTracker: configure done!"; + ldmx_log(info) << "\nClustering threshold: " << minThr_ + << "\nExpected pad hit time: " << padTime_ + << "\nMax hit time delay: " << timeTolerance_ + << "\ndigis1 collection: " << digis1_collection_ + << "\ndigis2 collection: " << digis2_collection_ + << "\ndigis3 collection: " << digis3_collection_ + << "\nInput pass name: " << passName_ + << "\nOutput collection: " << output_collection_ + << "\nVerbosity: " << verbose_; + } + + return; +} + +void TrigScintFirmwareTracker::produce(framework::Event &event) { + // This processor takes in TS digis and outputs a track collection. It does so + // using clusterproducer_sw and trackproducer_hw, which are validated pieces + // of HLS code (though clusterproducer_sw has had its instances of pragmas + // excluded. I will comment on how clusterproducer and trackproducer work more + // thouroughly in them respectively, but generally the clusterproducer makes + // only two hit clusters (as ready that was all that was made from the + // original sw) and does so by making a digi map and running along channels + // numerically and pairing if possible. The trackproducer takes a LOOKUP array + // as a LUT and does track pattern mathcing. This depends on alignment through + // the A vector below. + + if (verbose_) { + ldmx_log(debug) + << "TrigScintFirmwareTracker: produce() starts! Event number: " + << event.getEventHeader().getEventNumber(); + } + ap_int<12> A[3] = {0, 0, 0}; + ap_int<12> LOOKUP[NCENT][COMBO][2]; + + // This line fills in the LOOKUP table used for patter matching latter. The + // array takes in as its first argument the centroid of a first pad cluster, + // then the next two take on which track pattern (of ~9) we are matching to + // and the last if we are matching to a cluster with two hits + for (int i = 0; i < NCENT; i++) { + for (int j = 0; j < COMBO; j++) { + LOOKUP[i][j][0] = (i - A[1] + A[0]); + LOOKUP[i][j][1] = (i - A[2] + A[0]); + if (j / 3 == 0) { + LOOKUP[i][j][0] -= 1; + } else if (j / 3 == 2) { + LOOKUP[i][j][0] += 1; + } + if (j % 3 == 0) { + LOOKUP[i][j][1] -= 1; + } else if (j % 3 == 2) { + LOOKUP[i][j][1] += 1; + } + if (not((LOOKUP[i][j][0] >= 0) and (LOOKUP[i][j][1] >= 0) and + (LOOKUP[i][j][0] < NCENT) and (LOOKUP[i][j][1] < NCENT))) { + LOOKUP[i][j][0] = -1; + LOOKUP[i][j][1] = -1; + } + } + } + // Here we instantiate arrays necessary to do the rest of it. + Hit HPad1[NHITS]; + Hit HPad2[NHITS]; + Hit HPad3[NHITS]; + + Cluster Pad1[NCLUS]; + Cluster Pad2[NCLUS]; + Cluster Pad3[NCLUS]; + Track outTrk[NTRK]; + + for (int j = 0; j < NHITS; j++) { + clearHit(HPad1[j]); + clearHit(HPad2[j]); + clearHit(HPad3[j]); + } + for (int j = 0; j < NCLUS; j++) { + if (j < NTRK) { + clearClus(Pad1[j]); + } + clearClus(Pad2[j]); + clearClus(Pad3[j]); + } + for (int j = 0; j < NTRK; j++) { + clearTrack(outTrk[j]); + } + // I am reading in the three digi collections + const auto digis1_{ + event.getCollection(digis1_collection_, passName_)}; + const auto digis3_{ + event.getCollection(digis2_collection_, passName_)}; + const auto digis2_{ + event.getCollection(digis3_collection_, passName_)}; + + if (verbose_) { + ldmx_log(debug) << "Got digi collection " << digis1_collection_ << "_" + << passName_ << " with " << digis1_.size() << " entries "; + } + + // The next collection of things fill in the firmware hit objects from reading + // in the digi collections the necessary information. The firmware hit objects + // only keep bID,mID,Time, and PE count. + int occupied[NCHAN]; + for (int i = 0; i < NCHAN; i++) { + occupied[i] = -1; + } + int count = 0; + for (const auto &digi : digis1_) { + if ((digi.getPE() > minThr_) and (digi.getBarID() <= NCHAN) and + (digi.getBarID() >= 0)) { + ap_int<12> bID = (ap_int<12>)(digi.getBarID()); + ap_int<12> Amp = (ap_int<12>)(digi.getPE()); + int index = count; + if (occupied[(int)digi.getBarID()] >= 0) { + if (HPad1[(int)occupied[(int)digi.getBarID()]].Amp < digi.getPE()) { + HPad1[(int)occupied[(int)digi.getBarID()]].bID = + (ap_int<12>)(digi.getBarID()); + HPad1[(int)occupied[(int)digi.getBarID()]].mID = + (ap_int<12>)(digi.getModuleID()); + HPad1[(int)occupied[(int)digi.getBarID()]].Amp = + (ap_int<12>)(digi.getPE()); + HPad1[(int)occupied[(int)digi.getBarID()]].Time = + (ap_int<12>)(digi.getTime()); + } + } else { + HPad1[count].bID = (ap_int<12>)(digi.getBarID()); + HPad1[count].mID = (ap_int<12>)(digi.getModuleID()); + HPad1[count].Amp = (ap_int<12>)(digi.getPE()); + HPad1[count].Time = (ap_int<12>)(digi.getTime()); + occupied[digi.getBarID()] = count; + count++; + } + } + } + + for (int i = 0; i < NCHAN; i++) { + occupied[i] = -1; + } + count = 0; + for (const auto &digi : digis2_) { + if ((digi.getPE() > minThr_) and (digi.getBarID() <= NCHAN) and + (digi.getBarID() >= 0)) { + ap_int<12> bID = (ap_int<12>)(digi.getBarID()); + ap_int<12> Amp = (ap_int<12>)(digi.getPE()); + int index = count; + if (occupied[(int)digi.getBarID()] >= 0) { + if (HPad2[(int)occupied[(int)digi.getBarID()]].Amp < digi.getPE()) { + HPad2[(int)occupied[(int)digi.getBarID()]].bID = + (ap_int<12>)(digi.getBarID()); + HPad2[(int)occupied[(int)digi.getBarID()]].mID = + (ap_int<12>)(digi.getModuleID()); + HPad2[(int)occupied[(int)digi.getBarID()]].Amp = + (ap_int<12>)(digi.getPE()); + HPad2[(int)occupied[(int)digi.getBarID()]].Time = + (ap_int<12>)(digi.getTime()); + } + } else { + HPad2[count].bID = (ap_int<12>)(digi.getBarID()); + HPad2[count].mID = (ap_int<12>)(digi.getModuleID()); + HPad2[count].Amp = (ap_int<12>)(digi.getPE()); + HPad2[count].Time = (ap_int<12>)(digi.getTime()); + occupied[digi.getBarID()] = count; + count++; + } + } + } + for (int i = 0; i < NCHAN; i++) { + occupied[i] = -1; + } + count = 0; + for (const auto &digi : digis3_) { + if ((digi.getPE() > minThr_) and (digi.getBarID() <= NCHAN) and + (digi.getBarID() >= 0)) { + ap_int<12> bID = (ap_int<12>)(digi.getBarID()); + ap_int<12> Amp = (ap_int<12>)(digi.getPE()); + int index = count; + if (occupied[(int)digi.getBarID()] >= 0) { + if (HPad3[(int)occupied[(int)digi.getBarID()]].Amp < digi.getPE()) { + HPad3[(int)occupied[(int)digi.getBarID()]].bID = + (ap_int<12>)(digi.getBarID()); + HPad3[(int)occupied[(int)digi.getBarID()]].mID = + (ap_int<12>)(digi.getModuleID()); + HPad3[(int)occupied[(int)digi.getBarID()]].Amp = + (ap_int<12>)(digi.getPE()); + HPad3[(int)occupied[(int)digi.getBarID()]].Time = + (ap_int<12>)(digi.getTime()); + } + } else { + HPad3[count].bID = (ap_int<12>)(digi.getBarID()); + HPad3[count].mID = (ap_int<12>)(digi.getModuleID()); + HPad3[count].Amp = (ap_int<12>)(digi.getPE()); + HPad3[count].Time = (ap_int<12>)(digi.getTime()); + occupied[digi.getBarID()] = count; + count++; + } + } + } + count = 0; + // These next lines here calls clusterproducer_sw(HPad1), which is just the + // validated firmware module. Since ap_* class is messy, I had to do some + // post-call cleanup before looping over the clusters and putting them into + // Point i which is feed into track producer + int counterN = 0; + std::array Point1 = clusterproducer_sw(HPad1); + int topSeed = 0; + for (int i = 0; i < NCLUS; i++) { + if ((Point1[i].Seed.Amp < 450) and (Point1[i].Seed.Amp > 30) and + (Point1[i].Seed.bID < (NCHAN + 1)) and (Point1[i].Seed.bID >= 0) and + (Point1[i].Sec.Amp < 450) and (counterN < NTRK)) { + if (Point1[i].Seed.bID >= topSeed) { + cpyHit(Pad1[counterN].Seed, Point1[i].Seed); + cpyHit(Pad1[counterN].Sec, Point1[i].Sec); + calcCent(Pad1[counterN]); + counterN++; + topSeed = Point1[i].Seed.bID; + } + } + } + std::array Point2 = clusterproducer_sw(HPad2); + topSeed = 0; + for (int i = 0; i < NCLUS; i++) { + if ((Point2[i].Seed.Amp < 450) and (Point2[i].Seed.Amp > 30) and + (Point2[i].Seed.bID < (NCHAN + 1)) and (Point2[i].Seed.bID >= 0) and + (Point2[i].Sec.Amp < 450)) { + if (Point2[i].Seed.bID >= topSeed) { + cpyHit(Pad2[i].Seed, Point2[i].Seed); + cpyHit(Pad2[i].Sec, Point2[i].Sec); + calcCent(Pad2[i]); + topSeed = Point2[i].Seed.bID; + } + } + } + std::array Point3 = clusterproducer_sw(HPad3); + topSeed = 0; + for (int i = 0; i < NCLUS; i++) { + if ((Point3[i].Seed.Amp < 450) and (Point3[i].Seed.Amp > 30) and + (Point3[i].Seed.bID < (NCHAN + 1)) and (Point3[i].Seed.bID >= 0) and + (Point3[i].Sec.Amp < 450)) { + if (Point3[i].Seed.bID >= topSeed) { + cpyHit(Pad3[i].Seed, Point3[i].Seed); + cpyHit(Pad3[i].Sec, Point3[i].Sec); + calcCent(Pad3[i]); + topSeed = Point3[i].Seed.bID; + } + } + } + // I have stagged the digis into firmware digi objects and paired them into + // firmware cluster objects, so at this point I can insert them and the LUT + // into the trackproducer_hw to create the track collection I use makeTrack to + // revert the firmware track object back into a regular track object for + // analysis purposes + // + // NOTE: Pad1 has NTRK instead of NCLUS clusters for a reason: the firmware + // cannot facilitate NCLUS many tracks within its alloted bandwidth , we have + // to put a cut on them which is facilitated by a cut on the number of + // clusters in Pad1. Do not change this. + trackproducer_hw(Pad1, Pad2, Pad3, outTrk, LOOKUP); + for (int I = 0; I < NTRK; I++) { + if (outTrk[I].Pad1.Seed.Amp > 0) { + ldmx::TrigScintTrack trk = makeTrack(outTrk[I]); + tracks_.push_back(trk); + } + } + event.add(output_collection_, tracks_); + tracks_.resize(0); + + return; +} + +ldmx::TrigScintTrack TrigScintFirmwareTracker::makeTrack(Track outTrk) { + // This takes a firmware track object and reverts it into an ldmx track + // object, unfortunately only retaining that information of the track that is + // retained in the firmware track. + ldmx::TrigScintTrack tr; + float pe = outTrk.Pad1.Seed.Amp + outTrk.Pad1.Sec.Amp; + pe += outTrk.Pad2.Seed.Amp + outTrk.Pad2.Sec.Amp; + pe += outTrk.Pad3.Seed.Amp + outTrk.Pad3.Sec.Amp; + tr.setCentroid(calcTCent(outTrk)); + calcResid(outTrk); + tr.setPE(pe); + return tr; +} + +} // namespace trigscint + +DECLARE_PRODUCER_NS(trigscint, TrigScintFirmwareTracker);