diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index 424bcceaa63..2f3d9a67d11 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -211,12 +211,12 @@ int jetTrackFromHFShower(T const& jet, U const& /*tracks*/, V const& particles, hasMcParticle = true; auto const& particle = track.template mcParticle_as(); origin = RecoDecay::getParticleOrigin(particles, particle, searchUpToQuark); - if (origin == 1 || origin == 2) { // 1=charm , 2=beauty + if (origin == RecoDecay::OriginType::Prompt || origin == RecoDecay::OriginType::NonPrompt) { // 1=charm , 2=beauty hftrack = track; - if (origin == 1) { + if (origin == RecoDecay::OriginType::Prompt) { return JetTaggingSpecies::charm; } - if (origin == 2) { + if (origin == RecoDecay::OriginType::NonPrompt) { return JetTaggingSpecies::beauty; } } @@ -244,12 +244,12 @@ int jetParticleFromHFShower(T const& jet, U const& particles, typename U::iterat for (const auto& particle : jet.template tracks_as()) { hfparticle = particle; // for init if origin is 1 or 2, the particle is not hfparticle origin = RecoDecay::getParticleOrigin(particles, particle, searchUpToQuark); - if (origin == 1 || origin == 2) { // 1=charm , 2=beauty + if (origin == RecoDecay::OriginType::Prompt || origin == RecoDecay::OriginType::NonPrompt) { // 1=charm , 2=beauty hfparticle = particle; - if (origin == 1) { + if (origin == RecoDecay::OriginType::Prompt) { return JetTaggingSpecies::charm; } - if (origin == 2) { + if (origin == RecoDecay::OriginType::NonPrompt) { return JetTaggingSpecies::beauty; } } @@ -993,11 +993,11 @@ void analyzeJetSVInfo4ML(AnalysisJet const& myJet, AnyTracks const& /*allTracks* // Looping over the track info and putting them in the input vector template -void analyzeJetTrackInfo4ML(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, SecondaryVertices const& /*allSVs*/, std::vector& tracksParams, float trackPtMin = 0.5) +void analyzeJetTrackInfo4ML(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, SecondaryVertices const& /*allSVs*/, std::vector& tracksParams, float trackPtMin = 0.5, float trackDcaXYMax = 10.0, float trackDcaZMax = 10.0) { for (const auto& constituent : analysisJet.template tracks_as()) { - if (constituent.pt() < trackPtMin) { + if (constituent.pt() < trackPtMin || !trackAcceptanceWithDca(constituent, trackDcaXYMax, trackDcaZMax)) { continue; } @@ -1026,11 +1026,11 @@ void analyzeJetTrackInfo4ML(AnalysisJet const& analysisJet, AnyTracks const& /*a // Looping over the track info and putting them in the input vector without using any SV info template -void analyzeJetTrackInfo4MLnoSV(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, std::vector& tracksParams, float trackPtMin = 0.5) +void analyzeJetTrackInfo4MLnoSV(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, std::vector& tracksParams, float trackPtMin = 0.5, float trackDcaXYMax = 10.0, float trackDcaZMax = 10.0) { for (const auto& constituent : analysisJet.template tracks_as()) { - if (constituent.pt() < trackPtMin) { + if (constituent.pt() < trackPtMin || !trackAcceptanceWithDca(constituent, trackDcaXYMax, trackDcaZMax)) { continue; } diff --git a/PWGJE/Core/MlResponseHfTagging.h b/PWGJE/Core/MlResponseHfTagging.h index 462cb661505..2580c5296df 100644 --- a/PWGJE/Core/MlResponseHfTagging.h +++ b/PWGJE/Core/MlResponseHfTagging.h @@ -113,6 +113,11 @@ class MlResponseHfTagging : public MlResponse /// @return A vector of input shapes std::vector> getInputShape() const { return this->mModels[0].getInputShapes(); } + /// @brief Method to get the output shape of a model + /// \param imod is the index of the model + /// @return number of output nodes + int getOutputNodes(int imod = 0) const { return this->mModels[imod].getNumOutputNodes(); } + /// Method to fill the inputs of jet, tracks and secondary vertices /// \param jet is the b-jet candidate /// \param tracks is the vector of tracks associated to the jet diff --git a/PWGJE/TableProducer/jetTaggerHF.cxx b/PWGJE/TableProducer/jetTaggerHF.cxx index 0ab545fe414..2a37871c970 100644 --- a/PWGJE/TableProducer/jetTaggerHF.cxx +++ b/PWGJE/TableProducer/jetTaggerHF.cxx @@ -94,6 +94,7 @@ struct JetTaggerHFTask { Configurable> cutsMl{"cutsMl", {DefaultCutsMl[0], 1, 2, {"pT bin 0"}, {"score for default b-jet tagging", "uncer 1"}}, "ML selections per pT bin"}; Configurable nClassesMl{"nClassesMl", 2, "Number of classes in ML model"}; Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"feature1", "feature2"}, "Names of ML model input features"}; + Configurable useDb{"useDb", false, "Flag to use DB for ML model instead of the score"}; Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable> modelPathsCCDB{"modelPathsCCDB", std::vector{"Users/h/hahassan"}, "Paths of models on CCDB"}; @@ -297,16 +298,16 @@ struct JetTaggerHFTask { std::map metadata; resoFuncMatch = resoFuncMatching; - const int IPmethodResolutionFunctionSize = 7; + const int mIPmethodResolutionFunctionSize = 7; auto loadCCDBforIP = [&](const std::vector& paths, std::vector& targetVec, const std::string& name) { - if (paths.size() != IPmethodResolutionFunctionSize) { + if (paths.size() != mIPmethodResolutionFunctionSize) { usepTcategorize.value = false; LOG(info) << name << " does not have 7 entries. Disabling pT categorization (usepTcategorize = false)."; resoFuncMatch = 0; return; } - for (int i = 0; i < IPmethodResolutionFunctionSize; i++) { + for (int i = 0; i < mIPmethodResolutionFunctionSize; i++) { targetVec.push_back(ccdbApi.retrieveFromTFileAny(paths[i], metadata, -1)); } }; @@ -335,7 +336,7 @@ struct JetTaggerHFTask { } maxOrder = numCount + 1; // 0: untagged, >1 : N ordering - const int IPmethodNumOfParameters = 9; + const int mIPmethodNumOfParameters = 9; // Set up the resolution function switch (resoFuncMatch) { @@ -383,7 +384,7 @@ struct JetTaggerHFTask { for (size_t j = 0; j < resoFuncIncCCDB.size(); j++) { std::vector params; if (resoFuncIncCCDB[j]) { - for (int i = 0; i < IPmethodNumOfParameters; i++) { + for (int i = 0; i < mIPmethodNumOfParameters; i++) { params.emplace_back(resoFuncIncCCDB[j]->GetParameter(i)); } } @@ -397,7 +398,7 @@ struct JetTaggerHFTask { for (size_t j = 0; j < resoFuncBeautyCCDB.size(); j++) { std::vector params; if (resoFuncBeautyCCDB[j]) { - for (int i = 0; i < IPmethodNumOfParameters; i++) { + for (int i = 0; i < mIPmethodNumOfParameters; i++) { params.emplace_back(resoFuncBeautyCCDB[j]->GetParameter(i)); } } @@ -406,7 +407,7 @@ struct JetTaggerHFTask { for (size_t j = 0; j < resoFuncCharmCCDB.size(); j++) { std::vector params; if (resoFuncCharmCCDB[j]) { - for (int i = 0; i < IPmethodNumOfParameters; i++) { + for (int i = 0; i < mIPmethodNumOfParameters; i++) { params.emplace_back(resoFuncCharmCCDB[j]->GetParameter(i)); } } @@ -415,7 +416,7 @@ struct JetTaggerHFTask { for (size_t j = 0; j < resoFuncLfCCDB.size(); j++) { std::vector params; if (resoFuncLfCCDB[j]) { - for (int i = 0; i < IPmethodNumOfParameters; i++) { + for (int i = 0; i < mIPmethodNumOfParameters; i++) { params.emplace_back(resoFuncLfCCDB[j]->GetParameter(i)); } } @@ -427,7 +428,7 @@ struct JetTaggerHFTask { for (size_t j = 0; j < resoFuncDataCCDB.size(); j++) { std::vector params; if (resoFuncDataCCDB[j]) { - for (int i = 0; i < IPmethodNumOfParameters; i++) { + for (int i = 0; i < mIPmethodNumOfParameters; i++) { params.emplace_back(resoFuncDataCCDB[j]->GetParameter(i)); } } @@ -507,7 +508,7 @@ struct JetTaggerHFTask { std::vector svsParams; jettaggingutilities::analyzeJetSVInfo4ML(analysisJet, allTracks, allSVs, svsParams, svPtMin, svReductionFactor); - jettaggingutilities::analyzeJetTrackInfo4ML(analysisJet, allTracks, allSVs, tracksParams, trackPtMin); + jettaggingutilities::analyzeJetTrackInfo4ML(analysisJet, allTracks, allSVs, tracksParams, trackPtMin, trackDcaXYMax, trackDcaZMax); int nSVs = analysisJet.template secondaryVertices_as().size(); @@ -525,7 +526,15 @@ struct JetTaggerHFTask { bMlResponse.isSelectedMl(inputML, analysisJet.pt(), output); } - scoreML[analysisJet.globalIndex()] = output[0]; + if (bMlResponse.getOutputNodes() > 1) { + auto mDb = [](std::vector scores, float fC) { + return std::log(scores[2] / (fC * scores[1] + (1 - fC) * scores[0])); + }; + + scoreML[analysisJet.globalIndex()] = useDb ? mDb(output, fC) : output[2]; // 2 is the b-jet index + } else { + scoreML[analysisJet.globalIndex()] = output[0]; + } } } @@ -537,7 +546,7 @@ struct JetTaggerHFTask { std::vector tracksParams; std::vector svsParams; - jettaggingutilities::analyzeJetTrackInfo4MLnoSV(analysisJet, allTracks, tracksParams, trackPtMin); + jettaggingutilities::analyzeJetTrackInfo4MLnoSV(analysisJet, allTracks, tracksParams, trackPtMin, trackDcaXYMax, trackDcaZMax); jettaggingutilities::BJetParams jetparam = {analysisJet.pt(), analysisJet.eta(), analysisJet.phi(), static_cast(tracksParams.size()), 0, analysisJet.mass()}; tracksParams.resize(nJetConst); // resize to the number of inputs of the ML diff --git a/PWGJE/TableProducer/secondaryVertexReconstruction.cxx b/PWGJE/TableProducer/secondaryVertexReconstruction.cxx index de1b47cb830..c1f7598d176 100644 --- a/PWGJE/TableProducer/secondaryVertexReconstruction.cxx +++ b/PWGJE/TableProducer/secondaryVertexReconstruction.cxx @@ -72,6 +72,8 @@ struct SecondaryVertexReconstruction { Configurable ptMinTrack{"ptMinTrack", -1., "min. track pT"}; Configurable etaMinTrack{"etaMinTrack", -99999., "min. pseudorapidity"}; Configurable etaMaxTrack{"etaMaxTrack", 4., "max. pseudorapidity"}; + Configurable maxIPxy{"maxIPxy", 10, "maximum track DCA in xy plane"}; + Configurable maxIPz{"maxIPz", 10, "maximum track DCA in z direction"}; Configurable fillHistograms{"fillHistograms", true, "do validation plots"}; Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; @@ -288,7 +290,7 @@ struct SecondaryVertexReconstruction { for (size_t iprong = prongIndex; iprong < particles.size(); ++iprong) { const auto& testTrack = particles[iprong].template track_as(); - if (testTrack.pt() < ptMinTrack || testTrack.eta() < etaMinTrack || testTrack.eta() > etaMaxTrack) { + if (testTrack.pt() < ptMinTrack || testTrack.eta() < etaMinTrack || testTrack.eta() > etaMaxTrack || std::abs(testTrack.dcaXY()) > maxIPxy || std::abs(testTrack.dcaZ()) > maxIPz) { continue; } diff --git a/PWGJE/Tasks/bjetTaggingML.cxx b/PWGJE/Tasks/bjetTaggingML.cxx index 7b2638f37a5..afa34a243e6 100644 --- a/PWGJE/Tasks/bjetTaggingML.cxx +++ b/PWGJE/Tasks/bjetTaggingML.cxx @@ -58,6 +58,8 @@ struct BJetTaggingML { Configurable trackPtMax{"trackPtMax", 1000.0, "maximum track pT"}; Configurable trackEtaMin{"trackEtaMin", -0.9, "minimum track eta"}; Configurable trackEtaMax{"trackEtaMax", 0.9, "maximum track eta"}; + Configurable maxIPxy{"maxIPxy", 10, "maximum track DCA in xy plane"}; + Configurable maxIPz{"maxIPz", 10, "maximum track DCA in z direction"}; // track level configurables Configurable svPtMin{"svPtMin", 0.5, "minimum SV pT"}; @@ -68,6 +70,7 @@ struct BJetTaggingML { Configurable jetEtaMin{"jetEtaMin", -99.0, "minimum jet pseudorapidity"}; Configurable jetEtaMax{"jetEtaMax", 99.0, "maximum jet pseudorapidity"}; Configurable nJetConst{"nJetConst", 10, "maximum number of jet consistuents to be used for ML evaluation"}; + Configurable useDb{"useDb", false, "Flag whether to use the Db instead of the score for tagging"}; Configurable doDataDriven{"doDataDriven", false, "Flag whether to use fill THnSpase for data driven methods"}; @@ -86,10 +89,14 @@ struct BJetTaggingML { eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); + const AxisSpec axisDb{300, -10., 20., "#it{D}_{b}"}; + const AxisSpec axisScore{120, -0.1, 1.1, "Score"}; + const AxisSpec axisLogScore{120, 0., 30, "-log(1-score)"}; + registry.add("h_vertexZ", "Vertex Z;#it{Z} (cm)", {HistType::kTH1F, {{40, -20.0, 20.0}}}); - registry.add("h2_score_jetpT", "ML scores for inclusive jets;#it{p}_{T,jet} (GeV/#it{c});Score", {HistType::kTH2F, {{200, 0., 200.}, {120, -0.1, 1.1}}}); - registry.add("h2_logscore_jetpT", "ML scores for inclusive jets;#it{p}_{T,jet} (GeV/#it{c});- log(1 - Score)", {HistType::kTH2F, {{200, 0., 200.}, {120, 0, 30}}}); + registry.add("h2_score_jetpT", "ML scores for inclusive jets;#it{p}_{T,jet} (GeV/#it{c});Score", {HistType::kTH2F, {{200, 0., 200.}, useDb ? axisDb : axisScore}}); + registry.add("h2_logscore_jetpT", "ML scores for inclusive jets;#it{p}_{T,jet} (GeV/#it{c});- log(1 - Score)", {HistType::kTH2F, {{200, 0., 200.}, axisLogScore}}); registry.add("h2_nTracks_jetpT", "Number of tracks;#it{p}_{T,jet} (GeV/#it{c});nTracks", {HistType::kTH2F, {{200, 0., 200.}, {100, 0, 100.0}}}); registry.add("h2_nSV_jetpT", "Number of secondary vertices;#it{p}_{T,jet} (GeV/#it{c});nSVs", {HistType::kTH2F, {{200, 0., 200.}, {250, 0, 250.0}}}); @@ -102,18 +109,18 @@ struct BJetTaggingML { registry.add("h2_SVMass_jetpT", "Secondary vertex mass;#it{p}_{T,jet} (GeV/#it{c});#it{m}_{SV} (GeV/#it{c}^{2})", {HistType::kTH2F, {{200, 0., 200.}, {50, 0, 10}}}); if (doDataDriven) { - registry.add("hSparse_Incljets", "Inclusive jets Info;#it{p}_{T,jet} (GeV/#it{c});Score;-log(1-score);#it{m}_{jet} (GeV/#it{c}^{2});-log(JP);#it{m}_{SV} (GeV/#it{c}^{2});SVfE;", {HistType::kTHnSparseF, {{200, 0., 200.}, {120, -0.1, 1.1}, {120, 0, 30}, {50, 0, 50}, {375, 0, 30}, {50, 0, 10}, {50, 0, 1}}}); + registry.add("hSparse_Incljets", "Inclusive jets Info;#it{p}_{T,jet} (GeV/#it{c});Score;-log(1-score);#it{m}_{jet} (GeV/#it{c}^{2});-log(JP);#it{m}_{SV} (GeV/#it{c}^{2});SVfE;", {HistType::kTHnSparseF, {{200, 0., 200.}, useDb ? axisDb : axisScore, axisLogScore, {50, 0, 50}, {375, 0, 30}, {50, 0, 10}, {50, 0, 1}}}); if (doprocessMCJets) { - registry.add("hSparse_bjets", "Tagged b-jets Info;#it{p}_{T,jet} (GeV/#it{c});Score;-log(1-score);#it{m}_{jet} (GeV/#it{c}^{2});-log(JP);#it{m}_{SV} (GeV/#it{c}^{2});SVfE;", {HistType::kTHnSparseF, {{200, 0., 200.}, {120, -0.1, 1.1}, {120, 0, 30}, {50, 0, 50}, {375, 0, 30}, {50, 0, 10}, {50, 0, 1}}}); - registry.add("hSparse_cjets", "Tagged c-jets Info;#it{p}_{T,jet} (GeV/#it{c});Score;-log(1-score);#it{m}_{jet} (GeV/#it{c}^{2});-log(JP);#it{m}_{SV} (GeV/#it{c}^{2});SVfE;", {HistType::kTHnSparseF, {{200, 0., 200.}, {120, -0.1, 1.1}, {120, 0, 30}, {50, 0, 50}, {375, 0, 30}, {50, 0, 10}, {50, 0, 1}}}); - registry.add("hSparse_lfjets", "Tagged lf-jets Info;#it{p}_{T,jet} (GeV/#it{c});Score;-log(1-score);#it{m}_{jet} (GeV/#it{c}^{2});-log(JP);#it{m}_{SV} (GeV/#it{c}^{2});SVfE;", {HistType::kTHnSparseF, {{200, 0., 200.}, {120, -0.1, 1.1}, {120, 0, 30}, {50, 0, 50}, {375, 0, 30}, {50, 0, 10}, {50, 0, 1}}}); + registry.add("hSparse_bjets", "Tagged b-jets Info;#it{p}_{T,jet} (GeV/#it{c});Score;-log(1-score);#it{m}_{jet} (GeV/#it{c}^{2});-log(JP);#it{m}_{SV} (GeV/#it{c}^{2});SVfE;", {HistType::kTHnSparseF, {{200, 0., 200.}, useDb ? axisDb : axisScore, axisLogScore, {50, 0, 50}, {375, 0, 30}, {50, 0, 10}, {50, 0, 1}}}); + registry.add("hSparse_cjets", "Tagged c-jets Info;#it{p}_{T,jet} (GeV/#it{c});Score;-log(1-score);#it{m}_{jet} (GeV/#it{c}^{2});-log(JP);#it{m}_{SV} (GeV/#it{c}^{2});SVfE;", {HistType::kTHnSparseF, {{200, 0., 200.}, useDb ? axisDb : axisScore, axisLogScore, {50, 0, 50}, {375, 0, 30}, {50, 0, 10}, {50, 0, 1}}}); + registry.add("hSparse_lfjets", "Tagged lf-jets Info;#it{p}_{T,jet} (GeV/#it{c});Score;-log(1-score);#it{m}_{jet} (GeV/#it{c}^{2});-log(JP);#it{m}_{SV} (GeV/#it{c}^{2});SVfE;", {HistType::kTHnSparseF, {{200, 0., 200.}, useDb ? axisDb : axisScore, axisLogScore, {50, 0, 50}, {375, 0, 30}, {50, 0, 10}, {50, 0, 1}}}); } } if (doprocessMCJets) { - registry.add("h2_score_jetpT_bjet", "ML scores for b-jets;#it{p}_{T,jet} (GeV/#it{c});Score", {HistType::kTH2F, {{200, 0., 200.}, {120, -0.1, 1.1}}}); - registry.add("h2_logscore_jetpT_bjet", "ML scores for b-jets;#it{p}_{T,jet} (GeV/#it{c});- log(1 - Score)", {HistType::kTH2F, {{200, 0., 200.}, {120, 0, 30}}}); + registry.add("h2_score_jetpT_bjet", "ML scores for b-jets;#it{p}_{T,jet} (GeV/#it{c});Score", {HistType::kTH2F, {{200, 0., 200.}, useDb ? axisDb : axisScore}}); + registry.add("h2_logscore_jetpT_bjet", "ML scores for b-jets;#it{p}_{T,jet} (GeV/#it{c});- log(1 - Score)", {HistType::kTH2F, {{200, 0., 200.}, axisLogScore}}); registry.add("h2_SIPs2D_jetpT_bjet", "2D IP significance b-jets;#it{p}_{T,jet} (GeV/#it{c});IPs", {HistType::kTH2F, {{200, 0., 200.}, {100, -50.0, 50.0}}}); registry.add("h2_SIPs3D_jetpT_bjet", "3D IP significance b-jets;#it{p}_{T,jet} (GeV/#it{c});IPs", {HistType::kTH2F, {{200, 0., 200.}, {100, -50.0, 50.0}}}); registry.add("h2_LxyS_jetpT_bjet", "Decay length in XY b-jets;#it{p}_{T,jet} (GeV/#it{c});S#it{L}_{xy}", {HistType::kTH2F, {{200, 0., 200.}, {100, 0., 100.0}}}); @@ -121,8 +128,8 @@ struct BJetTaggingML { registry.add("h2_jetMass_jetpT_bjet", "Jet mass b-jets;#it{p}_{T,jet} (GeV/#it{c});#it{m}_{jet} (GeV/#it{c}^{2})", {HistType::kTH2F, {{200, 0., 200.}, {50, 0, 50.0}}}); registry.add("h2_SVMass_jetpT_bjet", "Secondary vertex mass b-jets;#it{p}_{T,jet} (GeV/#it{c});#it{m}_{SV} (GeV/#it{c}^{2})", {HistType::kTH2F, {{200, 0., 200.}, {50, 0, 10.0}}}); - registry.add("h2_score_jetpT_cjet", "ML scores for c-jets;#it{p}_{T,jet} (GeV/#it{c});Score", {HistType::kTH2F, {{200, 0., 200.}, {120, -0.1, 1.1}}}); - registry.add("h2_logscore_jetpT_cjet", "ML scores for c-jets;#it{p}_{T,jet} (GeV/#it{c});- log(1 - Score)", {HistType::kTH2F, {{200, 0., 200.}, {120, 0, 30}}}); + registry.add("h2_score_jetpT_cjet", "ML scores for c-jets;#it{p}_{T,jet} (GeV/#it{c});Score", {HistType::kTH2F, {{200, 0., 200.}, useDb ? axisDb : axisScore}}); + registry.add("h2_logscore_jetpT_cjet", "ML scores for c-jets;#it{p}_{T,jet} (GeV/#it{c});- log(1 - Score)", {HistType::kTH2F, {{200, 0., 200.}, axisLogScore}}); registry.add("h2_SIPs2D_jetpT_cjet", "2D IP significance c-jets;#it{p}_{T,jet} (GeV/#it{c});IPs", {HistType::kTH2F, {{200, 0., 200.}, {100, -50.0, 50.0}}}); registry.add("h2_SIPs3D_jetpT_cjet", "3D IP significance c-jets;#it{p}_{T,jet} (GeV/#it{c});IPs", {HistType::kTH2F, {{200, 0., 200.}, {100, -50.0, 50.0}}}); registry.add("h2_LxyS_jetpT_cjet", "Decay length in XY c-jets;#it{p}_{T,jet} (GeV/#it{c});S#it{L}_{xy}", {HistType::kTH2F, {{200, 0., 200.}, {100, 0., 100.0}}}); @@ -130,8 +137,8 @@ struct BJetTaggingML { registry.add("h2_jetMass_jetpT_cjet", "Jet mass c-jets;#it{p}_{T,jet} (GeV/#it{c});#it{m}_{jet} (GeV/#it{c}^{2})", {HistType::kTH2F, {{200, 0., 200.}, {50, 0, 50.0}}}); registry.add("h2_SVMass_jetpT_cjet", "Secondary vertex mass c-jets;#it{p}_{T,jet} (GeV/#it{c});#it{m}_{SV} (GeV/#it{c}^{2})", {HistType::kTH2F, {{200, 0., 200.}, {50, 0, 10.0}}}); - registry.add("h2_score_jetpT_lfjet", "ML scores for lf-jets;#it{p}_{T,jet} (GeV/#it{c});Score", {HistType::kTH2F, {{200, 0., 200.}, {120, -0.1, 1.1}}}); - registry.add("h2_logscore_jetpT_lfjet", "ML scores for lf-jets;#it{p}_{T,jet} (GeV/#it{c});- log(1 - Score)", {HistType::kTH2F, {{200, 0., 200.}, {120, 0, 30}}}); + registry.add("h2_score_jetpT_lfjet", "ML scores for lf-jets;#it{p}_{T,jet} (GeV/#it{c});Score", {HistType::kTH2F, {{200, 0., 200.}, useDb ? axisDb : axisScore}}); + registry.add("h2_logscore_jetpT_lfjet", "ML scores for lf-jets;#it{p}_{T,jet} (GeV/#it{c});- log(1 - Score)", {HistType::kTH2F, {{200, 0., 200.}, axisLogScore}}); registry.add("h2_SIPs2D_jetpT_lfjet", "2D IP significance lf-jet;#it{p}_{T,jet} (GeV/#it{c});IPs", {HistType::kTH2F, {{200, 0., 200.}, {100, -50.0, 50.0}}}); registry.add("h2_SIPs3D_jetpT_lfjet", "3D IP significance lf-jet;#it{p}_{T,jet} (GeV/#it{c});IPs", {HistType::kTH2F, {{200, 0., 200.}, {100, -50.0, 50.0}}}); registry.add("h2_LxyS_jetpT_lfjet", "Decay length in XY lf-jet;#it{p}_{T,jet} (GeV/#it{c});S#it{L}_{xy}", {HistType::kTH2F, {{200, 0., 200.}, {100, 0., 100.0}}}); @@ -205,11 +212,11 @@ struct BJetTaggingML { registry.fill(HIST("h2_SVMass_jetpT"), myJet.pt(), massSV, eventweight); if (doprocessMCJets) { - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h2_LxyS_jetpT_bjet"), myJet.pt(), candSV.decayLengthXY() / candSV.errorDecayLengthXY(), eventweight); registry.fill(HIST("h2_Dispersion_jetpT_bjet"), myJet.pt(), candSV.dispersion(), eventweight); registry.fill(HIST("h2_SVMass_jetpT_bjet"), myJet.pt(), massSV, eventweight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h2_LxyS_jetpT_cjet"), myJet.pt(), candSV.decayLengthXY() / candSV.errorDecayLengthXY(), eventweight); registry.fill(HIST("h2_Dispersion_jetpT_cjet"), myJet.pt(), candSV.dispersion(), eventweight); registry.fill(HIST("h2_SVMass_jetpT_cjet"), myJet.pt(), massSV, eventweight); @@ -229,7 +236,7 @@ struct BJetTaggingML { for (const auto& constituent : analysisJet.template tracks_as()) { - if (constituent.pt() < trackPtMin) { + if (constituent.pt() < trackPtMin || !jettaggingutilities::trackAcceptanceWithDca(constituent, maxIPxy, maxIPz)) { continue; } @@ -243,10 +250,10 @@ struct BJetTaggingML { registry.fill(HIST("h2_SIPs3D_jetpT"), analysisJet.pt(), sign * std::abs(constituent.dcaXYZ()) / constituent.sigmadcaXYZ(), eventweight); if (doprocessMCJets) { - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h2_SIPs2D_jetpT_bjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXY()) / constituent.sigmadcaXY(), eventweight); registry.fill(HIST("h2_SIPs3D_jetpT_bjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXYZ()) / constituent.sigmadcaXYZ(), eventweight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h2_SIPs2D_jetpT_cjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXY()) / constituent.sigmadcaXY(), eventweight); registry.fill(HIST("h2_SIPs3D_jetpT_cjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXYZ()) / constituent.sigmadcaXYZ(), eventweight); } else { @@ -307,12 +314,14 @@ struct BJetTaggingML { registry.fill(HIST("h2_nSV_jetpT"), analysisJet.pt(), nSVs < 250 ? nSVs : 249); registry.fill(HIST("h2_score_jetpT"), analysisJet.pt(), analysisJet.scoreML()); - registry.fill(HIST("h2_logscore_jetpT"), analysisJet.pt(), -1 * std::log(1 - analysisJet.scoreML())); + if (!useDb) { + registry.fill(HIST("h2_logscore_jetpT"), analysisJet.pt(), -1 * std::log(1 - analysisJet.scoreML())); + } registry.fill(HIST("h2_jetMass_jetpT"), analysisJet.pt(), analysisJet.mass()); if (doDataDriven) { - registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), -1 * std::log(1 - analysisJet.scoreML()), analysisJet.mass(), -1 * std::log(analysisJet.jetProb()), svsParams[0].svMass, svsParams[0].svfE); + registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), useDb ? 0 : -1 * std::log(1 - analysisJet.scoreML()), analysisJet.mass(), -1 * std::log(analysisJet.jetProb()), svsParams[0].svMass, svsParams[0].svfE); } } } @@ -369,33 +378,41 @@ struct BJetTaggingML { registry.fill(HIST("h2_nSV_jetpT"), analysisJet.pt(), nSVs < 250 ? nSVs : 249); registry.fill(HIST("h2_score_jetpT"), analysisJet.pt(), analysisJet.scoreML(), eventWeight); - registry.fill(HIST("h2_logscore_jetpT"), analysisJet.pt(), -1 * std::log(1 - analysisJet.scoreML()), eventWeight); + if (!useDb) { + registry.fill(HIST("h2_logscore_jetpT"), analysisJet.pt(), -1 * std::log(1 - analysisJet.scoreML()), eventWeight); + } registry.fill(HIST("h2_jetMass_jetpT"), analysisJet.pt(), analysisJet.mass(), eventWeight); if (doDataDriven) { - registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), -1 * std::log(1 - analysisJet.scoreML()), analysisJet.mass(), -1 * std::log(analysisJet.jetProb()), svsParams[0].svMass, svsParams[0].svfE, eventWeight); - if (jetFlavor == 2) { - registry.fill(HIST("hSparse_bjets"), analysisJet.pt(), analysisJet.scoreML(), -1 * std::log(1 - analysisJet.scoreML()), analysisJet.mass(), -1 * std::log(analysisJet.jetProb()), svsParams[0].svMass, svsParams[0].svfE, eventWeight); - } else if (jetFlavor == 1) { - registry.fill(HIST("hSparse_cjets"), analysisJet.pt(), analysisJet.scoreML(), -1 * std::log(1 - analysisJet.scoreML()), analysisJet.mass(), -1 * std::log(analysisJet.jetProb()), svsParams[0].svMass, svsParams[0].svfE, eventWeight); + registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), useDb ? 0 : -1 * std::log(1 - analysisJet.scoreML()), analysisJet.mass(), -1 * std::log(analysisJet.jetProb()), svsParams[0].svMass, svsParams[0].svfE, eventWeight); + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("hSparse_bjets"), analysisJet.pt(), analysisJet.scoreML(), useDb ? 0 : -1 * std::log(1 - analysisJet.scoreML()), analysisJet.mass(), -1 * std::log(analysisJet.jetProb()), svsParams[0].svMass, svsParams[0].svfE, eventWeight); + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("hSparse_cjets"), analysisJet.pt(), analysisJet.scoreML(), useDb ? 0 : -1 * std::log(1 - analysisJet.scoreML()), analysisJet.mass(), -1 * std::log(analysisJet.jetProb()), svsParams[0].svMass, svsParams[0].svfE, eventWeight); } else { - registry.fill(HIST("hSparse_lfjets"), analysisJet.pt(), analysisJet.scoreML(), -1 * std::log(1 - analysisJet.scoreML()), analysisJet.mass(), -1 * std::log(analysisJet.jetProb()), svsParams[0].svMass, svsParams[0].svfE, eventWeight); + registry.fill(HIST("hSparse_lfjets"), analysisJet.pt(), analysisJet.scoreML(), useDb ? 0 : -1 * std::log(1 - analysisJet.scoreML()), analysisJet.mass(), -1 * std::log(analysisJet.jetProb()), svsParams[0].svMass, svsParams[0].svfE, eventWeight); } } - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h2_score_jetpT_bjet"), analysisJet.pt(), analysisJet.scoreML(), eventWeight); - registry.fill(HIST("h2_logscore_jetpT_bjet"), analysisJet.pt(), -1 * std::log(1 - analysisJet.scoreML()), eventWeight); + if (!useDb) { + registry.fill(HIST("h2_logscore_jetpT_bjet"), analysisJet.pt(), -1 * std::log(1 - analysisJet.scoreML()), eventWeight); + } registry.fill(HIST("h2_jetMass_jetpT_bjet"), analysisJet.pt(), analysisJet.mass(), eventWeight); registry.fill(HIST("h_jetpT_detector_bjet"), analysisJet.pt(), eventWeight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h2_score_jetpT_cjet"), analysisJet.pt(), analysisJet.scoreML(), eventWeight); - registry.fill(HIST("h2_logscore_jetpT_cjet"), analysisJet.pt(), -1 * std::log(1 - analysisJet.scoreML()), eventWeight); + if (!useDb) { + registry.fill(HIST("h2_logscore_jetpT_cjet"), analysisJet.pt(), -1 * std::log(1 - analysisJet.scoreML()), eventWeight); + } registry.fill(HIST("h2_jetMass_jetpT_cjet"), analysisJet.pt(), analysisJet.mass(), eventWeight); registry.fill(HIST("h_jetpT_detector_cjet"), analysisJet.pt(), eventWeight); } else { registry.fill(HIST("h2_score_jetpT_lfjet"), analysisJet.pt(), analysisJet.scoreML(), eventWeight); - registry.fill(HIST("h2_logscore_jetpT_lfjet"), analysisJet.pt(), -1 * std::log(1 - analysisJet.scoreML()), eventWeight); + if (!useDb) { + registry.fill(HIST("h2_logscore_jetpT_lfjet"), analysisJet.pt(), -1 * std::log(1 - analysisJet.scoreML()), eventWeight); + } registry.fill(HIST("h2_jetMass_jetpT_lfjet"), analysisJet.pt(), analysisJet.mass(), eventWeight); registry.fill(HIST("h_jetpT_detector_lfjet"), analysisJet.pt(), eventWeight); } @@ -405,9 +422,9 @@ struct BJetTaggingML { continue; } - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_bjet"), analysisJet.pt(), mcpjet.pt(), eventWeight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_cjet"), analysisJet.pt(), mcpjet.pt(), eventWeight); } else { registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_lfjet"), analysisJet.pt(), mcpjet.pt(), eventWeight); @@ -440,9 +457,9 @@ struct BJetTaggingML { registry.fill(HIST("h_jetpT_particle_DetColl"), mcpjet.pt(), eventWeight); - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h_jetpT_particle_DetColl_bjet"), mcpjet.pt(), eventWeight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h_jetpT_particle_DetColl_cjet"), mcpjet.pt(), eventWeight); } else { registry.fill(HIST("h_jetpT_particle_DetColl_lfjet"), mcpjet.pt(), eventWeight); @@ -479,9 +496,9 @@ struct BJetTaggingML { int8_t jetFlavor = mcpjet.origin(); - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h_jetpT_particle_bjet"), mcpjet.pt(), eventWeight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h_jetpT_particle_cjet"), mcpjet.pt(), eventWeight); } else { registry.fill(HIST("h_jetpT_particle_lfjet"), mcpjet.pt(), eventWeight); diff --git a/PWGJE/Tasks/bjetTreeCreator.cxx b/PWGJE/Tasks/bjetTreeCreator.cxx index 4508600962c..ba7c8e49c6d 100644 --- a/PWGJE/Tasks/bjetTreeCreator.cxx +++ b/PWGJE/Tasks/bjetTreeCreator.cxx @@ -207,6 +207,9 @@ struct BJetTreeCreator { Configurable trackEtaMin{"trackEtaMin", -0.9, "minimum track eta"}; Configurable trackEtaMax{"trackEtaMax", 0.9, "maximum track eta"}; + Configurable maxIPxy{"maxIPxy", 10, "maximum track DCA in xy plane"}; + Configurable maxIPz{"maxIPz", 10, "maximum track DCA in z direction"}; + Configurable useQuarkDef{"useQuarkDef", true, "Flag whether to use quarks or hadrons for determining the jet flavor"}; // track level configurables @@ -297,20 +300,20 @@ struct BJetTreeCreator { //+jet registry.add("h_jet_pt", "jet_pt;#it{p}_{T}^{ch jet} (GeV/#it{c});Entries", {HistType::kTH1F, {{200, 0., 200.}}}); registry.add("h_jet_eta", "jet_eta;#it{#eta}_{ch jet};Entries", {HistType::kTH1F, {{200, -2., 2.}}}); - registry.add("h_jet_phi", "jet_phi;#it{#phi}_{ch jet};Entries", {HistType::kTH1F, {{200, 0., 2. * M_PI}}}); + registry.add("h_jet_phi", "jet_phi;#it{#phi}_{ch jet};Entries", {HistType::kTH1F, {{200, 0., o2::constants::math::TwoPI}}}); registry.add("h_jet_flav", "jet_flav;jet flavor;Entries", {HistType::kTH1F, {{4, 0., 4.}}}); registry.add("h_n_trks", "n_trks;#it{n}_{tracks};Entries", {HistType::kTH1F, {{50, 0., 50.}}}); registry.add("h_jet_mass", "jet_mass;#it{m}_{jet} (GeV/#it{c}^2);Entries", {HistType::kTH1F, {{200, 0., 50.}}}); - auto h_jet_flav = registry.get(HIST("h_jet_flav")); - h_jet_flav->GetXaxis()->SetBinLabel(1, "no mcparticle"); // 0 - h_jet_flav->GetXaxis()->SetBinLabel(2, "c-jet"); // 1 - h_jet_flav->GetXaxis()->SetBinLabel(3, "b-jet"); // 2 - h_jet_flav->GetXaxis()->SetBinLabel(4, "lf-jet"); // 3 + auto hJetFlavor = registry.get(HIST("h_jet_flav")); + hJetFlavor->GetXaxis()->SetBinLabel(1, "no mcparticle"); // 0 + hJetFlavor->GetXaxis()->SetBinLabel(2, "c-jet"); // 1 + hJetFlavor->GetXaxis()->SetBinLabel(3, "b-jet"); // 2 + hJetFlavor->GetXaxis()->SetBinLabel(4, "lf-jet"); // 3 registry.add("h_n_vertices", "n_vertices;#it{n}_{vertex};Entries", {HistType::kTH1F, {{50, 0., 50.}}}); //+trk registry.add("h_trk_pt", "trk_pt;#it{p}_{T} (GeV/#it{c});Entries", {HistType::kTH1F, {{200, 0., 100.}}}); registry.add("h_trk_eta", "trk_eta;#it{#eta};Entries", {HistType::kTH1F, {{200, -2., 2.}}}); - registry.add("h_trk_phi", "trk_phi;#it{#phi};Entries", {HistType::kTH1F, {{200, 0., 2. * M_PI}}}); + registry.add("h_trk_phi", "trk_phi;#it{#phi};Entries", {HistType::kTH1F, {{200, 0., o2::constants::math::TwoPI}}}); registry.add("h_trk_charge", "trk_charge;#it{q};Entries", {HistType::kTH1F, {{3, -1.5, 1.5}}}); registry.add("h_trk_dcaxy", "trk_dcaxy;#it{DCA}_{xy} (cm);Entries", {HistType::kTH1F, {{200, -0.1, 0.1}}}); registry.add("h_trk_dcaxyz", "trk_dcaxyz;#it{DCA}_{xyz} (cm);Entries", {HistType::kTH1F, {{200, -0.1, 0.1}}}); @@ -324,12 +327,12 @@ struct BJetTreeCreator { registry.add("h2_trk_jtrackpt_vs_origtrackpt", "JTracks::pt vs Tracks::pt", {HistType::kTH2F, {{200, 0., 100.}, {200, 0., 100.}}}); registry.add("h_trk_vtx_index", "trk_vtx_index;Vertex index;Entries", {HistType::kTH1F, {{20, 0., 20.}}}); registry.add("h_trk_origin", "trk_origin;Track origin;Entries", {HistType::kTH1F, {{5, 0., 5.}}}); - auto h_trk_origin = registry.get(HIST("h_trk_origin")); - h_trk_origin->GetXaxis()->SetBinLabel(1, "NotPhysPrim"); - h_trk_origin->GetXaxis()->SetBinLabel(2, "Charm"); - h_trk_origin->GetXaxis()->SetBinLabel(3, "Beauty"); - h_trk_origin->GetXaxis()->SetBinLabel(4, "Primary"); - h_trk_origin->GetXaxis()->SetBinLabel(5, "OtherSecondary"); + auto hTrackOrigin = registry.get(HIST("h_trk_origin")); + hTrackOrigin->GetXaxis()->SetBinLabel(1, "NotPhysPrim"); + hTrackOrigin->GetXaxis()->SetBinLabel(2, "Charm"); + hTrackOrigin->GetXaxis()->SetBinLabel(3, "Beauty"); + hTrackOrigin->GetXaxis()->SetBinLabel(4, "Primary"); + hTrackOrigin->GetXaxis()->SetBinLabel(5, "OtherSecondary"); } } @@ -404,11 +407,11 @@ struct BJetTreeCreator { registry.fill(HIST("h2_SVMass_jetpT"), myJet.pt(), massSV, eventweight); if (doprocessMCJets) { - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h2_LxyS_jetpT_bjet"), myJet.pt(), candSV.decayLengthXY() / candSV.errorDecayLengthXY(), eventweight); registry.fill(HIST("h2_Dispersion_jetpT_bjet"), myJet.pt(), candSV.chi2PCA(), eventweight); registry.fill(HIST("h2_SVMass_jetpT_bjet"), myJet.pt(), massSV, eventweight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h2_LxyS_jetpT_cjet"), myJet.pt(), candSV.decayLengthXY() / candSV.errorDecayLengthXY(), eventweight); registry.fill(HIST("h2_Dispersion_jetpT_cjet"), myJet.pt(), candSV.chi2PCA(), eventweight); registry.fill(HIST("h2_SVMass_jetpT_cjet"), myJet.pt(), massSV, eventweight); @@ -427,9 +430,9 @@ struct BJetTreeCreator { void analyzeJetTrackInfo(AnyCollision const& /*collision*/, AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, SecondaryVertices const& /*allSVs*/, std::vector& trackIndices, int jetFlavor = 0, double eventweight = 1.0) { - for (auto& constituent : analysisJet.template tracks_as()) { + for (const auto& constituent : analysisJet.template tracks_as()) { - if (constituent.pt() < trackPtMin) { + if (constituent.pt() < trackPtMin || !jettaggingutilities::trackAcceptanceWithDca(constituent, maxIPxy, maxIPz)) { continue; } @@ -437,11 +440,11 @@ struct BJetTreeCreator { double dotProduct = RecoDecay::dotProd(std::array{analysisJet.px(), analysisJet.py(), analysisJet.pz()}, std::array{constituent.px(), constituent.py(), constituent.pz()}); int sign = jettaggingutilities::getGeoSign(analysisJet, constituent); - float RClosestSV = 10.; + float dRClosestSV = 10.; for (const auto& candSV : analysisJet.template secondaryVertices_as()) { double deltaRTrackSV = jetutilities::deltaR(constituent, candSV); - if (deltaRTrackSV < RClosestSV) { - RClosestSV = deltaRTrackSV; + if (deltaRTrackSV < dRClosestSV) { + dRClosestSV = deltaRTrackSV; } } @@ -449,10 +452,10 @@ struct BJetTreeCreator { registry.fill(HIST("h2_SIPs3D_jetpT"), analysisJet.pt(), sign * std::abs(constituent.dcaXYZ()) / constituent.sigmadcaXYZ(), eventweight); if (doprocessMCJets) { - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h2_SIPs2D_jetpT_bjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXY()) / constituent.sigmadcaXY(), eventweight); registry.fill(HIST("h2_SIPs3D_jetpT_bjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXYZ()) / constituent.sigmadcaXYZ(), eventweight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h2_SIPs2D_jetpT_cjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXY()) / constituent.sigmadcaXY(), eventweight); registry.fill(HIST("h2_SIPs3D_jetpT_cjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXYZ()) / constituent.sigmadcaXYZ(), eventweight); } else { @@ -462,7 +465,7 @@ struct BJetTreeCreator { } if (produceTree) { - bjetTracksParamsTable(bjetParamsTable.lastIndex() + 1, constituent.pt(), constituent.eta(), dotProduct, dotProduct / analysisJet.p(), deltaRJetTrack, std::abs(constituent.dcaXY()) * sign, constituent.sigmadcaXY(), std::abs(constituent.dcaZ()) * sign, constituent.sigmadcaZ(), constituent.sigmadcaXYZ(), constituent.p() / analysisJet.p(), RClosestSV); + bjetTracksParamsTable(bjetParamsTable.lastIndex() + 1, constituent.pt(), constituent.eta(), dotProduct, dotProduct / analysisJet.p(), deltaRJetTrack, std::abs(constituent.dcaXY()) * sign, constituent.sigmadcaXY(), std::abs(constituent.dcaZ()) * sign, constituent.sigmadcaZ(), constituent.sigmadcaXYZ(), constituent.p() / analysisJet.p(), dRClosestSV); } trackIndices.push_back(bjetTracksParamsTable.lastIndex()); } @@ -472,7 +475,7 @@ struct BJetTreeCreator { void analyzeJetTrackInfoForGNN(AnyCollision const& /*collision*/, AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, AnyOriginalTracks const&, std::vector& trackIndices, int jetFlavor = 0, double eventweight = 1.0, TrackLabelMap* trkLabels = nullptr) { int trkIdx = -1; - for (auto& constituent : analysisJet.template tracks_as()) { + for (const auto& constituent : analysisJet.template tracks_as()) { trkIdx++; @@ -488,10 +491,10 @@ struct BJetTreeCreator { registry.fill(HIST("h2_SIPs3D_jetpT"), analysisJet.pt(), sign * std::abs(constituent.dcaXYZ()) / constituent.sigmadcaXYZ(), eventweight); if (doprocessMCJetsForGNN) { - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h2_SIPs2D_jetpT_bjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXY()) / constituent.sigmadcaXY(), eventweight); registry.fill(HIST("h2_SIPs3D_jetpT_bjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXYZ()) / constituent.sigmadcaXYZ(), eventweight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h2_SIPs2D_jetpT_cjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXY()) / constituent.sigmadcaXY(), eventweight); registry.fill(HIST("h2_SIPs3D_jetpT_cjet"), analysisJet.pt(), sign * std::abs(constituent.dcaXYZ()) / constituent.sigmadcaXYZ(), eventweight); } else { @@ -555,7 +558,7 @@ struct BJetTreeCreator { for (const auto& analysisJet : alljets) { bool jetIncluded = false; - for (auto jetR : jetRadiiValues) { + for (const auto& jetR : jetRadiiValues) { if (analysisJet.r() == static_cast(jetR * 100)) { jetIncluded = true; break; @@ -570,22 +573,22 @@ struct BJetTreeCreator { continue; } - std::vector tracksIndices; - std::vector SVsIndices; + std::vector indicesTracks; + std::vector indicesSVs; - analyzeJetSVInfo(analysisJet, allTracks, allSVs, SVsIndices); - analyzeJetTrackInfo(collision, analysisJet, allTracks, allSVs, tracksIndices); + analyzeJetSVInfo(analysisJet, allTracks, allSVs, indicesSVs); + analyzeJetTrackInfo(collision, analysisJet, allTracks, allSVs, indicesTracks); registry.fill(HIST("h2_jetMass_jetpT"), analysisJet.pt(), analysisJet.mass()); int nSVs = analysisJet.template secondaryVertices_as().size(); - registry.fill(HIST("h2_nTracks_jetpT"), analysisJet.pt(), tracksIndices.size()); + registry.fill(HIST("h2_nTracks_jetpT"), analysisJet.pt(), indicesTracks.size()); registry.fill(HIST("h2_nSV_jetpT"), analysisJet.pt(), nSVs < 250 ? nSVs : 249); if (produceTree) { - bjetConstituentsTable(bjetParamsTable.lastIndex() + 1, tracksIndices, SVsIndices); - bjetParamsTable(analysisJet.pt(), analysisJet.eta(), analysisJet.phi(), tracksIndices.size(), nSVs, analysisJet.mass(), 0, analysisJet.r()); + bjetConstituentsTable(bjetParamsTable.lastIndex() + 1, indicesTracks, indicesSVs); + bjetParamsTable(analysisJet.pt(), analysisJet.eta(), analysisJet.phi(), indicesTracks.size(), nSVs, analysisJet.mass(), 0, analysisJet.r()); } } } @@ -595,8 +598,8 @@ struct BJetTreeCreator { using MCPJetTable = soa::Filtered>; using FilteredCollisionMCD = soa::Filtered>; - Preslice McParticlesPerCollision = aod::jmcparticle::mcCollisionId; - Preslice McPJetsPerCollision = aod::jet::mcCollisionId; + Preslice mcParticlesPerCollision = aod::jmcparticle::mcCollisionId; + Preslice mcpJetsPerCollision = aod::jet::mcCollisionId; void processMCJets(FilteredCollisionMCD::iterator const& collision, MCDJetTable const& MCDjets, MCPJetTable const& MCPjets, JetTracksMCDwID const& allTracks, aod::JetParticles const& MCParticles, aod::MCDSecondaryVertex3Prongs const& allSVs) { @@ -606,13 +609,13 @@ struct BJetTreeCreator { registry.fill(HIST("h_vertexZ"), collision.posZ()); - auto const mcParticlesPerColl = MCParticles.sliceBy(McParticlesPerCollision, collision.mcCollisionId()); - auto const mcPJetsPerColl = MCPjets.sliceBy(McPJetsPerCollision, collision.mcCollisionId()); + auto const mcParticlesPerColl = MCParticles.sliceBy(mcParticlesPerCollision, collision.mcCollisionId()); + auto const mcPJetsPerColl = MCPjets.sliceBy(mcpJetsPerCollision, collision.mcCollisionId()); for (const auto& analysisJet : MCDjets) { bool jetIncluded = false; - for (auto jetR : jetRadiiValues) { + for (const auto& jetR : jetRadiiValues) { if (analysisJet.r() == static_cast(jetR * 100)) { jetIncluded = true; break; @@ -629,8 +632,8 @@ struct BJetTreeCreator { continue; } - std::vector tracksIndices; - std::vector SVsIndices; + std::vector indicesTracks; + std::vector indicesSVs; int16_t jetFlavor = 0; @@ -638,7 +641,7 @@ struct BJetTreeCreator { // jetFlavor = jettaggingutilities::mcdJetFromHFShower(analysisJet, allTracks, mcParticlesPerColl, (float)(analysisJet.r() / 100.)); // jetFlavor = jettaggingutilities::jetTrackFromHFShower(analysisJet, nonFilteredTracks, mcParticlesPerColl, hftrack); - for (auto& mcpjet : analysisJet.template matchedJetGeo_as()) { + for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { if (useQuarkDef) { jetFlavor = jettaggingutilities::getJetFlavor(mcpjet, mcParticlesPerColl); } else { @@ -651,20 +654,20 @@ struct BJetTreeCreator { continue; } - analyzeJetSVInfo(analysisJet, allTracks, allSVs, SVsIndices, jetFlavor, eventWeight); - analyzeJetTrackInfo(collision, analysisJet, allTracks, allSVs, tracksIndices, jetFlavor, eventWeight); + analyzeJetSVInfo(analysisJet, allTracks, allSVs, indicesSVs, jetFlavor, eventWeight); + analyzeJetTrackInfo(collision, analysisJet, allTracks, allSVs, indicesTracks, jetFlavor, eventWeight); int nSVs = analysisJet.template secondaryVertices_as().size(); registry.fill(HIST("h2_jetMass_jetpT"), analysisJet.pt(), analysisJet.mass(), eventWeight); - registry.fill(HIST("h2_nTracks_jetpT"), analysisJet.pt(), tracksIndices.size()); + registry.fill(HIST("h2_nTracks_jetpT"), analysisJet.pt(), indicesTracks.size()); registry.fill(HIST("h2_nSV_jetpT"), analysisJet.pt(), nSVs < 250 ? nSVs : 249); - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h2_jetMass_jetpT_bjet"), analysisJet.pt(), analysisJet.mass(), eventWeight); registry.fill(HIST("h_jetpT_detector_bjet"), analysisJet.pt(), eventWeight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h2_jetMass_jetpT_cjet"), analysisJet.pt(), analysisJet.mass(), eventWeight); registry.fill(HIST("h_jetpT_detector_cjet"), analysisJet.pt(), eventWeight); } else { @@ -672,15 +675,15 @@ struct BJetTreeCreator { registry.fill(HIST("h_jetpT_detector_lfjet"), analysisJet.pt(), eventWeight); } - for (auto& mcpjet : analysisJet.template matchedJetGeo_as()) { + for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { if (mcpjet.pt() > pTHatMaxMCP * pTHat) { continue; } - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_bjet"), analysisJet.pt(), mcpjet.pt(), eventWeight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_cjet"), analysisJet.pt(), mcpjet.pt(), eventWeight); } else { registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_lfjet"), analysisJet.pt(), mcpjet.pt(), eventWeight); @@ -688,8 +691,8 @@ struct BJetTreeCreator { } if (produceTree) { - bjetConstituentsTable(bjetParamsTable.lastIndex() + 1, tracksIndices, SVsIndices); - bjetParamsTable(analysisJet.pt(), analysisJet.eta(), analysisJet.phi(), tracksIndices.size(), nSVs, analysisJet.mass(), jetFlavor, analysisJet.r()); + bjetConstituentsTable(bjetParamsTable.lastIndex() + 1, indicesTracks, indicesSVs); + bjetParamsTable(analysisJet.pt(), analysisJet.eta(), analysisJet.phi(), indicesTracks.size(), nSVs, analysisJet.mass(), jetFlavor, analysisJet.r()); } } } @@ -706,13 +709,13 @@ struct BJetTreeCreator { registry.fill(HIST("h_vertexZ"), collision.posZ()); - auto const mcParticlesPerColl = MCParticles.sliceBy(McParticlesPerCollision, collision.mcCollisionId()); - auto const mcPJetsPerColl = MCPjets.sliceBy(McPJetsPerCollision, collision.mcCollisionId()); + auto const mcParticlesPerColl = MCParticles.sliceBy(mcParticlesPerCollision, collision.mcCollisionId()); + auto const mcPJetsPerColl = MCPjets.sliceBy(mcpJetsPerCollision, collision.mcCollisionId()); for (const auto& analysisJet : MCDjets) { bool jetIncluded = false; - for (auto jetR : jetRadiiValues) { + for (const auto& jetR : jetRadiiValues) { if (analysisJet.r() == static_cast(jetR * 100)) { jetIncluded = true; break; @@ -723,12 +726,12 @@ struct BJetTreeCreator { continue; } - std::vector tracksIndices; - std::vector SVsIndices; + std::vector indicesTracks; + std::vector indicesSVs; int16_t jetFlavor = 0; - for (auto& mcpjet : analysisJet.template matchedJetGeo_as()) { + for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { if (useQuarkDef) { jetFlavor = jettaggingutilities::getJetFlavor(mcpjet, mcParticlesPerColl); } else { @@ -745,24 +748,24 @@ struct BJetTreeCreator { //+ TrackLabelMap trkLabels{{"trkVtxIndex", {}}, {"trkOrigin", {}}}; int nVertices = jettaggingutilities::vertexClustering(collision.template mcCollision_as(), analysisJet, allTracks, MCParticles, origParticles, trkLabels, true, vtxRes, trackPtMin); - analyzeJetTrackInfoForGNN(collision, analysisJet, allTracks, origTracks, tracksIndices, jetFlavor, eventWeight, &trkLabels); + analyzeJetTrackInfoForGNN(collision, analysisJet, allTracks, origTracks, indicesTracks, jetFlavor, eventWeight, &trkLabels); registry.fill(HIST("h2_jetMass_jetpT"), analysisJet.pt(), analysisJet.mass(), eventWeight); - registry.fill(HIST("h2_nTracks_jetpT"), analysisJet.pt(), tracksIndices.size()); + registry.fill(HIST("h2_nTracks_jetpT"), analysisJet.pt(), indicesTracks.size()); //+jet registry.fill(HIST("h_jet_pt"), analysisJet.pt()); registry.fill(HIST("h_jet_eta"), analysisJet.eta()); registry.fill(HIST("h_jet_phi"), analysisJet.phi()); registry.fill(HIST("h_jet_flav"), jetFlavor); - registry.fill(HIST("h_n_trks"), tracksIndices.size()); + registry.fill(HIST("h_n_trks"), indicesTracks.size()); registry.fill(HIST("h_jet_mass"), analysisJet.mass()); registry.fill(HIST("h_n_vertices"), nVertices); - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h2_jetMass_jetpT_bjet"), analysisJet.pt(), analysisJet.mass(), eventWeight); registry.fill(HIST("h_jetpT_detector_bjet"), analysisJet.pt(), eventWeight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h2_jetMass_jetpT_cjet"), analysisJet.pt(), analysisJet.mass(), eventWeight); registry.fill(HIST("h_jetpT_detector_cjet"), analysisJet.pt(), eventWeight); } else { @@ -771,8 +774,8 @@ struct BJetTreeCreator { } if (produceTree) { - bjetConstituentsTable(bjetParamsTable.lastIndex() + 1, tracksIndices, SVsIndices); - bjetParamsTable(analysisJet.pt(), analysisJet.eta(), analysisJet.phi(), tracksIndices.size(), nVertices, analysisJet.mass(), jetFlavor, analysisJet.r()); + bjetConstituentsTable(bjetParamsTable.lastIndex() + 1, indicesTracks, indicesSVs); + bjetParamsTable(analysisJet.pt(), analysisJet.eta(), analysisJet.phi(), indicesTracks.size(), nVertices, analysisJet.mass(), jetFlavor, analysisJet.r()); } } } @@ -787,7 +790,7 @@ struct BJetTreeCreator { for (const auto& mcpjet : MCPjets) { bool jetIncluded = false; - for (auto jetR : jetRadiiValues) { + for (const auto& jetR : jetRadiiValues) { if (mcpjet.r() == static_cast(jetR * 100)) { jetIncluded = true; break; @@ -812,9 +815,9 @@ struct BJetTreeCreator { // jetFlavor = jettaggingutilities::mcpJetFromHFShower(mcpjet, MCParticles, (float)(mcpjet.r() / 100.)); } - if (jetFlavor == 2) { + if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h_jetpT_particle_bjet"), mcpjet.pt(), eventWeight); - } else if (jetFlavor == 1) { + } else if (jetFlavor == JetTaggingSpecies::charm) { registry.fill(HIST("h_jetpT_particle_cjet"), mcpjet.pt(), eventWeight); } else { registry.fill(HIST("h_jetpT_particle_lfjet"), mcpjet.pt(), eventWeight);