From 8d01fa35e944d5317cad2036c06df07972c4432d Mon Sep 17 00:00:00 2001
From: Aleksandra Poreba <aleksandra.poreba@cern.ch>
Date: Mon, 4 Mar 2024 14:49:35 +0100
Subject: [PATCH] Add triplet confirmation
---
.../src/TrigITkAccelerationTool.cxx | 5 +-
.../src/DoubletHelperFunctionsCuda_ITk.cuh | 9 +++
.../src/DoubletMatchingKernelCuda_ITk.cuh | 72 +++++++++++++++++--
.../src/SeedMakingDataStructures_ITk.h | 2 +-
.../ITkFastTrackFinderStandaloneConfig.py | 6 +-
5 files changed, 85 insertions(+), 9 deletions(-)
diff --git a/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationTool/src/TrigITkAccelerationTool.cxx b/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationTool/src/TrigITkAccelerationTool.cxx
index 9eb3c04a23a2..655172e3487e 100644
--- a/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationTool/src/TrigITkAccelerationTool.cxx
+++ b/Trigger/TrigAccel/TrigInDetAccel/TrigInDetAccelerationTool/src/TrigITkAccelerationTool.cxx
@@ -155,7 +155,7 @@ size_t TrigITkAccelerationTool::exportSeedMakingJob(const TrigCombinatorialSetti
for(int slice = 0;slice<nSlices;slice++) {
for(int layer = 0;layer<nLayers;layer++) {
int layerStart = spIdx;
- bool isBarrel = (layerTypes[layer] == 0);
+ bool isBarrel = (layerTypes[layer] == 0); // barrel = 0, ec != 0
std::vector<std::pair<int, const TrigSiSpacePointBase*> >& v = phiLArray[layer + slice*nLayers];
for(std::vector<std::pair<int, const TrigSiSpacePointBase*> >::iterator it = v.begin();it!=v.end();++it) {
const TrigSiSpacePointBase* sp = (*it).second;
@@ -196,6 +196,9 @@ int TrigITkAccelerationTool::extractTripletsFromOutput(std::shared_ptr<TrigAccel
output.clear();
for(int k=0;k<nTriplets;k++) {
+ // Check if a valid triplet was returned
+ if ((pOutput->m_innerIndex[k] == pOutput->m_outerIndex[k])) continue;
+
const TrigSiSpacePointBase& SPi = vsp[pOutput->m_innerIndex[k]];
const TrigSiSpacePointBase& SPm = vsp[pOutput->m_middleIndex[k]];
const TrigSiSpacePointBase& SPo = vsp[pOutput->m_outerIndex[k]];
diff --git a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletHelperFunctionsCuda_ITk.cuh b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletHelperFunctionsCuda_ITk.cuh
index b9c7eb1f57fe..73e5adb35e0c 100644
--- a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletHelperFunctionsCuda_ITk.cuh
+++ b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletHelperFunctionsCuda_ITk.cuh
@@ -32,6 +32,15 @@ __device__ static float getMaxDeltaLEta (float eta) {
else return eta*eta*eta*eta*1.7582417 + eta*eta*-129.67033 + 3324.61538;
}
+__device__ static int getInnerDoubletIdx (int pairIdx, int nOuter) {
+ return nOuter > 0 ? pairIdx/nOuter : 0;
+}
+
+__device__ static int getOuterDoubletIdx (int pairIdx, int nOuter, int startOfOuter) {
+ return startOfOuter + pairIdx % nOuter;
+}
+
+
}
#endif
diff --git a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMatchingKernelCuda_ITk.cuh b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMatchingKernelCuda_ITk.cuh
index dae141504f6c..a9ebbaf922cc 100644
--- a/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMatchingKernelCuda_ITk.cuh
+++ b/Trigger/TrigAccel/TrigInDetCUDA/src/DoubletMatchingKernelCuda_ITk.cuh
@@ -47,7 +47,9 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
__shared__ int PairIdx_array[MAX_TRIPLETS_ITk];
__shared__ float Q_array[MAX_TRIPLETS_ITk]; // Quality score for a triplet Q=d0*d0
__shared__ int sortedIdx[MAX_TRIPLETS_ITk];
+ __shared__ float pt_array[MAX_TRIPLETS_ITk]; //Store curvature of triplets
+ __shared__ int innerDoubletTriplets[MAX_NUMBER_DOUBLETS_ITk];
__shared__ int iDoublet;
__shared__ int startOfOuter;
@@ -69,9 +71,14 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
const float phiMinus = dSettings->m_phiMinus;
const bool isFullscan = (dSettings->m_isFullScan == 1);
-
for(int itemIdx = blockIdx.x;itemIdx<maxItem;itemIdx += gridDim.x) {
+ for(int innerIdx = threadIdx.x; innerIdx<MAX_NUMBER_DOUBLETS_ITk;innerIdx+=blockDim.x) {
+ innerDoubletTriplets[innerIdx] = 0;
+ }
+
+ __syncthreads();
+
if(threadIdx.x==0) {
nTriplets = 0;
@@ -170,9 +177,9 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
__syncthreads();
for(int pairIdx = threadIdx.x;pairIdx<nPairs;pairIdx += blockDim.x) {
-
- int doublet_i = pairIdx / nOuter; // inner doublet
- int doublet_j = startOfOuter + pairIdx % nOuter; //outer doublet
+ // Decode pair index into inner and outer doublet indices
+ int doublet_i = GPUTrackSeedingItkHelpers::getInnerDoubletIdx(pairIdx, nOuter);
+ int doublet_j = GPUTrackSeedingItkHelpers::getOuterDoubletIdx(pairIdx, nOuter, startOfOuter);
if(doublet_i >= MAX_NUMBER_DOUBLETS_ITk || doublet_j >=MAX_NUMBER_DOUBLETS_ITk ) continue;
@@ -241,18 +248,74 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
//Calculate Quality
float Q = d0*d0;
+ float pt = ptCoeff*std::sqrt(1+A*A)/(B);
int l = atomicAdd(&nTriplets, 1);
if(l<MAX_TRIPLETS_ITk) {
PairIdx_array[l] = pairIdx;
Q_array[l] = Q;
sortedIdx[l] = 0;
+ pt_array[l] = pt;
+
+ // Count the outer doublets for the inner doublet that passed the selection
+ atomicAdd(&innerDoubletTriplets[doublet_i], 1);
+
}
}
__syncthreads();
+ // Confirm the triplets by looking for duplicates
+ if (nOuter == 0) continue;
+
+ for(int doublet_i = threadIdx.x; doublet_i<iDoublet; doublet_i += blockDim.x) {
+ if (innerDoubletTriplets[doublet_i] == 0) continue;
+ int spiIdx = spIdx_array[doublet_i]; // index of inner spacepoint of the doublet
+
+ // Find pairs with the inner doublet
+ for (int l=0; l < nTriplets && l < MAX_TRIPLETS_ITk; ++l) {
+ // Check if this doublet belongs to this triplet pair
+ if (doublet_i != GPUTrackSeedingItkHelpers::getInnerDoubletIdx(PairIdx_array[l], nOuter)) continue;
+ int doublet_j = GPUTrackSeedingItkHelpers::getOuterDoubletIdx(PairIdx_array[l], nOuter, startOfOuter); //outer doublet
+ int spoIdx = spIdx_array[doublet_j]; // index of outer spacepoint of the doublet
+
+ int nDupes = 0;
+ // Look for duplicates for this triplet
+ for (int l2=0; l2 < nTriplets && l2 < MAX_TRIPLETS_ITk; ++l2) {
+ if (l == l2) continue;
+
+ // Check if they share the inner spacepoints
+ if (doublet_i != GPUTrackSeedingItkHelpers::getInnerDoubletIdx(PairIdx_array[l2], nOuter)) continue;
+ int other_doublet_j = GPUTrackSeedingItkHelpers::getOuterDoubletIdx(PairIdx_array[l2], nOuter, startOfOuter);
+ int otherSpoIdx = spIdx_array[other_doublet_j];
+
+ // Triplet duplicates (from the same track) will not lay on the same layer
+ bool isBarrel1 = (dSpacepoints->m_type[spoIdx] == 0); // barrel = 0, ec != 0
+ bool isBarrel2 = (dSpacepoints->m_type[otherSpoIdx] == 0);
+ if ( isBarrel1 && isBarrel2 && std::abs(dSpacepoints->m_r[spoIdx]-dSpacepoints->m_r[otherSpoIdx]) < 20 ) {
+ continue;
+ } else if ( !isBarrel1 && !isBarrel2 && std::abs(dSpacepoints->m_z[spoIdx]-dSpacepoints->m_z[otherSpoIdx]) < 20 ) {
+ continue;
+ }
+
+ // Triplet duplicates (from the same track) will have the same curvature direction
+ if (pt_array[l] * pt_array[l2] < 0) continue;
+
+ // Triplet duplicates (from the same track) will have the same pt within stddev (based on 1GeV single muon)
+ float dPt = std::abs(1./pt_array[l] - 1./pt_array[l2]);
+ if (dPt > 1*0.00015243) continue;
+
+ ++nDupes;
+ }
+
+ // Reject the seeds without duplicates - will not have a track extension
+ if (nDupes < 1) Q_array[l] += 10000;
+
+ }
+ }
+
+ __syncthreads();
if(nTriplets>TRIPLET_BUFFER_DEPTH_ITk) {//sorting
@@ -281,6 +344,7 @@ __global__ static void doubletMatchingKernel_ITk(TrigAccel::ITk::SEED_FINDER_SET
int k = atomicAdd(&d_Out->m_nSeeds, nT);
int nStored=0;
for(int tIdx=0;tIdx<nTriplets;tIdx++) {
+ if (Q_array[tIdx] > 10000) continue;
if(sortedIdx[tIdx]<TRIPLET_BUFFER_DEPTH_ITk) {//store this triplet
int pairIdx = PairIdx_array[tIdx];
diff --git a/Trigger/TrigAccel/TrigInDetCUDA/src/SeedMakingDataStructures_ITk.h b/Trigger/TrigAccel/TrigInDetCUDA/src/SeedMakingDataStructures_ITk.h
index c5a85f7a542d..7772ab2e249c 100644
--- a/Trigger/TrigAccel/TrigInDetCUDA/src/SeedMakingDataStructures_ITk.h
+++ b/Trigger/TrigAccel/TrigInDetCUDA/src/SeedMakingDataStructures_ITk.h
@@ -16,7 +16,7 @@ static constexpr unsigned int NUM_TRIPLET_BLOCKS_ITk = 1024;
static constexpr unsigned int NUM_TRIPLET_THREADS_ITk = 1024;
static constexpr unsigned int NUM_DOUBLET_THREADS_ITk = 16;
static constexpr unsigned int MAX_TRIPLETS_ITk = 300;
-static constexpr unsigned int TRIPLET_BUFFER_DEPTH_ITk = 3;
+static constexpr unsigned int TRIPLET_BUFFER_DEPTH_ITk = 2;
typedef struct doubletInfoITk {
public:
diff --git a/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/ITkFastTrackFinderStandaloneConfig.py b/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/ITkFastTrackFinderStandaloneConfig.py
index 49d5885581ae..e770f41e1776 100644
--- a/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/ITkFastTrackFinderStandaloneConfig.py
+++ b/Trigger/TrigAlgorithms/TrigFastTrackFinder/python/ITkFastTrackFinderStandaloneConfig.py
@@ -23,9 +23,9 @@ def ITkFastTrackFinderStandaloneCfg(flags, SiSPSeededTrackCollectionKey = None):
inDetAccelSvc = CompFactory.TrigInDetAccelerationSvc("TrigInDetAccelerationSvc")
inDetAccelSvc.useITkGeometry = True # Allows to read and export the ITk geometry
inDetAccelSvc.MiddleSpacePointLayers = [81000, 82000,
- 91005, 90014, 92000, 92001, 92002, 92003, 92004, 92005, 92006, 92007, 92008, 92009, 92010,
+ 91005, 91004, 90014, 90013, 92000, 92001, 92002, 92003, 92004, 92005, 92006, 92007, 92008, 92009, 92010,
92011, 92012, 92013, 92014, 92015, 92016, 92017, 92018, 92019, 92020, 92021, 92022,
- 71005, 70014, 72000, 72001, 72002, 72003, 72004, 72005, 72006, 72007, 72008, 72009, 72010,
+ 71005, 71004, 70014, 70013, 72000, 72001, 72002, 72003, 72004, 72005, 72006, 72007, 72008, 72009, 72010,
72011, 72012, 72013, 72014, 72015, 72016, 72017, 72018, 72019, 72020, 72021, 72022
]
acc.addService(inDetAccelSvc)
@@ -72,7 +72,7 @@ def ITkFastTrackFinderStandaloneCfg(flags, SiSPSeededTrackCollectionKey = None):
useNewLayerNumberScheme = True,
MinHits = 3,
ITkMode = True, # Allows ftf to use the new TrigTrackSeedGenerator for ITk
- useGPU = False,
+ useGPU = flags.Trigger.InDetTracking.doGPU,
StandaloneMode = True, # Allows ftf to be run as an offline algorithm with reco_tf
doTrackRefit = False,
FreeClustersCut = 1,
--
GitLab