Skip to content

Commit cc75848

Browse files
choich08365Changhwan Choi
choich08365
and
Changhwan Choi
authored
[PWGJE] Added THnSparse::Sumw2() call option for error calculation, Fixed bugs. (#10223)
Co-authored-by: Changhwan Choi <changhwan.choi@cern.ch>
1 parent 4bb0979 commit cc75848

File tree

2 files changed

+42
-59
lines changed

2 files changed

+42
-59
lines changed

‎PWGJE/Tasks/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -262,7 +262,7 @@ if(FastJet_FOUND)
262262
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
263263
COMPONENT_NAME Analysis)
264264
o2physics_add_dpl_workflow(bjet-tagging-gnn
265-
SOURCES bjetTaggingGNN.cxx
265+
SOURCES bjetTaggingGnn.cxx
266266
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
267267
COMPONENT_NAME Analysis)
268268
o2physics_add_dpl_workflow(jet-shape

PWGJE/Tasks/bjetTaggingGNN.cxx renamed to PWGJE/Tasks/bjetTaggingGnn.cxx

Lines changed: 41 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@
99
// granted to it by virtue of its status as an Intergovernmental Organization
1010
// or submit itself to any jurisdiction.
1111

12-
/// \file bjetTaggingGNN.cxx
12+
/// \file bjetTaggingGnn.cxx
1313
/// \brief b-jet tagging using GNN
1414
///
1515
/// \author Changhwan Choi <changhwan.choi@cern.ch>, Pusan National University
@@ -43,7 +43,7 @@ using namespace o2;
4343
using namespace o2::framework;
4444
using namespace o2::framework::expressions;
4545

46-
struct BjetTaggingGNN {
46+
struct BjetTaggingGnn {
4747

4848
HistogramRegistry registry;
4949

@@ -67,13 +67,6 @@ struct BjetTaggingGNN {
6767
// track level configurables
6868
Configurable<float> svPtMin{"svPtMin", 0.5, "minimum SV pT"};
6969

70-
// Configurable<float> prongsigmaLxyMax{"prongsigmaLxyMax", 100, "maximum sigma of decay length of prongs on xy plane"};
71-
// // Configurable<float> prongsigmaLxyzMax{"prongsigmaLxyzMax", 100, "maximum sigma of decay length of prongs on xyz plane"};
72-
// Configurable<float> prongIPxyMin{"prongIPxyMin", 0.008, "minimum impact paramter of prongs on xy plane [cm]"};
73-
// Configurable<float> prongIPxyMax{"prongIPxyMax", 1, "maximum impact parmeter of prongs on xy plane [cm]"};
74-
// Configurable<float> prongChi2PCAMin{"prongChi2PCAMin", 4, "minimum Chi2 PCA of decay length of prongs"};
75-
// Configurable<float> prongChi2PCAMax{"prongChi2PCAMax", 100, "maximum Chi2 PCA of decay length of prongs"};
76-
7770
// jet level configurables
7871
Configurable<float> jetPtMin{"jetPtMin", 5.0, "minimum jet pT"};
7972
Configurable<float> jetPtMax{"jetPtMax", 1000.0, "maximum jet pT"};
@@ -83,6 +76,7 @@ struct BjetTaggingGNN {
8376
Configurable<std::vector<double>> jetRadii{"jetRadii", std::vector<double>{0.4}, "jet resolution parameters"};
8477

8578
Configurable<bool> doDataDriven{"doDataDriven", false, "Flag whether to use fill THnSpase for data driven methods"};
79+
Configurable<bool> callSumw2{"callSumw2", false, "Flag whether to call THnSparse::Sumw2() for error calculation"};
8680

8781
std::vector<int> eventSelectionBits;
8882

@@ -110,15 +104,10 @@ struct BjetTaggingGNN {
110104
registry.add("h_Db", "", {HistType::kTH1F, {axisDbFine}});
111105
registry.add("h2_jetpT_Db", "", {HistType::kTH2F, {axisJetpT, axisDb}});
112106
registry.add("h2_jetpT_SVMass", "", {HistType::kTH2F, {axisJetpT, axisSVMass}});
113-
registry.add("h2_jetpT_SVEnergy", "", {HistType::kTH2F, {axisJetpT, axisSVEnergy}});
114-
registry.add("h2_jetpT_SLxy", "", {HistType::kTH2F, {axisJetpT, axisSLxy}});
115107
registry.add("h2_jetpT_jetMass", "", {HistType::kTH2F, {axisJetpT, axisJetMass}});
116108
registry.add("h2_jetpT_jetProb", "", {HistType::kTH2F, {axisJetpT, axisJetProb}});
117109
registry.add("h2_jetpT_nTracks", "", {HistType::kTH2F, {axisJetpT, axisNTracks}});
118110

119-
registry.add("h3_jetpT_nTracks_Db", "", {HistType::kTH3F, {axisJetpT, axisNTracks, axisDb}});
120-
registry.add("h3_mSV_eSV_slXY", ";m_{SV};E_{SV};SL_{XY}", {HistType::kTH3F, {{110, -2., 20.}, {110, -2., 20.}, {102, -2., 100.}}});
121-
122111
if (doprocessMCJets) {
123112
registry.add("h_jetpT_b", "b-jet", {HistType::kTH1F, {axisJetpT}});
124113
registry.add("h_jetpT_c", "c-jet", {HistType::kTH1F, {axisJetpT}});
@@ -132,12 +121,6 @@ struct BjetTaggingGNN {
132121
registry.add("h2_jetpT_SVMass_b", "b-jet", {HistType::kTH2F, {axisJetpT, axisSVMass}});
133122
registry.add("h2_jetpT_SVMass_c", "c-jet", {HistType::kTH2F, {axisJetpT, axisSVMass}});
134123
registry.add("h2_jetpT_SVMass_lf", "lf-jet", {HistType::kTH2F, {axisJetpT, axisSVMass}});
135-
registry.add("h2_jetpT_SVEnergy_b", "b-jet", {HistType::kTH2F, {axisJetpT, axisSVEnergy}});
136-
registry.add("h2_jetpT_SVEnergy_c", "c-jet", {HistType::kTH2F, {axisJetpT, axisSVEnergy}});
137-
registry.add("h2_jetpT_SVEnergy_lf", "lf-jet", {HistType::kTH2F, {axisJetpT, axisSVEnergy}});
138-
registry.add("h2_jetpT_SLxy_b", "b-jet", {HistType::kTH2F, {axisJetpT, axisSLxy}});
139-
registry.add("h2_jetpT_SLxy_c", "c-jet", {HistType::kTH2F, {axisJetpT, axisSLxy}});
140-
registry.add("h2_jetpT_SLxy_lf", "lf-jet", {HistType::kTH2F, {axisJetpT, axisSLxy}});
141124
registry.add("h2_jetpT_jetMass_b", "b-jet", {HistType::kTH2F, {axisJetpT, axisJetMass}});
142125
registry.add("h2_jetpT_jetMass_c", "c-jet", {HistType::kTH2F, {axisJetpT, axisJetMass}});
143126
registry.add("h2_jetpT_jetMass_lf", "lf-jet", {HistType::kTH2F, {axisJetpT, axisJetMass}});
@@ -161,6 +144,10 @@ struct BjetTaggingGNN {
161144
registry.add("h_Db_npp_b", "NotPhysPrim b-jet", {HistType::kTH1F, {axisDbFine}});
162145
registry.add("h_Db_npp_c", "NotPhysPrim c-jet", {HistType::kTH1F, {axisDbFine}});
163146
registry.add("h_Db_npp_lf", "NotPhysPrim lf-jet", {HistType::kTH1F, {axisDbFine}});
147+
// registry.add("h2_pT_dcaXY_pp", "tracks", {HistType::kTH2F, {axisJetpT, {200, 0., 1.}}});
148+
// registry.add("h2_pT_dcaXY_npp", "NotPhysPrim tracks", {HistType::kTH2F, {axisJetpT, {200, 0., 1.}}});
149+
// registry.add("h2_pT_dcaZ_pp", "tracks", {HistType::kTH2F, {axisJetpT, {200, 0., 2.}}});
150+
// registry.add("h2_pT_dcaZ_npp", "NotPhysPrim tracks", {HistType::kTH2F, {axisJetpT, {200, 0., 2.}}});
164151
}
165152

166153
if (doprocessMCTruthJets) {
@@ -171,13 +158,13 @@ struct BjetTaggingGNN {
171158
}
172159

173160
if (doDataDriven) {
174-
registry.add("hSparse_Incljets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}});
161+
registry.add("hSparse_Incljets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2);
175162
if (doprocessMCJets) {
176-
registry.add("hSparse_bjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}});
177-
registry.add("hSparse_cjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}});
178-
registry.add("hSparse_lfjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}});
179-
registry.add("hSparse_lfjets_none", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}});
180-
registry.add("hSparse_lfjets_matched", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}});
163+
registry.add("hSparse_bjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2);
164+
registry.add("hSparse_cjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2);
165+
registry.add("hSparse_lfjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2);
166+
registry.add("hSparse_lfjets_none", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2);
167+
registry.add("hSparse_lfjets_matched", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2);
181168
}
182169
}
183170
}
@@ -197,7 +184,7 @@ struct BjetTaggingGNN {
197184
int analyzeJetTrackInfo(AnyCollision const& /*collision*/, AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/ /*, int8_t jetFlavor = 0, double weight = 1.0*/)
198185
{
199186
int nTracks = 0;
200-
for (auto& constituent : analysisJet.template tracks_as<AnyTracks>()) {
187+
for (const auto& constituent : analysisJet.template tracks_as<AnyTracks>()) {
201188

202189
if (constituent.pt() < trackPtMin) {
203190
continue;
@@ -241,7 +228,7 @@ struct BjetTaggingGNN {
241228
void processDummy(FilteredCollision::iterator const& /*collision*/)
242229
{
243230
}
244-
PROCESS_SWITCH(BjetTaggingGNN, processDummy, "Dummy process function turned on by default", true);
231+
PROCESS_SWITCH(BjetTaggingGnn, processDummy, "Dummy process function turned on by default", true);
245232

246233
void processDataJets(FilteredCollision::iterator const& collision, DataJets const& alljets, JetTrackswID const& allTracks, SVTable const& allSVs)
247234
{
@@ -268,25 +255,23 @@ struct BjetTaggingGNN {
268255
int nTracks = analyzeJetTrackInfo(collision, analysisJet, allTracks);
269256

270257
float mSV = -1.f;
271-
float eSV = -1.f;
272-
float slXY = -1.f;
258+
// float eSV = -1.f;
259+
// float slXY = -1.f;
273260

274261
bool checkSV;
275262
// auto sv = jettaggingutilities::jetFromProngMaxDecayLength<MCDSVTable>(analysisJet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, false, &checkSV);
276263
auto sv = analyzeJetSVInfo(analysisJet, allSVs, checkSV);
277264

278265
if (checkSV) {
279266
mSV = sv.m();
280-
eSV = sv.e();
281-
slXY = sv.decayLengthXY() / sv.errorDecayLengthXY();
267+
// eSV = sv.e();
268+
// slXY = sv.decayLengthXY() / sv.errorDecayLengthXY();
282269
}
283270

284271
registry.fill(HIST("h_jetpT"), analysisJet.pt());
285272
registry.fill(HIST("h_Db"), analysisJet.scoreML());
286273
registry.fill(HIST("h2_jetpT_Db"), analysisJet.pt(), analysisJet.scoreML());
287274
registry.fill(HIST("h2_jetpT_SVMass"), analysisJet.pt(), mSV);
288-
registry.fill(HIST("h2_jetpT_SVEnergy"), analysisJet.pt(), eSV);
289-
registry.fill(HIST("h2_jetpT_SLxy"), analysisJet.pt(), slXY);
290275
registry.fill(HIST("h2_jetpT_jetMass"), analysisJet.pt(), analysisJet.mass());
291276
registry.fill(HIST("h2_jetpT_jetProb"), analysisJet.pt(), analysisJet.jetProb());
292277
registry.fill(HIST("h2_jetpT_nTracks"), analysisJet.pt(), nTracks);
@@ -296,13 +281,13 @@ struct BjetTaggingGNN {
296281
}
297282
}
298283
}
299-
PROCESS_SWITCH(BjetTaggingGNN, processDataJets, "jet information in Data", false);
284+
PROCESS_SWITCH(BjetTaggingGnn, processDataJets, "jet information in Data", false);
300285

301286
using MCDJetTable = soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets, aod::ChargedMCDetectorLevelJetFlavourDef, aod::ChargedMCDetectorLevelJetTags, aod::ChargedMCDetectorLevelJetEventWeights, aod::MCDSecondaryVertex3ProngIndices>>;
302287
using MCPJetTable = soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents, aod::ChargedMCParticleLevelJetsMatchedToChargedMCDetectorLevelJets, aod::ChargedMCParticleLevelJetFlavourDef, aod::ChargedMCParticleLevelJetEventWeights>>;
303288
using FilteredCollisionMCD = soa::Filtered<soa::Join<aod::JetCollisions, aod::JCollisionPIs, aod::JMcCollisionLbs>>;
304289

305-
void processMCJets(FilteredCollisionMCD::iterator const& collision, MCDJetTable const& MCDjets, MCPJetTable const& /*MCPjets*/, JetTracksMCDwID const& allTracks, MCDSVTable const& allSVs, aod::JetParticles const& /*MCParticles*/)
290+
void processMCJets(FilteredCollisionMCD::iterator const& collision, MCDJetTable const& MCDjets, MCPJetTable const& /*MCPjets*/, JetTracksMCDwID const& /*allTracks*/, MCDSVTable const& allSVs, aod::JetParticles const& /*MCParticles*/)
306291
{
307292
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
308293
return;
@@ -332,49 +317,51 @@ struct BjetTaggingGNN {
332317

333318
int8_t jetFlavor = analysisJet.origin();
334319

335-
int nTracks = analyzeJetTrackInfo(collision, analysisJet, allTracks /*, jetFlavor, weight*/);
320+
// int nTracks = analyzeJetTrackInfo(collision, analysisJet, allTracks /*, jetFlavor, weight*/);
321+
int nTracks = 0;
336322

337323
int nNppTracks = 0;
338-
for (auto& constituent : analysisJet.template tracks_as<JetTracksMCDwID>()) {
339-
if (!constituent.has_mcParticle() || !constituent.template mcParticle_as<aod::JetParticles>().isPhysicalPrimary() || constituent.pt() < trackPtMin) {
324+
for (const auto& constituent : analysisJet.template tracks_as<JetTracksMCDwID>()) {
325+
if (constituent.pt() < trackPtMin) {
326+
continue;
327+
}
328+
if (!constituent.has_mcParticle() || !constituent.template mcParticle_as<aod::JetParticles>().isPhysicalPrimary()) {
329+
// registry.fill(HIST("h2_pT_dcaXY_npp"), constituent.pt(), constituent.dcaXY());
330+
// registry.fill(HIST("h2_pT_dcaZ_npp"), constituent.pt(), constituent.dcaZ());
340331
++nNppTracks;
332+
} else {
333+
// registry.fill(HIST("h2_pT_dcaXY_pp"), constituent.pt(), constituent.dcaXY());
334+
// registry.fill(HIST("h2_pT_dcaZ_pp"), constituent.pt(), constituent.dcaZ());
341335
}
336+
++nTracks;
342337
}
343338

344339
float mSV = -1.f;
345-
float eSV = -1.f;
346-
float slXY = -1.f;
340+
// float eSV = -1.f;
341+
// float slXY = -1.f;
347342

348343
bool checkSV;
349-
// auto sv = jettaggingutilities::jetFromProngMaxDecayLength<MCDSVTable>(analysisJet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, false, &checkSV);
350344
auto sv = analyzeJetSVInfo(analysisJet, allSVs, checkSV /*, jetFlavor, weight*/);
351345

352346
if (checkSV) {
353347
mSV = sv.m();
354-
eSV = sv.e();
355-
slXY = sv.decayLengthXY() / sv.errorDecayLengthXY();
356-
registry.fill(HIST("h3_mSV_eSV_slXY"), mSV, eSV, slXY, weight);
348+
// eSV = sv.e();
349+
// slXY = sv.decayLengthXY() / sv.errorDecayLengthXY();
357350
}
358351

359352
registry.fill(HIST("h_jetpT"), analysisJet.pt(), weight);
360353
registry.fill(HIST("h_Db"), analysisJet.scoreML(), weight);
361354
registry.fill(HIST("h2_jetpT_Db"), analysisJet.pt(), analysisJet.scoreML(), weight);
362355
registry.fill(HIST("h2_jetpT_SVMass"), analysisJet.pt(), mSV, weight);
363-
registry.fill(HIST("h2_jetpT_SVEnergy"), analysisJet.pt(), eSV, weight);
364-
registry.fill(HIST("h2_jetpT_SLxy"), analysisJet.pt(), slXY, weight);
365356
registry.fill(HIST("h2_jetpT_jetMass"), analysisJet.pt(), analysisJet.mass(), weight);
366357
registry.fill(HIST("h2_jetpT_jetProb"), analysisJet.pt(), analysisJet.jetProb(), weight);
367358
registry.fill(HIST("h2_jetpT_nTracks"), analysisJet.pt(), nTracks, weight);
368359

369-
registry.fill(HIST("h3_jetpT_nTracks_Db"), analysisJet.pt(), nTracks, analysisJet.scoreML(), weight);
370-
371360
if (jetFlavor == JetTaggingSpecies::beauty) {
372361
registry.fill(HIST("h_jetpT_b"), analysisJet.pt(), weight);
373362
registry.fill(HIST("h_Db_b"), analysisJet.scoreML(), weight);
374363
registry.fill(HIST("h2_jetpT_Db_b"), analysisJet.pt(), analysisJet.scoreML(), weight);
375364
registry.fill(HIST("h2_jetpT_SVMass_b"), analysisJet.pt(), mSV, weight);
376-
registry.fill(HIST("h2_jetpT_SVEnergy_b"), analysisJet.pt(), eSV, weight);
377-
registry.fill(HIST("h2_jetpT_SLxy_b"), analysisJet.pt(), slXY, weight);
378365
registry.fill(HIST("h2_jetpT_jetMass_b"), analysisJet.pt(), analysisJet.mass(), weight);
379366
registry.fill(HIST("h2_jetpT_jetProb_b"), analysisJet.pt(), analysisJet.jetProb(), weight);
380367
registry.fill(HIST("h2_jetpT_nTracks_b"), analysisJet.pt(), nTracks, weight);
@@ -383,8 +370,6 @@ struct BjetTaggingGNN {
383370
registry.fill(HIST("h_Db_c"), analysisJet.scoreML(), weight);
384371
registry.fill(HIST("h2_jetpT_Db_c"), analysisJet.pt(), analysisJet.scoreML(), weight);
385372
registry.fill(HIST("h2_jetpT_SVMass_c"), analysisJet.pt(), mSV, weight);
386-
registry.fill(HIST("h2_jetpT_SVEnergy_c"), analysisJet.pt(), eSV, weight);
387-
registry.fill(HIST("h2_jetpT_SLxy_c"), analysisJet.pt(), slXY, weight);
388373
registry.fill(HIST("h2_jetpT_jetMass_c"), analysisJet.pt(), analysisJet.mass(), weight);
389374
registry.fill(HIST("h2_jetpT_jetProb_c"), analysisJet.pt(), analysisJet.jetProb(), weight);
390375
registry.fill(HIST("h2_jetpT_nTracks_c"), analysisJet.pt(), nTracks, weight);
@@ -393,8 +378,6 @@ struct BjetTaggingGNN {
393378
registry.fill(HIST("h_Db_lf"), analysisJet.scoreML(), weight);
394379
registry.fill(HIST("h2_jetpT_Db_lf"), analysisJet.pt(), analysisJet.scoreML(), weight);
395380
registry.fill(HIST("h2_jetpT_SVMass_lf"), analysisJet.pt(), mSV, weight);
396-
registry.fill(HIST("h2_jetpT_SVEnergy_lf"), analysisJet.pt(), eSV, weight);
397-
registry.fill(HIST("h2_jetpT_SLxy_lf"), analysisJet.pt(), slXY, weight);
398381
registry.fill(HIST("h2_jetpT_jetMass_lf"), analysisJet.pt(), analysisJet.mass(), weight);
399382
registry.fill(HIST("h2_jetpT_jetProb_lf"), analysisJet.pt(), analysisJet.jetProb(), weight);
400383
registry.fill(HIST("h2_jetpT_nTracks_lf"), analysisJet.pt(), nTracks, weight);
@@ -452,7 +435,7 @@ struct BjetTaggingGNN {
452435
}
453436
}
454437
}
455-
PROCESS_SWITCH(BjetTaggingGNN, processMCJets, "jet information in MC", false);
438+
PROCESS_SWITCH(BjetTaggingGnn, processMCJets, "jet information in MC", false);
456439

457440
Filter mccollisionFilter = nabs(aod::jmccollision::posZ) < vertexZCut;
458441
using FilteredCollisionMCP = soa::Filtered<aod::JMcCollisions>;
@@ -493,11 +476,11 @@ struct BjetTaggingGNN {
493476
}
494477
}
495478
}
496-
PROCESS_SWITCH(BjetTaggingGNN, processMCTruthJets, "truth jet information", false);
479+
PROCESS_SWITCH(BjetTaggingGnn, processMCTruthJets, "truth jet information", false);
497480
};
498481

499482
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
500483
{
501484
return WorkflowSpec{
502-
adaptAnalysisTask<BjetTaggingGNN>(cfgc, TaskName{"bjet-tagging-gnn"})};
485+
adaptAnalysisTask<BjetTaggingGnn>(cfgc)};
503486
}

0 commit comments

Comments
 (0)