From efa1dc16d91f2ecfb26ee7276358b0d11ea1a137 Mon Sep 17 00:00:00 2001 From: Changhwan Choi Date: Thu, 16 Oct 2025 21:36:01 +0900 Subject: [PATCH] Updated GNN b-jet analysis histograms, applied proper event selections --- PWGJE/TableProducer/jetTaggerHF.cxx | 6 +- PWGJE/Tasks/bjetTaggingGnn.cxx | 230 +++++++++++++++++++++------- 2 files changed, 179 insertions(+), 57 deletions(-) diff --git a/PWGJE/TableProducer/jetTaggerHF.cxx b/PWGJE/TableProducer/jetTaggerHF.cxx index 70d88de3e45..44899e4c702 100644 --- a/PWGJE/TableProducer/jetTaggerHF.cxx +++ b/PWGJE/TableProducer/jetTaggerHF.cxx @@ -280,10 +280,10 @@ struct JetTaggerHFTask { if (jet.pt() >= jetpTMin) { if (scoreML[jet.globalIndex()] < dbMin) { dbRange = 0.5; // underflow - } else if (scoreML[jet.globalIndex()] >= dbMax) { - dbRange = 2.5; // overflow - } else { + } else if (scoreML[jet.globalIndex()] < dbMax) { dbRange = 1.5; // in range + } else { + dbRange = 2.5; // overflow } registry.fill(HIST("h2_count_db"), 3.5, dbRange); // incl jet if constexpr (isMC) { diff --git a/PWGJE/Tasks/bjetTaggingGnn.cxx b/PWGJE/Tasks/bjetTaggingGnn.cxx index a71c0d7f75b..9d2c82b47c9 100644 --- a/PWGJE/Tasks/bjetTaggingGnn.cxx +++ b/PWGJE/Tasks/bjetTaggingGnn.cxx @@ -45,6 +45,7 @@ struct BjetTaggingGnn { // event level configurables Configurable vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"}; Configurable eventSelections{"eventSelections", "sel8", "choose event selection"}; + Configurable eventSelectionsSel{"eventSelectionsSel", "sel8", "choose event selection"}; Configurable useEventWeight{"useEventWeight", true, "Flag whether to scale histograms with the event weight"}; Configurable pTHatMaxMCD{"pTHatMaxMCD", 999.0, "maximum fraction of hard scattering for jet acceptance in detector MC"}; @@ -92,7 +93,7 @@ struct BjetTaggingGnn { jetRadiiValues = (std::vector)jetRadii; eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); - eventSelectionBitsSel = jetderiveddatautilities::initialiseEventSelectionBits("sel8"); + eventSelectionBitsSel = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelectionsSel)); registry.add("h_vertexZ", "Vertex Z;#it{Z} (cm)", {HistType::kTH1F, {{100, -20.0, 20.0}}}, callSumw2); // registry.add("h_vertexZ_truth", "Vertex Z truth;#it{Z} (cm)", {HistType::kTH1F, {{100, -20.0, 20.0}}}, callSumw2); @@ -160,6 +161,11 @@ struct BjetTaggingGnn { registry.add("h_Db_npp_b", "NotPhysPrim b-jet", {HistType::kTH1F, {axisDbFine}}); registry.add("h_Db_npp_c", "NotPhysPrim c-jet", {HistType::kTH1F, {axisDbFine}}); registry.add("h_Db_npp_lf", "NotPhysPrim lf-jet", {HistType::kTH1F, {axisDbFine}}); + registry.add("h_jetpT_matched", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_matched", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_b_matched", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_matched", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("hSparse_overflow_b", "", {HistType::kTHnSparseF, {axisJetpT, axisJetpT, axisNTracks}}); } if (doprocessMCJetsSel) { @@ -182,14 +188,32 @@ struct BjetTaggingGnn { registry.add("h_jetpT_particle_b_tvx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); } + if (doprocessMCTruthJetsOld) { + registry.add("h_jetpT_particle_old", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_old", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_old", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_lf_old", "particle lf-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + } + + if (doprocessMCCollision) { + registry.add("h_jetpT_particle2", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b2", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c2", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_lf2", "particle lf-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_sel2", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_sel2", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_tvx2", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_tvx2", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + } + if (doDataDriven) { - registry.add("hSparse_Incljets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2); + registry.add("hSparse_Incljets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); if (doprocessMCJets) { - registry.add("hSparse_bjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2); - registry.add("hSparse_cjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2); - registry.add("hSparse_lfjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2); - registry.add("hSparse_lfjets_none", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2); - registry.add("hSparse_lfjets_matched", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisSVMass, axisJetMass, axisNTracks}}, callSumw2); + registry.add("hSparse_bjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); + registry.add("hSparse_cjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); + registry.add("hSparse_lfjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); + registry.add("hSparse_lfjets_none", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); + registry.add("hSparse_lfjets_matched", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); } } } @@ -210,7 +234,11 @@ struct BjetTaggingGnn { using AnalysisCollisionsMCD = soa::Join; using FilteredCollisionsMCD = soa::Filtered; + + Filter mccollisionFilter = nabs(aod::jmccollision::posZ) < vertexZCut; + using FilteredCollisionMCP = soa::Filtered; using MCPJets = soa::Join; + using FilteredMCPJets = soa::Filtered; using SVTable = aod::DataSecondaryVertex3Prongs; using MCDSVTable = aod::MCDSecondaryVertex3Prongs; @@ -232,40 +260,40 @@ struct BjetTaggingGnn { return nTracks; } - template - SecondaryVertices::iterator analyzeJetSVInfo(AnalysisJet const& analysisJet, SecondaryVertices const& allSVs, bool& checkSV /*, int8_t jetFlavor = 0, double weight = 1.0*/) - { - using SVType = typename SecondaryVertices::iterator; + // template + // SecondaryVertices::iterator analyzeJetSVInfo(AnalysisJet const& analysisJet, SecondaryVertices const& allSVs, bool& checkSV /*, int8_t jetFlavor = 0, double weight = 1.0*/) + // { + // using SVType = typename SecondaryVertices::iterator; - auto compare = [](SVType& sv1, SVType& sv2) { - return (sv1.decayLengthXY() / sv1.errorDecayLengthXY()) > (sv2.decayLengthXY() / sv2.errorDecayLengthXY()); - }; + // auto compare = [](SVType& sv1, SVType& sv2) { + // return (sv1.decayLengthXY() / sv1.errorDecayLengthXY()) > (sv2.decayLengthXY() / sv2.errorDecayLengthXY()); + // }; - auto svs = analysisJet.template secondaryVertices_as(); + // auto svs = analysisJet.template secondaryVertices_as(); - std::sort(svs.begin(), svs.end(), compare); + // std::sort(svs.begin(), svs.end(), compare); - checkSV = false; - for (const auto& candSV : svs) { + // checkSV = false; + // for (const auto& candSV : svs) { - if (candSV.pt() < svPtMin) { - continue; - } + // if (candSV.pt() < svPtMin) { + // continue; + // } - checkSV = true; - return candSV; - } + // checkSV = true; + // return candSV; + // } - // No SV found - return *allSVs.begin(); - } + // // No SV found + // return *allSVs.begin(); + // } void processDummy(FilteredCollisions::iterator const& /*collision*/) { } PROCESS_SWITCH(BjetTaggingGnn, processDummy, "Dummy process function turned on by default", true); - void processDataJets(FilteredCollisions::iterator const& collision, FilteredDataJets const& alljets, JetTrackswID const& allTracks, SVTable const& allSVs) + void processDataJets(FilteredCollisions::iterator const& collision, FilteredDataJets const& alljets, JetTrackswID const& allTracks) { if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { return; @@ -290,21 +318,21 @@ struct BjetTaggingGnn { int nTracks = analyzeJetTrackInfo(collision, analysisJet, allTracks); - float mSV = -1.f; + // float mSV = -1.f; - bool checkSV; - auto sv = analyzeJetSVInfo(analysisJet, allSVs, checkSV); + // bool checkSV; + // auto sv = analyzeJetSVInfo(analysisJet, allSVs, checkSV); - if (checkSV) { - mSV = sv.m(); - } + // if (checkSV) { + // mSV = sv.m(); + // } registry.fill(HIST("h_jetpT"), analysisJet.pt()); registry.fill(HIST("h_Db"), analysisJet.scoreML()); registry.fill(HIST("h2_jetpT_Db"), analysisJet.pt(), analysisJet.scoreML()); if (doDataDriven) { - registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), mSV, analysisJet.mass(), nTracks); + registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), nTracks); } } } @@ -355,7 +383,7 @@ struct BjetTaggingGnn { } PROCESS_SWITCH(BjetTaggingGnn, processDataJetsSel, "jet information in Data (sel8)", false); - void processMCJets(FilteredCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, MCPJets const& /*MCPjets*/, JetTracksMCDwID const& /*allTracks*/, MCDSVTable const& allSVs, aod::JetParticles const& /*MCParticles*/) + void processMCJets(FilteredCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, JetTracksMCDwID const& /*allTracks*/, FilteredMCPJets const& /*MCPjets*/, aod::JetParticles const& /*MCParticles*/) { float weightEvt = useEventWeight ? collision.weight() : 1.f; if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { @@ -406,14 +434,14 @@ struct BjetTaggingGnn { ++nTracks; } - float mSV = -1.f; + // float mSV = -1.f; - bool checkSV; - auto sv = analyzeJetSVInfo(analysisJet, allSVs, checkSV /*, jetFlavor, weight*/); + // bool checkSV; + // auto sv = analyzeJetSVInfo(analysisJet, allSVs, checkSV /*, jetFlavor, weight*/); - if (checkSV) { - mSV = sv.m(); - } + // if (checkSV) { + // mSV = sv.m(); + // } registry.fill(HIST("h_jetpT"), analysisJet.pt(), weight); registry.fill(HIST("h_Db"), analysisJet.scoreML(), weight); @@ -454,29 +482,33 @@ struct BjetTaggingGnn { } if (doDataDriven) { - registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), mSV, analysisJet.mass(), nTracks, weight); + registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weight); if (jetFlavor == JetTaggingSpecies::beauty) { - registry.fill(HIST("hSparse_bjets"), analysisJet.pt(), analysisJet.scoreML(), mSV, analysisJet.mass(), nTracks, weight); + registry.fill(HIST("hSparse_bjets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weight); } else if (jetFlavor == JetTaggingSpecies::charm) { - registry.fill(HIST("hSparse_cjets"), analysisJet.pt(), analysisJet.scoreML(), mSV, analysisJet.mass(), nTracks, weight); + registry.fill(HIST("hSparse_cjets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weight); } else { - registry.fill(HIST("hSparse_lfjets"), analysisJet.pt(), analysisJet.scoreML(), mSV, analysisJet.mass(), nTracks, weight); + registry.fill(HIST("hSparse_lfjets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weight); if (jetFlavor == JetTaggingSpecies::none) { - registry.fill(HIST("hSparse_lfjets_none"), analysisJet.pt(), analysisJet.scoreML(), mSV, analysisJet.mass(), nTracks, weight); + registry.fill(HIST("hSparse_lfjets_none"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weight); } else { - registry.fill(HIST("hSparse_lfjets_matched"), analysisJet.pt(), analysisJet.scoreML(), mSV, analysisJet.mass(), nTracks, weight); + registry.fill(HIST("hSparse_lfjets_matched"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weight); } } } - for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { + for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { if (mcpjet.pt() > pTHatMaxMCP * pTHat) { continue; } registry.fill(HIST("h2_Response_DetjetpT_PartjetpT"), analysisJet.pt(), mcpjet.pt(), weight); + registry.fill(HIST("h_jetpT_matched"), analysisJet.pt(), weight); + registry.fill(HIST("h_jetpT_particle_matched"), mcpjet.pt(), weight); if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_b"), analysisJet.pt(), mcpjet.pt(), weight); + registry.fill(HIST("h_jetpT_b_matched"), analysisJet.pt(), weight); + registry.fill(HIST("h_jetpT_particle_b_matched"), mcpjet.pt(), weight); } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_c"), analysisJet.pt(), mcpjet.pt(), weight); } else { @@ -487,7 +519,7 @@ struct BjetTaggingGnn { } PROCESS_SWITCH(BjetTaggingGnn, processMCJets, "jet information in MC", false); - void processMCJetsSel(AnalysisCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, MCPJets const& /*MCPjets*/) + void processMCJetsSel(AnalysisCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, FilteredMCPJets const& /*MCPjets*/) { float weightEvt = useEventWeight ? collision.weight() : 1.f; registry.fill(HIST("h_event_counter"), 0.5, weightEvt); @@ -520,7 +552,7 @@ struct BjetTaggingGnn { registry.fill(HIST("h_jetpT_b_tvx"), analysisJet.pt(), weight); } - for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { + for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { if (mcpjet.pt() > pTHatMaxMCP * pTHat) { continue; } @@ -565,7 +597,7 @@ struct BjetTaggingGnn { registry.fill(HIST("h_jetpT_b_sel"), analysisJet.pt(), weight); } - for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { + for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { if (mcpjet.pt() > pTHatMaxMCP * pTHat) { continue; } @@ -581,7 +613,7 @@ struct BjetTaggingGnn { PresliceUnsorted collisionsPerMCPCollision = aod::jmccollisionlb::mcCollisionId; - void processMCTruthJets(MCPJets::iterator const& mcpjet, aod::JetParticles const& /*MCParticles*/, aod::JetMcCollisions const& /*mcCollisions*/, AnalysisCollisionsMCD const& collisions) + void processMCTruthJets(FilteredMCPJets::iterator const& mcpjet, aod::JetParticles const& /*MCParticles*/, aod::JetMcCollisions const& /*mcCollisions*/, AnalysisCollisionsMCD const& collisions) { bool jetIncluded = false; for (const auto& jetR : jetRadiiValues) { @@ -635,7 +667,9 @@ struct BjetTaggingGnn { } PROCESS_SWITCH(BjetTaggingGnn, processMCTruthJets, "truth jet information", false); - void processMCCollision(aod::McCollisions::iterator const& mcCollision, AnalysisCollisionsMCD const& collisions) + Preslice mcpjetsPerMCPCollision = aod::jmccollisionlb::mcCollisionId; + + void processMCCollision(aod::McCollisions::iterator const& mcCollision, FilteredMCPJets const& mcpjets, AnalysisCollisionsMCD const& collisions, aod::JetParticles const& /*MCParticles*/) { float weightEvt = useEventWeight ? mcCollision.weight() : 1.f; registry.fill(HIST("h_event_counter"), 3.5, weightEvt); // McColl(INEL) @@ -648,9 +682,97 @@ struct BjetTaggingGnn { registry.fill(HIST("h_event_counter"), 5.5, weightEvt); // McColl(-> Coll+TVX+Sel8+...) } } + auto mcpjetspermcpcollision = mcpjets.sliceBy(mcpjetsPerMCPCollision, mcCollision.globalIndex()); + for (const auto& mcpjet : mcpjetspermcpcollision) { + bool jetIncluded = false; + for (const auto& jetR : jetRadiiValues) { + if (mcpjet.r() == static_cast(jetR * 100)) { + jetIncluded = true; + break; + } + } + + if (!jetIncluded) { + continue; + } + + float weight = useEventWeight ? mcpjet.eventWeight() : 1.0; + float pTHat = 10. / (std::pow(mcpjet.eventWeight(), 1.0 / pTHatExponent)); + if (mcpjet.pt() > pTHatMaxMCP * pTHat) { + continue; + } + + int8_t jetFlavor = mcpjet.origin(); + + registry.fill(HIST("h_jetpT_particle_tvx2"), mcpjet.pt(), weight); + + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("h_jetpT_particle_b_tvx2"), mcpjet.pt(), weight); + } + + if (collisionspermccollision.size() >= 1) { + if (jetderiveddatautilities::selectCollision(collisionspermccollision.begin(), eventSelectionBitsSel)) { + registry.fill(HIST("h_jetpT_particle_sel2"), mcpjet.pt(), weight); + + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("h_jetpT_particle_b_sel2"), mcpjet.pt(), weight); + } + } + + if (jetderiveddatautilities::selectCollision(collisionspermccollision.begin(), eventSelectionBits) && std::fabs(collisionspermccollision.begin().posZ()) < vertexZCut) { + registry.fill(HIST("h_jetpT_particle2"), mcpjet.pt(), weight); + + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("h_jetpT_particle_b2"), mcpjet.pt(), weight); + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("h_jetpT_particle_c2"), mcpjet.pt(), weight); + } else { + registry.fill(HIST("h_jetpT_particle_lf2"), mcpjet.pt(), weight); + } + } + } + } } PROCESS_SWITCH(BjetTaggingGnn, processMCCollision, "mc collision information", false); + void processMCTruthJetsOld(FilteredCollisionMCP::iterator const& /*collision*/, FilteredMCPJets const& MCPjets, aod::JetParticles const& /*MCParticles*/) + { + + for (const auto& mcpjet : MCPjets) { + + bool jetIncluded = false; + for (const auto& jetR : jetRadiiValues) { + if (mcpjet.r() == static_cast(jetR * 100)) { + jetIncluded = true; + break; + } + } + + if (!jetIncluded) { + continue; + } + + float weight = useEventWeight ? mcpjet.eventWeight() : 1.0; + float pTHat = 10. / (std::pow(mcpjet.eventWeight(), 1.0 / pTHatExponent)); + if (mcpjet.pt() > pTHatMaxMCP * pTHat) { + continue; + } + + int8_t jetFlavor = mcpjet.origin(); + + registry.fill(HIST("h_jetpT_particle_old"), mcpjet.pt(), weight); + + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("h_jetpT_particle_b_old"), mcpjet.pt(), weight); + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("h_jetpT_particle_c_old"), mcpjet.pt(), weight); + } else { + registry.fill(HIST("h_jetpT_particle_lf_old"), mcpjet.pt(), weight); + } + } + } + PROCESS_SWITCH(BjetTaggingGnn, processMCTruthJetsOld, "truth jet information", false); + PresliceUnsorted> perFoundBC = aod::evsel::foundBCId; void processBCs(soa::Join const& bcs, soa::Join const& collisions)