diff --git a/DataFormats/Common/src/classes_def.xml b/DataFormats/Common/src/classes_def.xml
index 87f037e09d5c0..b14858cc8a163 100644
--- a/DataFormats/Common/src/classes_def.xml
+++ b/DataFormats/Common/src/classes_def.xml
@@ -205,11 +205,19 @@
+
+
-
+
+
+
+
+
+
+
diff --git a/RecoTracker/LST/plugins/LSTOutputConverter.cc b/RecoTracker/LST/plugins/LSTOutputConverter.cc
index 200db3cc5c9e6..d20026c89a012 100644
--- a/RecoTracker/LST/plugins/LSTOutputConverter.cc
+++ b/RecoTracker/LST/plugins/LSTOutputConverter.cc
@@ -165,28 +165,15 @@ void LSTOutputConverter::produce(edm::Event& iEvent, const edm::EventSetup& iSet
recHits.push_back(hit.clone());
}
- // pixel-seeded TCs from LST always have 4 pixel hits
- unsigned int const nPixelHits = (iType == lst::LSTObjType::T5 || iType == lst::LSTObjType::T4) ? 0 : 4;
- unsigned int nHits = 0;
- switch (iType) {
- case lst::LSTObjType::T5:
- nHits = lst::Params_T5::kHits;
- break;
- case lst::LSTObjType::pT3:
- nHits = lst::Params_pT3::kHits;
- break;
- case lst::LSTObjType::pT5:
- nHits = lst::Params_pT5::kHits;
- break;
- case lst::LSTObjType::pLS:
- nHits = lst::Params_pLS::kHits;
- break;
- case lst::LSTObjType::T4:
- nHits = lst::Params_T4::kHits;
- break;
+ // The pixel hits are packed into first kPixelLayerSlots layer slots.
+ for (unsigned int layerSlot = lst::Params_TC::kPixelLayerSlots; layerSlot < lst::Params_TC::kLayers; ++layerSlot) {
+ for (unsigned int hitSlot = 0; hitSlot < lst::Params_TC::kHitsPerLayer; ++hitSlot) {
+ unsigned int hitIdx = lstOutput_view.hitIndices()[i][layerSlot][hitSlot];
+ if (hitIdx == lst::kTCEmptyHitIdx)
+ continue;
+ recHits.push_back(OTHits[hitIdx]->clone());
+ }
}
- for (unsigned int j = nPixelHits; j < nHits; j++)
- recHits.push_back(OTHits[lstOutput_view.hitIndices()[i][j]]->clone());
recHits.sort([](const auto& a, const auto& b) {
const auto asub = a.det()->subDetector();
@@ -215,7 +202,7 @@ void LSTOutputConverter::produce(edm::Event& iEvent, const edm::EventSetup& iSet
if (includeNonpLSTSs_ || (iType == lst::LSTObjType::T5 || iType == lst::LSTObjType::T4)) {
using Hit = SeedingHitSet::ConstRecHitPointer;
std::vector hitsForSeed;
- hitsForSeed.reserve(nHits);
+ hitsForSeed.reserve(recHits.size());
int n = 0;
unsigned int firstLayer;
for (auto const& hit : recHits) {
diff --git a/RecoTracker/LSTCore/interface/Common.h b/RecoTracker/LSTCore/interface/Common.h
index 0fa84d28aee54..78c1b90bad4e7 100644
--- a/RecoTracker/LSTCore/interface/Common.h
+++ b/RecoTracker/LSTCore/interface/Common.h
@@ -35,6 +35,9 @@ namespace lst {
constexpr unsigned int size_superbins = 45000;
+ constexpr uint16_t kTCEmptyLowerModule = 0xFFFF; // Sentinel for empty lowerModule index
+ constexpr unsigned int kTCEmptyHitIdx = 0xFFFFFFFF; // Sentinel for empty hit slots
+
// Half precision wrapper functions.
#if defined(FP16_Base)
#define __F2H __float2half
@@ -94,10 +97,15 @@ namespace lst {
using ArrayUxHits = edm::StdArray;
};
struct Params_TC {
- static constexpr int kLayers = 7, kHits = 14;
+ static constexpr int kLayers = 13;
+ static constexpr int kHitsPerLayer = 2;
+ // Number of layers resevered for pixel hits.
+ static constexpr int kPixelLayerSlots = 2;
+ static constexpr int kHits = kLayers * kHitsPerLayer;
using ArrayU8xLayers = edm::StdArray;
using ArrayU16xLayers = edm::StdArray;
- using ArrayUxHits = edm::StdArray;
+ using ArrayUHitsPerLayer = edm::StdArray;
+ using ArrayUxHits = edm::StdArray;
};
using ArrayIx2 = edm::StdArray;
diff --git a/RecoTracker/LSTCore/interface/QuintupletsSoA.h b/RecoTracker/LSTCore/interface/QuintupletsSoA.h
index 1d3c6bd59408e..21fa680463a73 100644
--- a/RecoTracker/LSTCore/interface/QuintupletsSoA.h
+++ b/RecoTracker/LSTCore/interface/QuintupletsSoA.h
@@ -32,6 +32,7 @@ namespace lst {
SOA_COLUMN(float, rzChiSquared), // r-z only chi2
SOA_COLUMN(float, chiSquared),
SOA_COLUMN(float, nonAnchorChiSquared),
+ SOA_COLUMN(float, dnnScore),
SOA_COLUMN(float, dBeta1),
SOA_COLUMN(float, dBeta2));
diff --git a/RecoTracker/LSTCore/interface/alpaka/Common.h b/RecoTracker/LSTCore/interface/alpaka/Common.h
index ddef1abc9c7c5..5ed0546351490 100644
--- a/RecoTracker/LSTCore/interface/alpaka/Common.h
+++ b/RecoTracker/LSTCore/interface/alpaka/Common.h
@@ -45,6 +45,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
HOST_DEVICE_CONSTANT int kNTripletThreshold = 1000;
// To be updated with std::numeric_limits::infinity() in the code and data files
HOST_DEVICE_CONSTANT float kVerticalModuleSlope = 123456789.0;
+ HOST_DEVICE_CONSTANT int kLogicalOTLayers = 11; // logical OT layers are 1..11
HOST_DEVICE_CONSTANT float kMiniDeltaTilted[3] = {0.26f, 0.26f, 0.26f};
HOST_DEVICE_CONSTANT float kMiniDeltaFlat[6] = {0.26f, 0.16f, 0.16f, 0.18f, 0.18f, 0.18f};
diff --git a/RecoTracker/LSTCore/src/alpaka/LSTEvent.dev.cc b/RecoTracker/LSTCore/src/alpaka/LSTEvent.dev.cc
index a683bdbd03a28..bd0be0b948628 100644
--- a/RecoTracker/LSTCore/src/alpaka/LSTEvent.dev.cc
+++ b/RecoTracker/LSTCore/src/alpaka/LSTEvent.dev.cc
@@ -644,6 +644,26 @@ void LSTEvent::createTrackCandidates(bool no_pls_dupclean, bool tc_pls_triplets)
pixelSegmentsDC_->const_view(),
tc_pls_triplets);
+ // Get number of TCs to configure grid
+ auto nTrackCandidates_buf_h = cms::alpakatools::make_host_buffer(queue_);
+ auto nTrackCandidates_buf_d =
+ cms::alpakatools::make_device_view(queue_, (*trackCandidatesBaseDC_)->nTrackCandidates());
+ alpaka::memcpy(queue_, nTrackCandidates_buf_h, nTrackCandidates_buf_d);
+ alpaka::wait(queue_);
+ unsigned int nTC = *nTrackCandidates_buf_h.data();
+
+ auto const wd = cms::alpakatools::make_workdiv(nTC, 128);
+
+ alpaka::exec(queue_,
+ wd,
+ ExtendTrackCandidatesFromDupT5{},
+ modules_.const_view(),
+ rangesDC_->const_view(),
+ quintupletsDC_->const_view(),
+ quintupletsDC_->const_view(),
+ trackCandidatesBaseDC_->view(),
+ trackCandidatesExtendedDC_->view());
+
// Check if either n_max_pixel_track_candidates or n_max_nonpixel_track_candidates was reached
auto nTrackCanpT5Host_buf = cms::alpakatools::make_host_buffer(queue_);
auto nTrackCanpT3Host_buf = cms::alpakatools::make_host_buffer(queue_);
diff --git a/RecoTracker/LSTCore/src/alpaka/NeuralNetwork.h b/RecoTracker/LSTCore/src/alpaka/NeuralNetwork.h
index 9e1ff995188d2..dc401c71cbbac 100644
--- a/RecoTracker/LSTCore/src/alpaka/NeuralNetwork.h
+++ b/RecoTracker/LSTCore/src/alpaka/NeuralNetwork.h
@@ -202,7 +202,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
const unsigned int mdIndex5,
const float innerRadius,
const float outerRadius,
- const float bridgeRadius) {
+ const float bridgeRadius,
+ float& dnnScore) {
// Constants
constexpr unsigned int kInputFeatures = 23;
constexpr unsigned int kHiddenFeatures = 32;
@@ -281,7 +282,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
// Layer 3: Linear + Sigmoid
linear_layer(x_2, x_3, dnn::t5dnn::wgtT_output_layer, dnn::t5dnn::bias_output_layer);
- float x_5 = sigmoid_activation(acc, x_3[0]);
+ dnnScore = sigmoid_activation(acc, x_3[0]);
// Get the bin index based on abs(eta) of first hit and t5_pt
float t5_pt = innerRadius * lst::k2Rinv1GeVf * 2;
@@ -289,8 +290,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
uint8_t pt_index = (t5_pt > 5.0f);
uint8_t bin_index = (eta1 > 2.5f) ? (dnn::kEtaBins - 1) : static_cast(eta1 / dnn::kEtaSize);
- // Compare x_5 to the cut value for the relevant bin
- return x_5 > dnn::t5dnn::kWp[pt_index][bin_index];
+ // Compare output to the cut value for the relevant bin
+ return dnnScore > dnn::t5dnn::kWp[pt_index][bin_index];
}
} // namespace t5dnn
diff --git a/RecoTracker/LSTCore/src/alpaka/Quintuplet.h b/RecoTracker/LSTCore/src/alpaka/Quintuplet.h
index c4585c62f9f79..a943facc28be4 100644
--- a/RecoTracker/LSTCore/src/alpaka/Quintuplet.h
+++ b/RecoTracker/LSTCore/src/alpaka/Quintuplet.h
@@ -46,7 +46,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
uint8_t layer,
unsigned int quintupletIndex,
const float (&t5Embed)[Params_T5::kEmbed],
- bool tightCutFlag) {
+ bool tightCutFlag,
+ float dnnScore) {
quintuplets.tripletIndices()[quintupletIndex][0] = innerTripletIndex;
quintuplets.tripletIndices()[quintupletIndex][1] = outerTripletIndex;
@@ -88,6 +89,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
quintuplets.nonAnchorChiSquared()[quintupletIndex] = nonAnchorChiSquared;
quintuplets.dBeta1()[quintupletIndex] = dBeta1;
quintuplets.dBeta2()[quintupletIndex] = dBeta2;
+ quintuplets.dnnScore()[quintupletIndex] = dnnScore;
CMS_UNROLL_LOOP
for (unsigned int i = 0; i < Params_T5::kEmbed; ++i) {
@@ -1491,6 +1493,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
float& nonAnchorChiSquared,
float& dBeta1,
float& dBeta2,
+ float& dnnScore,
bool& tightCutFlag,
float (&t5Embed)[Params_T5::kEmbed],
const float ptCut) {
@@ -1531,7 +1534,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
fifthMDIndex,
innerRadius,
outerRadius,
- bridgeRadius);
+ bridgeRadius,
+ dnnScore);
tightCutFlag = tightCutFlag and inference; // T5-in-TC cut
if (!inference) // T5-building cut
return false;
@@ -1754,7 +1758,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
uint16_t lowerModule5 = triplets.lowerModuleIndices()[outerTripletIndex][2];
float innerRadius, outerRadius, bridgeRadius, regressionCenterX, regressionCenterY, regressionRadius,
- rzChiSquared, chiSquared, nonAnchorChiSquared, dBeta1, dBeta2; //required for making distributions
+ rzChiSquared, chiSquared, nonAnchorChiSquared, dBeta1, dBeta2,
+ dnnScore; //required for making distributions
float t5Embed[Params_T5::kEmbed] = {0.f};
@@ -1783,6 +1788,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
nonAnchorChiSquared,
dBeta1,
dBeta2,
+ dnnScore,
tightCutFlag,
t5Embed,
ptCut);
@@ -1841,7 +1847,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
layer,
quintupletIndex,
t5Embed,
- tightCutFlag);
+ tightCutFlag,
+ dnnScore);
triplets.partOfT5()[quintuplets.tripletIndices()[quintupletIndex][0]] = true;
triplets.partOfT5()[quintuplets.tripletIndices()[quintupletIndex][1]] = true;
@@ -1890,7 +1897,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
uint16_t lowerModule5 = triplets.lowerModuleIndices()[outerTripletIndex][2];
float innerRadius, outerRadius, bridgeRadius, regressionCenterX, regressionCenterY, regressionRadius,
- rzChiSquared, chiSquared, nonAnchorChiSquared, dBeta1, dBeta2; //required for making distributions
+ rzChiSquared, chiSquared, nonAnchorChiSquared, dBeta1, dBeta2,
+ dnnScore; //required for making distributions
float t5Embed[Params_T5::kEmbed] = {0.f};
@@ -1919,6 +1927,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
nonAnchorChiSquared,
dBeta1,
dBeta2,
+ dnnScore,
tightCutFlag,
t5Embed,
ptCut);
@@ -1974,7 +1983,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
layer,
quintupletIndex,
t5Embed,
- tightCutFlag);
+ tightCutFlag,
+ dnnScore);
triplets.partOfT5()[quintuplets.tripletIndices()[quintupletIndex][0]] = true;
triplets.partOfT5()[quintuplets.tripletIndices()[quintupletIndex][1]] = true;
@@ -2045,7 +2055,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
float innerRadius, outerRadius, bridgeRadius;
float regCx, regCy, regR;
- float rzChi2, chi2, nonAnchorChi2, dBeta1, dBeta2;
+ float rzChi2, chi2, nonAnchorChi2, dBeta1, dBeta2, dnnScore;
float t5Embed[Params_T5::kEmbed] = {0.f};
bool tightFlag = false;
@@ -2072,6 +2082,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
nonAnchorChi2,
dBeta1,
dBeta2,
+ dnnScore,
tightFlag,
t5Embed,
ptCut);
diff --git a/RecoTracker/LSTCore/src/alpaka/TrackCandidate.h b/RecoTracker/LSTCore/src/alpaka/TrackCandidate.h
index 8cfe00506588b..dbb67c79805f5 100644
--- a/RecoTracker/LSTCore/src/alpaka/TrackCandidate.h
+++ b/RecoTracker/LSTCore/src/alpaka/TrackCandidate.h
@@ -1,6 +1,8 @@
#ifndef RecoTracker_LSTCore_src_alpaka_TrackCandidate_h
#define RecoTracker_LSTCore_src_alpaka_TrackCandidate_h
+#include
+
#include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
#include "FWCore/Utilities/interface/CMSUnrollLoop.h"
#include "HeterogeneousCore/AlpakaMath/interface/deltaPhi.h"
@@ -35,11 +37,40 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
candsExtended.objectIndices()[trackCandidateIndex][0] = trackletIndex;
candsExtended.objectIndices()[trackCandidateIndex][1] = trackletIndex;
+ // Initialize all slots to empty
+ auto& tcHits = candsBase.hitIndices()[trackCandidateIndex];
+ auto& tcModules = candsExtended.lowerModuleIndices()[trackCandidateIndex];
+ auto& tcLayers = candsExtended.logicalLayers()[trackCandidateIndex];
+ CMS_UNROLL_LOOP for (int layerSlot = 0; layerSlot < Params_TC::kLayers; ++layerSlot) {
+ tcLayers[layerSlot] = 0;
+ tcModules[layerSlot] = lst::kTCEmptyLowerModule;
+ tcHits[layerSlot][0] = lst::kTCEmptyHitIdx;
+ tcHits[layerSlot][1] = lst::kTCEmptyHitIdx;
+ }
+
// Order explanation in https://github.com/SegmentLinking/TrackLooper/issues/267
- candsBase.hitIndices()[trackCandidateIndex][0] = hitIndices[0];
- candsBase.hitIndices()[trackCandidateIndex][1] = hitIndices[2];
- candsBase.hitIndices()[trackCandidateIndex][2] = hitIndices[1];
- candsBase.hitIndices()[trackCandidateIndex][3] = hitIndices[3];
+ tcHits[0][0] = hitIndices[0];
+ tcHits[0][1] = hitIndices[2];
+ tcHits[1][0] = hitIndices[1];
+ tcHits[1][1] = hitIndices[3];
+ }
+
+ ALPAKA_FN_ACC ALPAKA_FN_INLINE void addTrackCandidateLayerHits(
+ TrackCandidatesBase& candsBase,
+ TrackCandidatesExtended& candsExtended,
+ unsigned int trackCandidateIndex,
+ int layerSlot, // 0..12 (0/1 = pixel, 2..12 = OT logical layers 1..11)
+ uint8_t logicalLayer, // 0 for pixel, 1..11 for OT
+ uint16_t lowerModule,
+ unsigned int hitIndex0,
+ unsigned int hitIndex1) {
+ auto& tcModules = candsExtended.lowerModuleIndices()[trackCandidateIndex];
+ auto& tcLayers = candsExtended.logicalLayers()[trackCandidateIndex];
+ auto& tcHits = candsBase.hitIndices()[trackCandidateIndex];
+ tcLayers[layerSlot] = logicalLayer;
+ tcModules[layerSlot] = lowerModule;
+ tcHits[layerSlot][0] = hitIndex0;
+ tcHits[layerSlot][1] = hitIndex1;
}
ALPAKA_FN_ACC ALPAKA_FN_INLINE void addTrackCandidateToMemory(TrackCandidatesBase& candsBase,
@@ -63,50 +94,58 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
candsExtended.objectIndices()[trackCandidateIndex][0] = innerTrackletIndex;
candsExtended.objectIndices()[trackCandidateIndex][1] = outerTrackletIndex;
- size_t limits = trackCandidateType == LSTObjType::pT5 ? Params_pT5::kLayers : Params_pT3::kLayers;
-
- //send the starting pointer to the logicalLayer and hitIndices
- for (size_t i = 0; i < limits; i++) {
- candsExtended.logicalLayers()[trackCandidateIndex][i] = logicalLayerIndices[i];
- candsExtended.lowerModuleIndices()[trackCandidateIndex][i] = lowerModuleIndices[i];
+ // Initialize all slots to empty
+ auto& tcHits = candsBase.hitIndices()[trackCandidateIndex];
+ auto& tcModules = candsExtended.lowerModuleIndices()[trackCandidateIndex];
+ auto& tcLayers = candsExtended.logicalLayers()[trackCandidateIndex];
+ CMS_UNROLL_LOOP for (int layerSlot = 0; layerSlot < Params_TC::kLayers; ++layerSlot) {
+ tcLayers[layerSlot] = 0; // 0 is "pixel" when filled
+ tcModules[layerSlot] = lst::kTCEmptyLowerModule;
+ tcHits[layerSlot][0] = lst::kTCEmptyHitIdx;
+ tcHits[layerSlot][1] = lst::kTCEmptyHitIdx;
}
- for (size_t i = 0; i < 2 * limits; i++) {
- candsBase.hitIndices()[trackCandidateIndex][i] = hitIndices[i];
- }
- candsExtended.centerX()[trackCandidateIndex] = __F2H(centerX);
- candsExtended.centerY()[trackCandidateIndex] = __F2H(centerY);
- candsExtended.radius()[trackCandidateIndex] = __F2H(radius);
- }
- ALPAKA_FN_ACC ALPAKA_FN_INLINE void addT4TrackCandidateToMemory(TrackCandidatesBase& candsBase,
- TrackCandidatesExtended& candsExtended,
- LSTObjType trackCandidateType,
- unsigned int innerTrackletIndex,
- unsigned int outerTrackletIndex,
- uint8_t* logicalLayerIndices,
- uint16_t* lowerModuleIndices,
- unsigned int* hitIndices,
- int pixelSeedIndex,
- float centerX,
- float centerY,
- float radius,
- unsigned int trackCandidateIndex,
- unsigned int directObjectIndex) {
- candsBase.trackCandidateType()[trackCandidateIndex] = trackCandidateType;
- candsExtended.directObjectIndices()[trackCandidateIndex] = directObjectIndex;
- candsBase.pixelSeedIndex()[trackCandidateIndex] = pixelSeedIndex;
+ // Configuration based on Type
+ int nLayersToProcess = 0;
+ int nPixelLayers = 0;
+
+ if (trackCandidateType == LSTObjType::T5) {
+ nLayersToProcess = Params_T5::kLayers;
+ } else if (trackCandidateType == LSTObjType::pT5) {
+ nLayersToProcess = Params_pT5::kLayers;
+ nPixelLayers = Params_TC::kPixelLayerSlots;
+ } else if (trackCandidateType == LSTObjType::T4) {
+ nLayersToProcess = Params_T4::kLayers;
+ } else if (trackCandidateType == LSTObjType::pT3) {
+ nLayersToProcess = Params_pT3::kLayers;
+ nPixelLayers = Params_TC::kPixelLayerSlots;
+ }
- candsExtended.objectIndices()[trackCandidateIndex][0] = innerTrackletIndex;
- candsExtended.objectIndices()[trackCandidateIndex][1] = outerTrackletIndex;
+ CMS_UNROLL_LOOP
+ for (int i = 0; i < Params_TC::kLayers; ++i) {
+ if (i >= nLayersToProcess)
+ break;
+
+ uint8_t logicalLayer = logicalLayerIndices[i];
+ uint16_t lowerModule = lowerModuleIndices[i];
+ unsigned int hit0 = hitIndices[2 * i];
+ unsigned int hit1 = hitIndices[2 * i + 1];
+
+ int layerSlot;
+
+ if (i < nPixelLayers) {
+ // Pixel layers occupy slots 0 and 1 strictly
+ layerSlot = i;
+ logicalLayer = 0;
+ } else {
+ // OT layers are mapped: (LogicalLayer - 1) + kPixelLayerSlots
+ layerSlot = (logicalLayer - 1) + Params_TC::kPixelLayerSlots;
+ }
- //send the starting pointer to the logicalLayer and hitIndices
- for (int i = 0; i < Params_T4::kLayers; i++) {
- candsExtended.logicalLayers()[trackCandidateIndex][i] = logicalLayerIndices[i];
- candsExtended.lowerModuleIndices()[trackCandidateIndex][i] = lowerModuleIndices[i];
- }
- for (int i = 0; i < 2 * Params_T4::kLayers; i++) {
- candsBase.hitIndices()[trackCandidateIndex][i] = hitIndices[i];
+ addTrackCandidateLayerHits(
+ candsBase, candsExtended, trackCandidateIndex, layerSlot, logicalLayer, lowerModule, hit0, hit1);
}
+
candsExtended.centerX()[trackCandidateIndex] = __F2H(centerX);
candsExtended.centerY()[trackCandidateIndex] = __F2H(centerY);
candsExtended.radius()[trackCandidateIndex] = __F2H(radius);
@@ -609,6 +648,227 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
}
};
+ ALPAKA_FN_ACC ALPAKA_FN_INLINE void countSharedT5HitsAndFindUnmatched(QuintupletsConst quintuplets,
+ unsigned int testT5Index,
+ unsigned int refT5Index,
+ int& sharedHitCount,
+ int& unmatchedLayerSlot) {
+ sharedHitCount = 0;
+ unmatchedLayerSlot = -1;
+
+ CMS_UNROLL_LOOP
+ for (int layerIndex = 0; layerIndex < Params_T5::kLayers; ++layerIndex) {
+ const unsigned int candidateHit0 = quintuplets.hitIndices()[testT5Index][2 * layerIndex + 0];
+ const unsigned int candidateHit1 = quintuplets.hitIndices()[testT5Index][2 * layerIndex + 1];
+
+ bool hit0InBase = false;
+ bool hit1InBase = false;
+
+ CMS_UNROLL_LOOP
+ for (int baseHitIndex = 0; baseHitIndex < Params_T5::kHits; ++baseHitIndex) {
+ const unsigned int baseHit = quintuplets.hitIndices()[refT5Index][baseHitIndex];
+ if (candidateHit0 == baseHit)
+ hit0InBase = true;
+ if (candidateHit1 == baseHit)
+ hit1InBase = true;
+ if (hit0InBase && hit1InBase)
+ break;
+ }
+
+ if (hit0InBase)
+ ++sharedHitCount;
+ if (hit1InBase)
+ ++sharedHitCount;
+
+ if (!hit0InBase && !hit1InBase) {
+ unmatchedLayerSlot = layerIndex;
+ }
+ }
+ }
+
+ struct ExtendTrackCandidatesFromDupT5 {
+ static constexpr int kPackedScoreShift = 32; // Bit shift for packing the score (upper 32 bits).
+ static constexpr int kPackedIndexShift = 4; // Bit shift for packing the T5 index.
+ static constexpr unsigned int kPackedIndexMask = 0xFFFFFFF; // Mask for the 28-bit T5 index.
+ static constexpr unsigned int kPackedSlotMask = 0xF; // Mask for the 4-bit layer slot.
+ static constexpr int kT5DuplicateMinSharedHits = 8; // Minimum shared hits required to consider a T5 for merging.
+
+ ALPAKA_FN_ACC void operator()(Acc1D const& acc,
+ ModulesConst modules,
+ ObjectRangesConst ranges,
+ QuintupletsConst quintuplets,
+ QuintupletsOccupancyConst quintupletsOccupancy,
+ TrackCandidatesBase candsBase,
+ TrackCandidatesExtended candsExtended) const {
+ // Assert that the number of blocks equals nTC
+ ALPAKA_ASSERT_ACC((alpaka::getWorkDiv(acc)[0u] == candsBase.nTrackCandidates()));
+
+ // Shared memory: Packed best candidate info per logical OT layer (1..11)
+ // Format: [Score (32 bits) | Quintuplet Index (28 bits) | Layer Slot (4 bits)]
+ uint64_t* sharedBestPacked = alpaka::declareSharedVar(acc);
+
+ // In 1D, the Grid index [0] is the Block index (which is our TC index)
+ const unsigned int trackCandidateIndex = alpaka::getIdx(acc)[0u];
+
+ // Initialize shared memory once per block
+ if (cms::alpakatools::once_per_block(acc)) {
+ for (int logicalLayerBin = 0; logicalLayerBin < lst::kLogicalOTLayers; ++logicalLayerBin) {
+ sharedBestPacked[logicalLayerBin] = 0; // 0.0f score, 0 index, 0 slot
+ }
+ }
+ alpaka::syncBlockThreads(acc);
+
+ // Only consider merging T5s onto T5/pT5 TCs for now
+ const LSTObjType trackCandidateType = candsBase.trackCandidateType()[trackCandidateIndex];
+ if (!(trackCandidateType == LSTObjType::T5 || trackCandidateType == LSTObjType::pT5))
+ return;
+
+ // Base quintuplet (for T5: objectIndices[0], for pT5: objectIndices[1])
+ const unsigned int refT5Index = (trackCandidateType == LSTObjType::T5)
+ ? candsExtended.objectIndices()[trackCandidateIndex][0]
+ : candsExtended.objectIndices()[trackCandidateIndex][1];
+
+ const float baseEta = __H2F(quintuplets.eta()[refT5Index]);
+ const float basePhi = __H2F(quintuplets.phi()[refT5Index]);
+ const uint8_t baseStartLogicalLayer = quintuplets.logicalLayers()[refT5Index][0];
+
+ // Module range to scan
+ int lowerModuleBegin = 0;
+ int lowerModuleEnd = static_cast(modules.nLowerModules());
+
+ // If starting at layer 1, restrict to the module of the second hit
+ if (baseStartLogicalLayer == 1) {
+ lowerModuleBegin = quintuplets.lowerModuleIndices()[refT5Index][1];
+ lowerModuleEnd = lowerModuleBegin + 1;
+ }
+
+ // Linear thread index and block dimension in 1D
+ const int threadIndexFlat = alpaka::getIdx(acc)[0u];
+ const int blockDimFlat = alpaka::getWorkDiv(acc)[0u];
+
+ // Scan over lower modules
+ for (int lowerModuleIndex = lowerModuleBegin + threadIndexFlat; lowerModuleIndex < lowerModuleEnd;
+ lowerModuleIndex += blockDimFlat) {
+ const int firstQuintupletInModule = ranges.quintupletModuleIndices()[lowerModuleIndex];
+ if (firstQuintupletInModule == -1)
+ continue;
+
+ const unsigned int nQuintupletsInModule = quintupletsOccupancy.nQuintuplets()[lowerModuleIndex];
+
+ // Scan over quintuplets in this module
+ for (unsigned int quintupletOffset = 0; quintupletOffset < nQuintupletsInModule; ++quintupletOffset) {
+ const unsigned int testT5Index = firstQuintupletInModule + quintupletOffset;
+ if (testT5Index == refT5Index)
+ continue;
+
+ // Require different starting layer for now
+ if (quintuplets.logicalLayers()[testT5Index][0] == baseStartLogicalLayer)
+ continue;
+
+ // Quick eta / phi window selection
+ const float candidateEta = __H2F(quintuplets.eta()[testT5Index]);
+ if (alpaka::math::abs(acc, baseEta - candidateEta) > 0.1f)
+ continue;
+
+ const float candidatePhi = __H2F(quintuplets.phi()[testT5Index]);
+ if (alpaka::math::abs(acc, cms::alpakatools::deltaPhi(acc, basePhi, candidatePhi)) > 0.1f)
+ continue;
+
+ // Embedding distance
+ float embedDistance2 = 0.f;
+ CMS_UNROLL_LOOP
+ for (unsigned int embedIndex = 0; embedIndex < Params_T5::kEmbed; ++embedIndex) {
+ const float diff =
+ quintuplets.t5Embed()[refT5Index][embedIndex] - quintuplets.t5Embed()[testT5Index][embedIndex];
+ embedDistance2 += diff * diff;
+ }
+ if (embedDistance2 > 1.0f)
+ continue;
+
+ // Hit matching against base T5 hits
+ int sharedHitCount = 0;
+ int unmatchedLayerSlot = -1;
+
+ countSharedT5HitsAndFindUnmatched(quintuplets, testT5Index, refT5Index, sharedHitCount, unmatchedLayerSlot);
+
+ if (sharedHitCount < kT5DuplicateMinSharedHits)
+ continue;
+ if (unmatchedLayerSlot < 0)
+ continue;
+
+ // Candidate score is the T5 DNN output
+ const float candidateScore = quintuplets.dnnScore()[testT5Index];
+
+ // New OT logical layer from the candidate
+ const uint8_t newLogicalLayer = quintuplets.logicalLayers()[testT5Index][unmatchedLayerSlot];
+ const int logicalLayerBin = static_cast(newLogicalLayer) - 1; // 0..kLogicalOTLayers-1
+
+ // Pack the Score, Index, and Slot into 64 bits for atomic update
+ uint64_t scoreBits = (uint64_t)std::bit_cast(candidateScore);
+ uint64_t newPacked = (scoreBits << kPackedScoreShift) |
+ ((uint64_t)(testT5Index & kPackedIndexMask) << kPackedIndexShift) |
+ (unmatchedLayerSlot & kPackedSlotMask);
+
+ // Atomic CAS loop
+ uint64_t oldPacked = sharedBestPacked[logicalLayerBin];
+ while (true) {
+ const float oldScore = std::bit_cast((int)(oldPacked >> kPackedScoreShift));
+
+ // If we aren't strictly better, stop trying
+ if (candidateScore <= oldScore) {
+ break;
+ }
+
+ // Try to swap. atomicCas returns the value that was previously there
+ uint64_t assumedOld = alpaka::atomicCas(
+ acc, &sharedBestPacked[logicalLayerBin], oldPacked, newPacked, alpaka::hierarchy::Threads{});
+
+ if (assumedOld == oldPacked) {
+ break; // Success
+ } else {
+ oldPacked = assumedOld; // Failed, update view and retry
+ }
+ }
+ }
+ }
+
+ // One thread per block finalizes by actually extending the TC
+ alpaka::syncBlockThreads(acc);
+
+ if (cms::alpakatools::once_per_block(acc)) {
+ CMS_UNROLL_LOOP
+ for (int logicalLayerBin = 0; logicalLayerBin < lst::kLogicalOTLayers; ++logicalLayerBin) {
+ uint64_t bestPacked = sharedBestPacked[logicalLayerBin];
+
+ // Check if valid update occurred (non-zero score)
+ int bestScoreInt = (int)(bestPacked >> kPackedScoreShift);
+ if (bestScoreInt == 0)
+ continue;
+
+ // Unpack Index (bits 4-31) and Slot (bits 0-3)
+ const int bestT5Index = (int)((bestPacked >> kPackedIndexShift) & kPackedIndexMask);
+ const int bestT5LayerSlot = (int)(bestPacked & kPackedSlotMask);
+
+ const uint8_t logicalLayer = static_cast(logicalLayerBin + 1); // 1..11
+ const int layerSlot = (logicalLayer - 1) + Params_TC::kPixelLayerSlots; // OT slots for TC
+
+ const uint16_t lowerModuleIndex = quintuplets.lowerModuleIndices()[bestT5Index][bestT5LayerSlot];
+ const unsigned int hitIndex0 = quintuplets.hitIndices()[bestT5Index][2 * bestT5LayerSlot + 0];
+ const unsigned int hitIndex1 = quintuplets.hitIndices()[bestT5Index][2 * bestT5LayerSlot + 1];
+
+ addTrackCandidateLayerHits(candsBase,
+ candsExtended,
+ trackCandidateIndex,
+ layerSlot,
+ logicalLayer,
+ lowerModuleIndex,
+ hitIndex0,
+ hitIndex1);
+ }
+ }
+ }
+ };
+
struct AddpT5asTrackCandidate {
ALPAKA_FN_ACC void operator()(Acc1D const& acc,
uint16_t nLowerModules,
@@ -694,20 +954,20 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
break;
} else {
alpaka::atomicAdd(acc, &candsExtended.nTrackCandidatesT4(), 1u, alpaka::hierarchy::Threads{});
- addT4TrackCandidateToMemory(candsBase,
- candsExtended,
- LSTObjType::T4,
- quadrupletIndex,
- quadrupletIndex,
- quadruplets.logicalLayers()[quadrupletIndex].data(),
- quadruplets.lowerModuleIndices()[quadrupletIndex].data(),
- quadruplets.hitIndices()[quadrupletIndex].data(),
- -1 /*no pixel seed index for T4s*/,
- quadruplets.regressionCenterX()[quadrupletIndex],
- quadruplets.regressionCenterY()[quadrupletIndex],
- quadruplets.regressionRadius()[quadrupletIndex],
- trackCandidateIdx,
- quadrupletIndex);
+ addTrackCandidateToMemory(candsBase,
+ candsExtended,
+ LSTObjType::T4,
+ quadrupletIndex,
+ quadrupletIndex,
+ quadruplets.logicalLayers()[quadrupletIndex].data(),
+ quadruplets.lowerModuleIndices()[quadrupletIndex].data(),
+ quadruplets.hitIndices()[quadrupletIndex].data(),
+ -1 /*no pixel seed index for T4s*/,
+ quadruplets.regressionCenterX()[quadrupletIndex],
+ quadruplets.regressionCenterY()[quadrupletIndex],
+ quadruplets.regressionRadius()[quadrupletIndex],
+ trackCandidateIdx,
+ quadrupletIndex);
quadruplets.partOfTC()[quadrupletIndex] = true;
}
}
diff --git a/RecoTracker/LSTCore/standalone/code/core/AccessHelper.cc b/RecoTracker/LSTCore/standalone/code/core/AccessHelper.cc
index aac6ef9fc6b55..b6419ef97045c 100644
--- a/RecoTracker/LSTCore/standalone/code/core/AccessHelper.cc
+++ b/RecoTracker/LSTCore/standalone/code/core/AccessHelper.cc
@@ -543,3 +543,38 @@ std::tuple, std::vector> getHitIdxsAndHi
break;
}
}
+
+std::pair, std::vector> getHitIdxsAndTypesFromTC(LSTEvent* event,
+ unsigned int tc_idx) {
+ auto const& base = event->getTrackCandidatesBase();
+ auto const& ext = event->getTrackCandidatesExtended();
+ auto const& hitsBase = event->getInput();
+
+ std::vector hitIdx;
+ hitIdx.reserve(Params_TC::kHits);
+ std::vector hitType;
+ hitType.reserve(Params_TC::kHits);
+
+ for (int layerSlot = 0; layerSlot < Params_TC::kLayers; ++layerSlot) {
+ if (ext.lowerModuleIndices()[tc_idx][layerSlot] == lst::kTCEmptyLowerModule)
+ continue;
+
+ for (unsigned int hitSlot = 0; hitSlot < Params_TC::kHitsPerLayer; ++hitSlot) {
+ const unsigned int hitLocal = base.hitIndices()[tc_idx][layerSlot][hitSlot];
+
+ if (hitLocal == lst::kTCEmptyHitIdx)
+ continue;
+
+ // Get the GLOBAL ntuple indices
+ const unsigned int hitGlobal = hitsBase.idxs()[hitLocal];
+
+ // Determine the type from the hit's detid (1 = pixel, otherwise OT)
+ const unsigned int type = (hitsBase.detid()[hitLocal] == 1) ? 0u : 4u;
+
+ // Push the GLOBAL index and type
+ hitIdx.push_back(hitGlobal);
+ hitType.push_back(type);
+ }
+ }
+ return {hitIdx, hitType};
+}
diff --git a/RecoTracker/LSTCore/standalone/code/core/AccessHelper.h b/RecoTracker/LSTCore/standalone/code/core/AccessHelper.h
index 12b70ffb30b56..4145ccd965ad2 100644
--- a/RecoTracker/LSTCore/standalone/code/core/AccessHelper.h
+++ b/RecoTracker/LSTCore/standalone/code/core/AccessHelper.h
@@ -94,5 +94,7 @@ std::vector getLSsFromTC(LSTEvent* event, unsigned int TC);
std::vector getHitsFromTC(LSTEvent* event, unsigned int TC);
std::tuple, std::vector> getHitIdxsAndHitTypesFromTC(LSTEvent* event,
unsigned int TC);
+std::pair, std::vector> getHitIdxsAndTypesFromTC(LSTEvent* event,
+ unsigned int tc_idx);
#endif
diff --git a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc
index 2e8e428a65147..6428196e2c096 100644
--- a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc
+++ b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.cc
@@ -310,13 +310,17 @@ void createTrackCandidateBranches() {
//
// The container will hold per entry a track candidate built by LST in the event.
//
- ana.tx->createBranch>("tc_pt"); // pt
- ana.tx->createBranch>("tc_eta"); // eta
- ana.tx->createBranch>("tc_phi"); // phi
+ ana.tx->createBranch>("tc_pt"); // pt
+ ana.tx->createBranch>("tc_eta"); // eta
+ ana.tx->createBranch>("tc_phi"); // phi
+ ana.tx->createBranch>("tc_pMatched");
ana.tx->createBranch>("tc_type"); // type = 7 (pT5), 5 (pT3), 4 (T5), 8 (pLS), 9 (T4)
ana.tx->createBranch>("tc_isFake"); // 1 if tc is fake 0 other if not
ana.tx->createBranch>("tc_isDuplicate"); // 1 if tc is duplicate 0 other if not
ana.tx->createBranch>("tc_simIdx"); // idx of best matched (highest nhit and > 75%) simulated track
+ ana.tx->createBranch>("tc_nhitOT");
+ ana.tx->createBranch>("tc_nhits");
+ ana.tx->createBranch>("tc_nlayers");
// list of idx of all matched (> 0%) simulated track
ana.tx->createBranch>>("tc_simIdxAll");
// list of idx of all matched (> 0%) simulated track
@@ -469,6 +473,8 @@ void createQuintupletBranches() {
// The container will hold per entry a quintuplet built by LST in the event.
//
// pt (computed based on average of the 4 circles formed by, (1, 2, 3), (2, 3, 4), (3, 4, 5), (1, 3, 5)
+ ana.tx->createBranch>>("t5_embed");
+ ana.tx->createBranch>("t5_dnnScore");
ana.tx->createBranch>("t5_pt");
ana.tx->createBranch>("t5_eta"); // eta (computed based on last anchor hit's eta)
ana.tx->createBranch>("t5_phi"); // phi (computed based on first anchor hit's phi)
@@ -1581,6 +1587,14 @@ std::map setQuintupletBranches(LSTEvent* event,
ana.tx->pushbackToBranch("t5_bridgeRadius", __H2F(quintuplets.bridgeRadius()[t5Idx]));
ana.tx->pushbackToBranch("t5_outerRadius", __H2F(quintuplets.outerRadius()[t5Idx]));
ana.tx->pushbackToBranch("t5_pMatched", percent_matched);
+
+ std::vector current_t5_embed;
+ for (unsigned int i_embed = 0; i_embed < Params_T5::kEmbed; ++i_embed) {
+ current_t5_embed.push_back(quintuplets.t5Embed()[t5Idx][i_embed]);
+ }
+ ana.tx->pushbackToBranch>("t5_embed", current_t5_embed);
+ ana.tx->pushbackToBranch("t5_dnnScore", quintuplets.dnnScore()[t5Idx]);
+
bool isfake = true;
for (size_t isim = 0; isim < simidx.size(); ++isim) {
if (simidxfrac[isim] > matchfrac) {
@@ -2195,6 +2209,7 @@ void setTrackCandidateBranches(LSTEvent* event,
std::vector simidxfrac; // list of match fraction for each matched sim idx
// The following function reads off and computes the matched sim track indices
+ float percent_matched;
std::tie(type, pt, eta, phi, isFake, simidx, simidxfrac) = parseTrackCandidateAllMatch(event,
tc_idx,
trk_ph2_x,
@@ -2203,8 +2218,32 @@ void setTrackCandidateBranches(LSTEvent* event,
trk_simhit_simTrkIdx,
trk_ph2_simHitIdx,
trk_pix_simHitIdx,
+ percent_matched,
matchfrac);
+ int nPixHits = 0, nOtHits = 0, nLayers = 0;
+ for (int layerSlot = 0; layerSlot < Params_TC::kLayers; ++layerSlot) {
+ if (trackCandidatesExtended.lowerModuleIndices()[tc_idx][layerSlot] == lst::kTCEmptyLowerModule)
+ continue;
+
+ ++nLayers;
+ const bool isPixel = (trackCandidatesExtended.logicalLayers()[tc_idx][layerSlot] == 0);
+
+ for (unsigned int hitSlot = 0; hitSlot < Params_TC::kHitsPerLayer; ++hitSlot) {
+ if (trackCandidatesBase.hitIndices()[tc_idx][layerSlot][hitSlot] == lst::kTCEmptyHitIdx)
+ continue;
+
+ if (isPixel)
+ nPixHits++;
+ else
+ nOtHits++;
+ }
+ }
+
+ ana.tx->pushbackToBranch("tc_nhitOT", nOtHits);
+ ana.tx->pushbackToBranch("tc_nhits", nPixHits + nOtHits);
+ ana.tx->pushbackToBranch("tc_nlayers", nLayers);
+
// Fill some branches for this track candidate
ana.tx->pushbackToBranch("tc_pt", pt);
ana.tx->pushbackToBranch("tc_eta", eta);
@@ -2278,6 +2317,7 @@ void setTrackCandidateBranches(LSTEvent* event,
}
ana.tx->pushbackToBranch("tc_isFake", isFake);
+ ana.tx->pushbackToBranch("tc_pMatched", percent_matched);
// For this tc, keep track of all the simidx that are matched
tc_simIdxAll.push_back(simidx);
@@ -2838,6 +2878,10 @@ std::tuple> parseTrackCandidate(
break;
}
+ if (type == LSTObjType::T5 || type == LSTObjType::pT5) {
+ std::tie(hit_idx, hit_type) = getHitIdxsAndTypesFromTC(event, idx);
+ }
+
// Perform matching
std::vector simidx = matchedSimTrkIdxs(
hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx, false, matchfrac);
@@ -2856,6 +2900,7 @@ std::tuple, std::vector>
std::vector const& trk_simhit_simTrkIdx,
std::vector> const& trk_ph2_simHitIdx,
std::vector> const& trk_pix_simHitIdx,
+ float& percent_matched,
float matchfrac) {
// Get the type of the track candidate
auto const& trackCandidatesBase = event->getTrackCandidatesBase();
@@ -2882,11 +2927,15 @@ std::tuple, std::vector>
break;
}
+ if (type == LSTObjType::T5 || type == LSTObjType::pT5) {
+ std::tie(hit_idx, hit_type) = getHitIdxsAndTypesFromTC(event, idx);
+ }
+
// Perform matching
std::vector simidx;
std::vector simidxfrac;
std::tie(simidx, simidxfrac) = matchedSimTrkIdxsAndFracs(
- hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx, false, matchfrac);
+ hit_idx, hit_type, trk_simhit_simTrkIdx, trk_ph2_simHitIdx, trk_pix_simHitIdx, false, matchfrac, &percent_matched);
int isFake = simidx.size() == 0;
return {type, pt, eta, phi, isFake, simidx, simidxfrac};
diff --git a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h
index cb47f97ba7764..ee3786782200c 100644
--- a/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h
+++ b/RecoTracker/LSTCore/standalone/code/core/write_lst_ntuple.h
@@ -101,6 +101,7 @@ std::tuple, std::vector>
std::vector const& trk_simhit_simTrkIdx,
std::vector> const& trk_ph2_simHitIdx,
std::vector> const& trk_pix_simHitIdx,
+ float& percent_matched,
float matchfrac = 0.75);
std::tuple, std::vector> parsepT5(LSTEvent* event,
unsigned int);
diff --git a/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py b/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py
index 8faeb1d8dd978..d7b6aee0da42d 100755
--- a/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py
+++ b/RecoTracker/LSTCore/standalone/efficiency/python/lst_plot_performance.py
@@ -8,7 +8,7 @@
from math import sqrt
sel_choices = ["base", "loweta", "xtr", "vtr", "none"]
-metric_choices = ["eff", "fakerate", "duplrate", "fakeorduplrate"]
+metric_choices = ["eff", "fakerate", "duplrate", "fakeorduplrate", "avgOTlen"]
variable_choices = ["pt", "ptmtv", "ptlow", "eta", "phi", "dxy", "dz", "vxy", "deltaEta", "deltaPhi", "deltaR", "jet_eta", "jet_phi", "jet_pt"]
objecttype_choices = ["TC", "pT5", "T5", "pT3", "pLS", "T4", "pT5_lower", "pT3_lower", "T5_lower", "pLS_lower"]
#lowerObjectType = ["pT5_lower", "pT3_lower", "T5_lower"]
@@ -82,6 +82,10 @@ def plot(args):
params = process_arguments_into_params(args)
+ if params["metric"] == "avgOTlen":
+ params["objecttype"] = "TC"
+ params["breakdown"] = False
+
if params["metric"] == "eff":
params["output_name"] = "{objecttype}_{selection}_{pdgid}_{charge}_{metric}_{variable}".format(**params)
else:
@@ -123,7 +127,12 @@ def plot(args):
# Numerator histograms
numer = []
- numer.append(params["input_file"].Get(params["numer"]).Clone())
+ if params["metric"] == "avgOTlen":
+ h_num = params["input_file"].Get(params["numer"]).Clone()
+ h_numSq = params["input_file"].Get(params["numer"].replace("numer", "numerSq")).Clone()
+ numer.append((h_num, h_numSq))
+ else:
+ numer.append(params["input_file"].Get(params["numer"]).Clone())
breakdown_hist_types = ["pT5", "pT3", "T5", "pLS", "T4"]
print("breakdown = ", params["breakdown"])
@@ -136,8 +145,13 @@ def plot(args):
if params["compare"]:
for f in params["additional_input_files"]:
- hist = f.Get(params["numer"])
- numer.append(hist.Clone())
+ if params["metric"] == "avgOTlen":
+ h_num = f.Get(params["numer"]).Clone()
+ h_numSq = f.Get(params["numer"].replace("numer", "numerSq")).Clone()
+ numer.append((h_num, h_numSq))
+ else:
+ hist = f.Get(params["numer"])
+ numer.append(hist.Clone())
hist = f.Get(params["denom"])
denom.append(hist.Clone())
@@ -243,6 +257,7 @@ def process_arguments_into_params(args):
if params["metric"] == "duplrate": params["metricsuffix"] = "dr"
if params["metric"] == "fakerate": params["metricsuffix"] = "fr"
if params["metric"] == "fakeorduplrate": params["metricsuffix"] = "fdr"
+ if params["metric"] == "avgOTlen": params["metricsuffix"] = "ol"
# If breakdown it must be object type of TC
params["breakdown"] = args.breakdown
@@ -302,24 +317,34 @@ def draw_ratio(nums, dens, params):
# Rebin if necessary
if "scalar" in params["output_name"] and "ptscalar" not in params["output_name"]:
for num in nums:
- num.Rebin(180)
+ if isinstance(num, tuple):
+ num[0].Rebin(180)
+ num[1].Rebin(180)
+ else:
+ num.Rebin(180)
for den in dens:
den.Rebin(180)
# Rebin if necessary
if "coarse" in params["output_name"] and "ptcoarse" not in params["output_name"]:
for num in nums:
- num.Rebin(6)
+ if isinstance(num, tuple):
+ num[0].Rebin(6)
+ num[1].Rebin(6)
+ else:
+ num.Rebin(6)
for den in dens:
den.Rebin(6)
# Deal with overflow bins for pt plots
if "pt" in params["output_name"] or "vxy" in params["output_name"]:
for num in nums:
- overFlowBin = num.GetBinContent(num.GetNbinsX() + 1)
- lastBin = num.GetBinContent(num.GetNbinsX())
- num.SetBinContent(num.GetNbinsX(), lastBin + overFlowBin)
- num.SetBinError(num.GetNbinsX(), sqrt(lastBin + overFlowBin))
+ h_list = num if isinstance(num, tuple) else [num]
+ for h in h_list:
+ overFlowBin = h.GetBinContent(h.GetNbinsX() + 1)
+ lastBin = h.GetBinContent(h.GetNbinsX())
+ h.SetBinContent(h.GetNbinsX(), lastBin + overFlowBin)
+ h.SetBinError(h.GetNbinsX(), sqrt(lastBin + overFlowBin))
for den in dens:
overFlowBin = den.GetBinContent(den.GetNbinsX() + 1)
lastBin = den.GetBinContent(den.GetNbinsX())
@@ -329,15 +354,39 @@ def draw_ratio(nums, dens, params):
# Create efficiency graphs
teffs = []
effs = []
- for num, den in zip(nums, dens):
- teff = r.TEfficiency(num, den)
- eff = teff.CreateGraph()
- teffs.append(teff)
- effs.append(eff)
-
- # print(effs)
+ if params["metric"] == "avgOTlen":
+ for (num, num2), den in zip(nums, dens):
+ g = r.TGraphErrors()
+ ip = 0
+ nb = den.GetNbinsX()
+ for ib in range(1, nb + 1):
+ cnt = den.GetBinContent(ib)
+ S = num.GetBinContent(ib)
+ S2 = num2.GetBinContent(ib)
+ if cnt > 0:
+ x = den.GetXaxis().GetBinCenter(ib)
+ ex = den.GetXaxis().GetBinWidth(ib) / 2.0
+ mean = S / cnt
+ y = mean
+ mean2 = S2 / cnt
+ variance = mean2 - (mean * mean)
+ if variance < 0: variance = 0.0
+ stddev = sqrt(variance)
+ ey = stddev / sqrt(cnt)
+ g.SetPoint(ip, x, y)
+ g.SetPointError(ip, ex, ey)
+ ip += 1
+ if g.GetN() == 0:
+ g.SetPoint(0, 0.0, 0.0)
+ g.SetPointError(0, 0.0, 0.0)
+ effs.append(g)
+ else:
+ for num, den in zip(nums, dens):
+ teff = r.TEfficiency(num, den)
+ eff = teff.CreateGraph()
+ teffs.append(teff)
+ effs.append(eff)
- #
hist_name_suffix = ""
if params["xcoarse"]:
hist_name_suffix += "coarse"
@@ -348,14 +397,12 @@ def draw_ratio(nums, dens, params):
outputname = params["output_name"]
for den in dens:
den.Write(den.GetName() + hist_name_suffix, r.TObject.kOverwrite)
- # eff_den = r.TGraphAsymmErrors(den)
- # eff_den.SetName(outputname+"_den")
- # eff_den.Write("", r.TObject.kOverwrite)
for num in nums:
- num.Write(num.GetName() + hist_name_suffix, r.TObject.kOverwrite)
- # eff_num = r.TGraphAsymmErrors(num)
- # eff_num.SetName(outputname+"_num")
- # eff_num.Write("", r.TObject.kOverwrite)
+ if isinstance(num, tuple):
+ num[0].Write(num[0].GetName() + hist_name_suffix, r.TObject.kOverwrite)
+ num[1].Write(num[1].GetName() + hist_name_suffix, r.TObject.kOverwrite)
+ else:
+ num.Write(num.GetName() + hist_name_suffix, r.TObject.kOverwrite)
for eff in effs:
eff.SetName(outputname)
eff.Write("", r.TObject.kOverwrite)
@@ -369,6 +416,8 @@ def parse_plot_name(output_name):
rtnstr = ["Fake-or-duplicate Rate of"]
elif "fake" in output_name:
rtnstr = ["Fake Rate of"]
+ elif "avgOTlen" in output_name:
+ rtnstr = ["Average OT length of"]
elif "dup" in output_name:
rtnstr = ["Duplicate Rate of"]
elif "inefficiency" in output_name:
@@ -452,6 +501,8 @@ def set_label(eff, output_name, raw_number):
eff.GetXaxis().SetTitle(title)
if "fakeorduplrate" in output_name:
eff.GetYaxis().SetTitle("Fake-or-duplicate Rate")
+ elif "avgOTlen" in output_name:
+ eff.GetYaxis().SetTitle("")
elif "fakerate" in output_name:
eff.GetYaxis().SetTitle("Fake Rate")
elif "duplrate" in output_name:
@@ -518,7 +569,7 @@ def draw_label(params):
particleselection = ((", Particle:" + pdgidstr) if pdgidstr else "" ) + ((", Charge:" + chargestr) if chargestr else "" )
fiducial_label += particleselection
# If fake rate or duplicate rate plot follow the following fiducial label rule
- elif "fakerate" in output_name or "duplrate" in output_name or "fakeorduplrate" in output_name:
+ elif "fakerate" in output_name or "duplrate" in output_name or "fakeorduplrate" in output_name or "avgOTlen" in output_name:
if "_pt" in output_name:
fiducial_label = "|#eta| < {eta}".format(eta=etacut)
elif "_eta" in output_name:
@@ -595,18 +646,22 @@ def draw_plot(effs, nums, dens, params):
yaxis_min = effs[0].GetY()[i]
# Set Yaxis range
- effs[0].GetYaxis().SetRangeUser(0, 1.02)
- if "zoom" not in output_name:
- effs[0].GetYaxis().SetRangeUser(0, 1.02)
- else:
- if "fakeorduplrate" in output_name:
- effs[0].GetYaxis().SetRangeUser(0.0, yaxis_max * 1.1)
- elif "fakerate" in output_name:
- effs[0].GetYaxis().SetRangeUser(0.0, yaxis_max * 1.1)
- elif "duplrate" in output_name:
+ if "avgOTlen" in output_name:
+ if "zoom" not in output_name:
effs[0].GetYaxis().SetRangeUser(0.0, yaxis_max * 1.1)
else:
- effs[0].GetYaxis().SetRangeUser(0.6, 1.02)
+ effs[0].GetYaxis().SetRangeUser(8.0, 12.0)
+ else:
+ effs[0].GetYaxis().SetRangeUser(0, 1.02)
+ if "zoom" not in output_name:
+ effs[0].GetYaxis().SetRangeUser(0, 1.02)
+ else:
+ if "fakerate" in output_name or "fakeorduplrate" in output_name:
+ effs[0].GetYaxis().SetRangeUser(0.0, yaxis_max * 1.1)
+ elif "duplrate" in output_name:
+ effs[0].GetYaxis().SetRangeUser(0.0, yaxis_max * 1.1)
+ else:
+ effs[0].GetYaxis().SetRangeUser(0.6, 1.02)
# Set xaxis range
if "_eta" in output_name:
@@ -625,8 +680,10 @@ def draw_plot(effs, nums, dens, params):
effs[0].SetName(output_name)
for i, num in enumerate(nums):
- set_label(num, output_name, raw_number=True)
- num.Draw("hist")
+ # Extract the Histogram if it's a tuple (sum, sumSq)
+ h_to_draw = num[0] if isinstance(num, tuple) else num
+ set_label(h_to_draw, output_name, raw_number=True)
+ h_to_draw.Draw("hist")
c1.SaveAs("{}".format(output_fullpath.replace("/mtv/var/", "/mtv/num/").replace(".pdf", "_num{}.pdf".format(i))))
c1.SaveAs("{}".format(output_fullpath.replace("/mtv/var/", "/mtv/num/").replace(".pdf", "_num{}.png".format(i))))
@@ -639,11 +696,13 @@ def draw_plot(effs, nums, dens, params):
# Double ratio if more than one nums are provided
# Take the first num as the base
if len(nums) > 1:
- base = nums[0].Clone()
- base.Divide(nums[0], dens[0], 1, 1, "B") #Binomial
+ base_obj = nums[0][0] if isinstance(nums[0], tuple) else nums[0]
+ base = base_obj.Clone()
+ base.Divide(base, dens[0], 1, 1, "B") #Binomial
others = []
for num, den in zip(nums[1:], dens[1:]):
- other = num.Clone()
+ other_obj = num[0] if isinstance(num, tuple) else num
+ other = other_obj.Clone()
other.Divide(other, den, 1, 1, "B")
others.append(other)
@@ -674,6 +733,7 @@ def plot_standard_performance_plots(args):
"fakerate": ["pt", "ptlow", "ptmtv", "eta", "phi"],
"duplrate": ["pt", "ptlow", "ptmtv", "eta", "phi"],
"fakeorduplrate": ["pt", "ptlow", "ptmtv", "eta", "phi"],
+ "avgOTlen": ["eta"],
}
if (args.jet_branches): variables["eff"] = ["pt", "ptlow", "ptmtv", "eta", "phi", "dxy", "dz", "vxy", "deltaEta", "deltaPhi", "deltaR", "jet_eta", "jet_phi", "jet_pt"]
sels = {
@@ -681,6 +741,7 @@ def plot_standard_performance_plots(args):
"fakerate": ["none"],
"duplrate": ["none"],
"fakeorduplrate": ["none"],
+ "avgOTlen": ["none"],
}
xcoarses = {
"pt": [False],
@@ -750,18 +811,31 @@ def plot_standard_performance_plots(args):
"T5_lower":[False],
"pLS_lower":[False],
},
+ "avgOTlen":{
+ "TC":[False],
+ "pT5":[False],
+ "pT3":[False],
+ "T5":[False],
+ "pLS":[False],
+ "pT5_lower":[False],
+ "pT3_lower":[False],
+ "T5_lower":[False],
+ "pLS_lower":[False],
+ },
}
pdgids = {
"eff": [0, 11, 13, 211, 321],
"fakerate": [0],
"duplrate": [0],
"fakeorduplrate": [0],
+ "avgOTlen": [0],
}
charges = {
"eff":[0, 1, -1],
"fakerate":[0],
"duplrate":[0],
"fakeorduplrate":[0],
+ "avgOTlen":[0],
}
if args.metric:
@@ -806,15 +880,15 @@ def plot_standard_performance_plots(args):
breakdowns = {"eff":{"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]},
"fakerate": {"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]},
"duplrate": {"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]},
- "fakeorduplrate": {"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]}}
-
-
+ "fakeorduplrate": {"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]},
+ "avgOTlen":{"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]}}
else:
# Only eff / TC matters here
breakdowns = {"eff":{"TC":[True], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]},
"fakerate": {"TC":[True], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]},
"duplrate": {"TC":[True], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]},
- "fakeorduplrate": {"TC":[True], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]}}
+ "fakeorduplrate": {"TC":[True], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]},
+ "avgOTlen":{"TC":[False], "pT5_lower":[False], "pT3_lower":[False], "T5_lower":[False], "pLS_lower":[False]}}
if args.yzoom:
args.yzooms = [args.yzoom]
@@ -826,7 +900,7 @@ def plot_standard_performance_plots(args):
for variable in variables[metric]:
for yzoom in yzooms:
for xcoarse in xcoarses[variable]:
- for typ in types:
+ for typ in (["TC"] if metric == "avgOTlen" else types):
print("type = ",typ)
for breakdown in breakdowns[metric][typ]:
for pdgid in pdgids[metric]:
diff --git a/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc b/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc
index 56bac3b6de395..af1ab9d3f9edc 100644
--- a/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc
+++ b/RecoTracker/LSTCore/standalone/efficiency/src/performance.cc
@@ -176,6 +176,20 @@ int main(int argc, char** argv) {
}
}
+ // ---------------- OT track length sets ----------------
+ std::vector list_OLSetDef;
+ list_OLSetDef.push_back(RecoTrackSetDefinition(
+ /* name */
+ "TC",
+ /* pass */ [&](unsigned int itc) { return 1; },
+ /* sel */ [&](unsigned int itc) { return 1; },
+ /* pt */ [&]() { return lstEff.getVF("tc_pt"); },
+ /* eta */ [&]() { return lstEff.getVF("tc_eta"); },
+ /* phi */ [&]() { return lstEff.getVF("tc_phi"); },
+ /* type */ [&]() { return lstEff.getVI("tc_type"); }));
+ bookOTLengthSets(list_OLSetDef);
+ // --------------------------------------------------------
+
bookEfficiencySets(list_effSetDef);
// creating a set of fake rate plots
@@ -555,6 +569,7 @@ int main(int argc, char** argv) {
fillFakeRateSets(list_FRSetDef);
fillDuplicateRateSets(list_DRSetDef);
fillFakeOrDuplicateRateSets(list_FDRSetDef);
+ fillOTLengthSets(list_OLSetDef);
// Reset all temporary variables necessary for histogramming
ana.cutflow.fill();
@@ -585,6 +600,73 @@ int main(int argc, char** argv) {
// ---------------------------------------------------------=============================================-------------------------------------------------------------------
// ---------------------------------------------------------=============================================-------------------------------------------------------------------
+// ---------------------------------------------------------=============================================-------------------------------------------------------------------
+void bookOTLengthSets(std::vector& OLsets) {
+ for (auto& OLset : OLsets) {
+ bookOTLengthSet(OLset);
+ }
+}
+
+//__________________________________________________________________________________________________________________________________________________________________________
+void bookOTLengthSet(RecoTrackSetDefinition& OLset) {
+ TString category_name = OLset.set_name;
+
+ ana.tx.createBranch>(category_name + "_ol_denom_eta");
+ ana.tx.createBranch>(category_name + "_ol_numer_eta");
+ ana.tx.createBranch>(category_name + "_ol_numerSq_eta");
+
+ ana.histograms.addVecHistogram(category_name + "_ol_denom_eta", 180, -4.5, 4.5, [&, category_name]() {
+ return ana.tx.getBranchLazy>(category_name + "_ol_denom_eta");
+ });
+ ana.histograms.addVecHistogram(category_name + "_ol_numer_eta", 180, -4.5, 4.5, [&, category_name]() {
+ return ana.tx.getBranchLazy>(category_name + "_ol_numer_eta");
+ });
+ ana.histograms.addVecHistogram(category_name + "_ol_numerSq_eta", 180, -4.5, 4.5, [&, category_name]() {
+ return ana.tx.getBranchLazy>(category_name + "_ol_numerSq_eta");
+ });
+}
+
+//__________________________________________________________________________________________________________________________________________________________________________
+void fillOTLengthSets(std::vector& OLsets) {
+ for (auto& OLset : OLsets) {
+ const std::vector& pt = OLset.pt();
+ const std::vector& eta = OLset.eta();
+ for (unsigned int itc = 0; itc < pt.size(); ++itc) {
+ fillOTLengthSet(static_cast(itc), OLset, pt[itc], eta[itc]);
+ }
+ }
+}
+
+//__________________________________________________________________________________________________________________________________________________________________________
+void fillOTLengthSet(int itc, RecoTrackSetDefinition& OLset, float pt, float eta) {
+ TString category_name = OLset.set_name;
+ const bool sel = OLset.sel(itc); // mirror FR/DR pattern (selection gate)
+ if (!sel)
+ return;
+
+ // Pull OT length from the relevant TC branch
+ const auto& tc_nhitOT = lstEff.getVI("tc_nhitOT");
+ if (itc < 0 || itc >= static_cast(tc_nhitOT.size()))
+ return;
+ const int otlen = tc_nhitOT[itc];
+
+ if (pt > ana.pt_cut) {
+ // Denominator: one entry per TC
+ ana.tx.pushbackToBranch(category_name + "_ol_denom_eta", eta);
+
+ // Numerator: 'otlen' entries at the same eta bin sum, ignore pLSs
+ if (otlen > 0) {
+ for (int k = 0; k < otlen; ++k) {
+ ana.tx.pushbackToBranch(category_name + "_ol_numer_eta", eta);
+ }
+ int otlenSq = otlen * otlen;
+ for (int k = 0; k < otlenSq; ++k) {
+ ana.tx.pushbackToBranch(category_name + "_ol_numerSq_eta", eta);
+ }
+ }
+ }
+}
+
//__________________________________________________________________________________________________________________________________________________________________________
void bookEfficiencySets(std::vector& effsets) {
for (auto& effset : effsets)
diff --git a/RecoTracker/LSTCore/standalone/efficiency/src/performance.h b/RecoTracker/LSTCore/standalone/efficiency/src/performance.h
index a46686014fb83..97e64cb751456 100644
--- a/RecoTracker/LSTCore/standalone/efficiency/src/performance.h
+++ b/RecoTracker/LSTCore/standalone/efficiency/src/performance.h
@@ -14,6 +14,8 @@ void bookDuplicateRateSets(std::vector& DRset);
void bookDuplicateRateSet(RecoTrackSetDefinition& DRset);
void bookFakeOrDuplicateRateSets(std::vector& FDRset);
void bookFakeOrDuplicateRateSet(RecoTrackSetDefinition& FDRset);
+void bookOTLengthSets(std::vector& OLsets);
+void bookOTLengthSet(RecoTrackSetDefinition& OLset);
void fillEfficiencySets(std::vector& effset);
@@ -48,6 +50,8 @@ void fillEfficiencySet(int isimtrk,
float vtx_y,
float vtx_z);
+void fillOTLengthSets(std::vector& OLsets);
+void fillOTLengthSet(int itc, RecoTrackSetDefinition& OLset, float pt, float eta);
void fillFakeRateSets(std::vector& FRset);
void fillFakeRateSet(int isimtrk, RecoTrackSetDefinition& FRset, float pt, float eta, float phi);
void fillDuplicateRateSets(std::vector& DRset);