Skip to content

Making jet prob a number and common #22

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
hanseopark merged 1 commit into from
Dec 23, 2024
Merged

Making jet prob a number and common #22

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 32 additions & 35 deletions PWGJE/Core/JetTaggingUtilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,32 +116,6 @@ bool isCHadron(int pc)
return (std::find(bPdG.begin(), bPdG.end(), std::abs(pc)) != bPdG.end());
}

template <typename T, typename U, typename V = float>
uint8_t setTaggingIPBit(T const& jet, U const& jtracks, V const& trackDcaXYMax, V const& trackDcaZMax, V const& tagPointForIP, int const& minIPCount)
{
uint8_t bit = 0;
if (isGreatherTanTaggingPoint(jet, jtracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount, false)) {
SETBIT(bit, TaggingMethodNonML::IPs);
}
if (isGreatherTanTaggingPoint(jet, jtracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount, true)) {
SETBIT(bit, TaggingMethodNonML::IPs3D);
}
return bit;
}

template <typename T, typename U, typename V = float>
uint8_t setTaggingSVBit(T const& jet, U const& prongs, V const& prongChi2PCAMax, V const& prongsigmaLxyMax, V const& svDispersionMax, V const& tagPointForSV)
{
uint8_t bit = 0;
if(isTaggedJetSV(jet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, svDispersionMax, false, tagPointForSV) {
SETBIT(bit, TaggingMethodNonML::SV);
}
if(isTaggedJetSV(jet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, svDispersionMax, true, tagPointForSV) {
SETBIT(bit, TaggingMethodNonML::SV3D);
}
return bit;
}

/**
* returns the globalIndex of the earliest mother of a particle in the shower. returns -1 if a suitable mother is not found
*
Expand Down Expand Up @@ -522,7 +496,7 @@ int getGeoSign(T const& jet, U const& jtrack)
* in a vector in descending order.
*/
template <typename T, typename U, typename Vec = std::vector<float>>
void orderForIPJetTracks(T const& jet, U const& /*jtracks*/, float const& trackDcaXYMax, float const& trackDcaZMax, Vec& vecSignImpSig, bool useIPxyz)
void orderForIPJetTracks(T const& jet, U const& /*jtracks*/, float trackDcaXYMax, float trackDcaZMax, Vec& vecSignImpSig, bool useIPxyz)
{
for (auto const& jtrack : jet.template tracks_as<U>()) {
if (!trackAcceptanceWithDca(jtrack, trackDcaXYMax, trackDcaZMax))
Expand All @@ -543,7 +517,7 @@ void orderForIPJetTracks(T const& jet, U const& /*jtracks*/, float const& trackD
* Checks if a jet is greater than the given tagging working point based on the signed impact parameter significances
*/
template <typename T, typename U>
bool isGreaterThanTaggingPoint(T const& jet, U const& jtracks, float const& trackDcaXYMax, float const& trackDcaZMax, float const& taggingPoint = 1.0, int const& cnt = 1, bool useIPxyz = false)
bool isGreaterThanTaggingPoint(T const& jet, U const& jtracks, float trackDcaXYMax, float trackDcaZMax, float taggingPoint = 1.0, int cnt = 1, bool useIPxyz = false)
{
if (cnt == 0) {
return true; // untagged
Expand Down Expand Up @@ -595,7 +569,7 @@ std::unique_ptr<TF1> setResolutionFunction(T const& vecParams)
* impact parameter significance.
*/
template <typename T, typename U>
float getTrackProbability(T const& fResoFuncjet, U const& track, const float& minSignImpXYSig = -40)
float getTrackProbability(T const& fResoFuncjet, U const& track, float minSignImpXYSig = -40)
{
float probTrack = 0;
auto varSignImpXYSig = std::abs(track.dcaXY()) / track.sigmadcaXY();
Expand Down Expand Up @@ -625,11 +599,8 @@ float getTrackProbability(T const& fResoFuncjet, U const& track, const float& mi
* geometric sign.
*/
template <typename T, typename U, typename V>
float getJetProbability(T const& fResoFuncjet, U const& jet, V const& jtracks, float const& trackDcaXYMax, float const& trackDcaZMax, const int& cnt, const float& tagPoint = 1.0, const float& minSignImpXYSig = -10, bool useIPxy = true)
float getJetProbability(T const& fResoFuncjet, U const& jet, V const& /*jtracks*/, float trackDcaXYMax, float trackDcaZMax, float minSignImpXYSig = -10)
{
if (!(isGreaterThanTaggingPoint(jet, jtracks, trackDcaXYMax, trackDcaZMax, tagPoint, cnt, useIPxy)))
return -1;

std::vector<float> jetTracksPt;
float trackjetProb = 1.;

Expand Down Expand Up @@ -661,7 +632,7 @@ float getJetProbability(T const& fResoFuncjet, U const& jet, V const& jtracks, f

// For secaondy vertex method utilites
template <typename ProngType, typename JetType>
typename ProngType::iterator jetFromProngMaxDecayLength(const JetType& jet, float const& prongChi2PCAMin, float const& prongChi2PCAMax, float const& prongsigmaLxyMax, float const& prongIPxyMin, float const& prongIPxyMax, const bool& doXYZ = false, bool* checkSv = nullptr)
typename ProngType::iterator jetFromProngMaxDecayLength(const JetType& jet, float const& prongChi2PCAMin, float prongChi2PCAMax, float prongsigmaLxyMax, float prongIPxyMin, float prongIPxyMax, bool doXYZ = false, bool* checkSv = nullptr)
{
if (checkSv)
*checkSv = false;
Expand All @@ -685,7 +656,7 @@ typename ProngType::iterator jetFromProngMaxDecayLength(const JetType& jet, floa
}

template <typename T, typename U>
bool isTaggedJetSV(T const jet, U const& /*prongs*/, float const& prongChi2PCAMin, float const& prongChi2PCAMax, float const& prongsigmaLxyMax, float const& prongIPxyMin, float const& prongIPxyMax, float svDispersionMax, float const& doXYZ = false, float const& tagPointForSV = 15.)
bool isTaggedJetSV(T const jet, U const& /*prongs*/, float prongChi2PCAMin, float prongChi2PCAMax, float prongsigmaLxyMax, float prongIPxyMin, float prongIPxyMax, float svDispersionMax, float doXYZ = false, float tagPointForSV = 15.)
{
bool checkSv = false;
auto bjetCand = jetFromProngMaxDecayLength<U>(jet, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, doXYZ, &checkSv);
Expand All @@ -703,6 +674,32 @@ bool isTaggedJetSV(T const jet, U const& /*prongs*/, float const& prongChi2PCAMi
return true;
}

template <typename T, typename U, typename V = float>
uint8_t setTaggingIPBit(T const& jet, U const& jtracks, V trackDcaXYMax, V trackDcaZMax, V tagPointForIP, int minIPCount)
{
uint8_t bit = 0;
if (isGreaterThanTaggingPoint(jet, jtracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount, false)) {
SETBIT(bit, TaggingMethodNonML::IPs);
}
if (isGreaterThanTaggingPoint(jet, jtracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount, true)) {
SETBIT(bit, TaggingMethodNonML::IPs3D);
}
return bit;
}

template <typename T, typename U, typename V = float>
uint8_t setTaggingSVBit(T const& jet, U const& prongs, V prongChi2PCAMin, V prongChi2PCAMax, V prongsigmaLxyMax, float prongIPxyMin, float prongIPxyMax, V svDispersionMax, V tagPointForSV)
{
uint8_t bit = 0;
if (isTaggedJetSV(jet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, svDispersionMax, false, tagPointForSV)) {
SETBIT(bit, TaggingMethodNonML::SV);
}
if (isTaggedJetSV(jet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, svDispersionMax, true, tagPointForSV)) {
SETBIT(bit, TaggingMethodNonML::SV3D);
}
return bit;
}

/**
* Clusters jet constituent tracks into groups of tracks originating from same mcParticle position (trkVtxIndex), and finds each track origin (trkOrigin). (for GNN b-jet tagging)
* @param trkLabels Track labels for GNN vertex and track origin predictions. trkVtxIndex: The index value of each vertex (cluster) which is determined by the function. trkOrigin: The category of the track origin (0: not physical primary, 1: charm, 2: beauty, 3: primary vertex, 4: other secondary vertex).
Expand Down
4 changes: 2 additions & 2 deletions PWGJE/DataModel/JetTagging.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,8 +131,8 @@ JETSV_TABLES_DEF(Charged, SecondaryVertex2Prong, "2PRONG");
#define JETTAGGING_TABLE_DEF(_jet_type_, _name_, _description_) \
namespace _name_##tagging \
{ \
DECLARE_SOA_COLUMN(bitTaggedjetNonML, BitTaggedjetNonML, uint8_t); \
DECLARE_SOA_COLUMN(JetProb, jetProb, std::vector<float>); \
DECLARE_SOA_COLUMN(BitTaggedjetNonML, bitTaggedjetNonML, uint8_t); \
DECLARE_SOA_COLUMN(JetProb, jetProb, float); \
DECLARE_SOA_COLUMN(ScoreJetML, scoreML, float); \
} \
DECLARE_SOA_TABLE(_jet_type_##Tags, "AOD", _description_ "Tags", _name_##tagging::BitTaggedjetNonML, _name_##tagging::JetProb, _name_##tagging::ScoreJetML);
Expand Down
2 changes: 1 addition & 1 deletion PWGJE/TableProducer/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@ o2physics_add_dpl_workflow(jet-hf-definition

o2physics_add_dpl_workflow(jet-taggerhf
SOURCES jetTaggerHF.cxx
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore O2Physics::MLCore
COMPONENT_NAME Analysis)

o2physics_add_dpl_workflow(estimator-rho
Expand Down
10 changes: 3 additions & 7 deletions PWGJE/TableProducer/heavyFlavourDefinition.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,6 @@
#include <memory>
#include <vector>

#include <TF1.h>
#include <TH1.h>

#include "Framework/AnalysisTask.h"
#include "Framework/AnalysisDataModel.h"
#include "Framework/ASoA.h"
Expand All @@ -29,7 +26,6 @@
#include "PWGJE/DataModel/Jet.h"
#include "PWGJE/DataModel/JetTagging.h"
#include "PWGJE/Core/JetTaggingUtilities.h"
#include "PWGJE/Core/JetDerivedDataUtilities.h"

using namespace o2;
using namespace o2::framework;
Expand All @@ -44,7 +40,7 @@ struct HeavyFlavourDefinitionTask {
Configurable<float> maxDeltaR{"maxDeltaR", 0.25, "maximum distance of jet axis from flavour initiating parton"};
Configurable<bool> searchUpToQuark{"searchUpToQuark", true, "Finding first mother in particles to quark"};

using JetTracksMCD = soa::Join<aod::JetTracksMCD, aod::JTrackExtras, aod::JTrackPIs>;
using JetTracksMCD = soa::Join<aod::JetTracksMCD, aod::JTrackPIs>;
Preslice<aod::JetParticles> particlesPerCollision = aod::jmcparticle::mcCollisionId;

void init(InitContext const&)
Expand Down Expand Up @@ -114,8 +110,8 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

std::vector<o2::framework::DataProcessorSpec> tasks;

tasks.emplace_back(adaptAnalysisTask<JetHfDefinitionCharged>(cfgc));
tasks.emplace_back(adaptAnalysisTask<JetHfDefinitionFull>(cfgc));
tasks.emplace_back(adaptAnalysisTask<JetHfDefinitionCharged>(cfgc, SetDefaultProcesses{}, TaskName{"jet-hf-definition-charged"}));
tasks.emplace_back(adaptAnalysisTask<JetHfDefinitionFull>(cfgc, SetDefaultProcesses{}, TaskName{"jet-hf-definition-full"}));

return WorkflowSpec{tasks};
}
103 changes: 45 additions & 58 deletions PWGJE/TableProducer/jetTaggerHF.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,36 +108,28 @@ struct JetTaggerHFTask {
std::unique_ptr<TF1> fSignImpXYSigLfJetMC = nullptr;

std::vector<uint8_t> decisionNonML;
std::vector<std::vector<float>> decisionJP;
std::vector<float> scoreML;

template <typename T, typename U>
void calculateJetProbability(int origin, T const& jet, U const& jtracks, std::vector<float>& jetProb, bool const& isMC = true)
float calculateJetProbability(int origin, T const& jet, U const& jtracks, bool const& isMC = true)
{
jetProb.clear();
jetProb.reserve(maxOrder);
for (int order = 0; order < maxOrder; order++) {
if (!isMC) {
jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigData, jet, jtracks, trackDcaXYMax, trackDcaZMax, order, tagPointForIP, minSignImpXYSig));
float jetProb = -1.0;
if (!isMC) {
jetProb = jettaggingutilities::getJetProbability(fSignImpXYSigData, jet, jtracks, trackDcaXYMax, trackDcaZMax, minSignImpXYSig);
} else {
if (useResoFuncFromIncJet) {
jetProb = jettaggingutilities::getJetProbability(fSignImpXYSigIncJetMC, jet, jtracks, trackDcaXYMax, trackDcaZMax, minSignImpXYSig);
} else {
if (useResoFuncFromIncJet) {
jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigIncJetMC, jet, jtracks, trackDcaXYMax, trackDcaZMax, order, tagPointForIP, minSignImpXYSig));
if (origin == JetTaggingSpecies::charm) {
jetProb = jettaggingutilities::getJetProbability(fSignImpXYSigCharmJetMC, jet, jtracks, trackDcaXYMax, trackDcaZMax, minSignImpXYSig);
} else if (origin == JetTaggingSpecies::beauty) {
jetProb = jettaggingutilities::getJetProbability(fSignImpXYSigBeautyJetMC, jet, jtracks, trackDcaXYMax, trackDcaZMax, minSignImpXYSig);
} else {
if (origin == JetTaggingSpecies::charm) {
jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigCharmJetMC, jet, jtracks, trackDcaXYMax, trackDcaZMax, order, tagPointForIP, minSignImpXYSig));
}
if (origin == JetTaggingSpecies::beauty) {
jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigBeautyJetMC, jet, jtracks, trackDcaXYMax, trackDcaZMax, order, tagPointForIP, minSignImpXYSig));
}
if (origin == JetTaggingSpecies::lightflavour) {
jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigLfJetMC, jet, jtracks, trackDcaXYMax, trackDcaZMax, order, tagPointForIP, minSignImpXYSig));
}
if (origin != JetTaggingSpecies::charm && origin != JetTaggingSpecies::beauty && origin != JetTaggingSpecies::lightflavour) {
jetProb.push_back(-1);
}
jetProb = jettaggingutilities::getJetProbability(fSignImpXYSigLfJetMC, jet, jtracks, trackDcaXYMax, trackDcaZMax, minSignImpXYSig);
}
}
}
return jetProb;
}

template <typename T, typename U>
Expand Down Expand Up @@ -244,11 +236,11 @@ struct JetTaggerHFTask {
if (trackProbQA) {
AxisSpec trackProbabilityAxis = {binTrackProbability, "Track proability"};
AxisSpec jetFlavourAxis = {binJetFlavour, "Jet flavour"};
if (doprocessDataIP) {
if (doprocessFillTablesData) {
registry.add("h_pos_track_probability", "positive track probability", {HistType::kTH1F, {{trackProbabilityAxis}}});
registry.add("h_neg_track_probability", "negative track probability", {HistType::kTH1F, {{trackProbabilityAxis}}});
}
if (doprocessMCDIP) {
if (doprocessFillTablesMCD) {
registry.add("h2_pos_track_probability_flavour", "positive track probability", {HistType::kTH2F, {{trackProbabilityAxis}, {jetFlavourAxis}}});
registry.add("h2_neg_track_probability_flavour", "negative track probability", {HistType::kTH2F, {{trackProbabilityAxis}, {jetFlavourAxis}}});
}
Expand Down Expand Up @@ -299,70 +291,65 @@ struct JetTaggerHFTask {
}
PROCESS_SWITCH(JetTaggerHFTask, processDummy, "Dummy process", true);

void processSetUp(JetTable const& jets)
void processSetup(JetTable const& jets)
{
decisionNonML.clear();
decisionJP.clear();
scoreML.clear();
decisionNonML.resize(jets.size());
decisionJP.resize(jets.size());
scoreML.resize(jets.size());
}
PROCESS_SWITCH(JetTaggerHFTask, processSetup, "Setup initialization and size of jets for filling table", true);

void processDataIP(aod::JetCollision const& /*collision*/, JetTable const& jets, JetTracksExt const& jtracks)
void processIP(JetTable const& jets, JetTracksExt const& jtracks)
{
for (auto& jet : jets) {
if (useJetProb) {
calculateJetProbability(0, jet, jtracks, jetProb, false);
if (trackProbQA) {
evaluateTrackProbQA(0, jet, jtracks, false);
}
}
decisionJP[jet.globalIndex()] = jetProb;
uint8_t bit = jettaggingutilities::setTaggingIPBit(jet, jtracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount);
decisionNonML[jet.globalIndex()] |= bit;
}
}
PROCESS_SWITCH(JetTaggerHFTask, processDataIP, "Fill tagging decision for data jets", false);
PROCESS_SWITCH(JetTaggerHFTask, processIP, "Fill tagging decision for jet with the IP algorithm", false);

void processMCDIP(aod::JetCollision const& /*collision*/, soa::Join<JetTable, aod::ChargedMCDetectorLevelJetFlavourDef> const& jets, JetTracksExt const& jtracks)
void processSV(soa::Join<JetTable, SVIndicesTable> const& jets, SVTable const& prongs)
{
for (auto& jet : jets) {
int origin = jet.origin();
if (useJetProb) {
calculateJetProbability(origin, jet, jtracks, jetProb);
if (trackProbQA) {
evaluateTrackProbQA(origin, jet, jtracks);
}
}
decisionJP[jet.globalIndex()] = jetProb;
uint8_t bit = jettaggingutilities::setTaggingIPBit(jet, jtracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount);
uint8_t bit = jettaggingutilities::setTaggingSVBit(jet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, prongIPxyMin, prongIPxyMax, svDispersionMax, tagPointForSV);
decisionNonML[jet.globalIndex()] |= bit;
}
}
PROCESS_SWITCH(JetTaggerHFTask, processMCDIP, "Fill tagging decision for mcd jets", false);

void processSV(aod::JetCollision const& /*collision*/, soa::Join<JetTable, SVIndicesTable> const& jets, JetTracksExt const& /*jtracks*/, SVTable const& prongs)
{
for (auto& jet : jets) {
uint8_t bit = jettaggingutilities::setTaggingSVBit(jet, prongs, prongChi2PCAMin, prongChi2PCAMax, prongsigmaLxyMax, svDispersionMax, tagPointForSV);
decisionNonML[jet.globalIndex()] |= bit;
}
}
PROCESS_SWITCH(JetTaggerHFTask, processSV, "Fill tagging decision for charged jets", false);
PROCESS_SWITCH(JetTaggerHFTask, processSV, "Fill tagging decision for jets with the SV", false);

void processAlgorithmML(soa::Join<JetTable, SVIndicesTable> const& allJets, JetTracksExt const& allTracks, SVTable const& allSVs)
{
analyzeJetAlgorithmML(allJets, allTracks, allSVs);
}
PROCESS_SWITCH(JetTaggerHFTask, processAlgorithmML, "Fill ML evaluation score for charged jets", false);

void processFillTables(JetTable::iterator const& jet)
void processFillTablesData(JetTable::iterator const& jet, JetTracksExt const& jtracks)
{
taggingTableMCD(decisionNonML[jet.globalIndex()], decisionJP[jet.globalIndex()], scoreML[jet.globalIndex()]);
float jetProb = -1.0;
if (useJetProb) {
jetProb = calculateJetProbability(0, jet, jtracks, false);
if (trackProbQA) {
evaluateTrackProbQA(0, jet, jtracks, false);
}
}
taggingTable(decisionNonML[jet.globalIndex()], jetProb, scoreML[jet.globalIndex()]);
}
PROCESS_SWITCH(JetTaggerHFTask, processFillTablesData, "Fill Tables for tagging decision, jet probability, and ML score on charged jets in data", false);

void processFillTablesMCD(soa::Join<JetTable, aod::ChargedMCDetectorLevelJetFlavourDef>::iterator const& jet, JetTracksExt const& jtracks)
{
float jetProb = -1.0;
int origin = jet.origin();
if (useJetProb) {
jetProb = calculateJetProbability(origin, jet, jtracks);
if (trackProbQA) {
evaluateTrackProbQA(origin, jet, jtracks);
}
}
taggingTable(decisionNonML[jet.globalIndex()], jetProb, scoreML[jet.globalIndex()]);
}
PROCESS_SWITCH(JetTaggerHFTask, processFillTables, "Fill Tables for tagging decision and ML score on charged jets", false);
PROCESS_SWITCH(JetTaggerHFTask, processFillTablesMCD, "Fill Tables for tagging decision, jet probability, and ML score on charged jets in MC", false);
};

using JetTaggerHFDataCharged = JetTaggerHFTask<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents>, aod::DataSecondaryVertex3ProngIndices, aod::DataSecondaryVertex3Prongs, aod::ChargedJetTags>;
Expand Down
Loading
Loading