Skip to content

Commit 3513c51

Browse files
choich08365Changhwan Choi
andauthored
[Tools,PWGJE] Implemented full inference of GNN b-jet tagger (#18682)
Co-authored-by: Changhwan Choi <changhwan.choi@cern.ch>
1 parent 53d7827 commit 3513c51

10 files changed

Lines changed: 1108 additions & 97 deletions

File tree

‎PWGJE/Core/JetTaggingUtilities.h‎

Lines changed: 50 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1130,9 +1130,13 @@ void analyzeJetTrackInfo4GNNwExtra(AnalysisJet const& analysisJet, AnyTracks con
11301130

11311131
// Looping over the track info and putting them in the input vector (for GNN b-jet tagging)
11321132
template <typename AnalysisJet, typename AnyTracks>
1133-
void analyzeJetTrackInfo4GNN(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, std::vector<std::vector<float>>& tracksParams, float trackPtMin = 0.5, float trackDcaXYMax = 10.0, float trackDcaZMax = 10.0, int64_t nMaxConstit = 40)
1133+
void analyzeJetTrackInfo4GNN(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, std::vector<std::vector<float>>& tracksParams, float trackPtMin = 0.5, float trackDcaXYMax = 10.0, float trackDcaZMax = 10.0, int64_t nMaxConstit = 40, std::vector<std::size_t>* constituentIndices = nullptr)
11341134
{
1135-
for (const auto& constituent : analysisJet.template tracks_as<AnyTracks>()) {
1135+
// Optionally retain constituent-array positions in exactly the model-node
1136+
// order, including replacements when the maximum track count is exceeded.
1137+
auto const tracks = analysisJet.template tracks_as<AnyTracks>();
1138+
for (std::size_t constituentIndex = 0; constituentIndex < tracks.size(); ++constituentIndex) {
1139+
auto const& constituent = tracks[constituentIndex];
11361140

11371141
if (constituent.pt() < trackPtMin || !trackAcceptanceWithDca(constituent, trackDcaXYMax, trackDcaZMax)) {
11381142
continue;
@@ -1142,15 +1146,22 @@ void analyzeJetTrackInfo4GNN(AnalysisJet const& analysisJet, AnyTracks const& /*
11421146

11431147
if (static_cast<int64_t>(tracksParams.size()) < nMaxConstit) {
11441148
tracksParams.emplace_back(std::vector<float>{constituent.pt(), constituent.phi(), constituent.eta(), static_cast<float>(constituent.sign()), std::abs(constituent.dcaXY()) * sign, constituent.sigmadcaXY(), std::abs(constituent.dcaZ()) * sign, constituent.sigmadcaZ()});
1149+
if (constituentIndices) {
1150+
constituentIndices->push_back(constituentIndex);
1151+
}
11451152
} else {
11461153
// If there are more than nMaxConstit constituents in the jet, select only nMaxConstit constituents with the highest DCA_XY significance.
11471154
size_t minIdx = 0;
11481155
for (size_t i = 0; i < tracksParams.size(); ++i) {
11491156
if (tracksParams[i][4] / tracksParams[i][5] < tracksParams[minIdx][4] / tracksParams[minIdx][5])
11501157
minIdx = i;
11511158
}
1152-
if (std::abs(constituent.dcaXY()) * sign / constituent.sigmadcaXY() > tracksParams[minIdx][4] / tracksParams[minIdx][5])
1159+
if (std::abs(constituent.dcaXY()) * sign / constituent.sigmadcaXY() > tracksParams[minIdx][4] / tracksParams[minIdx][5]) {
11531160
tracksParams[minIdx] = std::vector<float>{constituent.pt(), constituent.phi(), constituent.eta(), static_cast<float>(constituent.sign()), std::abs(constituent.dcaXY()) * sign, constituent.sigmadcaXY(), std::abs(constituent.dcaZ()) * sign, constituent.sigmadcaZ()};
1161+
if (constituentIndices) {
1162+
(*constituentIndices)[minIdx] = constituentIndex;
1163+
}
1164+
}
11541165
}
11551166
}
11561167
}
@@ -1218,6 +1229,9 @@ T getDbFromProb(const std::vector<T>& probs, double fC = 0.018)
12181229
template <typename T>
12191230
T getDb(const std::vector<T>& logits, double fC = 0.018)
12201231
{
1232+
if (logits.size() < 3) {
1233+
LOG(fatal) << "getDb expects at least 3 class logits, received " << logits.size();
1234+
}
12211235
auto softmax = [](const std::vector<T>& logits) {
12221236
std::vector<T> res;
12231237
T maxLogit = *std::max_element(logits.begin(), logits.end());
@@ -1234,11 +1248,41 @@ T getDb(const std::vector<T>& logits, double fC = 0.018)
12341248

12351249
std::vector<T> softmaxLogits = softmax(logits);
12361250

1237-
if (softmaxLogits[1] == 0. && softmaxLogits[2] == 0.) {
1238-
LOG(debug) << "jettaggingutilities::Db, Divide by zero: softmaxLogits = (" << softmaxLogits[0] << ", " << softmaxLogits[1] << ", " << softmaxLogits[2] << ")";
1251+
auto denominator = fC * softmaxLogits[1] + (1. - fC) * softmaxLogits[2];
1252+
if (!std::isfinite(softmaxLogits[0]) || !std::isfinite(denominator) || softmaxLogits[0] <= 0. || denominator <= 0.) {
1253+
LOG(debug) << "jettaggingutilities::getDb, non-positive or non-finite numerator/denominator";
1254+
return std::numeric_limits<T>::quiet_NaN();
1255+
}
1256+
1257+
return std::log(softmaxLogits[0] / denominator);
1258+
}
1259+
1260+
// Track-origin discriminant. Logit order: [background, charm, beauty, PV, other SV].
1261+
template <typename T>
1262+
T getTrackDb(const std::vector<T>& trackLogits, double trackFC)
1263+
{
1264+
constexpr std::size_t NTrackOrigins = 5;
1265+
if (trackLogits.size() != NTrackOrigins) {
1266+
LOG(fatal) << "getTrackDb expects 5 track-origin logits, received " << trackLogits.size();
12391267
}
12401268

1241-
return std::log(softmaxLogits[0] / (fC * softmaxLogits[1] + (1. - fC) * softmaxLogits[2]));
1269+
std::array<T, NTrackOrigins> softmaxLogits;
1270+
T maxLogit = *std::max_element(trackLogits.begin(), trackLogits.end());
1271+
T sumLogit = 0.;
1272+
for (std::size_t i = 0; i < NTrackOrigins; ++i) {
1273+
softmaxLogits[i] = std::exp(trackLogits[i] - maxLogit);
1274+
sumLogit += softmaxLogits[i];
1275+
}
1276+
for (auto& probability : softmaxLogits) {
1277+
probability /= sumLogit;
1278+
}
1279+
1280+
T denominator = trackFC * softmaxLogits[1] + (1. - trackFC) * (softmaxLogits[0] + softmaxLogits[3] + softmaxLogits[4]);
1281+
if (!std::isfinite(softmaxLogits[2]) || !std::isfinite(denominator) || softmaxLogits[2] <= 0. || denominator <= 0.) {
1282+
LOG(debug) << "jettaggingutilities::getTrackDb, non-positive or non-finite numerator/denominator";
1283+
return std::numeric_limits<T>::quiet_NaN();
1284+
}
1285+
return std::log(softmaxLogits[2] / denominator);
12421286
}
12431287

12441288
}; // namespace jettaggingutilities

‎PWGJE/DataModel/JetTagging.h‎

Lines changed: 106 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,8 @@
3030
#include <array>
3131
#include <cmath>
3232
#include <cstdint>
33+
#include <utility>
34+
#include <vector>
3335

3436
namespace o2::aod
3537
{
@@ -165,6 +167,110 @@ DECLARE_SOA_DYNAMIC_COLUMN(IsTagged, isTagged, [](uint16_t bit, BJetTaggingMetho
165167
JETTAGGING_TABLES_DEF(Charged, "C");
166168
JETTAGGING_TABLES_DEF(Full, "F");
167169

170+
// Tables for full inference of GNN-based b-jet tagger (including track origin, vertex predictions)
171+
namespace jetgnnindex
172+
{
173+
DECLARE_SOA_INDEX_COLUMN(JTrack, track); //! index in JTracks of this GNN input track (not the parent AOD Tracks index)
174+
DECLARE_SOA_COLUMN(TrackIndex, trackIndex, uint8_t); //! node index (0..255) in the GNN input after track selection/truncation; need not match the constituent-array position
175+
// DECLARE_SOA_COLUMN(EdgeIndex, edgeIndex, uint8_t); //! unique in-jet edge index for the (i,j) constituent-track pair, i<j, triangular order i*(2N-i-1)/2+(j-i-1)
176+
// DECLARE_SOA_DYNAMIC_COLUMN(TrackIds, trackIds, //! inverse of the triangular edge indexing above: for n constituent tracks, recovers the (i,j) in-jet track index pair (i<j) this edge connects
177+
// [](uint8_t edgeIndex, int n) -> std::pair<uint8_t, uint8_t> {
178+
// int i = 0;
179+
// int offset = 0;
180+
// while (offset + (n - 1 - i) <= static_cast<int>(edgeIndex)) {
181+
// offset += (n - 1 - i);
182+
// i++;
183+
// }
184+
// int j = i + 1 + (static_cast<int>(edgeIndex) - offset);
185+
// return std::pair<uint8_t, uint8_t>{static_cast<uint8_t>(i), static_cast<uint8_t>(j)};
186+
// });
187+
} // namespace jetgnnindex
188+
189+
namespace jetgnnlabel
190+
{
191+
DECLARE_SOA_COLUMN(Origin, origin, uint8_t); //! trkOrigin from vertexClustering() (0: non-physical primary, 1: charm, 2: beauty, 3: primary vertex, 4: other secondary)
192+
DECLARE_SOA_COLUMN(VertexIndex, vertexIndex, int8_t); //! trkVtxIndex from vertexClustering() (-1: unassigned, 0: primary vertex, >0: secondary vertex, ranked by distance to PV)
193+
// DECLARE_SOA_COLUMN(SameVertex, sameVertex, bool); //! true if the (i,j) constituent-track pair of this edge share the same in-jet vertex cluster (same trkVtxIndex from vertexClustering())
194+
} // namespace jetgnnlabel
195+
196+
namespace jetgnnscore
197+
{
198+
DECLARE_SOA_COLUMN(ScoreML, scoreML, float); //! GNN discriminant: one track-origin D_b per GnnTracks row
199+
} // namespace jetgnnscore
200+
201+
// GnnTracks has one row per evaluated (jet, track) pair. Its jet and track indices
202+
// identify the pair independently of the ordering in JetConstituents.
203+
// GnnTrackIds is a positional jet extension: write one row for every jet in
204+
// Jets order, including an empty array for jets without evaluated tracks. Store
205+
// GnnTracks row indices (e.g. the backbone producer's lastIndex()), not JTrackIds.
206+
#define DECLARE_JET_GNNTRACKINDEX_TABLE(_jet_type_, _description_) \
207+
namespace _jet_type_##gnntrackindextable \
208+
{ \
209+
DECLARE_SOA_INDEX_COLUMN(_jet_type_, jet); \
210+
} \
211+
DECLARE_SOA_TABLE(_jet_type_##GnnTracks, "AOD", _description_, \
212+
o2::soa::Index<>, \
213+
_jet_type_##gnntrackindextable::_jet_type_##Id, \
214+
jetgnnindex::JTrackId, \
215+
jetgnnindex::TrackIndex); \
216+
namespace _jet_type_##gnntrackindices \
217+
{ \
218+
DECLARE_SOA_ARRAY_INDEX_COLUMN_FULL(_jet_type_##GnnTrack, tracksGnn, int32_t, _jet_type_##GnnTracks, ""); \
219+
} \
220+
DECLARE_SOA_TABLE(_jet_type_##GnnTrackIds, "AOD", _description_ "I", \
221+
_jet_type_##gnntrackindices::_jet_type_##GnnTrackIds);
222+
223+
/*
224+
#define DECLARE_JET_GNNEDGEINDEX_TABLE(_jet_type_, _description_) \
225+
namespace _jet_type_##gnnedgeindextable \
226+
{ \
227+
DECLARE_SOA_INDEX_COLUMN(_jet_type_, jet); \
228+
} \
229+
DECLARE_SOA_TABLE(_jet_type_##GnnEdges, "AOD", _description_, \
230+
o2::soa::Index<>, \
231+
_jet_type_##gnnedgeindextable::_jet_type_##Id, \
232+
jetgnnindex::EdgeIndex, \
233+
jetgnnindex::TrackIds<jetgnnindex::EdgeIndex>); \
234+
namespace _jet_type_##gnnedgeindices \
235+
{ \
236+
DECLARE_SOA_ARRAY_INDEX_COLUMN(_jet_type_##GnnEdge, edges); \
237+
} \
238+
DECLARE_SOA_TABLE(_jet_type_##GnnEdgeIds, "AOD", _description_ "I", \
239+
_jet_type_##gnnedgeindices::_jet_type_##GnnEdgeIds);
240+
*/
241+
242+
// In-jet track/edge index backbone tables (own the in-jet ordering + the back-reference to the original track)
243+
DECLARE_JET_GNNTRACKINDEX_TABLE(ChargedJet, "CJGNT");
244+
DECLARE_JET_GNNTRACKINDEX_TABLE(ChargedMCDetectorLevelJet, "CDJGNT");
245+
// DECLARE_JET_GNNEDGEINDEX_TABLE(ChargedJet, "CJGNE");
246+
// DECLARE_JET_GNNEDGEINDEX_TABLE(ChargedMCDetectorLevelJet, "CDJGNE");
247+
248+
#undef DECLARE_JET_GNNTRACKINDEX_TABLE
249+
// #undef DECLARE_JET_GNNEDGEINDEX_TABLE
250+
251+
// Truth labels from vertexClustering() (MC detector level only). Thin tables: join positionally with
252+
// ChargedMCDetectorLevelJetGnnTracks / ChargedMCDetectorLevelJetGnnEdges (same producer loop, same row order).
253+
// TrackOriginDef must have exactly one row per GnnTracks backbone row, aligned
254+
// with GnnTrackOriginTags after GNN input selection/truncation. It is not a JTracks
255+
// extension: a JTrack used by multiple jets has a truth row for each backbone row.
256+
// tracksGnnIds() addresses this shared row position through a join containing
257+
// the backbone, the origin definition and the prediction scores.
258+
// VertexIndexDef uses the same GnnTracks row order and keeps the full-jet vertex
259+
// labels, without renumbering after GNN track selection/truncation.
260+
DECLARE_SOA_TABLE(ChargedMCDetectorLevelJetTrackOriginDef, "AOD", "CDJTOD", jetgnnlabel::Origin);
261+
DECLARE_SOA_TABLE(ChargedMCDetectorLevelJetVertexIndexDef, "AOD", "CDJTVD", jetgnnlabel::VertexIndex);
262+
// DECLARE_SOA_TABLE(ChargedMCDetectorLevelJetEdgeDef, "AOD", "CDJEDF", jetgnnlabel::SameVertex);
263+
264+
// GNN track-origin outputs, data and MC detector level. Thin tables: write one
265+
// scalar score per GnnTracks row, in the same row order as the backbone table.
266+
DECLARE_SOA_TABLE(ChargedJetGnnTrackOriginTags, "AOD", "CJTRKT", jetgnnscore::ScoreML);
267+
DECLARE_SOA_TABLE(ChargedMCDetectorLevelJetGnnTrackOriginTags, "AOD", "CDJTRKT", jetgnnscore::ScoreML);
268+
269+
// GNN ScoreML outputs, edge level, data and MC detector level. Thin tables: join positionally with
270+
// ChargedJetGnnEdges / ChargedMCDetectorLevelJetGnnEdges.
271+
// DECLARE_SOA_TABLE(ChargedJetGnnEdgeTags, "AOD", "CJEDGT", jetgnnscore::ScoreML);
272+
// DECLARE_SOA_TABLE(ChargedMCDetectorLevelJetGnnEdgeTags, "AOD", "CDJEDGT", jetgnnscore::ScoreML);
273+
168274
} // namespace o2::aod
169275

170276
#endif // PWGJE_DATAMODEL_JETTAGGING_H_

‎PWGJE/TableProducer/CMakeLists.txt‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -86,6 +86,16 @@ o2physics_add_dpl_workflow(jet-taggerhf
8686
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore O2Physics::MLCore
8787
COMPONENT_NAME Analysis)
8888

89+
o2physics_add_dpl_workflow(jet-gnn-definition
90+
SOURCES jetGnnDefinition.cxx
91+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
92+
COMPONENT_NAME Analysis)
93+
94+
o2physics_add_dpl_workflow(jet-gnn-inference
95+
SOURCES jetGnnInference.cxx
96+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore O2Physics::MLCore
97+
COMPONENT_NAME Analysis)
98+
8999
o2physics_add_dpl_workflow(estimator-rho
90100
SOURCES rhoEstimator.cxx
91101
REUSE_FROM JetDataModelPCH

0 commit comments

Comments
 (0)