From 9f138a2977f7d32ff7769910e73339b6411c738e Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Mon, 3 Jun 2024 21:30:40 +0200 Subject: [PATCH 01/21] remote unused parameter --- PWGJE/Core/JetTaggingUtilities.h | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index bd33bdf1cb8..4c2017a21ef 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -368,20 +368,19 @@ std::unique_ptr setResolutionFunction(T const& vecParams) * * @param fResoFuncjet The resolution function for the jet, used to model the distribution of impact * parameter significances for tracks associated with the jet. - * @param collision The collision event data, providing context for calculating geometric signs. - * @param jet The specific jet being analyzed. * @param track The track for which the probability is being calculated. * @param minSignImpXYSig The minimum significance of the impact parameter in the XY plane, used as - * the lower limit for integration of the resolution function. Defaults to -10. + * the lower limit for integration of the resolution function. Defaults to -40. * @return The calculated probability of the track being associated with the jet, based on its * impact parameter significance. */ -template -float getTrackProbability(T const& fResoFuncjet, U const& collision, V const& jet, W const& track, const float& minSignImpXYSig = -10) +template +float getTrackProbability(T const& fResoFuncjet, U const& track, const float& minSignImpXYSig = -40) { float probTrack = 0.; - auto varSignImpXYSig = getGeoSign(collision, jet, track) * TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); - probTrack = fResoFuncjet->Integral(minSignImpXYSig, -1 * TMath::Abs(varSignImpXYSig)) / fResoFuncjet->Integral(minSignImpXYSig, 0); + auto varSignImpXYSig = TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); + if (-varSignImpXYSig < minSignImpXYSig) varSignImpXYSig = -minSignImpXYSig+0.01; // To avoid overflow for integral + probTrack = fResoFuncjet->Integral(minSignImpXYSig, -varSignImpXYSig) / fResoFuncjet->Integral(minSignImpXYSig, 0); return probTrack; } @@ -417,7 +416,7 @@ float getJetProbability(T const& fResoFuncjet, U const& collision, V const& jet, for (auto& jtrack : jet.template tracks_as()) { auto track = jtrack.template track_as(); - float probTrack = getTrackProbability(fResoFuncjet, collision, jet, track, minSignImpXYSig); + float probTrack = getTrackProbability(fResoFuncjet, track, minSignImpXYSig); auto geoSign = getGeoSign(collision, jet, track); if (geoSign > 0) { // only take positive sign track for JP calculation From 41df8927d25dceec16d280f2f44ab8afb87e7df4 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Mon, 3 Jun 2024 21:32:14 +0200 Subject: [PATCH 02/21] modiciation of track counting --- PWGJE/Tasks/jettaggerhfQA.cxx | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index f53452001a6..1e040267027 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -517,17 +517,26 @@ struct JetTaggerHFQA { std::sort(vecSignImpZSigTC.begin(), vecSignImpZSigTC.end(), sortImp); std::sort(vecSignImpXYZSigTC.begin(), vecSignImpXYZSigTC.end(), sortImp); - if (vecSignImpXYSigTC.empty() || vecSignImpZSigTC.empty() || vecSignImpXYZSigTC.empty()) + if (vecSignImpXYSigTC.empty()) continue; for (int order = 1; order <= numOrder; order++) { if (fillIPxy && order < vecSignImpXYSigTC.size()) { - registry.fill(HIST("h3_sign_impact_parameter_xy_significance_tc_flavour"), vecSignImpXYSigTC[order][0], order, jetflavour); + registry.fill(HIST("h3_sign_impact_parameter_xy_significance_tc_flavour"), vecSignImpXYSigTC[order - 1][0], order, jetflavour); } - if (fillIPz && order < vecSignImpZSigTC.size()) { - registry.fill(HIST("h3_sign_impact_parameter_z_significance_tc_flavour"), vecSignImpZSigTC[order][0], order, jetflavour); + } + if (vecSignImpZSigTC.empty()) + continue; + for (int order = 1; order <= numOrder; order++) { + if (fillIPxy && order < vecSignImpZSigTC.size()) { + registry.fill(HIST("h3_sign_impact_parameter_z_significance_tc_flavour"), vecSignImpZSigTC[order - 1][0], order, jetflavour); } - if (fillIPxyz && order < vecSignImpXYZSigTC.size()) { - registry.fill(HIST("h3_sign_impact_parameter_xyz_significance_tc_flavour"), vecSignImpXYZSigTC[order][0], order, jetflavour); + } + + if (vecSignImpXYZSigTC.empty()) + continue; + for (int order = 1; order <= numOrder; order++) { + if (fillIPxy && order < vecSignImpXYZSigTC.size()) { + registry.fill(HIST("h3_sign_impact_parameter_xyz_significance_tc_flavour"), vecSignImpXYZSigTC[order - 1][0], order, jetflavour); } } } From e295b7447138aeda5424f402397dcda0f7a58176 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Mon, 3 Jun 2024 21:33:13 +0200 Subject: [PATCH 03/21] modification of resolution function to add flavour --- PWGJE/TableProducer/jettaggerhf.cxx | 127 ++++++++++++++++++++++++---- 1 file changed, 109 insertions(+), 18 deletions(-) diff --git a/PWGJE/TableProducer/jettaggerhf.cxx b/PWGJE/TableProducer/jettaggerhf.cxx index b7249a74e20..db91736b711 100644 --- a/PWGJE/TableProducer/jettaggerhf.cxx +++ b/PWGJE/TableProducer/jettaggerhf.cxx @@ -45,37 +45,103 @@ struct JetTaggerHFTask { Configurable trackProbQA{"trackProbQA", false, "fill track probability histograms separately for geometric positive and negative tracks for QA"}; Configurable doSV{"doSV", false, "fill table for secondary vertex algorithm"}; Configurable numCount{"numCount", 3, "number of track counting"}; - Configurable> paramsResoFunc{"paramsResoFunc", std::vector{1306800, -0.1049, 0.861425, 13.7547, 0.977967, 8.96823, 0.151595, 6.94499, 0.0250301}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7))"}; - Configurable> paramsResoFuncMC{"paramsResoFuncMC", std::vector{61145.8, -0.082085, 0.706361, 10.0794, 1.15412, 5.81011, 0.188979, 3.8514, 0.032457}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7)))"}; + Configurable resoFuncMatching{"resoFuncMatching", 0, "matching parameters of resolution function as MC samble (0: custom, 1: custom & inc, 2: MB, 3: MB & inc, 4: JJ, 5: JJ & inc)"}; + Configurable> paramsResoFuncData{"paramsResoFuncData", std::vector{1306800, -0.1049, 0.861425, 13.7547, 0.977967, 8.96823, 0.151595, 6.94499, 0.0250301}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7))"}; + Configurable> paramsResoFuncIncJetMC{"paramsResoFuncIncJetMC", std::vector{1908803.027, -0.059, 0.895, 13.467, 1.005, 8.867, 0.098, 6.929, 0.011}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7)))"}; + Configurable> paramsResoFuncCharmJetMC{"paramsResoFuncCharmJetMC", std::vector{282119.753, -0.065, 0.893, 11.608, 0.945, 8.029, 0.131, 6.244, 0.027}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7)))"}; + Configurable> paramsResoFuncBeautyJetMC{"paramsResoFuncBeautyJetMC", std::vector{74901.583, -0.082, 0.874, 10.332, 0.941, 7.352, 0.097, 6.220, 0.022}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7)))"}; + Configurable> paramsResoFuncLfJetMC{"paramsResoFuncLfJetMC", std::vector{1539435.343, -0.061, 0.896, 13.272, 1.034, 5.884, 0.004, 7.843, 0.090}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7)))"}; Configurable minSignImpXYSig{"minsIPs", -40.0, "minimum of signed impact parameter significance"}; Configurable tagPoint{"tagPoint", 2.5, "tagging working point"}; ConfigurableAxis binTrackProbability{"binTrackProbability", {100, 0.f, 1.f}, ""}; + ConfigurableAxis binJetFlavour{"binJetFlavour", {6, -0.5, 5.5}, ""}; using JetTagTracksData = soa::Join; using JetTagTracksMCD = soa::Join; using OriTracksData = soa::Join; using OriTracksMCD = soa::Join; - std::unique_ptr fSignImpXYSig = nullptr; - std::vector vecParams; + std::vector vecParamsData; + std::vector vecParamsIncJetMC; + std::vector vecParamsCharmJetMC; + std::vector vecParamsBeautyJetMC; + std::vector vecParamsLfJetMC; std::vector jetProb; + bool useResoFuncFromIncJet = false; int maxOrder = -1; + int resoFuncMatch = 0; + std::unique_ptr fSignImpXYSigData = nullptr; + std::unique_ptr fSignImpXYSigIncJetMC = nullptr; + std::unique_ptr fSignImpXYSigCharmJetMC = nullptr; + std::unique_ptr fSignImpXYSigBeautyJetMC = nullptr; + std::unique_ptr fSignImpXYSigLfJetMC = nullptr; + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; void init(InitContext const&) { maxOrder = numCount + 1; // 0: untagged, >1 : N ordering - if (doprocessData) { - vecParams = (std::vector)paramsResoFunc; - } - if (doprocessMCD) { - vecParams = (std::vector)paramsResoFuncMC; + + // Set up the resolution function + resoFuncMatch = resoFuncMatching; + switch (resoFuncMatch) { + case 0: + vecParamsData = (std::vector) paramsResoFuncData; + vecParamsIncJetMC = (std::vector) paramsResoFuncIncJetMC; + LOG(info) << "defined parameters of resolution function: custom"; + break; + case 1: + vecParamsData = (std::vector) paramsResoFuncData; + vecParamsCharmJetMC = (std::vector) paramsResoFuncCharmJetMC; + vecParamsBeautyJetMC = (std::vector) paramsResoFuncBeautyJetMC; + vecParamsLfJetMC = (std::vector) paramsResoFuncLfJetMC; + useResoFuncFromIncJet = true; + LOG(info) << "defined parameters of resolution function: custom & use inclusive distribution"; + break; + case 2: // TODO + vecParamsData = (std::vector) paramsResoFuncData; + vecParamsCharmJetMC = {282119.753, -0.065, 0.893, 11.608, 0.945, 8.029, 0.131, 6.244, 0.027}; + vecParamsBeautyJetMC = {74901.583, -0.082, 0.874, 10.332, 0.941, 7.352, 0.097, 6.220, 0.022}; + vecParamsLfJetMC = {1539435.343, -0.061, 0.896, 13.272, 1.034, 5.884, 0.004, 7.843, 0.090}; + LOG(info) << "defined parameters of resolution function: PYTHIA8, MB, LHC23d1k"; + break; + case 3: // TODO + vecParamsData = (std::vector) paramsResoFuncData; + vecParamsIncJetMC = {1908803.027, -0.059, 0.895, 13.467, 1.005, 8.867, 0.098, 6.929, 0.011}; + LOG(info) << "defined parameters of resolution function: PYTHIA8, MB, LHC23d1k & use inclusive distribution"; + useResoFuncFromIncJet = true; + break; + case 4: // TODO + vecParamsData = (std::vector) paramsResoFuncData; + vecParamsCharmJetMC = {281446.003, -0.063, 0.894, 11.598, 0.943, 8.025, 0.130, 6.227, 0.027}; + vecParamsBeautyJetMC = {74839.065, -0.081, 0.875, 10.314, 0.939, 7.326, 0.101, 6.309, 0.024}; + vecParamsLfJetMC = {1531580.038, -0.062, -0.896, 13.267, 1.034, 5.866, 0.004, 7.836, 0.090}; + LOG(info) << "defined parameters of resolution function: PYTHIA8, JJ, LHC23d4"; + break; + case 5: // TODO + vecParamsData = (std::vector) paramsResoFuncData; + vecParamsIncJetMC = {1900387.527, -0.059, 0.895, 13.461, 1.004, 8.860, 0.098, 6.931, 0.011}; + LOG(info) << "defined parameters of resolution function: PYTHIA8, JJ, LHC23d4 & use inclusive distribution"; + useResoFuncFromIncJet = true; + break; + default: + LOG(fatal) << "undefined parameters of resolution function. Fix it!"; + break; } - fSignImpXYSig = jettaggingutilities::setResolutionFunction(vecParams); + + fSignImpXYSigData = jettaggingutilities::setResolutionFunction(vecParamsData); + fSignImpXYSigIncJetMC = jettaggingutilities::setResolutionFunction(vecParamsIncJetMC); + fSignImpXYSigCharmJetMC = jettaggingutilities::setResolutionFunction(vecParamsCharmJetMC); + fSignImpXYSigBeautyJetMC = jettaggingutilities::setResolutionFunction(vecParamsBeautyJetMC); + fSignImpXYSigLfJetMC = jettaggingutilities::setResolutionFunction(vecParamsLfJetMC); + + + // Use QA for effectivness of track probability if (trackProbQA) { - AxisSpec TrackProbabilityAxis = {binTrackProbability, "Track proability"}; - registry.add("h_pos_track_probability", "track probability", {HistType::kTH1F, {{TrackProbabilityAxis}}}); - registry.add("h_neg_track_probability", "track probability", {HistType::kTH1F, {{TrackProbabilityAxis}}}); + AxisSpec trackProbabilityAxis = {binTrackProbability, "Track proability"}; + AxisSpec jetFlavourAxis = {binJetFlavour, "Jet flavour"}; + 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}}}); } } @@ -93,7 +159,7 @@ struct JetTaggerHFTask { jetProb.clear(); jetProb.reserve(maxOrder); for (int order = 0; order < maxOrder; order++) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSig, collision, jet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigData, collision, jet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); } } // if (doSV) algorithm2 = jettaggingutilities::Algorithm2((mcdjet, tracks); @@ -117,17 +183,42 @@ struct JetTaggerHFTask { jetProb.clear(); jetProb.reserve(maxOrder); for (int order = 0; order < maxOrder; order++) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSig, collision, mcdjet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); + if (useResoFuncFromIncJet) { + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigIncJetMC, collision, mcdjet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); + } else { + if (origin == JetTaggingSpecies::charm) { + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigCharmJetMC, collision, mcdjet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); + } + if (origin == JetTaggingSpecies::beauty) { + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigBeautyJetMC, collision, mcdjet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); + } + if (origin == JetTaggingSpecies::lightflavour) { + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigLfJetMC, collision, mcdjet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); + } + } } if (trackProbQA) { for (auto& jtrack : mcdjet.template tracks_as()) { auto track = jtrack.template track_as(); auto geoSign = jettaggingutilities::getGeoSign(collision, mcdjet, track); - float probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSig, collision, mcdjet, track, minSignImpXYSig); + float probTrack = -1; + if (useResoFuncFromIncJet) { + probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigIncJetMC, track, minSignImpXYSig); + } else { + if (origin == JetTaggingSpecies::charm) { + probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigCharmJetMC, track, minSignImpXYSig); + } + if (origin == JetTaggingSpecies::beauty) { + probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigBeautyJetMC, track, minSignImpXYSig); + } + if (origin == JetTaggingSpecies::lightflavour) { + probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigLfJetMC, track, minSignImpXYSig); + } + } if (geoSign > 0) { - registry.fill(HIST("h_pos_track_probability"), probTrack); + registry.fill(HIST("h2_pos_track_probability_flavour"), probTrack, origin); } else { - registry.fill(HIST("h_neg_track_probability"), probTrack); + registry.fill(HIST("h2_neg_track_probability_flavour"), probTrack, origin); } } } From 8a8f7e7c93728175cd21ddf9bb672815c432adb7 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Mon, 3 Jun 2024 21:34:20 +0200 Subject: [PATCH 04/21] fix clang-format --- PWGJE/Core/JetTaggingUtilities.h | 3 ++- PWGJE/TableProducer/jettaggerhf.cxx | 23 +++++++++++------------ 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index 4c2017a21ef..abbb5fb7eb2 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -379,7 +379,8 @@ float getTrackProbability(T const& fResoFuncjet, U const& track, const float& mi { float probTrack = 0.; auto varSignImpXYSig = TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); - if (-varSignImpXYSig < minSignImpXYSig) varSignImpXYSig = -minSignImpXYSig+0.01; // To avoid overflow for integral + if (-varSignImpXYSig < minSignImpXYSig) + varSignImpXYSig = -minSignImpXYSig + 0.01; // To avoid overflow for integral probTrack = fResoFuncjet->Integral(minSignImpXYSig, -varSignImpXYSig) / fResoFuncjet->Integral(minSignImpXYSig, 0); return probTrack; diff --git a/PWGJE/TableProducer/jettaggerhf.cxx b/PWGJE/TableProducer/jettaggerhf.cxx index db91736b711..225f7919db6 100644 --- a/PWGJE/TableProducer/jettaggerhf.cxx +++ b/PWGJE/TableProducer/jettaggerhf.cxx @@ -81,45 +81,45 @@ struct JetTaggerHFTask { void init(InitContext const&) { maxOrder = numCount + 1; // 0: untagged, >1 : N ordering - + // Set up the resolution function resoFuncMatch = resoFuncMatching; switch (resoFuncMatch) { case 0: - vecParamsData = (std::vector) paramsResoFuncData; - vecParamsIncJetMC = (std::vector) paramsResoFuncIncJetMC; + vecParamsData = (std::vector)paramsResoFuncData; + vecParamsIncJetMC = (std::vector)paramsResoFuncIncJetMC; LOG(info) << "defined parameters of resolution function: custom"; break; case 1: - vecParamsData = (std::vector) paramsResoFuncData; - vecParamsCharmJetMC = (std::vector) paramsResoFuncCharmJetMC; - vecParamsBeautyJetMC = (std::vector) paramsResoFuncBeautyJetMC; - vecParamsLfJetMC = (std::vector) paramsResoFuncLfJetMC; + vecParamsData = (std::vector)paramsResoFuncData; + vecParamsCharmJetMC = (std::vector)paramsResoFuncCharmJetMC; + vecParamsBeautyJetMC = (std::vector)paramsResoFuncBeautyJetMC; + vecParamsLfJetMC = (std::vector)paramsResoFuncLfJetMC; useResoFuncFromIncJet = true; LOG(info) << "defined parameters of resolution function: custom & use inclusive distribution"; break; case 2: // TODO - vecParamsData = (std::vector) paramsResoFuncData; + vecParamsData = (std::vector)paramsResoFuncData; vecParamsCharmJetMC = {282119.753, -0.065, 0.893, 11.608, 0.945, 8.029, 0.131, 6.244, 0.027}; vecParamsBeautyJetMC = {74901.583, -0.082, 0.874, 10.332, 0.941, 7.352, 0.097, 6.220, 0.022}; vecParamsLfJetMC = {1539435.343, -0.061, 0.896, 13.272, 1.034, 5.884, 0.004, 7.843, 0.090}; LOG(info) << "defined parameters of resolution function: PYTHIA8, MB, LHC23d1k"; break; case 3: // TODO - vecParamsData = (std::vector) paramsResoFuncData; + vecParamsData = (std::vector)paramsResoFuncData; vecParamsIncJetMC = {1908803.027, -0.059, 0.895, 13.467, 1.005, 8.867, 0.098, 6.929, 0.011}; LOG(info) << "defined parameters of resolution function: PYTHIA8, MB, LHC23d1k & use inclusive distribution"; useResoFuncFromIncJet = true; break; case 4: // TODO - vecParamsData = (std::vector) paramsResoFuncData; + vecParamsData = (std::vector)paramsResoFuncData; vecParamsCharmJetMC = {281446.003, -0.063, 0.894, 11.598, 0.943, 8.025, 0.130, 6.227, 0.027}; vecParamsBeautyJetMC = {74839.065, -0.081, 0.875, 10.314, 0.939, 7.326, 0.101, 6.309, 0.024}; vecParamsLfJetMC = {1531580.038, -0.062, -0.896, 13.267, 1.034, 5.866, 0.004, 7.836, 0.090}; LOG(info) << "defined parameters of resolution function: PYTHIA8, JJ, LHC23d4"; break; case 5: // TODO - vecParamsData = (std::vector) paramsResoFuncData; + vecParamsData = (std::vector)paramsResoFuncData; vecParamsIncJetMC = {1900387.527, -0.059, 0.895, 13.461, 1.004, 8.860, 0.098, 6.931, 0.011}; LOG(info) << "defined parameters of resolution function: PYTHIA8, JJ, LHC23d4 & use inclusive distribution"; useResoFuncFromIncJet = true; @@ -135,7 +135,6 @@ struct JetTaggerHFTask { fSignImpXYSigBeautyJetMC = jettaggingutilities::setResolutionFunction(vecParamsBeautyJetMC); fSignImpXYSigLfJetMC = jettaggingutilities::setResolutionFunction(vecParamsLfJetMC); - // Use QA for effectivness of track probability if (trackProbQA) { AxisSpec trackProbabilityAxis = {binTrackProbability, "Track proability"}; From d26b819f4b75e7ba2f5d4cd7871688904624cdb7 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Mon, 3 Jun 2024 21:40:49 +0200 Subject: [PATCH 05/21] fix the calcualtion --- PWGJE/Core/JetTaggingUtilities.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index abbb5fb7eb2..9ea872269e0 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -377,10 +377,10 @@ std::unique_ptr setResolutionFunction(T const& vecParams) template float getTrackProbability(T const& fResoFuncjet, U const& track, const float& minSignImpXYSig = -40) { - float probTrack = 0.; + float probTrack = 0. auto varSignImpXYSig = TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); if (-varSignImpXYSig < minSignImpXYSig) - varSignImpXYSig = -minSignImpXYSig + 0.01; // To avoid overflow for integral + varSignImpXYSig = -minSignImpXYSig - 0.01; // To avoid overflow for integral probTrack = fResoFuncjet->Integral(minSignImpXYSig, -varSignImpXYSig) / fResoFuncjet->Integral(minSignImpXYSig, 0); return probTrack; From 887392ac9f5ac9c112ab7ff75f41ca7291cb19f5 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Mon, 3 Jun 2024 21:46:51 +0200 Subject: [PATCH 06/21] Fix clang-format --- PWGJE/Core/JetTaggingUtilities.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index 9ea872269e0..ab4bc658d1f 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -377,8 +377,7 @@ std::unique_ptr setResolutionFunction(T const& vecParams) template float getTrackProbability(T const& fResoFuncjet, U const& track, const float& minSignImpXYSig = -40) { - float probTrack = 0. - auto varSignImpXYSig = TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); + float probTrack = 0. auto varSignImpXYSig = TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); if (-varSignImpXYSig < minSignImpXYSig) varSignImpXYSig = -minSignImpXYSig - 0.01; // To avoid overflow for integral probTrack = fResoFuncjet->Integral(minSignImpXYSig, -varSignImpXYSig) / fResoFuncjet->Integral(minSignImpXYSig, 0); From d538be05b7907d5a674dc63c03e87afb9cbd14df Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Mon, 3 Jun 2024 22:42:04 +0200 Subject: [PATCH 07/21] fix mistake --- PWGJE/Core/JetTaggingUtilities.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index ab4bc658d1f..99a2691c80b 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -377,7 +377,8 @@ std::unique_ptr setResolutionFunction(T const& vecParams) template float getTrackProbability(T const& fResoFuncjet, U const& track, const float& minSignImpXYSig = -40) { - float probTrack = 0. auto varSignImpXYSig = TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); + float probTrack = 0; + auto varSignImpXYSig = TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); if (-varSignImpXYSig < minSignImpXYSig) varSignImpXYSig = -minSignImpXYSig - 0.01; // To avoid overflow for integral probTrack = fResoFuncjet->Integral(minSignImpXYSig, -varSignImpXYSig) / fResoFuncjet->Integral(minSignImpXYSig, 0); From 1dca25a40ca7a27d5dc47d029bb26b7398cf5a31 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Sun, 9 Jun 2024 23:32:55 +0200 Subject: [PATCH 08/21] Processing to update QA of sv --- PWGJE/Tasks/jettaggerhfQA.cxx | 80 +++++++++++++++++++++++------------ 1 file changed, 53 insertions(+), 27 deletions(-) diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index 1e040267027..4722246d45e 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -96,7 +96,6 @@ struct JetTaggerHFQA { void init(InitContext const&) { - // Axis AxisSpec jetFlavourAxis = {binJetFlavour, "Jet flavour"}; AxisSpec jetPtAxis = {binJetPt, "#it{p}_{T, jet}"}; @@ -236,12 +235,26 @@ struct JetTaggerHFQA { } if (doprocessSV2ProngMCD) { registry.add("h_2prong_nprongs", "", {HistType::kTH1F, {{nprongsAxis}}}); + registry.add("h2_jet_pt_2prong_Lxy", "", {HistType::kTH2F, {{jetPtAxis}, {LxyAxis}}}); + registry.add("h2_jet_pt_2prong_Sxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); + registry.add("h2_jet_pt_2prong_Lxyz", "", {HistType::kTH2F, {{jetPtAxis}, {LxyzAxis}}}); + registry.add("h2_jet_pt_2prong_Sxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); + registry.add("h2_jet_pt_2prong_Sxy_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); + registry.add("h2_jet_pt_2prong_Sxyz_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); + registry.add("h2_2prong_Sxy_sigmaLxy", "", {HistType::kTH2F, {{SxyAxis}, {sigmaLxyAxis}}}); + registry.add("h2_2prong_Sxyz_sigmaLxyz", "", {HistType::kTH2F, {{SxyzAxis}, {sigmaLxyzAxis}}}); + registry.add("h2_jet_pt_2prong_sigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyAxis}}}); + registry.add("h2_jet_pt_2prong_sigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyzAxis}}}); + } + + if (doprocessSV2ProngMCD) { + registry.add("h_2prong_nprongs", "", {HistType::kTH2F, {{nprongsAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Lxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {LxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Sxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Lxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {LxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Sxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_2prong_Sxy_flavour_N1", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_2prong_Sxyz_flavour_N1", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_2prong_Sxy_N1_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_2prong_Sxyz_N1_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_2prong_Sxy_sigmaLxy_flavour", "", {HistType::kTH3F, {{SxyAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); registry.add("h3_2prong_Sxyz_sigmaLxyz_flavour", "", {HistType::kTH3F, {{SxyzAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_sigmaLxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); @@ -253,8 +266,8 @@ struct JetTaggerHFQA { registry.add("h3_jet_pt_3prong_Sxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_Lxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {LxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_Sxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_3prong_Sxy_flavour_N1", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_3prong_Sxyz_flavour_N1", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_3prong_Sxy_N1_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_3prong_Sxyz_N1_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_3prong_Sxy_sigmaLxy_flavour", "", {HistType::kTH3F, {{SxyAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); registry.add("h3_3prong_Sxyz_sigmaLxyz_flavour", "", {HistType::kTH3F, {{SxyzAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_sigmaLxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); @@ -578,26 +591,32 @@ struct JetTaggerHFQA { } } - Preslice perColFor2Prong = aod::hf_cand::collisionId; - Preslice perColFor3Prong = aod::hf_cand::collisionId; + template + void fillHistogramSV2ProngData(T const& collision, U const& jets, V const& prongs) + { + LOG(info) << "2Prong"; + } + + template + void fillHistogramSV3ProngData(T const& collision, U const& jets, V const& prongs) + { + LOG(info) << "3Prong"; + } + template void fillHistogramSV2ProngMCD(T const& collision, U const& mcdjets, V const& prongs) { - int numOfSV = 0; int maxSxy = 0; int maxSxyz = 0; for (const auto& mcdjet : mcdjets) { + int numOfSV = 0; if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { continue; } auto origin = mcdjet.origin(); - auto prongsPerCol = prongs.sliceBy(perColFor2Prong, collision.globalIndex()); - for (const auto& prong : prongsPerCol) { - auto deltaR = jetutilities::deltaR(mcdjet, prong); - if (deltaR > mcdjet.r() / 100.0) { - continue; - } + registry.fill(HIST("h_2prong_nprongs"), mcdjet.template secondaryVertices_as().size()); + for (const auto& prong : mcdjet.template secondaryVertices_as()) { numOfSV++; auto Lxy = prong.decayLengthXY(); auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); @@ -616,10 +635,10 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_jet_pt_2prong_sigmaLxy_flavour"), mcdjet.pt(), prong.errorDecayLengthXY(), origin); registry.fill(HIST("h3_jet_pt_2prong_sigmaLxyz_flavour"), mcdjet.pt(), prong.errorDecayLength(), origin); } - registry.fill(HIST("h3_jet_pt_2prong_Sxy_flavour_N1"), mcdjet.pt(), maxSxy, origin); - registry.fill(HIST("h3_jet_pt_2prong_Sxyz_flavour_N1"), mcdjet.pt(), maxSxyz, origin); + registry.fill(HIST("h3_jet_pt_2prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); + registry.fill(HIST("h3_jet_pt_2prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin); + registry.fill(HIST("h_2prong_nprongs"), numOfSV); } - registry.fill(HIST("h_2prong_nprongs"), numOfSV); } template @@ -633,12 +652,7 @@ struct JetTaggerHFQA { continue; } auto origin = mcdjet.origin(); - auto prongsPerCol = prongs.sliceBy(perColFor3Prong, collision.globalIndex()); - for (const auto& prong : prongsPerCol) { - auto deltaR = jetutilities::deltaR(mcdjet, prong); - if (deltaR > mcdjet.r() / 100.0) { - continue; - } + for (const auto& prong : mcdjet.template secondaryVertices_as()) { numOfSV++; auto Lxy = prong.decayLengthXY(); auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); @@ -657,8 +671,8 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_jet_pt_3prong_sigmaLxy_flavour"), mcdjet.pt(), prong.errorDecayLengthXY(), origin); registry.fill(HIST("h3_jet_pt_3prong_sigmaLxyz_flavour"), mcdjet.pt(), prong.errorDecayLength(), origin); } - registry.fill(HIST("h3_jet_pt_2prong_Sxy_flavour_N1"), mcdjet.pt(), maxSxy, origin); - registry.fill(HIST("h3_jet_pt_2prong_Sxyz_flavour_N1"), mcdjet.pt(), maxSxyz, origin); + registry.fill(HIST("h3_jet_pt_2prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); + registry.fill(HIST("h3_jet_pt_2prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin); } registry.fill(HIST("h_3prong_nprongs"), numOfSV); } @@ -720,13 +734,25 @@ struct JetTaggerHFQA { } PROCESS_SWITCH(JetTaggerHFQA, processJPMCD, "Fill jet probability imformation for mcd jets", false); - void processSV2ProngMCD(soa::Filtered::iterator const& jcollision, JetTagTableMCD const& mcdjets, aod::HfCand2Prong const& prongs) + void processSV2ProngData(soa::Filtered::iterator const& jcollision, soa::Join const& jets, aod::DataSecondaryVertex2Prongs const& prongs) + { + fillHistogramSV2ProngMCD(jcollision, jets, prongs); + } + PROCESS_SWITCH(JetTaggerHFQA, processSV2ProngData, "Fill 2prong impormation for data jets", false); + + void processSV3ProngData(soa::Filtered::iterator const& jcollision, soa::Join const& jets, aod::DataSecondaryVertex3Prongs const& prongs) + { + fillHistogramSV3ProngMCD(jcollision, jets, prongs); + } + PROCESS_SWITCH(JetTaggerHFQA, processSV3ProngData, "Fill 3prong impormation for data jets", false); + + void processSV2ProngMCD(soa::Filtered::iterator const& jcollision, soa::Join const& mcdjets, aod::MCDSecondaryVertex2Prongs const& prongs) { fillHistogramSV2ProngMCD(jcollision, mcdjets, prongs); } PROCESS_SWITCH(JetTaggerHFQA, processSV2ProngMCD, "Fill 2prong impormation for mcd jets", false); - void processSV3ProngMCD(soa::Filtered::iterator const& jcollision, JetTagTableMCD const& mcdjets, aod::HfCand3Prong const& prongs) + void processSV3ProngMCD(soa::Filtered::iterator const& jcollision, soa::Join const& mcdjets, aod::MCDSecondaryVertex3Prongs const& prongs) { fillHistogramSV3ProngMCD(jcollision, mcdjets, prongs); } From 53d7017eb4f946224654110705f925b11eaf13dd Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Wed, 12 Jun 2024 14:32:06 +0200 Subject: [PATCH 09/21] Being updating SV QA --- PWGJE/Tasks/jettaggerhfQA.cxx | 148 +++++++++++++++++++++++++--------- 1 file changed, 109 insertions(+), 39 deletions(-) diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index 4722246d45e..bf4ae40865a 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -117,8 +117,8 @@ struct JetTaggerHFQA { AxisSpec SxyAxis = {binSxy, "S_{XY}"}; AxisSpec LxyzAxis = {binLxyz, "L_{XYZ} [cm]"}; AxisSpec SxyzAxis = {binSxyz, "S_{XYZ}"}; - AxisSpec sigmaLxyAxis = {binSigmaLxy, "#simga_{L_{XY}} [cm]"}; - AxisSpec sigmaLxyzAxis = {binSigmaLxyz, "#simga_{L_{XYZ}} [cm]"}; + AxisSpec sigmaLxyAxis = {binSigmaLxy, "#sigma_{L_{XY}} [cm]"}; + AxisSpec sigmaLxyzAxis = {binSigmaLxyz, "#sigma_{L_{XYZ}} [cm]"}; numberOfJetFlavourSpecies = static_cast(numFlavourSpecies); @@ -233,7 +233,7 @@ struct JetTaggerHFQA { registry.add("h3_jet_pt_JP_N3_flavour", "jet pt jet probability flavour N3", {HistType::kTH3F, {{jetPtAxis}, {JetProbabilityAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_neg_log_JP_N3_flavour", "jet pt log jet probability flavour N3", {HistType::kTH3F, {{jetPtAxis}, {JetProbabilityLogAxis}, {jetFlavourAxis}}}); } - if (doprocessSV2ProngMCD) { + if (doprocessSV2ProngData) { registry.add("h_2prong_nprongs", "", {HistType::kTH1F, {{nprongsAxis}}}); registry.add("h2_jet_pt_2prong_Lxy", "", {HistType::kTH2F, {{jetPtAxis}, {LxyAxis}}}); registry.add("h2_jet_pt_2prong_Sxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); @@ -246,9 +246,21 @@ struct JetTaggerHFQA { registry.add("h2_jet_pt_2prong_sigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyAxis}}}); registry.add("h2_jet_pt_2prong_sigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyzAxis}}}); } - + if (doprocessSV3ProngData) { + registry.add("h_3prong_nprongs", "", {HistType::kTH1F, {{nprongsAxis}}}); + registry.add("h2_jet_pt_3prong_Lxy", "", {HistType::kTH2F, {{jetPtAxis}, {LxyAxis}}}); + registry.add("h2_jet_pt_3prong_Sxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); + registry.add("h2_jet_pt_3prong_Lxyz", "", {HistType::kTH2F, {{jetPtAxis}, {LxyzAxis}}}); + registry.add("h2_jet_pt_3prong_Sxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); + registry.add("h2_jet_pt_3prong_Sxy_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); + registry.add("h2_jet_pt_3prong_Sxyz_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); + registry.add("h2_3prong_Sxy_sigmaLxy", "", {HistType::kTH2F, {{SxyAxis}, {sigmaLxyAxis}}}); + registry.add("h2_3prong_Sxyz_sigmaLxyz", "", {HistType::kTH2F, {{SxyzAxis}, {sigmaLxyzAxis}}}); + registry.add("h2_jet_pt_3prong_sigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyAxis}}}); + registry.add("h2_jet_pt_3prong_sigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyzAxis}}}); + } if (doprocessSV2ProngMCD) { - registry.add("h_2prong_nprongs", "", {HistType::kTH2F, {{nprongsAxis}, {jetFlavourAxis}}}); + registry.add("h2_2prong_nprongs_flavour", "", {HistType::kTH2F, {{nprongsAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Lxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {LxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Sxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Lxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {LxyzAxis}, {jetFlavourAxis}}}); @@ -261,7 +273,7 @@ struct JetTaggerHFQA { registry.add("h3_jet_pt_2prong_sigmaLxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); } if (doprocessSV3ProngMCD) { - registry.add("h_3prong_nprongs", "", {HistType::kTH1F, {{nprongsAxis}}}); + registry.add("h2_3prong_nprongs_flavour", "", {HistType::kTH2F, {{nprongsAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_Lxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {LxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_Sxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_Lxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {LxyzAxis}, {jetFlavourAxis}}}); @@ -307,6 +319,32 @@ struct JetTaggerHFQA { return 1; } + template + double getMaxSxyForJet(const JetType& mcdjet) + { + double maxSxy = 0; + for (const auto& prong : mcdjet.template secondaryVertices_as()) { + int Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); + if (maxSxy < Sxy) { + maxSxy = Sxy; + } + } + return maxSxy; + } + + template + double getMaxSxyzForJet(const JetType& mcdjet) + { + double maxSxyz = 0; + for (const auto& prong : mcdjet.template secondaryVertices_as()) { + int Sxyz = prong.decayLength() / prong.errorDecayLength(); + if (maxSxyz < Sxyz) { + maxSxyz = Sxyz; + } + } + return maxSxyz; + } + template void fillHistogramIPsData(T const& collision, U const& jets, V const& /*jtracks*/, W const& /*tracks*/) { @@ -594,38 +632,77 @@ struct JetTaggerHFQA { template void fillHistogramSV2ProngData(T const& collision, U const& jets, V const& prongs) { - LOG(info) << "2Prong"; + for (const auto& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + continue; + } + auto maxSxy = getMaxSxyForJet(jet); + auto maxSxyz = getMaxSxyzForJet(jet); + registry.fill(HIST("h_2prong_nprongs"), jet.template secondaryVertices_as().size()); + for (const auto& prong : jet.template secondaryVertices_as()) { + auto Lxy = prong.decayLengthXY(); + auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); + auto Lxyz = prong.decayLength(); + auto Sxyz = prong.decayLength() / prong.errorDecayLength(); + registry.fill(HIST("h2_jet_pt_2prong_Lxy"), jet.pt(), Lxy); + registry.fill(HIST("h2_jet_pt_2prong_Sxy"), jet.pt(), Sxy); + registry.fill(HIST("h2_jet_pt_2prong_Lxyz"), jet.pt(), Lxyz); + registry.fill(HIST("h2_jet_pt_2prong_Sxyz"), jet.pt(), Sxyz); + registry.fill(HIST("h2_2prong_Sxy_sigmaLxy"), Sxy, prong.errorDecayLengthXY()); + registry.fill(HIST("h2_2prong_Sxyz_sigmaLxyz"), Sxyz, prong.errorDecayLength()); + registry.fill(HIST("h2_jet_pt_2prong_sigmaLxy"), jet.pt(), prong.errorDecayLengthXY()); + registry.fill(HIST("h2_jet_pt_2prong_sigmaLxyz"), jet.pt(), prong.errorDecayLength()); + } + registry.fill(HIST("h2_jet_pt_2prong_Sxy_N1"), jet.pt(), maxSxy); + registry.fill(HIST("h2_jet_pt_2prong_Sxyz_N1"), jet.pt(), maxSxyz); + } } template void fillHistogramSV3ProngData(T const& collision, U const& jets, V const& prongs) { - LOG(info) << "3Prong"; + for (const auto& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + continue; + } + auto maxSxy = getMaxSxyForJet(jet); + auto maxSxyz = getMaxSxyzForJet(jet); + registry.fill(HIST("h_3prong_nprongs"), jet.template secondaryVertices_as().size()); + for (const auto& prong : jet.template secondaryVertices_as()) { + auto Lxy = prong.decayLengthXY(); + auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); + auto Lxyz = prong.decayLength(); + auto Sxyz = prong.decayLength() / prong.errorDecayLength(); + registry.fill(HIST("h2_jet_pt_3prong_Lxy"), jet.pt(), Lxy); + registry.fill(HIST("h2_jet_pt_3prong_Sxy"), jet.pt(), Sxy); + registry.fill(HIST("h2_jet_pt_3prong_Lxyz"), jet.pt(), Lxyz); + registry.fill(HIST("h2_jet_pt_3prong_Sxyz"), jet.pt(), Sxyz); + registry.fill(HIST("h2_3prong_Sxy_sigmaLxy"), Sxy, prong.errorDecayLengthXY()); + registry.fill(HIST("h2_3prong_Sxyz_sigmaLxyz"), Sxyz, prong.errorDecayLength()); + registry.fill(HIST("h2_jet_pt_3prong_sigmaLxy"), jet.pt(), prong.errorDecayLengthXY()); + registry.fill(HIST("h2_jet_pt_3prong_sigmaLxyz"), jet.pt(), prong.errorDecayLength()); + } + registry.fill(HIST("h2_jet_pt_3prong_Sxy_N1"), jet.pt(), maxSxy); + registry.fill(HIST("h2_jet_pt_3prong_Sxyz_N1"), jet.pt(), maxSxyz); + } } - template void fillHistogramSV2ProngMCD(T const& collision, U const& mcdjets, V const& prongs) { - int maxSxy = 0; - int maxSxyz = 0; for (const auto& mcdjet : mcdjets) { - int numOfSV = 0; if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { continue; } auto origin = mcdjet.origin(); - registry.fill(HIST("h_2prong_nprongs"), mcdjet.template secondaryVertices_as().size()); + auto maxSxy = getMaxSxyForJet(mcdjet); + auto maxSxyz = getMaxSxyzForJet(mcdjet); + registry.fill(HIST("h2_2prong_nprongs_flavour"), mcdjet.template secondaryVertices_as().size(), origin); for (const auto& prong : mcdjet.template secondaryVertices_as()) { - numOfSV++; auto Lxy = prong.decayLengthXY(); auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); auto Lxyz = prong.decayLength(); auto Sxyz = prong.decayLength() / prong.errorDecayLength(); - if (maxSxy < Sxy) - maxSxy = Sxy; - if (maxSxyz < Sxyz) - maxSxyz = Sxyz; registry.fill(HIST("h3_jet_pt_2prong_Lxy_flavour"), mcdjet.pt(), Lxy, origin); registry.fill(HIST("h3_jet_pt_2prong_Sxy_flavour"), mcdjet.pt(), Sxy, origin); registry.fill(HIST("h3_jet_pt_2prong_Lxyz_flavour"), mcdjet.pt(), Lxyz, origin); @@ -637,31 +714,25 @@ struct JetTaggerHFQA { } registry.fill(HIST("h3_jet_pt_2prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); registry.fill(HIST("h3_jet_pt_2prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin); - registry.fill(HIST("h_2prong_nprongs"), numOfSV); } } template void fillHistogramSV3ProngMCD(T const& collision, U const& mcdjets, V const& prongs) { - int numOfSV = 0; - int maxSxy = 0; - int maxSxyz = 0; for (const auto& mcdjet : mcdjets) { if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { continue; } auto origin = mcdjet.origin(); + auto maxSxy = getMaxSxyForJet(mcdjet); + auto maxSxyz = getMaxSxyzForJet(mcdjet); + registry.fill(HIST("h2_3prong_nprongs_flavour"), mcdjet.template secondaryVertices_as().size(), origin); for (const auto& prong : mcdjet.template secondaryVertices_as()) { - numOfSV++; auto Lxy = prong.decayLengthXY(); auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); auto Lxyz = prong.decayLength(); auto Sxyz = prong.decayLength() / prong.errorDecayLength(); - if (maxSxy < Sxy) - maxSxy = Sxy; - if (maxSxyz < Sxyz) - maxSxyz = Sxyz; registry.fill(HIST("h3_jet_pt_3prong_Lxy_flavour"), mcdjet.pt(), Lxy, origin); registry.fill(HIST("h3_jet_pt_3prong_Sxy_flavour"), mcdjet.pt(), Sxy, origin); registry.fill(HIST("h3_jet_pt_3prong_Lxyz_flavour"), mcdjet.pt(), Lxyz, origin); @@ -671,10 +742,9 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_jet_pt_3prong_sigmaLxy_flavour"), mcdjet.pt(), prong.errorDecayLengthXY(), origin); registry.fill(HIST("h3_jet_pt_3prong_sigmaLxyz_flavour"), mcdjet.pt(), prong.errorDecayLength(), origin); } - registry.fill(HIST("h3_jet_pt_2prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); - registry.fill(HIST("h3_jet_pt_2prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin); + registry.fill(HIST("h3_jet_pt_3prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); + registry.fill(HIST("h3_jet_pt_3prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin); } - registry.fill(HIST("h_3prong_nprongs"), numOfSV); } void processDummy(aod::Collision const&, aod::Tracks const&) @@ -708,19 +778,19 @@ struct JetTaggerHFQA { registry.fill(HIST("h_impact_parameter_xyz_significance"), varImpXYZSig); } } - PROCESS_SWITCH(JetTaggerHFQA, processTracksDca, "Fill inclusive tracks' impormation for data", false); + PROCESS_SWITCH(JetTaggerHFQA, processTracksDca, "Fill inclusive tracks' imformation for data", false); void processIPsData(soa::Filtered::iterator const& jcollision, JetTagTableData const& jets, JetTagTracksData const& jtracks, OriTracksData const& tracks) { fillHistogramIPsData(jcollision, jets, jtracks, tracks); } - PROCESS_SWITCH(JetTaggerHFQA, processIPsData, "Fill impact parameter impormation for data jets", false); + PROCESS_SWITCH(JetTaggerHFQA, processIPsData, "Fill impact parameter imformation for data jets", false); void processIPsMCD(soa::Filtered::iterator const& jcollision, JetTagTableMCD const& mcdjets, JetTagTracksMCD const& jtracks, OriTracksMCD const& tracks, JetParticles&) { fillHistogramIPsMCD(jcollision, mcdjets, jtracks, tracks); } - PROCESS_SWITCH(JetTaggerHFQA, processIPsMCD, "Fill impact parameter impormation for mcd jets", false); + PROCESS_SWITCH(JetTaggerHFQA, processIPsMCD, "Fill impact parameter imformation for mcd jets", false); void processJPData(soa::Filtered::iterator const& jcollision, JetTagTableData const& jets, JetTagTracksData const&) { @@ -736,27 +806,27 @@ struct JetTaggerHFQA { void processSV2ProngData(soa::Filtered::iterator const& jcollision, soa::Join const& jets, aod::DataSecondaryVertex2Prongs const& prongs) { - fillHistogramSV2ProngMCD(jcollision, jets, prongs); + fillHistogramSV2ProngData(jcollision, jets, prongs); } - PROCESS_SWITCH(JetTaggerHFQA, processSV2ProngData, "Fill 2prong impormation for data jets", false); + PROCESS_SWITCH(JetTaggerHFQA, processSV2ProngData, "Fill 2prong imformation for data jets", false); void processSV3ProngData(soa::Filtered::iterator const& jcollision, soa::Join const& jets, aod::DataSecondaryVertex3Prongs const& prongs) { - fillHistogramSV3ProngMCD(jcollision, jets, prongs); + fillHistogramSV3ProngData(jcollision, jets, prongs); } - PROCESS_SWITCH(JetTaggerHFQA, processSV3ProngData, "Fill 3prong impormation for data jets", false); + PROCESS_SWITCH(JetTaggerHFQA, processSV3ProngData, "Fill 2prong imformation for data jets", false); void processSV2ProngMCD(soa::Filtered::iterator const& jcollision, soa::Join const& mcdjets, aod::MCDSecondaryVertex2Prongs const& prongs) { fillHistogramSV2ProngMCD(jcollision, mcdjets, prongs); } - PROCESS_SWITCH(JetTaggerHFQA, processSV2ProngMCD, "Fill 2prong impormation for mcd jets", false); + PROCESS_SWITCH(JetTaggerHFQA, processSV2ProngMCD, "Fill 2prong imformation for mcd jets", false); void processSV3ProngMCD(soa::Filtered::iterator const& jcollision, soa::Join const& mcdjets, aod::MCDSecondaryVertex3Prongs const& prongs) { fillHistogramSV3ProngMCD(jcollision, mcdjets, prongs); } - PROCESS_SWITCH(JetTaggerHFQA, processSV3ProngMCD, "Fill 3prong impormation for mcd jets", false); + PROCESS_SWITCH(JetTaggerHFQA, processSV3ProngMCD, "Fill 3prong imformation for mcd jets", false); }; using JetTaggerQAChargedDataJets = soa::Join; From 0194505a0fe482d9fdfacde950bb5789020b82f3 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Fri, 21 Jun 2024 16:55:36 +0200 Subject: [PATCH 10/21] Add distribution of jet pt with flavour when removed by cut selection for efficiency and purity using sv --- PWGJE/Tasks/jettaggerhfQA.cxx | 100 ++++++++++++++++++++++++++-------- 1 file changed, 76 insertions(+), 24 deletions(-) diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index bf4ae40865a..6b6e9d8a026 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -53,13 +53,13 @@ struct JetTaggerHFQA { Configurable trackPtMax{"trackPtMax", 100.0, "maximum pT acceptance for tracks"}; Configurable jetEtaMin{"jetEtaMin", -99.0, "minimum jet pseudorapidity"}; Configurable jetEtaMax{"jetEtaMax", 99.0, "maximum jet pseudorapidity"}; - Configurable prong2LxySigmaMax{"prong2LxySigmaMax", 0.03, "maximum sigma of decay length of 2-prong on xy plane"}; + Configurable prong2sigmaLxyMax{"prong2simgaLxyMax", 0.03, "maximum sigma of decay length of 2-prong on xy plane"}; Configurable prong2SxyMin{"prong2SxyMin", 7, "minimum decay length significance of 2-prong on xy plane"}; - Configurable prong2LxyzSigmaMax{"prong2LxyzSigmazMax", 0.03, "maximum sigma of decay length of 2-prong on xyz plane"}; + Configurable prong2sigmaLxyzMax{"prong2sigmaLxyzMax", 0.03, "maximum sigma of decay length of 2-prong on xyz plane"}; Configurable prong2SxyzMin{"prong2SxyzMin", 7, "minimum decay length significance of 2-prong on xyz plane"}; - Configurable prong3LxySigmaMax{"prong3LxySigmaMax", 0.03, "maximum sigma of decay length of 3-prong on xy plane"}; + Configurable prong3sigmaLxyMax{"prong3sigmaLxyMax", 0.03, "maximum sigma of decay length of 3-prong on xy plane"}; Configurable prong3SxyMin{"prong3SxyMin", 7, "minimum decay length significance of 3-prong on xy plane"}; - Configurable prong3LxyzSigmaMax{"prong3LxyzSigmazMax", 0.03, "maximum sigma of decay length of 3-prong on xyz plane"}; + Configurable prong3sigmaLxyzMax{"prong3sigmaLxyzMax", 0.03, "maximum sigma of decay length of 3-prong on xyz plane"}; Configurable prong3SxyzMin{"prong3SxyzMin", 7, "minimum decay length significance of 3-prong on xyz plane"}; Configurable numFlavourSpecies{"numFlavourSpecies", 6, "number of jet flavour species"}; Configurable numOrder{"numOrder", 6, "number of ordering"}; @@ -91,6 +91,14 @@ struct JetTaggerHFQA { int numberOfJetFlavourSpecies = 6; int trackSelection = -1; + float maxSigmaLxy2Prong = 0.; + float minSxy2Prong = 0.; + float maxSigmaLxyz2Prong = 0.; + float minSxyz2Prong = 0.; + float maxSigmaLxy3Prong = 0.; + float minSxy3Prong = 0.; + float maxSigmaLxyz3Prong = 0.; + float minSxyz3Prong = 0.; HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -121,6 +129,14 @@ struct JetTaggerHFQA { AxisSpec sigmaLxyzAxis = {binSigmaLxyz, "#sigma_{L_{XYZ}} [cm]"}; numberOfJetFlavourSpecies = static_cast(numFlavourSpecies); + maxSigmaLxy2Prong = static_cast(prong2sigmaLxyMax); + minSxy2Prong = static_cast(prong2SxyMin); + maxSigmaLxyz2Prong = static_cast(prong2sigmaLxyzMax); + minSxyz2Prong = static_cast(prong2SxyzMin); + maxSigmaLxy3Prong = static_cast(prong3sigmaLxyMax); + minSxy3Prong = static_cast(prong3SxyMin); + maxSigmaLxyz3Prong = static_cast(prong3sigmaLxyzMax); + minSxyz3Prong = static_cast(prong3SxyzMin); trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); if (doprocessTracksDca) { @@ -271,6 +287,10 @@ struct JetTaggerHFQA { registry.add("h3_2prong_Sxyz_sigmaLxyz_flavour", "", {HistType::kTH3F, {{SxyzAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_sigmaLxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_sigmaLxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_2prong_Sxy_flavour_cutSxyAndsigmaLxy", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_2prong_Sxyz_flavour_cutSxyzAndsigmaLxyz", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_2prong_Sxy_N1_flavour_cutSxyAndsigmaLxy", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_2prong_Sxyz_N1_flavour_cutSxyzAndsigmaLxyz", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); } if (doprocessSV3ProngMCD) { registry.add("h2_3prong_nprongs_flavour", "", {HistType::kTH2F, {{nprongsAxis}, {jetFlavourAxis}}}); @@ -284,6 +304,10 @@ struct JetTaggerHFQA { registry.add("h3_3prong_Sxyz_sigmaLxyz_flavour", "", {HistType::kTH3F, {{SxyzAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_sigmaLxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_sigmaLxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_3prong_Sxy_flavour_cutSxyAndsigmaLxy", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_3prong_Sxyz_flavour_cutSxyzAndsigmaLxyz", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_3prong_Sxy_N1_flavour_cutSxyAndsigmaLxy", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_3prong_Sxyz_N1_flavour_cutSxyzAndsigmaLxyz", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); } } @@ -309,40 +333,44 @@ struct JetTaggerHFQA { return 1; } - template - bool prongAcceptance(T const& prong, U const& maxDecayLengthSigma, V const& minDecayLengthSig) + bool prongAcceptance(float sigmaDecayLength, float decayLengthSig, float maxSigmaDecayLength, float minDecayLengthSig) { - auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); - if (prong.errorDecayLengthXY() < maxDecayLengthSigma || Sxy > minDecayLengthSig) + if ((sigmaDecayLength < maxSigmaDecayLength) && (decayLengthSig > minDecayLengthSig)) return 0; return 1; } template - double getMaxSxyForJet(const JetType& mcdjet) + std::tuple getMaxSxyForJet(const JetType& mcdjet) { - double maxSxy = 0; + float maxSxy = 0; + float correspondingErrorDecayLengthXY = 0; + for (const auto& prong : mcdjet.template secondaryVertices_as()) { - int Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); + float Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); if (maxSxy < Sxy) { maxSxy = Sxy; + correspondingErrorDecayLengthXY = prong.errorDecayLengthXY(); } } - return maxSxy; + return std::make_tuple(maxSxy, correspondingErrorDecayLengthXY); } template - double getMaxSxyzForJet(const JetType& mcdjet) + std::tuple getMaxSxyzForJet(const JetType& mcdjet) { - double maxSxyz = 0; + float maxSxyz = 0; + float correspondingErrorDecayLength = 0; + for (const auto& prong : mcdjet.template secondaryVertices_as()) { - int Sxyz = prong.decayLength() / prong.errorDecayLength(); + float Sxyz = prong.decayLength() / prong.errorDecayLength(); if (maxSxyz < Sxyz) { maxSxyz = Sxyz; + correspondingErrorDecayLength = prong.errorDecayLength(); } } - return maxSxyz; + return std::make_tuple(maxSxyz, correspondingErrorDecayLength); } template @@ -636,8 +664,8 @@ struct JetTaggerHFQA { if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { continue; } - auto maxSxy = getMaxSxyForJet(jet); - auto maxSxyz = getMaxSxyzForJet(jet); + auto [maxSxy, sigmaLxy] = getMaxSxyForJet(jet); + auto [maxSxyz, sigmaLxyz] = getMaxSxyzForJet(jet); registry.fill(HIST("h_2prong_nprongs"), jet.template secondaryVertices_as().size()); for (const auto& prong : jet.template secondaryVertices_as()) { auto Lxy = prong.decayLengthXY(); @@ -665,8 +693,8 @@ struct JetTaggerHFQA { if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { continue; } - auto maxSxy = getMaxSxyForJet(jet); - auto maxSxyz = getMaxSxyzForJet(jet); + auto [maxSxy, sigmaLxy] = getMaxSxyForJet(jet); + auto [maxSxyz, sigmaLxyz] = getMaxSxyzForJet(jet); registry.fill(HIST("h_3prong_nprongs"), jet.template secondaryVertices_as().size()); for (const auto& prong : jet.template secondaryVertices_as()) { auto Lxy = prong.decayLengthXY(); @@ -695,8 +723,8 @@ struct JetTaggerHFQA { continue; } auto origin = mcdjet.origin(); - auto maxSxy = getMaxSxyForJet(mcdjet); - auto maxSxyz = getMaxSxyzForJet(mcdjet); + auto [maxSxy, sigmaLxy] = getMaxSxyForJet(mcdjet); + auto [maxSxyz, sigmaLxyz] = getMaxSxyzForJet(mcdjet); registry.fill(HIST("h2_2prong_nprongs_flavour"), mcdjet.template secondaryVertices_as().size(), origin); for (const auto& prong : mcdjet.template secondaryVertices_as()) { auto Lxy = prong.decayLengthXY(); @@ -711,9 +739,21 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_2prong_Sxyz_sigmaLxyz_flavour"), Sxyz, prong.errorDecayLength(), origin); registry.fill(HIST("h3_jet_pt_2prong_sigmaLxy_flavour"), mcdjet.pt(), prong.errorDecayLengthXY(), origin); registry.fill(HIST("h3_jet_pt_2prong_sigmaLxyz_flavour"), mcdjet.pt(), prong.errorDecayLength(), origin); + if (prongAcceptance(prong.errorDecayLengthXY(), Sxy, maxSigmaLxy2Prong, minSxy2Prong)) { + registry.fill(HIST("h3_jet_pt_2prong_Sxy_flavour_cutSxyAndsigmaLxy"), mcdjet.pt(), Sxy, origin); + } + if (prongAcceptance(prong.errorDecayLength(), Sxyz, maxSigmaLxyz2Prong, minSxyz2Prong)) { + registry.fill(HIST("h3_jet_pt_2prong_Sxyz_flavour_cutSxyzAndsigmaLxyz"), mcdjet.pt(), Sxyz, origin); + } } registry.fill(HIST("h3_jet_pt_2prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); registry.fill(HIST("h3_jet_pt_2prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin); + if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy2Prong, minSxy2Prong)) { + registry.fill(HIST("h3_jet_pt_2prong_Sxy_N1_flavour_cutSxyAndsigmaLxy"), mcdjet.pt(), maxSxy, origin); + } + if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy2Prong, minSxy2Prong)) { + registry.fill(HIST("h3_jet_pt_2prong_Sxyz_N1_flavour_cutSxyzAndsigmaLxyz"), mcdjet.pt(), maxSxyz, origin); + } } } @@ -725,8 +765,8 @@ struct JetTaggerHFQA { continue; } auto origin = mcdjet.origin(); - auto maxSxy = getMaxSxyForJet(mcdjet); - auto maxSxyz = getMaxSxyzForJet(mcdjet); + auto [maxSxy, sigmaLxy] = getMaxSxyForJet(mcdjet); + auto [maxSxyz, sigmaLxyz] = getMaxSxyzForJet(mcdjet); registry.fill(HIST("h2_3prong_nprongs_flavour"), mcdjet.template secondaryVertices_as().size(), origin); for (const auto& prong : mcdjet.template secondaryVertices_as()) { auto Lxy = prong.decayLengthXY(); @@ -741,9 +781,21 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_3prong_Sxyz_sigmaLxyz_flavour"), Sxyz, prong.errorDecayLength(), origin); registry.fill(HIST("h3_jet_pt_3prong_sigmaLxy_flavour"), mcdjet.pt(), prong.errorDecayLengthXY(), origin); registry.fill(HIST("h3_jet_pt_3prong_sigmaLxyz_flavour"), mcdjet.pt(), prong.errorDecayLength(), origin); + if (prongAcceptance(prong.errorDecayLengthXY(), Sxy, maxSigmaLxy3Prong, minSxy3Prong)) { + registry.fill(HIST("h3_jet_pt_3prong_Sxy_flavour_cutSxyAndsigmaLxy"), mcdjet.pt(), Sxy, origin); + } + if (prongAcceptance(prong.errorDecayLength(), Sxyz, maxSigmaLxyz3Prong, minSxyz3Prong)) { + registry.fill(HIST("h3_jet_pt_3prong_Sxyz_flavour_cutSxyzAndsigmaLxyz"), mcdjet.pt(), Sxyz, origin); + } } registry.fill(HIST("h3_jet_pt_3prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); registry.fill(HIST("h3_jet_pt_3prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin); + if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy3Prong, minSxy3Prong)) { + registry.fill(HIST("h3_jet_pt_3prong_Sxy_N1_flavour_cutSxyAndsigmaLxy"), mcdjet.pt(), maxSxy, origin); + } + if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy3Prong, minSxy3Prong)) { + registry.fill(HIST("h3_jet_pt_3prong_Sxyz_N1_flavour_cutSxyzAndsigmaLxyz"), mcdjet.pt(), maxSxyz, origin); + } } } From 98e3ea79cbdf9e3d6a03aef48e8ce584182075cb Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Fri, 28 Jun 2024 04:32:01 +0200 Subject: [PATCH 11/21] updating sv --- PWGJE/Tasks/jettaggerhfQA.cxx | 39 ++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index 2aa17f5d006..fdebdd20b0e 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -73,6 +73,7 @@ struct JetTaggerHFQA { ConfigurableAxis binNtracks{"binNtracks", {100, 0., 100.}, ""}; ConfigurableAxis binTrackPt{"binTrackPt", {200, 0.f, 100.f}, ""}; ConfigurableAxis binImpactParameterXY{"binImpactParameterXY", {800, -400.5f, 400.5f}, ""}; + ConfigurableAxis binSigmaImpactParameterXY{"binImpactSigmaParameterXY", {800, 0.f, 100.f}, ""}; ConfigurableAxis binImpactParameterXYSignificance{"binImpactParameterXYSignificance", {800, -40.5f, 40.5f}, ""}; ConfigurableAxis binImpactParameterZ{"binImpactParameterZ", {800, -400.5f, 400.5f}, ""}; ConfigurableAxis binImpactParameterZSignificance{"binImpactParameterZSignificance", {800, -40.5f, 40.5f}, ""}; @@ -112,6 +113,7 @@ struct JetTaggerHFQA { AxisSpec ntracksAxis = {binNtracks, "#it{N}_{tracks}"}; AxisSpec trackPtAxis = {binTrackPt, "#it{p}_{T}^{track}"}; AxisSpec impactParameterXYAxis = {binImpactParameterXY, "IP_{XY} [#mum]"}; + AxisSpec sigmaImpactParameterXYAxis = {binSigmaImpactParameterXY, "#sigma_{XY} [#mum]"}; AxisSpec impactParameterXYSignificanceAxis = {binImpactParameterXYSignificance, "IPs_{XY}"}; AxisSpec impactParameterZAxis = {binImpactParameterZ, "IP_{Z} [#mum]"}; AxisSpec impactParameterZSignificanceAxis = {binImpactParameterZSignificance, "IPs_{Z}"}; @@ -186,6 +188,7 @@ struct JetTaggerHFQA { registry.add("h3_jet_pt_track_phi_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {phiAxis}, {jetFlavourAxis}}}); if (fillIPxy) { registry.add("h3_jet_pt_impact_parameter_xy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {impactParameterXYAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_sigma_impact_parameter_xy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaImpactParameterXYAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_sign_impact_parameter_xy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {impactParameterXYAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_impact_parameter_xy_significance_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {impactParameterXYSignificanceAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_sign_impact_parameter_xy_significance_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {impactParameterXYSignificanceAxis}, {jetFlavourAxis}}}); @@ -261,6 +264,10 @@ struct JetTaggerHFQA { registry.add("h2_2prong_Sxyz_sigmaLxyz", "", {HistType::kTH2F, {{SxyzAxis}, {sigmaLxyzAxis}}}); registry.add("h2_jet_pt_2prong_sigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyAxis}}}); registry.add("h2_jet_pt_2prong_sigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyzAxis}}}); + registry.add("h2_jet_pt_2prong_Sxy_cutSxyAndsigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); + registry.add("h2_jet_pt_2prong_Sxyz_cutSxyzAndsigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); + registry.add("h2_jet_pt_2prong_Sxy_N1_cutSxyAndsigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); + registry.add("h2_jet_pt_2prong_Sxyz_N1_cutSxyzAndsigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); } if (doprocessSV3ProngData) { registry.add("h_3prong_nprongs", "", {HistType::kTH1F, {{nprongsAxis}}}); @@ -274,6 +281,10 @@ struct JetTaggerHFQA { registry.add("h2_3prong_Sxyz_sigmaLxyz", "", {HistType::kTH2F, {{SxyzAxis}, {sigmaLxyzAxis}}}); registry.add("h2_jet_pt_3prong_sigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyAxis}}}); registry.add("h2_jet_pt_3prong_sigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyzAxis}}}); + registry.add("h2_jet_pt_3prong_Sxy_cutSxyAndsigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); + registry.add("h2_jet_pt_3prong_Sxyz_cutSxyzAndsigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); + registry.add("h2_jet_pt_3prong_Sxy_N1_cutSxyAndsigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); + registry.add("h2_jet_pt_3prong_Sxyz_N1_cutSxyzAndsigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); } if (doprocessSV2ProngMCD) { registry.add("h2_2prong_nprongs_flavour", "", {HistType::kTH2F, {{nprongsAxis}, {jetFlavourAxis}}}); @@ -491,10 +502,12 @@ struct JetTaggerHFQA { if (fillIPxy) { float varImpXY, varSignImpXY, varImpXYSig, varSignImpXYSig; varImpXY = track.dcaXY() * jettaggingutilities::cmTomum; + float varSigmaImpXY = track.dcaXY() * jettaggingutilities::cmTomum; varSignImpXY = geoSign * TMath::Abs(track.dcaXY()) * jettaggingutilities::cmTomum; - varImpXYSig = track.dcaXY() / TMath::Sqrt(track.sigmaDcaXY2()); + varImpXYSig = TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); varSignImpXYSig = geoSign * TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); registry.fill(HIST("h3_jet_pt_impact_parameter_xy_flavour"), mcdjet.pt(), varImpXY, jetflavour); + registry.fill(HIST("h3_jet_pt_sigma_impact_parameter_xy_flavour"), mcdjet.pt(), varSigmaImpXY, jetflavour); registry.fill(HIST("h3_jet_pt_sign_impact_parameter_xy_flavour"), mcdjet.pt(), varSignImpXY, jetflavour); registry.fill(HIST("h3_jet_pt_impact_parameter_xy_significance_flavour"), mcdjet.pt(), varImpXYSig, jetflavour); registry.fill(HIST("h3_jet_pt_sign_impact_parameter_xy_significance_flavour"), mcdjet.pt(), varSignImpXYSig, jetflavour); @@ -680,9 +693,21 @@ struct JetTaggerHFQA { registry.fill(HIST("h2_2prong_Sxyz_sigmaLxyz"), Sxyz, prong.errorDecayLength()); registry.fill(HIST("h2_jet_pt_2prong_sigmaLxy"), jet.pt(), prong.errorDecayLengthXY()); registry.fill(HIST("h2_jet_pt_2prong_sigmaLxyz"), jet.pt(), prong.errorDecayLength()); + if (prongAcceptance(prong.errorDecayLengthXY(), Sxy, maxSigmaLxy2Prong, minSxy2Prong)) { + registry.fill(HIST("h2_jet_pt_2prong_Sxy_cutSxyAndsigmaLxy"), jet.pt(), Sxy); + } + if (prongAcceptance(prong.errorDecayLength(), Sxyz, maxSigmaLxyz2Prong, minSxyz2Prong)) { + registry.fill(HIST("h2_jet_pt_2prong_Sxyz_cutSxyzAndsigmaLxyz"), jet.pt(), Sxyz); + } } registry.fill(HIST("h2_jet_pt_2prong_Sxy_N1"), jet.pt(), maxSxy); registry.fill(HIST("h2_jet_pt_2prong_Sxyz_N1"), jet.pt(), maxSxyz); + if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy2Prong, minSxy2Prong)) { + registry.fill(HIST("h2_jet_pt_2prong_Sxy_N1_cutSxyAndsigmaLxy"), jet.pt(), maxSxy); + } + if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy2Prong, minSxy2Prong)) { + registry.fill(HIST("h2_jet_pt_2prong_Sxyz_N1_cutSxyzAndsigmaLxyz"), jet.pt(), maxSxyz); + } } } @@ -709,9 +734,21 @@ struct JetTaggerHFQA { registry.fill(HIST("h2_3prong_Sxyz_sigmaLxyz"), Sxyz, prong.errorDecayLength()); registry.fill(HIST("h2_jet_pt_3prong_sigmaLxy"), jet.pt(), prong.errorDecayLengthXY()); registry.fill(HIST("h2_jet_pt_3prong_sigmaLxyz"), jet.pt(), prong.errorDecayLength()); + if (prongAcceptance(prong.errorDecayLengthXY(), Sxy, maxSigmaLxy3Prong, minSxy3Prong)) { + registry.fill(HIST("h2_jet_pt_3prong_Sxy_cutSxyAndsigmaLxy"), jet.pt(), Sxy); + } + if (prongAcceptance(prong.errorDecayLength(), Sxyz, maxSigmaLxyz3Prong, minSxyz3Prong)) { + registry.fill(HIST("h2_jet_pt_3prong_Sxyz_cutSxyzAndsigmaLxyz"), jet.pt(), Sxyz); + } } registry.fill(HIST("h2_jet_pt_3prong_Sxy_N1"), jet.pt(), maxSxy); registry.fill(HIST("h2_jet_pt_3prong_Sxyz_N1"), jet.pt(), maxSxyz); + if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy3Prong, minSxy3Prong)) { + registry.fill(HIST("h2_jet_pt_3prong_Sxy_N1_cutSxyAndsigmaLxy"), jet.pt(), maxSxy); + } + if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy3Prong, minSxy3Prong)) { + registry.fill(HIST("h2_jet_pt_3prong_Sxyz_N1_cutSxyzAndsigmaLxyz"), jet.pt(), maxSxyz); + } } } From 11a131fc9169cc7537b12895a0e36336d7034ed1 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Tue, 2 Jul 2024 13:00:52 +0200 Subject: [PATCH 12/21] fix prong acceptance --- PWGJE/Tasks/jettaggerhfQA.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index fdebdd20b0e..726f2a8036a 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -347,9 +347,9 @@ struct JetTaggerHFQA { bool prongAcceptance(float sigmaDecayLength, float decayLengthSig, float maxSigmaDecayLength, float minDecayLengthSig) { if ((sigmaDecayLength < maxSigmaDecayLength) && (decayLengthSig > minDecayLengthSig)) - return 0; + return 1; - return 1; + return 0; } template From 01f9cd5ac01e60a872c16eb6e1b66e18e5fe753b Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Tue, 2 Jul 2024 17:14:23 +0200 Subject: [PATCH 13/21] fix bool and TMath value --- PWGJE/Tasks/jettaggerhfQA.cxx | 52 +++++++++++++++++------------------ 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index 726f2a8036a..0850eba10f1 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -339,17 +339,17 @@ struct JetTaggerHFQA { bool trackAcceptance(T const& track) { if (track.pt() < trackPtMin || track.pt() > trackPtMax) - return 0; + return false; - return 1; + return true; } bool prongAcceptance(float sigmaDecayLength, float decayLengthSig, float maxSigmaDecayLength, float minDecayLengthSig) { if ((sigmaDecayLength < maxSigmaDecayLength) && (decayLengthSig > minDecayLengthSig)) - return 1; + return true; - return 0; + return false; } template @@ -404,9 +404,9 @@ struct JetTaggerHFQA { if (fillIPxy) { float varImpXY, varSignImpXY, varImpXYSig, varSignImpXYSig; varImpXY = track.dcaXY() * jettaggingutilities::cmTomum; - varSignImpXY = geoSign * TMath::Abs(track.dcaXY()) * jettaggingutilities::cmTomum; - varImpXYSig = track.dcaXY() / TMath::Sqrt(track.sigmaDcaXY2()); - varSignImpXYSig = geoSign * TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); + varSignImpXY = geoSign * std::abs(track.dcaXY()) * jettaggingutilities::cmTomum; + varImpXYSig = track.dcaXY() / std::sqrt(track.sigmaDcaXY2()); + varSignImpXYSig = geoSign * std::abs(track.dcaXY()) / std::sqrt(track.sigmaDcaXY2()); registry.fill(HIST("h2_jet_pt_impact_parameter_xy"), jet.pt(), varImpXY); registry.fill(HIST("h2_jet_pt_sign_impact_parameter_xy"), jet.pt(), varSignImpXY); registry.fill(HIST("h2_jet_pt_impact_parameter_xy_significance"), jet.pt(), varImpXYSig); @@ -416,9 +416,9 @@ struct JetTaggerHFQA { if (fillIPz) { float varImpZ, varSignImpZ, varImpZSig, varSignImpZSig; varImpZ = track.dcaZ() * jettaggingutilities::cmTomum; - varSignImpZ = geoSign * TMath::Abs(track.dcaZ()) * jettaggingutilities::cmTomum; - varImpZSig = track.dcaZ() / TMath::Sqrt(track.sigmaDcaZ2()); - varSignImpZSig = geoSign * TMath::Abs(track.dcaZ()) / TMath::Sqrt(track.sigmaDcaZ2()); + varSignImpZ = geoSign * std::abs(track.dcaZ()) * jettaggingutilities::cmTomum; + varImpZSig = track.dcaZ() / std::sqrt(track.sigmaDcaZ2()); + varSignImpZSig = geoSign * std::abs(track.dcaZ()) / std::sqrt(track.sigmaDcaZ2()); registry.fill(HIST("h2_jet_pt_impact_parameter_z"), jet.pt(), varImpZ); registry.fill(HIST("h2_jet_pt_sign_impact_parameter_z"), jet.pt(), varSignImpZ); registry.fill(HIST("h2_jet_pt_impact_parameter_z_significance"), jet.pt(), varImpZSig); @@ -430,9 +430,9 @@ struct JetTaggerHFQA { float dcaXYZ = jtrack.dcaXYZ(); float sigmaDcaXYZ2 = jtrack.sigmaDcaXYZ2(); varImpXYZ = dcaXYZ * jettaggingutilities::cmTomum; - varSignImpXYZ = geoSign * TMath::Abs(dcaXYZ) * jettaggingutilities::cmTomum; - varImpXYZSig = dcaXYZ / TMath::Sqrt(sigmaDcaXYZ2); - varSignImpXYZSig = geoSign * TMath::Abs(dcaXYZ) / TMath::Sqrt(sigmaDcaXYZ2); + varSignImpXYZ = geoSign * std::abs(dcaXYZ) * jettaggingutilities::cmTomum; + varImpXYZSig = dcaXYZ / std::sqrt(sigmaDcaXYZ2); + varSignImpXYZSig = geoSign * std::abs(dcaXYZ) / std::sqrt(sigmaDcaXYZ2); registry.fill(HIST("h2_jet_pt_impact_parameter_xyz"), jet.pt(), varImpXYZ); registry.fill(HIST("h2_jet_pt_sign_impact_parameter_xyz"), jet.pt(), varSignImpXYZ); registry.fill(HIST("h2_jet_pt_impact_parameter_xyz_significance"), jet.pt(), varImpXYZSig); @@ -503,9 +503,9 @@ struct JetTaggerHFQA { float varImpXY, varSignImpXY, varImpXYSig, varSignImpXYSig; varImpXY = track.dcaXY() * jettaggingutilities::cmTomum; float varSigmaImpXY = track.dcaXY() * jettaggingutilities::cmTomum; - varSignImpXY = geoSign * TMath::Abs(track.dcaXY()) * jettaggingutilities::cmTomum; - varImpXYSig = TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); - varSignImpXYSig = geoSign * TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); + varSignImpXY = geoSign * std::abs(track.dcaXY()) * jettaggingutilities::cmTomum; + varImpXYSig = std::abs(track.dcaXY()) / std::sqrt(track.sigmaDcaXY2()); + varSignImpXYSig = geoSign * std::abs(track.dcaXY()) / std::sqrt(track.sigmaDcaXY2()); registry.fill(HIST("h3_jet_pt_impact_parameter_xy_flavour"), mcdjet.pt(), varImpXY, jetflavour); registry.fill(HIST("h3_jet_pt_sigma_impact_parameter_xy_flavour"), mcdjet.pt(), varSigmaImpXY, jetflavour); registry.fill(HIST("h3_jet_pt_sign_impact_parameter_xy_flavour"), mcdjet.pt(), varSignImpXY, jetflavour); @@ -524,9 +524,9 @@ struct JetTaggerHFQA { if (fillIPz) { float varImpZ, varSignImpZ, varImpZSig, varSignImpZSig; varImpZ = track.dcaZ() * jettaggingutilities::cmTomum; - varSignImpZ = geoSign * TMath::Abs(track.dcaZ()) * jettaggingutilities::cmTomum; - varImpZSig = track.dcaZ() / TMath::Sqrt(track.sigmaDcaZ2()); - varSignImpZSig = geoSign * TMath::Abs(track.dcaZ()) / TMath::Sqrt(track.sigmaDcaZ2()); + varSignImpZ = geoSign * std::abs(track.dcaZ()) * jettaggingutilities::cmTomum; + varImpZSig = track.dcaZ() / std::sqrt(track.sigmaDcaZ2()); + varSignImpZSig = geoSign * std::abs(track.dcaZ()) / std::sqrt(track.sigmaDcaZ2()); registry.fill(HIST("h3_jet_pt_impact_parameter_z_flavour"), mcdjet.pt(), varImpZ, jetflavour); registry.fill(HIST("h3_jet_pt_sign_impact_parameter_z_flavour"), mcdjet.pt(), varSignImpZ, jetflavour); registry.fill(HIST("h3_jet_pt_impact_parameter_z_significance_flavour"), mcdjet.pt(), varImpZSig, jetflavour); @@ -546,9 +546,9 @@ struct JetTaggerHFQA { float dcaXYZ = jtrack.dcaXYZ(); float sigmaDcaXYZ2 = jtrack.sigmaDcaXYZ2(); varImpXYZ = dcaXYZ * jettaggingutilities::cmTomum; - varSignImpXYZ = geoSign * TMath::Abs(dcaXYZ) * jettaggingutilities::cmTomum; - varImpXYZSig = dcaXYZ / TMath::Sqrt(sigmaDcaXYZ2); - varSignImpXYZSig = geoSign * TMath::Abs(dcaXYZ) / TMath::Sqrt(sigmaDcaXYZ2); + varSignImpXYZ = geoSign * std::abs(dcaXYZ) * jettaggingutilities::cmTomum; + varImpXYZSig = dcaXYZ / std::sqrt(sigmaDcaXYZ2); + varSignImpXYZSig = geoSign * std::abs(dcaXYZ) / std::sqrt(sigmaDcaXYZ2); registry.fill(HIST("h3_jet_pt_impact_parameter_xyz_flavour"), mcdjet.pt(), varImpXYZ, jetflavour); registry.fill(HIST("h3_jet_pt_sign_impact_parameter_xyz_flavour"), mcdjet.pt(), varSignImpXYZ, jetflavour); registry.fill(HIST("h3_jet_pt_impact_parameter_xyz_significance_flavour"), mcdjet.pt(), varImpXYZSig, jetflavour); @@ -642,7 +642,7 @@ struct JetTaggerHFQA { continue; } registry.fill(HIST("h2_jet_pt_JP"), jet.pt(), jet.jetProb()[0]); - registry.fill(HIST("h2_jet_pt_neg_log_JP"), jet.pt(), -1 * TMath::Log(jet.jetProb()[0])); + registry.fill(HIST("h2_jet_pt_neg_log_JP"), jet.pt(), -1 * std::log(jet.jetProb()[0])); registry.fill(HIST("h2_jet_pt_JP_N1"), jet.pt(), jet.jetProb()[1]); registry.fill(HIST("h2_jet_pt_neg_log_JP_N1"), jet.pt(), -1 * TMath::Log(jet.jetProb()[1])); registry.fill(HIST("h2_jet_pt_JP_N2"), jet.pt(), jet.jetProb()[2]); @@ -851,13 +851,13 @@ struct JetTaggerHFQA { float varImpXY, varImpXYSig, varImpZ, varImpZSig, varImpXYZ, varImpXYZSig; varImpXY = track.dcaXY() * jettaggingutilities::cmTomum; - varImpXYSig = track.dcaXY() / TMath::Sqrt(track.sigmaDcaXY2()); + varImpXYSig = track.dcaXY() / std::sqrt(track.sigmaDcaXY2()); varImpZ = track.dcaZ() * jettaggingutilities::cmTomum; - varImpZSig = track.dcaZ() / TMath::Sqrt(track.sigmaDcaZ2()); + varImpZSig = track.dcaZ() / std::sqrt(track.sigmaDcaZ2()); float dcaXYZ = jtrack.dcaXYZ(); float sigmaDcaXYZ2 = jtrack.sigmaDcaXYZ2(); varImpXYZ = dcaXYZ * jettaggingutilities::cmTomum; - varImpXYZSig = dcaXYZ / TMath::Sqrt(sigmaDcaXYZ2); + varImpXYZSig = dcaXYZ / std::sqrt(sigmaDcaXYZ2); registry.fill(HIST("h_impact_parameter_xy"), varImpXY); registry.fill(HIST("h_impact_parameter_xy_significance"), varImpXYSig); From 34c62790631f18405e3fb9bf7658ffa64be3695f Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Wed, 10 Jul 2024 17:56:04 +0200 Subject: [PATCH 14/21] Urgent fix tagger point for efficiency and purity --- PWGJE/Core/JetTaggingUtilities.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index e1d8a767d72..5e8793ed718 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -374,7 +374,7 @@ bool isGreaterThanTaggingPoint(T const& collision, U const& jet, V const& jtrack orderForIPJetTracks(collision, jet, jtracks, tracks, vecSignImpSig); if (vecSignImpSig.size() > static_cast::size_type>(cnt) - 1) { for (int i = 0; i < cnt; i++) { - if (0 < vecSignImpSig[i] && vecSignImpSig[i] < taggingPoint) { // tagger point set + if (vecSignImpSig[i] < taggingPoint) { // tagger point set return false; } } From 86510556e019c22d3e9aeecf279a4a2a5b1af830 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Wed, 7 Aug 2024 11:50:20 +0200 Subject: [PATCH 15/21] Add configuration about searchUpToQuark which is chossen between quark and hadron level for flavour definition --- PWGJE/Core/JetTaggingUtilities.h | 16 ++++++++-------- PWGJE/TableProducer/jettaggerhf.cxx | 5 +++-- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index 5e8793ed718..f06e8144b5b 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -109,7 +109,7 @@ int getOriginalHFMotherIndex(const typename T::iterator& hfparticle) * @param hftrack track passed as reference which is then replaced by the first track that originated from an HF shower */ template -int jetTrackFromHFShower(T const& jet, U const& /*tracks*/, V const& particles, typename U::iterator& hftrack) +int jetTrackFromHFShower(T const& jet, U const& /*tracks*/, V const& particles, typename U::iterator& hftrack, bool searchUpToQuark) { bool hasMcParticle = false; @@ -120,7 +120,7 @@ int jetTrackFromHFShower(T const& jet, U const& /*tracks*/, V const& particles, } hasMcParticle = true; auto const& particle = track.template mcParticle_as(); - origin = RecoDecay::getCharmHadronOrigin(particles, particle, true); + origin = RecoDecay::getCharmHadronOrigin(particles, particle, searchUpToQuark); if (origin == 1 || origin == 2) { // 1=charm , 2=beauty hftrack = track; if (origin == 1) { @@ -146,12 +146,12 @@ int jetTrackFromHFShower(T const& jet, U const& /*tracks*/, V const& particles, * @param hfparticle particle passed as reference which is then replaced by the first track that originated from an HF shower */ template -int jetParticleFromHFShower(T const& jet, U const& particles, typename U::iterator& hfparticle) +int jetParticleFromHFShower(T const& jet, U const& particles, typename U::iterator& hfparticle, bool searchUpToQuark) { int origin = -1; for (const auto& particle : jet.template tracks_as()) { - origin = RecoDecay::getCharmHadronOrigin(particles, particle, true); + origin = RecoDecay::getCharmHadronOrigin(particles, particle, searchUpToQuark); if (origin == 1 || origin == 2) { // 1=charm , 2=beauty hfparticle = particle; if (origin == 1) { @@ -174,11 +174,11 @@ int jetParticleFromHFShower(T const& jet, U const& particles, typename U::iterat */ template -int mcdJetFromHFShower(T const& jet, U const& tracks, V const& particles, float dRMax = 0.25) +int mcdJetFromHFShower(T const& jet, U const& tracks, V const& particles, float dRMax = 0.25, bool searchUpToQuark = false) { typename U::iterator hftrack; - int origin = jetTrackFromHFShower(jet, tracks, particles, hftrack); + int origin = jetTrackFromHFShower(jet, tracks, particles, hftrack, searchUpToQuark); if (origin == JetTaggingSpecies::charm || origin == JetTaggingSpecies::beauty) { if (!hftrack.has_mcParticle()) { return JetTaggingSpecies::none; @@ -215,11 +215,11 @@ int mcdJetFromHFShower(T const& jet, U const& tracks, V const& particles, float */ template -int mcpJetFromHFShower(T const& jet, U const& particles, float dRMax = 0.25) +int mcpJetFromHFShower(T const& jet, U const& particles, float dRMax = 0.25, bool searchUpToQuark = false) { typename U::iterator hfparticle; - int origin = jetParticleFromHFShower(jet, particles, hfparticle); + int origin = jetParticleFromHFShower(jet, particles, hfparticle, searchUpToQuark); if (origin == JetTaggingSpecies::charm || origin == JetTaggingSpecies::beauty) { int originalHFMotherIndex = getOriginalHFMotherIndex(hfparticle); diff --git a/PWGJE/TableProducer/jettaggerhf.cxx b/PWGJE/TableProducer/jettaggerhf.cxx index 225f7919db6..24991cc6fd0 100644 --- a/PWGJE/TableProducer/jettaggerhf.cxx +++ b/PWGJE/TableProducer/jettaggerhf.cxx @@ -41,6 +41,7 @@ struct JetTaggerHFTask { Configurable maxDeltaR{"maxDeltaR", 0.25, "maximum distance of jet axis from flavour initiating parton"}; Configurable removeGluonShower{"removeGluonShower", true, "find jet origin removed gluon spliting"}; // true:: remove gluon spliting + Configurable searchUpToQuark{"searchUpToQuark", true, "Finding first mother in particles to quark"}; Configurable useJetProb{"useJetProb", false, "fill table for track counting algorithm"}; Configurable trackProbQA{"trackProbQA", false, "fill track probability histograms separately for geometric positive and negative tracks for QA"}; Configurable doSV{"doSV", false, "fill table for secondary vertex algorithm"}; @@ -173,9 +174,9 @@ struct JetTaggerHFTask { typename JetTagTracksMCD::iterator hftrack; int origin = 0; if (removeGluonShower) - origin = jettaggingutilities::mcdJetFromHFShower(mcdjet, jtracks, particles, maxDeltaR); + origin = jettaggingutilities::mcdJetFromHFShower(mcdjet, jtracks, particles, maxDeltaR, searchUpToQuark); else - origin = jettaggingutilities::jetTrackFromHFShower(mcdjet, jtracks, particles, hftrack); + origin = jettaggingutilities::jetTrackFromHFShower(mcdjet, jtracks, particles, hftrack, searchUpToQuark); int algorithm2 = 0; int algorithm3 = 0; if (useJetProb) { From 8397675d3964010619eed482f98f84ed9a0c89c7 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Fri, 16 Aug 2024 13:18:57 +0200 Subject: [PATCH 16/21] Being developing tagging --- PWGJE/Tasks/jettaggerhfQA.cxx | 197 +++++++++++++++++++++++++++++++++- 1 file changed, 193 insertions(+), 4 deletions(-) diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index c12115ff89e..3fc17cde144 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -36,7 +36,7 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -template +template struct JetTaggerHFQA { // task on/off configuration @@ -84,9 +84,9 @@ struct JetTaggerHFQA { ConfigurableAxis binJetProbabilityLog{"binJetProbabilityLog", {100, 0.f, 10.f}, ""}; ConfigurableAxis binNprongs{"binNprongs", {100, 0., 100.}, ""}; ConfigurableAxis binLxy{"binLxy", {200, 0, 20.f}, ""}; - ConfigurableAxis binSxy{"binSxy", {200, 0, 200.f}, ""}; + ConfigurableAxis binSxy{"binSxy", {1000, 0, 1000.f}, ""}; ConfigurableAxis binLxyz{"binLxyz", {200, 0, 20.f}, ""}; - ConfigurableAxis binSxyz{"binSxyz", {200, 0, 200.f}, ""}; + ConfigurableAxis binSxyz{"binSxyz", {1000, 0, 1000.f}, ""}; ConfigurableAxis binSigmaLxy{"binSigmaLxy", {100, 0., 0.1}, ""}; ConfigurableAxis binSigmaLxyz{"binSigmaLxyz", {100, 0., 0.1}, ""}; @@ -101,6 +101,62 @@ struct JetTaggerHFQA { float maxSigmaLxyz3Prong = 0.; float minSxyz3Prong = 0.; + class bjetCandSV { + public: + bjetCandSV() = default; + + bjetCandSV(float xpv, float ypv, float zpv, float xsv, float ysv, float zsv, + float px, float py, float pz, float e, float m, float chi2PCA, + float errorDecayLength, float errorDecayLengthXY, + float rSecondaryVertex, float pt, float p, + std::array pVector, float eta, float phi, + float y, float decayLength, float decayLengthXY, + float decayLengthNormalised, float decayLengthXYNormalised, + float cpa, float impactParameterXY) + : xPVertex(xpv), yPVertex(ypv), zPVertex(zpv), + xSecondaryVertex(xsv), ySecondaryVertex(ysv), zSecondaryVertex(zsv), + px(px), py(py), pz(pz), e(e), m(m), chi2PCA(chi2PCA), + errorDecayLength(errorDecayLength), errorDecayLengthXY(errorDecayLengthXY), + rSecondaryVertex(rSecondaryVertex), pt(pt), p(p), + pVector(pVector), eta(eta), phi(phi), + y(y), decayLength(decayLength), decayLengthXY(decayLengthXY), + decayLengthNormalised(decayLengthNormalised), decayLengthXYNormalised(decayLengthXYNormalised), + cpa(cpa), impactParameterXY(impactParameterXY) + {} + + ~bjetCandSV() = default; + + float xPVertex = 0.0f; + float yPVertex = 0.0f; + float zPVertex = 0.0f; + float xSecondaryVertex = 0.0f; + float ySecondaryVertex = 0.0f; + float zSecondaryVertex = 0.0f; + float px = 0.0f; + float py = 0.0f; + float pz = 0.0f; + float e = 0.0f; + float m = 0.0f; + float chi2PCA = 0.0f; + float errorDecayLength = 0.0f; + float errorDecayLengthXY = 0.0f; + + float rSecondaryVertex = 0.0f; + float pt = 0.0f; + float p = 0.0f; + std::array pVector = {0.0f, 0.0f, 0.0f}; + float eta = 0.0f; + float phi = 0.0f; + float y = 0.0f; + float decayLength = 0.0f; + float decayLengthXY = 0.0f; + float decayLengthNormalised = 0.0f; + float decayLengthXYNormalised = 0.0f; + float cpa = 0.0f; + float impactParameterXY = 0.0f; + }; + + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; void init(InitContext const&) @@ -232,6 +288,10 @@ struct JetTaggerHFQA { registry.add("h3_sign_impact_parameter_xyz_significance_tc_flavour", "", {HistType::kTH3F, {{impactParameterXYZSignificanceAxis}, {numOrderAxis}, {jetFlavourAxis}}}); } } + if (doprocessIPsMCPMCDMatched) { + registry.add("h3_response_matrix_jet_pt_jet_pt_part_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {jetPtAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_flavour_flavour_run2", "", {HistType::kTH3F, {{jetPtAxis}, {jetFlavourAxis}, {jetFlavourAxis}}}); + } if (doprocessJPData) { registry.add("h2_jet_pt_JP", "jet pt jet probability untagged", {HistType::kTH2F, {{jetPtAxis}, {JetProbabilityAxis}}}); registry.add("h2_jet_pt_neg_log_JP", "jet pt jet probabilityun tagged", {HistType::kTH2F, {{jetPtAxis}, {JetProbabilityLogAxis}}}); @@ -329,6 +389,8 @@ struct JetTaggerHFQA { using JetTagTracksMCD = soa::Join; using OriTracksData = soa::Join; using OriTracksMCD = soa::Join; + using JetTagTableMCDMCPMatched = soa::Join; + using JetTagTableMCPMCDMatched = soa::Join; std::function&, const std::vector&)> sortImp = [](const std::vector& a, const std::vector& b) { @@ -360,6 +422,7 @@ struct JetTaggerHFQA { for (const auto& prong : mcdjet.template secondaryVertices_as()) { float Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); + if (prong.chi2PCA() > 10) continue; if (maxSxy < Sxy) { maxSxy = Sxy; correspondingErrorDecayLengthXY = prong.errorDecayLengthXY(); @@ -368,6 +431,90 @@ struct JetTaggerHFQA { return std::make_tuple(maxSxy, correspondingErrorDecayLengthXY); } + template + bjetCandSV getMaxSxyForJetTemp(const JetType& mcdjet) + { + float xPVertex = 0.0f; + float yPVertex = 0.0f; + float zPVertex = 0.0f; + float xSecondaryVertex = 0.0f; + float ySecondaryVertex = 0.0f; + float zSecondaryVertex = 0.0f; + float px = 0.0f; + float py = 0.0f; + float pz = 0.0f; + float e = 0.0f; + float m = 0.0f; + float chi2PCA = 0.0f; + float errorDecayLength = 0.0f; + float errorDecayLengthXY = 0.0f; + + float rSecondaryVertex = 0.0f; + float pt = 0.0f; + float p = 0.0f; + std::array pVector = {0.0f, 0.0f, 0.0f}; + float eta = 0.0f; + float phi = 0.0f; + float y = 0.0f; + float decayLength = 0.0f; + float decayLengthXY = 0.0f; + float decayLengthNormalised = 0.0f; + float decayLengthXYNormalised = 0.0f; + float cpa = 0.0f; + float impactParameterXY = 0.0f; + + float maxSxy = -1.0f; + + for (const auto& prong : mcdjet.template secondaryVertices_as()) { + float Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); + if (prong.chi2PCA() > 10) continue; + + if (maxSxy < Sxy) { + maxSxy = Sxy; + + xPVertex = prong.xPVertex(); + yPVertex = prong.yPVertex(); + zPVertex = prong.zPVertex(); + xSecondaryVertex = prong.xSecondaryVertex(); + ySecondaryVertex = prong.ySecondaryVertex(); + zSecondaryVertex = prong.zSecondaryVertex(); + px = prong.px(); + py = prong.py(); + pz = prong.pz(); + e = prong.e(); + m = prong.m(); + chi2PCA = prong.chi2PCA(); + errorDecayLength = prong.errorDecayLength(); + errorDecayLengthXY = prong.errorDecayLengthXY(); + rSecondaryVertex = prong.rSecondaryVertex(); + pt = prong.pt(); + p = prong.p(); + pVector = prong.pVector(); + eta = prong.eta(); + phi = prong.phi(); + y = prong.y(); + decayLength = prong.decayLength(); + decayLengthXY = prong.decayLengthXY(); + decayLengthNormalised = prong.decayLengthNormalised(); + decayLengthXYNormalised = prong.decayLengthXYNormalised(); + cpa = prong.cpa(); + impactParameterXY = prong.impactParameterXY(); + } + } + + return bjetCandSV( + xPVertex, yPVertex, zPVertex, + xSecondaryVertex, ySecondaryVertex, zSecondaryVertex, + px, py, pz, e, m, chi2PCA, + errorDecayLength, errorDecayLengthXY, + rSecondaryVertex, pt, p, + pVector, eta, phi, + y, decayLength, decayLengthXY, + decayLengthNormalised, decayLengthXYNormalised, + cpa, impactParameterXY + ); + } + template std::tuple getMaxSxyzForJet(const JetType& mcdjet) { @@ -634,6 +781,27 @@ struct JetTaggerHFQA { } } + Preslice particlesPerCollision = aod::jmcparticle::mcCollisionId; + template + void fillHistogramIPsMCPMCDMatched(T const& collision, U const& mcdjets, V const&, W const&, X const&, Y const& particles) + { + auto const particlesPerColl = particles.sliceBy(particlesPerCollision, collision.mcCollisionId()); + for (auto& mcdjet : mcdjets) { + if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + continue; + } + float eventWeight = mcdjet.eventWeight(); + int jetflavour = mcdjet.origin(); + int jetflavourRun2Def = -1; + //if (!mcdjet.has_matchedJetGeo()) continue; + for (auto& mcpjet : mcdjet.template matchedJetGeo_as()) { + jetflavourRun2Def = jettaggingutilities::getJetFlavor(mcpjet, particlesPerColl); + registry.fill(HIST("h3_response_matrix_jet_pt_jet_pt_part_flavour"), mcdjet.pt(), mcpjet.pt(), jetflavour, eventWeight); + } + registry.fill(HIST("h3_jet_pt_flavour_flavour_run2"), mcdjet.pt(), jetflavour, jetflavourRun2Def, eventWeight); + } + } + template void fillHistogramJPData(T const& /*collision*/, U const& jets) { @@ -803,9 +971,20 @@ struct JetTaggerHFQA { } auto origin = mcdjet.origin(); auto [maxSxy, sigmaLxy] = getMaxSxyForJet(mcdjet); + auto bjetCand = getMaxSxyForJetTemp(mcdjet); + //auto bjetCand = getMaxSxyForJetTemp(mcdjet); + + //auto maxSxy = bjetCand.decayLength(); + //auto sigmaLxy = bjetCand.errorDecayLength(); + auto [maxSxyz, sigmaLxyz] = getMaxSxyzForJet(mcdjet); registry.fill(HIST("h2_3prong_nprongs_flavour"), mcdjet.template secondaryVertices_as().size(), origin); + //if(origin==JetTaggingSpecies::beauty) LOG(info) << "nprongs: " << mcdjet.template secondaryVertices_as().size(); + if (mcdjet.template secondaryVertices_as().size() < 1) continue; for (const auto& prong : mcdjet.template secondaryVertices_as()) { + // + //auto dR = jetutilities::deltaR(mcdjet, prong); + // auto Lxy = prong.decayLengthXY(); auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); auto Lxyz = prong.decayLength(); @@ -824,7 +1003,9 @@ struct JetTaggerHFQA { if (prongAcceptance(prong.errorDecayLength(), Sxyz, maxSigmaLxyz3Prong, minSxyz3Prong)) { registry.fill(HIST("h3_jet_pt_3prong_Sxyz_flavour_cutSxyzAndsigmaLxyz"), mcdjet.pt(), Sxyz, origin); } + //if (origin==JetTaggingSpecies::beauty) LOG (info) << "deltaR: "<< dR << " prongs Sxy: " << Sxy; } + //if (origin==JetTaggingSpecies::beauty) LOG(info) << "flavour: " << origin << " maxSxy: " << maxSxy; registry.fill(HIST("h3_jet_pt_3prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); registry.fill(HIST("h3_jet_pt_3prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin); if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy3Prong, minSxy3Prong)) { @@ -881,6 +1062,12 @@ struct JetTaggerHFQA { } PROCESS_SWITCH(JetTaggerHFQA, processIPsMCD, "Fill impact parameter imformation for mcd jets", false); + void processIPsMCPMCDMatched(soa::Filtered>::iterator const& jcollision, JetTagTableMCDMCPMatched const& mcdjets, JetTagTableMCPMCDMatched const& mcpjets, JetTagTracksMCD const& jtracks, OriTracksMCD const& tracks, JetParticles& particles) + { + fillHistogramIPsMCPMCDMatched(jcollision, mcdjets, mcpjets, jtracks, tracks, particles); + } + PROCESS_SWITCH(JetTaggerHFQA, processIPsMCPMCDMatched, "Fill impact parameter imformation for mcp mcd mathced jets", false); + void processJPData(soa::Filtered::iterator const& jcollision, JetTagTableData const& jets, JetTagTracksData const&) { fillHistogramJPData(jcollision, jets); @@ -920,9 +1107,11 @@ struct JetTaggerHFQA { using JetTaggerQAChargedDataJets = soa::Join; using JetTaggerQAChargedMCDJets = soa::Join; +using JetTaggerQAChargedMCDMCPJets = soa::Join; using JetTaggerQAChargedMCPJets = soa::Join; +using JetTaggerQAChargedMCPMCDJets = soa::Join; -using JetTaggerQACharged = JetTaggerHFQA; +using JetTaggerQACharged = JetTaggerHFQA; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { From 99375b92ad24d02e160a7cc8f3affce31742ec29 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Fri, 16 Aug 2024 13:19:38 +0200 Subject: [PATCH 17/21] Being developing tagging --- PWGJE/Core/JetTaggingUtilities.h | 6 ++--- PWGJE/DataModel/JetTagging.h | 6 ++--- PWGJE/TableProducer/jettaggerhf.cxx | 38 +++++++++++++++++------------ 3 files changed, 29 insertions(+), 21 deletions(-) diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index f06e8144b5b..7388ea4f916 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -304,7 +304,7 @@ int16_t getJetFlavor(AnyJet const& jet, AllMCParticles const& mcparticles) if (dR < jet.r() / 100.f) { if (TMath::Abs(pdgcode) == 5) { - return 2; // Beauty jet + return JetTaggingSpecies::beauty; // Beauty jet } else { if (count > arraySize - 1) return 0; @@ -317,10 +317,10 @@ int16_t getJetFlavor(AnyJet const& jet, AllMCParticles const& mcparticles) for (int ij = 0; ij < count; ij++) { if (TMath::Abs(countpartcode[ij]) == 4) - return 1; // Charm jet + return JetTaggingSpecies::charm; // Charm jet } - return 0; // Light flavor jet + return JetTaggingSpecies::lightflavour; // Light flavor jet } /** diff --git a/PWGJE/DataModel/JetTagging.h b/PWGJE/DataModel/JetTagging.h index c5915396d07..b416d2789d9 100644 --- a/PWGJE/DataModel/JetTagging.h +++ b/PWGJE/DataModel/JetTagging.h @@ -123,10 +123,10 @@ JETSV_TABLES_DEF(Charged, SecondaryVertex2Prong, "2PRONG"); { \ DECLARE_SOA_COLUMN(Origin, origin, int); \ DECLARE_SOA_COLUMN(JetProb, jetProb, std::vector); \ - DECLARE_SOA_COLUMN(Algorithm2, algorithm2, int); \ - DECLARE_SOA_COLUMN(Algorithm3, algorithm3, int); \ + DECLARE_SOA_COLUMN(FlagbJetForIP, flagbjetForIP, int); \ + DECLARE_SOA_COLUMN(FlagbJetForSV, flagbjetForSV, int); \ } \ - DECLARE_SOA_TABLE(_jet_type_##Tags, "AOD", _description_ "Tags", _name_##tagging::Origin, _name_##tagging::JetProb, _name_##tagging::Algorithm2, _name_##tagging::Algorithm3); + DECLARE_SOA_TABLE(_jet_type_##Tags, "AOD", _description_ "Tags", _name_##tagging::Origin, _name_##tagging::JetProb, _name_##tagging::FlagbJetForIP, _name_##tagging::FlagbJetForSV); #define JETTAGGING_TABLES_DEF(_jet_type_, _description_) \ JETTAGGING_TABLE_DEF(_jet_type_##Jet, _jet_type_##jet, _description_) \ diff --git a/PWGJE/TableProducer/jettaggerhf.cxx b/PWGJE/TableProducer/jettaggerhf.cxx index 24991cc6fd0..34fa6ff4651 100644 --- a/PWGJE/TableProducer/jettaggerhf.cxx +++ b/PWGJE/TableProducer/jettaggerhf.cxx @@ -39,12 +39,13 @@ struct JetTaggerHFTask { Produces taggingTableData; Produces taggingTableMCD; + // jet flavour definition Configurable maxDeltaR{"maxDeltaR", 0.25, "maximum distance of jet axis from flavour initiating parton"}; Configurable removeGluonShower{"removeGluonShower", true, "find jet origin removed gluon spliting"}; // true:: remove gluon spliting Configurable searchUpToQuark{"searchUpToQuark", true, "Finding first mother in particles to quark"}; + // configuration about IP method Configurable useJetProb{"useJetProb", false, "fill table for track counting algorithm"}; Configurable trackProbQA{"trackProbQA", false, "fill track probability histograms separately for geometric positive and negative tracks for QA"}; - Configurable doSV{"doSV", false, "fill table for secondary vertex algorithm"}; Configurable numCount{"numCount", 3, "number of track counting"}; Configurable resoFuncMatching{"resoFuncMatching", 0, "matching parameters of resolution function as MC samble (0: custom, 1: custom & inc, 2: MB, 3: MB & inc, 4: JJ, 5: JJ & inc)"}; Configurable> paramsResoFuncData{"paramsResoFuncData", std::vector{1306800, -0.1049, 0.861425, 13.7547, 0.977967, 8.96823, 0.151595, 6.94499, 0.0250301}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7))"}; @@ -53,8 +54,13 @@ struct JetTaggerHFTask { Configurable> paramsResoFuncBeautyJetMC{"paramsResoFuncBeautyJetMC", std::vector{74901.583, -0.082, 0.874, 10.332, 0.941, 7.352, 0.097, 6.220, 0.022}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7)))"}; Configurable> paramsResoFuncLfJetMC{"paramsResoFuncLfJetMC", std::vector{1539435.343, -0.061, 0.896, 13.272, 1.034, 5.884, 0.004, 7.843, 0.090}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7)))"}; Configurable minSignImpXYSig{"minsIPs", -40.0, "minimum of signed impact parameter significance"}; - Configurable tagPoint{"tagPoint", 2.5, "tagging working point"}; + Configurable tagPointForIP{"tagPointForIP", 2.5, "tagging working point for IP"}; + Configurable minIPCount{"minSipCount", 2, "Select at least N signed impact parameter significance in jets"}; // default 2 + // configuration about SV method + Configurable doSV{"doSV", false, "fill table for secondary vertex algorithm"}; + Configurable minDecayLengthCount{"minDecayLengthCount", 1, "Select at least N decay length significance in jets"}; // default 1 + // axis spec ConfigurableAxis binTrackProbability{"binTrackProbability", {100, 0.f, 1.f}, ""}; ConfigurableAxis binJetFlavour{"binJetFlavour", {6, -0.5, 5.5}, ""}; @@ -153,17 +159,17 @@ struct JetTaggerHFTask { void processData(JetCollision const& collision, JetTableData const& jets, JetTagTracksData const& jtracks, OriTracksData const& tracks) { for (auto& jet : jets) { - int algorithm2 = 0; - int algorithm3 = 0; + int flagbJetForIP = 0; + int flagbJetForSV = 0; if (useJetProb) { jetProb.clear(); jetProb.reserve(maxOrder); for (int order = 0; order < maxOrder; order++) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigData, collision, jet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigData, collision, jet, jtracks, tracks, order, tagPointForIP, minSignImpXYSig)); } } - // if (doSV) algorithm2 = jettaggingutilities::Algorithm2((mcdjet, tracks); - taggingTableData(0, jetProb, algorithm2, algorithm3); + // if (doSV) flagbJetForSV = jettaggingutilities::FlagbJetForSV((mcdjet, tracks); + taggingTableData(0, jetProb, flagbJetForIP, flagbJetForSV); } } PROCESS_SWITCH(JetTaggerHFTask, processData, "Fill tagging decision for data jets", false); @@ -171,29 +177,31 @@ struct JetTaggerHFTask { void processMCD(JetCollision const& collision, JetTableMCD const& mcdjets, JetTagTracksMCD const& jtracks, OriTracksMCD const& tracks, JetParticles const& particles) { for (auto& mcdjet : mcdjets) { + int flagbJetForIP = 0; + int flagbJetForSV = 0; typename JetTagTracksMCD::iterator hftrack; int origin = 0; if (removeGluonShower) origin = jettaggingutilities::mcdJetFromHFShower(mcdjet, jtracks, particles, maxDeltaR, searchUpToQuark); else origin = jettaggingutilities::jetTrackFromHFShower(mcdjet, jtracks, particles, hftrack, searchUpToQuark); - int algorithm2 = 0; - int algorithm3 = 0; + //if (origin==JetTaggingSpecies::beauty || jettaggingutilities::isGreaterThanTaggingPoint(collision, mcdjet, jtracks, tracks, tagPointForIP, minIPCount)) + //flagbJetForIP = 1; if (useJetProb) { jetProb.clear(); jetProb.reserve(maxOrder); for (int order = 0; order < maxOrder; order++) { if (useResoFuncFromIncJet) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigIncJetMC, collision, mcdjet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigIncJetMC, collision, mcdjet, jtracks, tracks, order, tagPointForIP, minSignImpXYSig)); } else { if (origin == JetTaggingSpecies::charm) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigCharmJetMC, collision, mcdjet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigCharmJetMC, collision, mcdjet, jtracks, tracks, order, tagPointForIP, minSignImpXYSig)); } if (origin == JetTaggingSpecies::beauty) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigBeautyJetMC, collision, mcdjet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigBeautyJetMC, collision, mcdjet, jtracks, tracks, order, tagPointForIP, minSignImpXYSig)); } if (origin == JetTaggingSpecies::lightflavour) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigLfJetMC, collision, mcdjet, jtracks, tracks, order, tagPoint, minSignImpXYSig)); + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigLfJetMC, collision, mcdjet, jtracks, tracks, order, tagPointForIP, minSignImpXYSig)); } } } @@ -223,8 +231,8 @@ struct JetTaggerHFTask { } } } - // if (doSV) algorithm2 = jettaggingutilities::Algorithm2((mcdjet, tracks); - taggingTableMCD(origin, jetProb, algorithm2, algorithm3); + // if (doSV) flagbJetForSV = jettaggingutilities::FlagbJetForSV((mcdjet, tracks); + taggingTableMCD(origin, jetProb, flagbJetForIP, flagbJetForSV); } } PROCESS_SWITCH(JetTaggerHFTask, processMCD, "Fill tagging decision for mcd jets", false); From 18ee32cef48a4a7e3893f401ca14f8b975bde61e Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Fri, 23 Aug 2024 15:30:52 +0200 Subject: [PATCH 18/21] implement sv tagging jet and cut selection of ip method --- PWGJE/Core/JetTaggingUtilities.h | 216 +++++++++++++++- PWGJE/DataModel/JetTagging.h | 6 +- PWGJE/TableProducer/jettaggerhf.cxx | 178 ++++++++----- PWGJE/Tasks/jettaggerhfQA.cxx | 371 ++++------------------------ 4 files changed, 386 insertions(+), 385 deletions(-) diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index 7388ea4f916..2a8423dcf15 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -323,6 +323,32 @@ int16_t getJetFlavor(AnyJet const& jet, AllMCParticles const& mcparticles) return JetTaggingSpecies::lightflavour; // Light flavor jet } +/** + * return acceptance of track about DCA xy and z due to cut for QualityTracks + */ +template +bool trackAcceptanceWithDca(T const& track, float trackDcaXYMax, float trackDcaZMax) +{ + if (std::abs(track.dcaXY()) > trackDcaXYMax) + return false; + if (std::abs(track.dcaZ()) > trackDcaZMax) + return false; + return true; +} + +/** + * retrun acceptance of prong about chi2 and error of decay length due to cut for high quality secondary vertex + */ +template +bool prongAcceptance(T const& prong, float prongChi2PCAMax, float prongsigmaLxyMax) +{ + if (prong.chi2PCA() > prongChi2PCAMax) + return false; + if (prong.errorDecayLengthXY() > prongsigmaLxyMax) + return false; + return true; +} + /** * return geometric sign which is calculated scalar product between jet axis with DCA (track propagated to PV ) * positive and negative value are expected from primary vertex @@ -350,10 +376,11 @@ int getGeoSign(T const& collision, U const& jet, V const& track) * in a vector in descending order. */ template > -void orderForIPJetTracks(T const& collision, U const& jet, V const& /*jtracks*/, W const& /*tracks*/, Vec& vecSignImpSig) +void orderForIPJetTracks(T const& collision, U const& jet, V const& /*jtracks*/, W const& /*tracks*/, float const &trackDcaXYMax, float const &trackDcaZMax, Vec& vecSignImpSig) { for (auto& jtrack : jet.template tracks_as()) { auto track = jtrack.template track_as(); + if (!trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) continue; auto geoSign = getGeoSign(collision, jet, track); auto varSignImpXYSig = geoSign * TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); vecSignImpSig.push_back(varSignImpXYSig); @@ -364,14 +391,14 @@ void orderForIPJetTracks(T const& collision, U const& jet, V const& /*jtracks*/, /** * Checks if a jet is greater than the given tagging working point based on the signed impact parameter significances */ -template -bool isGreaterThanTaggingPoint(T const& collision, U const& jet, V const& jtracks, W const& tracks, X const& taggingPoint = 1.0, Y const& cnt = 1) +template +bool isGreaterThanTaggingPoint(T const& collision, U const& jet, V const& jtracks, W const& tracks, float const &trackDcaXYMax, float const &trackDcaZMax, float const& taggingPoint = 1.0, int const& cnt = 1) { if (cnt == 0) { return true; // untagged } std::vector vecSignImpSig; - orderForIPJetTracks(collision, jet, jtracks, tracks, vecSignImpSig); + orderForIPJetTracks(collision, jet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, vecSignImpSig); if (vecSignImpSig.size() > static_cast::size_type>(cnt) - 1) { for (int i = 0; i < cnt; i++) { if (vecSignImpSig[i] < taggingPoint) { // tagger point set @@ -448,9 +475,9 @@ float getTrackProbability(T const& fResoFuncjet, U const& track, const float& mi * geometric sign. */ template -float getJetProbability(T const& fResoFuncjet, U const& collision, V const& jet, W const& jtracks, X const& tracks, const int& cnt, const float& tagPoint = 1.0, const float& minSignImpXYSig = -10) +float getJetProbability(T const& fResoFuncjet, U const& collision, V const& jet, W const& jtracks, X const& tracks, float const& trackDcaXYMax, float const& trackDcaZMax, const int& cnt, const float& tagPoint = 1.0, const float& minSignImpXYSig = -10) { - if (!(isGreaterThanTaggingPoint(collision, jet, jtracks, tracks, tagPoint, cnt))) + if (!(isGreaterThanTaggingPoint(collision, jet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, tagPoint, cnt))) return -1; std::vector jetTracksPt; @@ -458,6 +485,7 @@ float getJetProbability(T const& fResoFuncjet, U const& collision, V const& jet, for (auto& jtrack : jet.template tracks_as()) { auto track = jtrack.template track_as(); + if (!trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) continue; float probTrack = getTrackProbability(fResoFuncjet, track, minSignImpXYSig); @@ -481,6 +509,182 @@ float getJetProbability(T const& fResoFuncjet, U const& collision, V const& jet, return JP; } +// For secaondy vertex method utilites + class bjetCandSV { + public: + bjetCandSV() = default; + + bjetCandSV(float xpv, float ypv, float zpv, float xsv, float ysv, float zsv, + float pxVal, float pyVal, float pzVal, float eVal, float mVal, float chi2Val, + float errDecayLength, float errDecayLengthXY, + float rSecVertex, float ptVal, float pVal, + std::array pVec, float etaVal, float phiVal, + float yVal, float decayLen, float decayLenXY, + float decayLenNorm, float decayLenXYNorm, + float cpaVal, float impParXY) + : m_xPVertex(xpv), m_yPVertex(ypv), m_zPVertex(zpv), + m_xSecondaryVertex(xsv), m_ySecondaryVertex(ysv), m_zSecondaryVertex(zsv), + m_px(pxVal), m_py(pyVal), m_pz(pzVal), m_e(eVal), m_m(mVal), m_chi2PCA(chi2Val), + m_errorDecayLength(errDecayLength), m_errorDecayLengthXY(errDecayLengthXY), + m_rSecondaryVertex(rSecVertex), m_pt(ptVal), m_p(pVal), + m_pVector(pVec), m_eta(etaVal), m_phi(phiVal), + m_y(yVal), m_decayLength(decayLen), m_decayLengthXY(decayLenXY), + m_decayLengthNormalised(decayLenNorm), m_decayLengthXYNormalised(decayLenXYNorm), + m_cpa(cpaVal), m_impactParameterXY(impParXY) + {} + + float xPVertex() const { return m_xPVertex; } + float yPVertex() const { return m_yPVertex; } + float zPVertex() const { return m_zPVertex; } + + float xSecondaryVertex() const { return m_xSecondaryVertex; } + float ySecondaryVertex() const { return m_ySecondaryVertex; } + float zSecondaryVertex() const { return m_zSecondaryVertex; } + + float px() const { return m_px; } + float py() const { return m_py; } + float pz() const { return m_pz; } + float e() const { return m_e; } + float m() const { return m_m; } + float chi2PCA() const { return m_chi2PCA; } + + float errorDecayLength() const { return m_errorDecayLength; } + float errorDecayLengthXY() const { return m_errorDecayLengthXY; } + + float rSecondaryVertex() const { return m_rSecondaryVertex; } + float pt() const { return m_pt; } + float p() const { return m_p; } + + std::array pVector() const { return m_pVector; } + + float eta() const { return m_eta; } + float phi() const { return m_phi; } + float y() const { return m_y; } + + float decayLength() const { return m_decayLength; } + float decayLengthXY() const { return m_decayLengthXY; } + float decayLengthNormalised() const { return m_decayLengthNormalised; } + float decayLengthXYNormalised() const { return m_decayLengthXYNormalised; } + + float cpa() const { return m_cpa; } + float impactParameterXY() const { return m_impactParameterXY; } + + private: + float m_xPVertex, m_yPVertex, m_zPVertex; + float m_xSecondaryVertex, m_ySecondaryVertex, m_zSecondaryVertex; + float m_px, m_py, m_pz, m_e, m_m, m_chi2PCA; + float m_errorDecayLength, m_errorDecayLengthXY; + float m_rSecondaryVertex, m_pt, m_p; + std::array m_pVector; + float m_eta, m_phi, m_y; + float m_decayLength, m_decayLengthXY, m_decayLengthNormalised, m_decayLengthXYNormalised; + float m_cpa, m_impactParameterXY; + }; + + template + bjetCandSV jetFromProngMaxDecayLength(const JetType& jet, bool doXYZ=false) + { + float xPVertex = 0.0f; + float yPVertex = 0.0f; + float zPVertex = 0.0f; + float xSecondaryVertex = 0.0f; + float ySecondaryVertex = 0.0f; + float zSecondaryVertex = 0.0f; + float px = 0.0f; + float py = 0.0f; + float pz = 0.0f; + float e = 0.0f; + float m = 0.0f; + float chi2PCA = 0.0f; + float errorDecayLength = 0.0f; + float errorDecayLengthXY = 0.0f; + + float rSecondaryVertex = 0.0f; + float pt = 0.0f; + float p = 0.0f; + std::array pVector = {0.0f, 0.0f, 0.0f}; + float eta = 0.0f; + float phi = 0.0f; + float y = 0.0f; + float decayLength = 0.0f; + float decayLengthXY = 0.0f; + float decayLengthNormalised = 0.0f; + float decayLengthXYNormalised = 0.0f; + float cpa = 0.0f; + float impactParameterXY = 0.0f; + + float maxSxy = -1.0f; + + for (const auto& prong : jet.template secondaryVertices_as()) { + float Sxy = -1.; + if (!doXYZ) { + Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); + } else { + Sxy = prong.decayLength() / prong.errorDecayLength(); + } + //if (!prongAcceptance(prong, doXYZ)) continue; + + if (maxSxy < Sxy) { + maxSxy = Sxy; + + xPVertex = prong.xPVertex(); + yPVertex = prong.yPVertex(); + zPVertex = prong.zPVertex(); + xSecondaryVertex = prong.xSecondaryVertex(); + ySecondaryVertex = prong.ySecondaryVertex(); + zSecondaryVertex = prong.zSecondaryVertex(); + px = prong.px(); + py = prong.py(); + pz = prong.pz(); + e = prong.e(); + m = prong.m(); + chi2PCA = prong.chi2PCA(); + errorDecayLength = prong.errorDecayLength(); + errorDecayLengthXY = prong.errorDecayLengthXY(); + rSecondaryVertex = prong.rSecondaryVertex(); + pt = prong.pt(); + p = prong.p(); + pVector = prong.pVector(); + eta = prong.eta(); + phi = prong.phi(); + y = prong.y(); + decayLength = prong.decayLength(); + decayLengthXY = prong.decayLengthXY(); + decayLengthNormalised = prong.decayLengthNormalised(); + decayLengthXYNormalised = prong.decayLengthXYNormalised(); + cpa = prong.cpa(); + impactParameterXY = prong.impactParameterXY(); + } + } + + return bjetCandSV( + xPVertex, yPVertex, zPVertex, + xSecondaryVertex, ySecondaryVertex, zSecondaryVertex, + px, py, pz, e, m, chi2PCA, + errorDecayLength, errorDecayLengthXY, + rSecondaryVertex, pt, p, + pVector, eta, phi, + y, decayLength, decayLengthXY, + decayLengthNormalised, decayLengthXYNormalised, + cpa, impactParameterXY + ); + } + + template + bool isTaggedJetSV (T const jet, U const& /*prongs*/, float const& doXYZ=false, float const& tagPointForSV = 15.) + { + auto bjetCand = jetFromProngMaxDecayLength(jet); + if (!doXYZ) { + auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); + if (maxSxy < tagPointForSV) return false; + } else { + auto maxSxyz = bjetCand.decayLength() / bjetCand.errorDecayLength(); + if (maxSxyz < tagPointForSV) return false; + } + return true; + } + + }; // namespace jettaggingutilities #endif // PWGJE_CORE_JETTAGGINGUTILITIES_H_ diff --git a/PWGJE/DataModel/JetTagging.h b/PWGJE/DataModel/JetTagging.h index b416d2789d9..860572af981 100644 --- a/PWGJE/DataModel/JetTagging.h +++ b/PWGJE/DataModel/JetTagging.h @@ -123,10 +123,10 @@ JETSV_TABLES_DEF(Charged, SecondaryVertex2Prong, "2PRONG"); { \ DECLARE_SOA_COLUMN(Origin, origin, int); \ DECLARE_SOA_COLUMN(JetProb, jetProb, std::vector); \ - DECLARE_SOA_COLUMN(FlagbJetForIP, flagbjetForIP, int); \ - DECLARE_SOA_COLUMN(FlagbJetForSV, flagbjetForSV, int); \ + DECLARE_SOA_COLUMN(FlagtaggedjetIP, flagtaggedjetIP, bool); \ + DECLARE_SOA_COLUMN(FlagtaggedjetSV, flagtaggedjetSV, bool); \ } \ - DECLARE_SOA_TABLE(_jet_type_##Tags, "AOD", _description_ "Tags", _name_##tagging::Origin, _name_##tagging::JetProb, _name_##tagging::FlagbJetForIP, _name_##tagging::FlagbJetForSV); + DECLARE_SOA_TABLE(_jet_type_##Tags, "AOD", _description_ "Tags", _name_##tagging::Origin, _name_##tagging::JetProb, _name_##tagging::FlagtaggedjetIP, _name_##tagging::FlagtaggedjetSV); #define JETTAGGING_TABLES_DEF(_jet_type_, _description_) \ JETTAGGING_TABLE_DEF(_jet_type_##Jet, _jet_type_##jet, _description_) \ diff --git a/PWGJE/TableProducer/jettaggerhf.cxx b/PWGJE/TableProducer/jettaggerhf.cxx index 34fa6ff4651..7f281dc639c 100644 --- a/PWGJE/TableProducer/jettaggerhf.cxx +++ b/PWGJE/TableProducer/jettaggerhf.cxx @@ -39,10 +39,18 @@ struct JetTaggerHFTask { Produces taggingTableData; Produces taggingTableMCD; + // configuration topological cut for track and sv + Configurable trackDcaXYMax{"trackDcaXYMax", 1, "minimum DCA xy acceptance for tracks [cm]"}; + Configurable trackDcaZMax{"trackDcaZMax", 2, "minimum DCA z acceptance for tracks [cm]"}; + Configurable prongsigmaLxyMax{"prongsigmaLxyMax", 100, "maximum sigma of decay length of prongs on xy plane"}; + Configurable prongsigmaLxyzMax{"prongsigmaLxyzMax", 100, "maximum sigma of decay length of prongs on xyz plane"}; + Configurable prongChi2PCAMax{"prongChi2PCAMax", 10, "maximum Chi2 PCA of decay length of prongs"}; + // jet flavour definition Configurable maxDeltaR{"maxDeltaR", 0.25, "maximum distance of jet axis from flavour initiating parton"}; Configurable removeGluonShower{"removeGluonShower", true, "find jet origin removed gluon spliting"}; // true:: remove gluon spliting Configurable searchUpToQuark{"searchUpToQuark", true, "Finding first mother in particles to quark"}; + // configuration about IP method Configurable useJetProb{"useJetProb", false, "fill table for track counting algorithm"}; Configurable trackProbQA{"trackProbQA", false, "fill track probability histograms separately for geometric positive and negative tracks for QA"}; @@ -55,17 +63,18 @@ struct JetTaggerHFTask { Configurable> paramsResoFuncLfJetMC{"paramsResoFuncLfJetMC", std::vector{1539435.343, -0.061, 0.896, 13.272, 1.034, 5.884, 0.004, 7.843, 0.090}, "parameters of gaus(0)+expo(3)+expo(5)+expo(7)))"}; Configurable minSignImpXYSig{"minsIPs", -40.0, "minimum of signed impact parameter significance"}; Configurable tagPointForIP{"tagPointForIP", 2.5, "tagging working point for IP"}; - Configurable minIPCount{"minSipCount", 2, "Select at least N signed impact parameter significance in jets"}; // default 2 + Configurable minIPCount{"minSipCount", 2, "Select at least N signed impact parameter significance in jets"}; // default 2 // configuration about SV method Configurable doSV{"doSV", false, "fill table for secondary vertex algorithm"}; - Configurable minDecayLengthCount{"minDecayLengthCount", 1, "Select at least N decay length significance in jets"}; // default 1 + Configurable useXYZForTagging{"useXYZForTagging", false, "Enable tagging decision using full XYZ DCA for secondary vertex algorithm"}; + Configurable tagPointForSV{"tagPointForSV", 15, "tagging working point for SV"}; // axis spec ConfigurableAxis binTrackProbability{"binTrackProbability", {100, 0.f, 1.f}, ""}; ConfigurableAxis binJetFlavour{"binJetFlavour", {6, -0.5, 5.5}, ""}; - using JetTagTracksData = soa::Join; - using JetTagTracksMCD = soa::Join; + using JetTagTracksData = soa::Join; + using JetTagTracksMCD = soa::Join; using OriTracksData = soa::Join; using OriTracksMCD = soa::Join; @@ -84,6 +93,58 @@ struct JetTaggerHFTask { std::unique_ptr fSignImpXYSigBeautyJetMC = nullptr; std::unique_ptr fSignImpXYSigLfJetMC = nullptr; + + template + void calculateJetProbabilityMCD(int origin, T const& collision, U const& mcdjet, JetTagTracksMCD const& jtracks, OriTracksMCD const& tracks, std::vector& jetProb) { + jetProb.clear(); + jetProb.reserve(maxOrder); + for (int order = 0; order < maxOrder; order++) { + if (useResoFuncFromIncJet) { + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigIncJetMC, collision, mcdjet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, order, tagPointForIP, minSignImpXYSig)); + } else { + if (origin == JetTaggingSpecies::charm) { + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigCharmJetMC, collision, mcdjet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, order, tagPointForIP, minSignImpXYSig)); + } + if (origin == JetTaggingSpecies::beauty) { + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigBeautyJetMC, collision, mcdjet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, order, tagPointForIP, minSignImpXYSig)); + } + if (origin == JetTaggingSpecies::lightflavour) { + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigLfJetMC, collision, mcdjet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, order, tagPointForIP, minSignImpXYSig)); + } + } + } + } + + template + void evaluateTrackProbQA(int origin, T const& collision, U const& mcdjet, JetTagTracksMCD const& /*jtracks*/, OriTracksMCD const& /*tracks*/) { + for (auto& jtrack : mcdjet.template tracks_as()) { + auto track = jtrack.template track_as(); + if (!track.has_mcParticle()) continue; + if (!jettaggingutilities::trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) continue; + auto geoSign = jettaggingutilities::getGeoSign(collision, mcdjet, track); + float probTrack = -1; + if (useResoFuncFromIncJet) { + probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigIncJetMC, track, minSignImpXYSig); + } else { + if (origin == JetTaggingSpecies::charm) { + probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigCharmJetMC, track, minSignImpXYSig); + } + if (origin == JetTaggingSpecies::beauty) { + probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigBeautyJetMC, track, minSignImpXYSig); + } + if (origin == JetTaggingSpecies::lightflavour) { + probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigLfJetMC, track, minSignImpXYSig); + } + } + if (geoSign > 0) { + registry.fill(HIST("h2_pos_track_probability_flavour"), probTrack, origin); + } else { + registry.fill(HIST("h2_neg_track_probability_flavour"), probTrack, origin); + } + } + } + + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; void init(InitContext const&) { @@ -159,84 +220,91 @@ struct JetTaggerHFTask { void processData(JetCollision const& collision, JetTableData const& jets, JetTagTracksData const& jtracks, OriTracksData const& tracks) { for (auto& jet : jets) { - int flagbJetForIP = 0; - int flagbJetForSV = 0; + bool flagtaggedjetIP = 0; + bool flagtaggedjetSV = 0; if (useJetProb) { jetProb.clear(); jetProb.reserve(maxOrder); for (int order = 0; order < maxOrder; order++) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigData, collision, jet, jtracks, tracks, order, tagPointForIP, minSignImpXYSig)); + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigData, collision, jet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, order, tagPointForIP, minSignImpXYSig)); } } - // if (doSV) flagbJetForSV = jettaggingutilities::FlagbJetForSV((mcdjet, tracks); - taggingTableData(0, jetProb, flagbJetForIP, flagbJetForSV); + if (jettaggingutilities::isGreaterThanTaggingPoint(collision, jet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount)) + flagtaggedjetIP = true; + taggingTableData(0, jetProb, flagtaggedjetIP, flagtaggedjetSV); } } PROCESS_SWITCH(JetTaggerHFTask, processData, "Fill tagging decision for data jets", false); + void processDataWithSV(JetCollision const& collision, soa::Join const& jets, JetTagTracksData const& jtracks, OriTracksData const& tracks, aod::DataSecondaryVertex3Prongs const& prongs) + { + for (auto& jet : jets) { + bool flagtaggedjetIP = 0; + bool flagtaggedjetSV = 0; + if (useJetProb) { + jetProb.clear(); + jetProb.reserve(maxOrder); + for (int order = 0; order < maxOrder; order++) { + jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigData, collision, jet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, order, tagPointForIP, minSignImpXYSig)); + } + } + if (jettaggingutilities::isGreaterThanTaggingPoint(collision, jet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount)) + flagtaggedjetIP = true; + flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(jet, prongs, useXYZForTagging, tagPointForSV); + taggingTableData(0, jetProb, flagtaggedjetIP, flagtaggedjetSV); + } + } + PROCESS_SWITCH(JetTaggerHFTask, processDataWithSV, "Fill tagging decision for data jets", false); + void processMCD(JetCollision const& collision, JetTableMCD const& mcdjets, JetTagTracksMCD const& jtracks, OriTracksMCD const& tracks, JetParticles const& particles) { for (auto& mcdjet : mcdjets) { - int flagbJetForIP = 0; - int flagbJetForSV = 0; + bool flagtaggedjetIP = 0; + bool flagtaggedjetSV = 0; typename JetTagTracksMCD::iterator hftrack; int origin = 0; if (removeGluonShower) origin = jettaggingutilities::mcdJetFromHFShower(mcdjet, jtracks, particles, maxDeltaR, searchUpToQuark); else origin = jettaggingutilities::jetTrackFromHFShower(mcdjet, jtracks, particles, hftrack, searchUpToQuark); - //if (origin==JetTaggingSpecies::beauty || jettaggingutilities::isGreaterThanTaggingPoint(collision, mcdjet, jtracks, tracks, tagPointForIP, minIPCount)) - //flagbJetForIP = 1; if (useJetProb) { - jetProb.clear(); - jetProb.reserve(maxOrder); - for (int order = 0; order < maxOrder; order++) { - if (useResoFuncFromIncJet) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigIncJetMC, collision, mcdjet, jtracks, tracks, order, tagPointForIP, minSignImpXYSig)); - } else { - if (origin == JetTaggingSpecies::charm) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigCharmJetMC, collision, mcdjet, jtracks, tracks, order, tagPointForIP, minSignImpXYSig)); - } - if (origin == JetTaggingSpecies::beauty) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigBeautyJetMC, collision, mcdjet, jtracks, tracks, order, tagPointForIP, minSignImpXYSig)); - } - if (origin == JetTaggingSpecies::lightflavour) { - jetProb.push_back(jettaggingutilities::getJetProbability(fSignImpXYSigLfJetMC, collision, mcdjet, jtracks, tracks, order, tagPointForIP, minSignImpXYSig)); - } - } - } + calculateJetProbabilityMCD(origin, collision, mcdjet, jtracks, tracks, jetProb); if (trackProbQA) { - for (auto& jtrack : mcdjet.template tracks_as()) { - auto track = jtrack.template track_as(); - auto geoSign = jettaggingutilities::getGeoSign(collision, mcdjet, track); - float probTrack = -1; - if (useResoFuncFromIncJet) { - probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigIncJetMC, track, minSignImpXYSig); - } else { - if (origin == JetTaggingSpecies::charm) { - probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigCharmJetMC, track, minSignImpXYSig); - } - if (origin == JetTaggingSpecies::beauty) { - probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigBeautyJetMC, track, minSignImpXYSig); - } - if (origin == JetTaggingSpecies::lightflavour) { - probTrack = jettaggingutilities::getTrackProbability(fSignImpXYSigLfJetMC, track, minSignImpXYSig); - } - } - if (geoSign > 0) { - registry.fill(HIST("h2_pos_track_probability_flavour"), probTrack, origin); - } else { - registry.fill(HIST("h2_neg_track_probability_flavour"), probTrack, origin); - } - } + evaluateTrackProbQA(origin, collision, mcdjet, jtracks, tracks); } } - // if (doSV) flagbJetForSV = jettaggingutilities::FlagbJetForSV((mcdjet, tracks); - taggingTableMCD(origin, jetProb, flagbJetForIP, flagbJetForSV); + if (jettaggingutilities::isGreaterThanTaggingPoint(collision, mcdjet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount)) + flagtaggedjetIP = true; + taggingTableMCD(origin, jetProb, flagtaggedjetIP, flagtaggedjetSV); } } PROCESS_SWITCH(JetTaggerHFTask, processMCD, "Fill tagging decision for mcd jets", false); + void processMCDWithSV(JetCollision const& collision, soa::Join const& mcdjets, JetTagTracksMCD const& jtracks, OriTracksMCD const& tracks, aod::MCDSecondaryVertex3Prongs const& prongs, JetParticles const& particles) + { + for (auto& mcdjet : mcdjets) { + bool flagtaggedjetIP = 0; + bool flagtaggedjetSV = 0; + typename JetTagTracksMCD::iterator hftrack; + int origin = 0; + if (removeGluonShower) + origin = jettaggingutilities::mcdJetFromHFShower(mcdjet, jtracks, particles, maxDeltaR, searchUpToQuark); + else + origin = jettaggingutilities::jetTrackFromHFShower(mcdjet, jtracks, particles, hftrack, searchUpToQuark); + if (useJetProb) { + calculateJetProbabilityMCD(origin, collision, mcdjet, jtracks, tracks, jetProb); + if (trackProbQA) { + evaluateTrackProbQA(origin, collision, mcdjet, jtracks, tracks); + } + } + if (jettaggingutilities::isGreaterThanTaggingPoint(collision, mcdjet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount)) + flagtaggedjetIP = true; + flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(mcdjet, prongs, useXYZForTagging, tagPointForSV); + taggingTableMCD(origin, jetProb, flagtaggedjetIP, flagtaggedjetSV); + } + } + PROCESS_SWITCH(JetTaggerHFTask, processMCDWithSV, "Fill tagging decision for mcd jets with sv", false); + void processTraining(JetCollision const& /*collision*/, JetTableMCD const& /*mcdjets*/, JetTagTracksMCD const& /*tracks*/) { // To create table for ML diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index 3fc17cde144..7d501a73e21 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -51,16 +51,13 @@ struct JetTaggerHFQA { Configurable trackEtaMax{"trackEtaMax", 0.9, "maximum eta acceptance for tracks"}; Configurable trackPtMin{"trackPtMin", 0.15, "minimum pT acceptance for tracks"}; Configurable trackPtMax{"trackPtMax", 100.0, "maximum pT acceptance for tracks"}; + Configurable trackDcaXYMax{"trackDcaXYMax", 1, "minimum DCA xy acceptance for tracks [cm]"}; + Configurable trackDcaZMax{"trackDcaZMax", 2, "minimum DCA z acceptance for tracks [cm]"}; Configurable jetEtaMin{"jetEtaMin", -99.0, "minimum jet pseudorapidity"}; Configurable jetEtaMax{"jetEtaMax", 99.0, "maximum jet pseudorapidity"}; - Configurable prong2sigmaLxyMax{"prong2simgaLxyMax", 0.03, "maximum sigma of decay length of 2-prong on xy plane"}; - Configurable prong2SxyMin{"prong2SxyMin", 7, "minimum decay length significance of 2-prong on xy plane"}; - Configurable prong2sigmaLxyzMax{"prong2sigmaLxyzMax", 0.03, "maximum sigma of decay length of 2-prong on xyz plane"}; - Configurable prong2SxyzMin{"prong2SxyzMin", 7, "minimum decay length significance of 2-prong on xyz plane"}; - Configurable prong3sigmaLxyMax{"prong3sigmaLxyMax", 0.03, "maximum sigma of decay length of 3-prong on xy plane"}; - Configurable prong3SxyMin{"prong3SxyMin", 7, "minimum decay length significance of 3-prong on xy plane"}; - Configurable prong3sigmaLxyzMax{"prong3sigmaLxyzMax", 0.03, "maximum sigma of decay length of 3-prong on xyz plane"}; - Configurable prong3SxyzMin{"prong3SxyzMin", 7, "minimum decay length significance of 3-prong on xyz plane"}; + Configurable prongChi2PCAMax{"prongChi2PCAMax", 10, "maximum Chi2 PCA of decay length of prongs"}; + Configurable prongsigmaLxyMax{"prongsigmaLxyMax", 100, "maximum sigma of decay length of prongs on xy plane"}; + Configurable prongsigmaLxyzMax{"prongsigmaLxyzMax", 100, "maximum sigma of decay length of prongs on xyz plane"}; Configurable numFlavourSpecies{"numFlavourSpecies", 6, "number of jet flavour species"}; Configurable numOrder{"numOrder", 6, "number of ordering"}; @@ -72,13 +69,13 @@ struct JetTaggerHFQA { ConfigurableAxis binPhi{"binPhi", {18 * 8, 0.f, 2. * TMath::Pi()}, ""}; ConfigurableAxis binNtracks{"binNtracks", {100, 0., 100.}, ""}; ConfigurableAxis binTrackPt{"binTrackPt", {200, 0.f, 100.f}, ""}; - ConfigurableAxis binImpactParameterXY{"binImpactParameterXY", {800, -400.5f, 400.5f}, ""}; + ConfigurableAxis binImpactParameterXY{"binImpactParameterXY", {801, -400.5f, 400.5f}, ""}; ConfigurableAxis binSigmaImpactParameterXY{"binImpactSigmaParameterXY", {800, 0.f, 100.f}, ""}; - ConfigurableAxis binImpactParameterXYSignificance{"binImpactParameterXYSignificance", {800, -40.5f, 40.5f}, ""}; - ConfigurableAxis binImpactParameterZ{"binImpactParameterZ", {800, -400.5f, 400.5f}, ""}; - ConfigurableAxis binImpactParameterZSignificance{"binImpactParameterZSignificance", {800, -40.5f, 40.5f}, ""}; - ConfigurableAxis binImpactParameterXYZ{"binImpactParameterXYZ", {2000, -1000.5f, 1000.5f}, ""}; - ConfigurableAxis binImpactParameterXYZSignificance{"binImpactParameterXYZSignificance", {2000, -100.5f, 100.5f}, ""}; + ConfigurableAxis binImpactParameterXYSignificance{"binImpactParameterXYSignificance", {801, -40.5f, 40.5f}, ""}; + ConfigurableAxis binImpactParameterZ{"binImpactParameterZ", {801, -400.5f, 400.5f}, ""}; + ConfigurableAxis binImpactParameterZSignificance{"binImpactParameterZSignificance", {801, -40.5f, 40.5f}, ""}; + ConfigurableAxis binImpactParameterXYZ{"binImpactParameterXYZ", {2001, -1000.5f, 1000.5f}, ""}; + ConfigurableAxis binImpactParameterXYZSignificance{"binImpactParameterXYZSignificance", {2001, -100.5f, 100.5f}, ""}; ConfigurableAxis binNumOrder{"binNumOrder", {6, 0.5, 6.5}, ""}; ConfigurableAxis binJetProbability{"binJetProbability", {100, 0.f, 1.f}, ""}; ConfigurableAxis binJetProbabilityLog{"binJetProbabilityLog", {100, 0.f, 10.f}, ""}; @@ -87,75 +84,12 @@ struct JetTaggerHFQA { ConfigurableAxis binSxy{"binSxy", {1000, 0, 1000.f}, ""}; ConfigurableAxis binLxyz{"binLxyz", {200, 0, 20.f}, ""}; ConfigurableAxis binSxyz{"binSxyz", {1000, 0, 1000.f}, ""}; + ConfigurableAxis binMass{"binMass", {50, 0, 10.f}, ""}; ConfigurableAxis binSigmaLxy{"binSigmaLxy", {100, 0., 0.1}, ""}; ConfigurableAxis binSigmaLxyz{"binSigmaLxyz", {100, 0., 0.1}, ""}; int numberOfJetFlavourSpecies = 6; int trackSelection = -1; - float maxSigmaLxy2Prong = 0.; - float minSxy2Prong = 0.; - float maxSigmaLxyz2Prong = 0.; - float minSxyz2Prong = 0.; - float maxSigmaLxy3Prong = 0.; - float minSxy3Prong = 0.; - float maxSigmaLxyz3Prong = 0.; - float minSxyz3Prong = 0.; - - class bjetCandSV { - public: - bjetCandSV() = default; - - bjetCandSV(float xpv, float ypv, float zpv, float xsv, float ysv, float zsv, - float px, float py, float pz, float e, float m, float chi2PCA, - float errorDecayLength, float errorDecayLengthXY, - float rSecondaryVertex, float pt, float p, - std::array pVector, float eta, float phi, - float y, float decayLength, float decayLengthXY, - float decayLengthNormalised, float decayLengthXYNormalised, - float cpa, float impactParameterXY) - : xPVertex(xpv), yPVertex(ypv), zPVertex(zpv), - xSecondaryVertex(xsv), ySecondaryVertex(ysv), zSecondaryVertex(zsv), - px(px), py(py), pz(pz), e(e), m(m), chi2PCA(chi2PCA), - errorDecayLength(errorDecayLength), errorDecayLengthXY(errorDecayLengthXY), - rSecondaryVertex(rSecondaryVertex), pt(pt), p(p), - pVector(pVector), eta(eta), phi(phi), - y(y), decayLength(decayLength), decayLengthXY(decayLengthXY), - decayLengthNormalised(decayLengthNormalised), decayLengthXYNormalised(decayLengthXYNormalised), - cpa(cpa), impactParameterXY(impactParameterXY) - {} - - ~bjetCandSV() = default; - - float xPVertex = 0.0f; - float yPVertex = 0.0f; - float zPVertex = 0.0f; - float xSecondaryVertex = 0.0f; - float ySecondaryVertex = 0.0f; - float zSecondaryVertex = 0.0f; - float px = 0.0f; - float py = 0.0f; - float pz = 0.0f; - float e = 0.0f; - float m = 0.0f; - float chi2PCA = 0.0f; - float errorDecayLength = 0.0f; - float errorDecayLengthXY = 0.0f; - - float rSecondaryVertex = 0.0f; - float pt = 0.0f; - float p = 0.0f; - std::array pVector = {0.0f, 0.0f, 0.0f}; - float eta = 0.0f; - float phi = 0.0f; - float y = 0.0f; - float decayLength = 0.0f; - float decayLengthXY = 0.0f; - float decayLengthNormalised = 0.0f; - float decayLengthXYNormalised = 0.0f; - float cpa = 0.0f; - float impactParameterXY = 0.0f; - }; - HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -183,18 +117,11 @@ struct JetTaggerHFQA { AxisSpec SxyAxis = {binSxy, "S_{XY}"}; AxisSpec LxyzAxis = {binLxyz, "L_{XYZ} [cm]"}; AxisSpec SxyzAxis = {binSxyz, "S_{XYZ}"}; + AxisSpec massAxis = {binMass, "#it{m}_{SV}"}; AxisSpec sigmaLxyAxis = {binSigmaLxy, "#sigma_{L_{XY}} [cm]"}; AxisSpec sigmaLxyzAxis = {binSigmaLxyz, "#sigma_{L_{XYZ}} [cm]"}; numberOfJetFlavourSpecies = static_cast(numFlavourSpecies); - maxSigmaLxy2Prong = static_cast(prong2sigmaLxyMax); - minSxy2Prong = static_cast(prong2SxyMin); - maxSigmaLxyz2Prong = static_cast(prong2sigmaLxyzMax); - minSxyz2Prong = static_cast(prong2SxyzMin); - maxSigmaLxy3Prong = static_cast(prong3sigmaLxyMax); - minSxy3Prong = static_cast(prong3SxyMin); - maxSigmaLxyz3Prong = static_cast(prong3sigmaLxyzMax); - minSxyz3Prong = static_cast(prong3SxyzMin); trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); if (doprocessTracksDca) { @@ -315,70 +242,50 @@ struct JetTaggerHFQA { if (doprocessSV2ProngData) { registry.add("h_2prong_nprongs", "", {HistType::kTH1F, {{nprongsAxis}}}); registry.add("h2_jet_pt_2prong_Lxy", "", {HistType::kTH2F, {{jetPtAxis}, {LxyAxis}}}); + registry.add("h2_jet_pt_2prong_sigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyAxis}}}); registry.add("h2_jet_pt_2prong_Sxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); registry.add("h2_jet_pt_2prong_Lxyz", "", {HistType::kTH2F, {{jetPtAxis}, {LxyzAxis}}}); + registry.add("h2_jet_pt_2prong_sigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyzAxis}}}); registry.add("h2_jet_pt_2prong_Sxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); registry.add("h2_jet_pt_2prong_Sxy_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); registry.add("h2_jet_pt_2prong_Sxyz_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); - registry.add("h2_2prong_Sxy_sigmaLxy", "", {HistType::kTH2F, {{SxyAxis}, {sigmaLxyAxis}}}); - registry.add("h2_2prong_Sxyz_sigmaLxyz", "", {HistType::kTH2F, {{SxyzAxis}, {sigmaLxyzAxis}}}); - registry.add("h2_jet_pt_2prong_sigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyAxis}}}); - registry.add("h2_jet_pt_2prong_sigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyzAxis}}}); - registry.add("h2_jet_pt_2prong_Sxy_cutSxyAndsigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); - registry.add("h2_jet_pt_2prong_Sxyz_cutSxyzAndsigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); - registry.add("h2_jet_pt_2prong_Sxy_N1_cutSxyAndsigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); - registry.add("h2_jet_pt_2prong_Sxyz_N1_cutSxyzAndsigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); + registry.add("h2_jet_pt_2prong_mass_N1", "", {HistType::kTH2F, {{jetPtAxis}, {massAxis}}}); } if (doprocessSV3ProngData) { registry.add("h_3prong_nprongs", "", {HistType::kTH1F, {{nprongsAxis}}}); registry.add("h2_jet_pt_3prong_Lxy", "", {HistType::kTH2F, {{jetPtAxis}, {LxyAxis}}}); + registry.add("h2_jet_pt_3prong_sigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyAxis}}}); registry.add("h2_jet_pt_3prong_Sxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); registry.add("h2_jet_pt_3prong_Lxyz", "", {HistType::kTH2F, {{jetPtAxis}, {LxyzAxis}}}); + registry.add("h2_jet_pt_3prong_sigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyzAxis}}}); registry.add("h2_jet_pt_3prong_Sxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); registry.add("h2_jet_pt_3prong_Sxy_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); registry.add("h2_jet_pt_3prong_Sxyz_N1", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); - registry.add("h2_3prong_Sxy_sigmaLxy", "", {HistType::kTH2F, {{SxyAxis}, {sigmaLxyAxis}}}); - registry.add("h2_3prong_Sxyz_sigmaLxyz", "", {HistType::kTH2F, {{SxyzAxis}, {sigmaLxyzAxis}}}); - registry.add("h2_jet_pt_3prong_sigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyAxis}}}); - registry.add("h2_jet_pt_3prong_sigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {sigmaLxyzAxis}}}); - registry.add("h2_jet_pt_3prong_Sxy_cutSxyAndsigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); - registry.add("h2_jet_pt_3prong_Sxyz_cutSxyzAndsigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); - registry.add("h2_jet_pt_3prong_Sxy_N1_cutSxyAndsigmaLxy", "", {HistType::kTH2F, {{jetPtAxis}, {SxyAxis}}}); - registry.add("h2_jet_pt_3prong_Sxyz_N1_cutSxyzAndsigmaLxyz", "", {HistType::kTH2F, {{jetPtAxis}, {SxyzAxis}}}); + registry.add("h2_jet_pt_3prong_mass_N1", "", {HistType::kTH2F, {{jetPtAxis}, {massAxis}}}); } if (doprocessSV2ProngMCD) { registry.add("h2_2prong_nprongs_flavour", "", {HistType::kTH2F, {{nprongsAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Lxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {LxyAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_2prong_sigmaLxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Sxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Lxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {LxyzAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_2prong_sigmaLxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Sxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Sxy_N1_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_2prong_Sxyz_N1_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); - registry.add("h3_2prong_Sxy_sigmaLxy_flavour", "", {HistType::kTH3F, {{SxyAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); - registry.add("h3_2prong_Sxyz_sigmaLxyz_flavour", "", {HistType::kTH3F, {{SxyzAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_2prong_sigmaLxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_2prong_sigmaLxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_2prong_Sxy_flavour_cutSxyAndsigmaLxy", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_2prong_Sxyz_flavour_cutSxyzAndsigmaLxyz", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_2prong_Sxy_N1_flavour_cutSxyAndsigmaLxy", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_2prong_Sxyz_N1_flavour_cutSxyzAndsigmaLxyz", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_2prong_mass_N1_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {massAxis}, {jetFlavourAxis}}}); } if (doprocessSV3ProngMCD) { registry.add("h2_3prong_nprongs_flavour", "", {HistType::kTH2F, {{nprongsAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_Lxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {LxyAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_3prong_sigmaLxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_Sxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_Lxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {LxyzAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_3prong_sigmaLxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_Sxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_Sxy_N1_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); registry.add("h3_jet_pt_3prong_Sxyz_N1_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); - registry.add("h3_3prong_Sxy_sigmaLxy_flavour", "", {HistType::kTH3F, {{SxyAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); - registry.add("h3_3prong_Sxyz_sigmaLxyz_flavour", "", {HistType::kTH3F, {{SxyzAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_3prong_sigmaLxy_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_3prong_sigmaLxyz_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {sigmaLxyzAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_3prong_Sxy_flavour_cutSxyAndsigmaLxy", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_3prong_Sxyz_flavour_cutSxyzAndsigmaLxyz", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_3prong_Sxy_N1_flavour_cutSxyAndsigmaLxy", "", {HistType::kTH3F, {{jetPtAxis}, {SxyAxis}, {jetFlavourAxis}}}); - registry.add("h3_jet_pt_3prong_Sxyz_N1_flavour_cutSxyzAndsigmaLxyz", "", {HistType::kTH3F, {{jetPtAxis}, {SxyzAxis}, {jetFlavourAxis}}}); + registry.add("h3_jet_pt_3prong_mass_N1_flavour", "", {HistType::kTH3F, {{jetPtAxis}, {massAxis}, {jetFlavourAxis}}}); } } @@ -406,131 +313,6 @@ struct JetTaggerHFQA { return true; } - bool prongAcceptance(float sigmaDecayLength, float decayLengthSig, float maxSigmaDecayLength, float minDecayLengthSig) - { - if ((sigmaDecayLength < maxSigmaDecayLength) && (decayLengthSig > minDecayLengthSig)) - return true; - - return false; - } - - template - std::tuple getMaxSxyForJet(const JetType& mcdjet) - { - float maxSxy = 0; - float correspondingErrorDecayLengthXY = 0; - - for (const auto& prong : mcdjet.template secondaryVertices_as()) { - float Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); - if (prong.chi2PCA() > 10) continue; - if (maxSxy < Sxy) { - maxSxy = Sxy; - correspondingErrorDecayLengthXY = prong.errorDecayLengthXY(); - } - } - return std::make_tuple(maxSxy, correspondingErrorDecayLengthXY); - } - - template - bjetCandSV getMaxSxyForJetTemp(const JetType& mcdjet) - { - float xPVertex = 0.0f; - float yPVertex = 0.0f; - float zPVertex = 0.0f; - float xSecondaryVertex = 0.0f; - float ySecondaryVertex = 0.0f; - float zSecondaryVertex = 0.0f; - float px = 0.0f; - float py = 0.0f; - float pz = 0.0f; - float e = 0.0f; - float m = 0.0f; - float chi2PCA = 0.0f; - float errorDecayLength = 0.0f; - float errorDecayLengthXY = 0.0f; - - float rSecondaryVertex = 0.0f; - float pt = 0.0f; - float p = 0.0f; - std::array pVector = {0.0f, 0.0f, 0.0f}; - float eta = 0.0f; - float phi = 0.0f; - float y = 0.0f; - float decayLength = 0.0f; - float decayLengthXY = 0.0f; - float decayLengthNormalised = 0.0f; - float decayLengthXYNormalised = 0.0f; - float cpa = 0.0f; - float impactParameterXY = 0.0f; - - float maxSxy = -1.0f; - - for (const auto& prong : mcdjet.template secondaryVertices_as()) { - float Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); - if (prong.chi2PCA() > 10) continue; - - if (maxSxy < Sxy) { - maxSxy = Sxy; - - xPVertex = prong.xPVertex(); - yPVertex = prong.yPVertex(); - zPVertex = prong.zPVertex(); - xSecondaryVertex = prong.xSecondaryVertex(); - ySecondaryVertex = prong.ySecondaryVertex(); - zSecondaryVertex = prong.zSecondaryVertex(); - px = prong.px(); - py = prong.py(); - pz = prong.pz(); - e = prong.e(); - m = prong.m(); - chi2PCA = prong.chi2PCA(); - errorDecayLength = prong.errorDecayLength(); - errorDecayLengthXY = prong.errorDecayLengthXY(); - rSecondaryVertex = prong.rSecondaryVertex(); - pt = prong.pt(); - p = prong.p(); - pVector = prong.pVector(); - eta = prong.eta(); - phi = prong.phi(); - y = prong.y(); - decayLength = prong.decayLength(); - decayLengthXY = prong.decayLengthXY(); - decayLengthNormalised = prong.decayLengthNormalised(); - decayLengthXYNormalised = prong.decayLengthXYNormalised(); - cpa = prong.cpa(); - impactParameterXY = prong.impactParameterXY(); - } - } - - return bjetCandSV( - xPVertex, yPVertex, zPVertex, - xSecondaryVertex, ySecondaryVertex, zSecondaryVertex, - px, py, pz, e, m, chi2PCA, - errorDecayLength, errorDecayLengthXY, - rSecondaryVertex, pt, p, - pVector, eta, phi, - y, decayLength, decayLengthXY, - decayLengthNormalised, decayLengthXYNormalised, - cpa, impactParameterXY - ); - } - - template - std::tuple getMaxSxyzForJet(const JetType& mcdjet) - { - float maxSxyz = 0; - float correspondingErrorDecayLength = 0; - - for (const auto& prong : mcdjet.template secondaryVertices_as()) { - float Sxyz = prong.decayLength() / prong.errorDecayLength(); - if (maxSxyz < Sxyz) { - maxSxyz = Sxyz; - correspondingErrorDecayLength = prong.errorDecayLength(); - } - } - return std::make_tuple(maxSxyz, correspondingErrorDecayLength); - } - template void fillHistogramIPsData(T const& collision, U const& jets, V const& /*jtracks*/, W const& /*tracks*/) { @@ -543,7 +325,7 @@ struct JetTaggerHFQA { auto track = jtrack.template track_as(); if (!trackAcceptance(track)) continue; - + if (!jettaggingutilities::trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) continue; // General parameters registry.fill(HIST("h3_jet_pt_track_pt_track_eta"), jet.pt(), track.pt(), track.eta()); registry.fill(HIST("h3_jet_pt_track_pt_track_phi"), jet.pt(), track.pt(), track.phi()); @@ -641,6 +423,8 @@ struct JetTaggerHFQA { auto track = jtrack.template track_as(); if (!trackAcceptance(track)) continue; + if (!jettaggingutilities::trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) continue; + // General parameters registry.fill(HIST("h3_jet_pt_track_pt_flavour"), mcdjet.pt(), track.pt(), jetflavour); registry.fill(HIST("h3_jet_pt_track_eta_flavour"), mcdjet.pt(), track.eta(), jetflavour); @@ -651,7 +435,7 @@ struct JetTaggerHFQA { varImpXY = track.dcaXY() * jettaggingutilities::cmTomum; float varSigmaImpXY = track.dcaXY() * jettaggingutilities::cmTomum; varSignImpXY = geoSign * std::abs(track.dcaXY()) * jettaggingutilities::cmTomum; - varImpXYSig = std::abs(track.dcaXY()) / std::sqrt(track.sigmaDcaXY2()); + varImpXYSig = track.dcaXY() / std::sqrt(track.sigmaDcaXY2()); varSignImpXYSig = geoSign * std::abs(track.dcaXY()) / std::sqrt(track.sigmaDcaXY2()); registry.fill(HIST("h3_jet_pt_impact_parameter_xy_flavour"), mcdjet.pt(), varImpXY, jetflavour); registry.fill(HIST("h3_jet_pt_sigma_impact_parameter_xy_flavour"), mcdjet.pt(), varSigmaImpXY, jetflavour); @@ -845,8 +629,6 @@ struct JetTaggerHFQA { if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { continue; } - auto [maxSxy, sigmaLxy] = getMaxSxyForJet(jet); - auto [maxSxyz, sigmaLxyz] = getMaxSxyzForJet(jet); registry.fill(HIST("h_2prong_nprongs"), jet.template secondaryVertices_as().size()); for (const auto& prong : jet.template secondaryVertices_as()) { auto Lxy = prong.decayLengthXY(); @@ -857,25 +639,15 @@ struct JetTaggerHFQA { registry.fill(HIST("h2_jet_pt_2prong_Sxy"), jet.pt(), Sxy); registry.fill(HIST("h2_jet_pt_2prong_Lxyz"), jet.pt(), Lxyz); registry.fill(HIST("h2_jet_pt_2prong_Sxyz"), jet.pt(), Sxyz); - registry.fill(HIST("h2_2prong_Sxy_sigmaLxy"), Sxy, prong.errorDecayLengthXY()); - registry.fill(HIST("h2_2prong_Sxyz_sigmaLxyz"), Sxyz, prong.errorDecayLength()); registry.fill(HIST("h2_jet_pt_2prong_sigmaLxy"), jet.pt(), prong.errorDecayLengthXY()); registry.fill(HIST("h2_jet_pt_2prong_sigmaLxyz"), jet.pt(), prong.errorDecayLength()); - if (prongAcceptance(prong.errorDecayLengthXY(), Sxy, maxSigmaLxy2Prong, minSxy2Prong)) { - registry.fill(HIST("h2_jet_pt_2prong_Sxy_cutSxyAndsigmaLxy"), jet.pt(), Sxy); - } - if (prongAcceptance(prong.errorDecayLength(), Sxyz, maxSigmaLxyz2Prong, minSxyz2Prong)) { - registry.fill(HIST("h2_jet_pt_2prong_Sxyz_cutSxyzAndsigmaLxyz"), jet.pt(), Sxyz); - } } + auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(jet); + auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); + auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(jet, true); + auto maxSxyz = bjetCandForXYZ.decayLength() / bjetCandForXYZ.errorDecayLength(); registry.fill(HIST("h2_jet_pt_2prong_Sxy_N1"), jet.pt(), maxSxy); registry.fill(HIST("h2_jet_pt_2prong_Sxyz_N1"), jet.pt(), maxSxyz); - if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy2Prong, minSxy2Prong)) { - registry.fill(HIST("h2_jet_pt_2prong_Sxy_N1_cutSxyAndsigmaLxy"), jet.pt(), maxSxy); - } - if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy2Prong, minSxy2Prong)) { - registry.fill(HIST("h2_jet_pt_2prong_Sxyz_N1_cutSxyzAndsigmaLxyz"), jet.pt(), maxSxyz); - } } } @@ -886,8 +658,6 @@ struct JetTaggerHFQA { if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { continue; } - auto [maxSxy, sigmaLxy] = getMaxSxyForJet(jet); - auto [maxSxyz, sigmaLxyz] = getMaxSxyzForJet(jet); registry.fill(HIST("h_3prong_nprongs"), jet.template secondaryVertices_as().size()); for (const auto& prong : jet.template secondaryVertices_as()) { auto Lxy = prong.decayLengthXY(); @@ -898,25 +668,15 @@ struct JetTaggerHFQA { registry.fill(HIST("h2_jet_pt_3prong_Sxy"), jet.pt(), Sxy); registry.fill(HIST("h2_jet_pt_3prong_Lxyz"), jet.pt(), Lxyz); registry.fill(HIST("h2_jet_pt_3prong_Sxyz"), jet.pt(), Sxyz); - registry.fill(HIST("h2_3prong_Sxy_sigmaLxy"), Sxy, prong.errorDecayLengthXY()); - registry.fill(HIST("h2_3prong_Sxyz_sigmaLxyz"), Sxyz, prong.errorDecayLength()); registry.fill(HIST("h2_jet_pt_3prong_sigmaLxy"), jet.pt(), prong.errorDecayLengthXY()); registry.fill(HIST("h2_jet_pt_3prong_sigmaLxyz"), jet.pt(), prong.errorDecayLength()); - if (prongAcceptance(prong.errorDecayLengthXY(), Sxy, maxSigmaLxy3Prong, minSxy3Prong)) { - registry.fill(HIST("h2_jet_pt_3prong_Sxy_cutSxyAndsigmaLxy"), jet.pt(), Sxy); - } - if (prongAcceptance(prong.errorDecayLength(), Sxyz, maxSigmaLxyz3Prong, minSxyz3Prong)) { - registry.fill(HIST("h2_jet_pt_3prong_Sxyz_cutSxyzAndsigmaLxyz"), jet.pt(), Sxyz); - } } + auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(jet); + auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); + auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(jet, true); + auto maxSxyz = bjetCandForXYZ.decayLength() / bjetCandForXYZ.errorDecayLength(); registry.fill(HIST("h2_jet_pt_3prong_Sxy_N1"), jet.pt(), maxSxy); registry.fill(HIST("h2_jet_pt_3prong_Sxyz_N1"), jet.pt(), maxSxyz); - if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy3Prong, minSxy3Prong)) { - registry.fill(HIST("h2_jet_pt_3prong_Sxy_N1_cutSxyAndsigmaLxy"), jet.pt(), maxSxy); - } - if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy3Prong, minSxy3Prong)) { - registry.fill(HIST("h2_jet_pt_3prong_Sxyz_N1_cutSxyzAndsigmaLxyz"), jet.pt(), maxSxyz); - } } } @@ -928,9 +688,8 @@ struct JetTaggerHFQA { continue; } auto origin = mcdjet.origin(); - auto [maxSxy, sigmaLxy] = getMaxSxyForJet(mcdjet); - auto [maxSxyz, sigmaLxyz] = getMaxSxyzForJet(mcdjet); registry.fill(HIST("h2_2prong_nprongs_flavour"), mcdjet.template secondaryVertices_as().size(), origin); + if (mcdjet.template secondaryVertices_as().size() < 1) continue; for (const auto& prong : mcdjet.template secondaryVertices_as()) { auto Lxy = prong.decayLengthXY(); auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); @@ -940,25 +699,17 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_jet_pt_2prong_Sxy_flavour"), mcdjet.pt(), Sxy, origin); registry.fill(HIST("h3_jet_pt_2prong_Lxyz_flavour"), mcdjet.pt(), Lxyz, origin); registry.fill(HIST("h3_jet_pt_2prong_Sxyz_flavour"), mcdjet.pt(), Sxyz, origin); - registry.fill(HIST("h3_2prong_Sxy_sigmaLxy_flavour"), Sxy, prong.errorDecayLengthXY(), origin); - registry.fill(HIST("h3_2prong_Sxyz_sigmaLxyz_flavour"), Sxyz, prong.errorDecayLength(), origin); registry.fill(HIST("h3_jet_pt_2prong_sigmaLxy_flavour"), mcdjet.pt(), prong.errorDecayLengthXY(), origin); registry.fill(HIST("h3_jet_pt_2prong_sigmaLxyz_flavour"), mcdjet.pt(), prong.errorDecayLength(), origin); - if (prongAcceptance(prong.errorDecayLengthXY(), Sxy, maxSigmaLxy2Prong, minSxy2Prong)) { - registry.fill(HIST("h3_jet_pt_2prong_Sxy_flavour_cutSxyAndsigmaLxy"), mcdjet.pt(), Sxy, origin); - } - if (prongAcceptance(prong.errorDecayLength(), Sxyz, maxSigmaLxyz2Prong, minSxyz2Prong)) { - registry.fill(HIST("h3_jet_pt_2prong_Sxyz_flavour_cutSxyzAndsigmaLxyz"), mcdjet.pt(), Sxyz, origin); - } } + auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet); + auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); + auto massSV = bjetCand.m(); + auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, true); + auto maxSxyz = bjetCandForXYZ.decayLength() / bjetCandForXYZ.errorDecayLength(); registry.fill(HIST("h3_jet_pt_2prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); registry.fill(HIST("h3_jet_pt_2prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin); - if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy2Prong, minSxy2Prong)) { - registry.fill(HIST("h3_jet_pt_2prong_Sxy_N1_flavour_cutSxyAndsigmaLxy"), mcdjet.pt(), maxSxy, origin); - } - if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy2Prong, minSxy2Prong)) { - registry.fill(HIST("h3_jet_pt_2prong_Sxyz_N1_flavour_cutSxyzAndsigmaLxyz"), mcdjet.pt(), maxSxyz, origin); - } + registry.fill(HIST("h3_jet_pt_2prong_mass_N1_flavour"), mcdjet.pt(), massSV, origin); } } @@ -970,21 +721,9 @@ struct JetTaggerHFQA { continue; } auto origin = mcdjet.origin(); - auto [maxSxy, sigmaLxy] = getMaxSxyForJet(mcdjet); - auto bjetCand = getMaxSxyForJetTemp(mcdjet); - //auto bjetCand = getMaxSxyForJetTemp(mcdjet); - - //auto maxSxy = bjetCand.decayLength(); - //auto sigmaLxy = bjetCand.errorDecayLength(); - - auto [maxSxyz, sigmaLxyz] = getMaxSxyzForJet(mcdjet); registry.fill(HIST("h2_3prong_nprongs_flavour"), mcdjet.template secondaryVertices_as().size(), origin); - //if(origin==JetTaggingSpecies::beauty) LOG(info) << "nprongs: " << mcdjet.template secondaryVertices_as().size(); if (mcdjet.template secondaryVertices_as().size() < 1) continue; for (const auto& prong : mcdjet.template secondaryVertices_as()) { - // - //auto dR = jetutilities::deltaR(mcdjet, prong); - // auto Lxy = prong.decayLengthXY(); auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); auto Lxyz = prong.decayLength(); @@ -993,27 +732,17 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_jet_pt_3prong_Sxy_flavour"), mcdjet.pt(), Sxy, origin); registry.fill(HIST("h3_jet_pt_3prong_Lxyz_flavour"), mcdjet.pt(), Lxyz, origin); registry.fill(HIST("h3_jet_pt_3prong_Sxyz_flavour"), mcdjet.pt(), Sxyz, origin); - registry.fill(HIST("h3_3prong_Sxy_sigmaLxy_flavour"), Sxy, prong.errorDecayLengthXY(), origin); - registry.fill(HIST("h3_3prong_Sxyz_sigmaLxyz_flavour"), Sxyz, prong.errorDecayLength(), origin); registry.fill(HIST("h3_jet_pt_3prong_sigmaLxy_flavour"), mcdjet.pt(), prong.errorDecayLengthXY(), origin); registry.fill(HIST("h3_jet_pt_3prong_sigmaLxyz_flavour"), mcdjet.pt(), prong.errorDecayLength(), origin); - if (prongAcceptance(prong.errorDecayLengthXY(), Sxy, maxSigmaLxy3Prong, minSxy3Prong)) { - registry.fill(HIST("h3_jet_pt_3prong_Sxy_flavour_cutSxyAndsigmaLxy"), mcdjet.pt(), Sxy, origin); - } - if (prongAcceptance(prong.errorDecayLength(), Sxyz, maxSigmaLxyz3Prong, minSxyz3Prong)) { - registry.fill(HIST("h3_jet_pt_3prong_Sxyz_flavour_cutSxyzAndsigmaLxyz"), mcdjet.pt(), Sxyz, origin); - } - //if (origin==JetTaggingSpecies::beauty) LOG (info) << "deltaR: "<< dR << " prongs Sxy: " << Sxy; } - //if (origin==JetTaggingSpecies::beauty) LOG(info) << "flavour: " << origin << " maxSxy: " << maxSxy; + auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet); + auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); + auto massSV = bjetCand.m(); + auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, true); + auto maxSxyz = bjetCandForXYZ.decayLength() / bjetCandForXYZ.errorDecayLength(); registry.fill(HIST("h3_jet_pt_3prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); registry.fill(HIST("h3_jet_pt_3prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin); - if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy3Prong, minSxy3Prong)) { - registry.fill(HIST("h3_jet_pt_3prong_Sxy_N1_flavour_cutSxyAndsigmaLxy"), mcdjet.pt(), maxSxy, origin); - } - if (prongAcceptance(sigmaLxy, maxSxy, maxSigmaLxy3Prong, minSxy3Prong)) { - registry.fill(HIST("h3_jet_pt_3prong_Sxyz_N1_flavour_cutSxyzAndsigmaLxyz"), mcdjet.pt(), maxSxyz, origin); - } + registry.fill(HIST("h3_jet_pt_3prong_mass_N1_flavour"), mcdjet.pt(), massSV, origin); } } From ef4c49024dc5129d37a9274f071c5980935fe1fb Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Fri, 23 Aug 2024 15:33:16 +0200 Subject: [PATCH 19/21] fix clang --- PWGJE/Core/JetTaggingUtilities.h | 340 ++++++++++++++-------------- PWGJE/TableProducer/jettaggerhf.cxx | 14 +- PWGJE/Tasks/jettaggerhfQA.cxx | 14 +- 3 files changed, 185 insertions(+), 183 deletions(-) diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index 2a8423dcf15..56c9fade744 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -342,9 +342,9 @@ bool trackAcceptanceWithDca(T const& track, float trackDcaXYMax, float trackDcaZ template bool prongAcceptance(T const& prong, float prongChi2PCAMax, float prongsigmaLxyMax) { - if (prong.chi2PCA() > prongChi2PCAMax) + if (prong.chi2PCA() > prongChi2PCAMax) return false; - if (prong.errorDecayLengthXY() > prongsigmaLxyMax) + if (prong.errorDecayLengthXY() > prongsigmaLxyMax) return false; return true; } @@ -376,11 +376,12 @@ int getGeoSign(T const& collision, U const& jet, V const& track) * in a vector in descending order. */ template > -void orderForIPJetTracks(T const& collision, U const& jet, V const& /*jtracks*/, W const& /*tracks*/, float const &trackDcaXYMax, float const &trackDcaZMax, Vec& vecSignImpSig) +void orderForIPJetTracks(T const& collision, U const& jet, V const& /*jtracks*/, W const& /*tracks*/, float const& trackDcaXYMax, float const& trackDcaZMax, Vec& vecSignImpSig) { for (auto& jtrack : jet.template tracks_as()) { auto track = jtrack.template track_as(); - if (!trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) continue; + if (!trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) + continue; auto geoSign = getGeoSign(collision, jet, track); auto varSignImpXYSig = geoSign * TMath::Abs(track.dcaXY()) / TMath::Sqrt(track.sigmaDcaXY2()); vecSignImpSig.push_back(varSignImpXYSig); @@ -392,7 +393,7 @@ void orderForIPJetTracks(T const& collision, U const& jet, V const& /*jtracks*/, * Checks if a jet is greater than the given tagging working point based on the signed impact parameter significances */ template -bool isGreaterThanTaggingPoint(T const& collision, U const& jet, V const& jtracks, W const& tracks, float const &trackDcaXYMax, float const &trackDcaZMax, float const& taggingPoint = 1.0, int const& cnt = 1) +bool isGreaterThanTaggingPoint(T const& collision, U const& jet, V const& jtracks, W const& tracks, float const& trackDcaXYMax, float const& trackDcaZMax, float const& taggingPoint = 1.0, int const& cnt = 1) { if (cnt == 0) { return true; // untagged @@ -485,7 +486,8 @@ float getJetProbability(T const& fResoFuncjet, U const& collision, V const& jet, for (auto& jtrack : jet.template tracks_as()) { auto track = jtrack.template track_as(); - if (!trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) continue; + if (!trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) + continue; float probTrack = getTrackProbability(fResoFuncjet, track, minSignImpXYSig); @@ -510,180 +512,174 @@ float getJetProbability(T const& fResoFuncjet, U const& collision, V const& jet, } // For secaondy vertex method utilites - class bjetCandSV { - public: - bjetCandSV() = default; - - bjetCandSV(float xpv, float ypv, float zpv, float xsv, float ysv, float zsv, - float pxVal, float pyVal, float pzVal, float eVal, float mVal, float chi2Val, - float errDecayLength, float errDecayLengthXY, - float rSecVertex, float ptVal, float pVal, - std::array pVec, float etaVal, float phiVal, - float yVal, float decayLen, float decayLenXY, - float decayLenNorm, float decayLenXYNorm, - float cpaVal, float impParXY) - : m_xPVertex(xpv), m_yPVertex(ypv), m_zPVertex(zpv), - m_xSecondaryVertex(xsv), m_ySecondaryVertex(ysv), m_zSecondaryVertex(zsv), - m_px(pxVal), m_py(pyVal), m_pz(pzVal), m_e(eVal), m_m(mVal), m_chi2PCA(chi2Val), - m_errorDecayLength(errDecayLength), m_errorDecayLengthXY(errDecayLengthXY), - m_rSecondaryVertex(rSecVertex), m_pt(ptVal), m_p(pVal), - m_pVector(pVec), m_eta(etaVal), m_phi(phiVal), - m_y(yVal), m_decayLength(decayLen), m_decayLengthXY(decayLenXY), - m_decayLengthNormalised(decayLenNorm), m_decayLengthXYNormalised(decayLenXYNorm), - m_cpa(cpaVal), m_impactParameterXY(impParXY) - {} - - float xPVertex() const { return m_xPVertex; } - float yPVertex() const { return m_yPVertex; } - float zPVertex() const { return m_zPVertex; } - - float xSecondaryVertex() const { return m_xSecondaryVertex; } - float ySecondaryVertex() const { return m_ySecondaryVertex; } - float zSecondaryVertex() const { return m_zSecondaryVertex; } - - float px() const { return m_px; } - float py() const { return m_py; } - float pz() const { return m_pz; } - float e() const { return m_e; } - float m() const { return m_m; } - float chi2PCA() const { return m_chi2PCA; } - - float errorDecayLength() const { return m_errorDecayLength; } - float errorDecayLengthXY() const { return m_errorDecayLengthXY; } - - float rSecondaryVertex() const { return m_rSecondaryVertex; } - float pt() const { return m_pt; } - float p() const { return m_p; } - - std::array pVector() const { return m_pVector; } - - float eta() const { return m_eta; } - float phi() const { return m_phi; } - float y() const { return m_y; } - - float decayLength() const { return m_decayLength; } - float decayLengthXY() const { return m_decayLengthXY; } - float decayLengthNormalised() const { return m_decayLengthNormalised; } - float decayLengthXYNormalised() const { return m_decayLengthXYNormalised; } - - float cpa() const { return m_cpa; } - float impactParameterXY() const { return m_impactParameterXY; } - - private: - float m_xPVertex, m_yPVertex, m_zPVertex; - float m_xSecondaryVertex, m_ySecondaryVertex, m_zSecondaryVertex; - float m_px, m_py, m_pz, m_e, m_m, m_chi2PCA; - float m_errorDecayLength, m_errorDecayLengthXY; - float m_rSecondaryVertex, m_pt, m_p; - std::array m_pVector; - float m_eta, m_phi, m_y; - float m_decayLength, m_decayLengthXY, m_decayLengthNormalised, m_decayLengthXYNormalised; - float m_cpa, m_impactParameterXY; - }; - - template - bjetCandSV jetFromProngMaxDecayLength(const JetType& jet, bool doXYZ=false) +class bjetCandSV +{ + public: + bjetCandSV() = default; + + bjetCandSV(float xpv, float ypv, float zpv, float xsv, float ysv, float zsv, + float pxVal, float pyVal, float pzVal, float eVal, float mVal, float chi2Val, + float errDecayLength, float errDecayLengthXY, + float rSecVertex, float ptVal, float pVal, + std::array pVec, float etaVal, float phiVal, + float yVal, float decayLen, float decayLenXY, + float decayLenNorm, float decayLenXYNorm, + float cpaVal, float impParXY) + : m_xPVertex(xpv), m_yPVertex(ypv), m_zPVertex(zpv), m_xSecondaryVertex(xsv), m_ySecondaryVertex(ysv), m_zSecondaryVertex(zsv), m_px(pxVal), m_py(pyVal), m_pz(pzVal), m_e(eVal), m_m(mVal), m_chi2PCA(chi2Val), m_errorDecayLength(errDecayLength), m_errorDecayLengthXY(errDecayLengthXY), m_rSecondaryVertex(rSecVertex), m_pt(ptVal), m_p(pVal), m_pVector(pVec), m_eta(etaVal), m_phi(phiVal), m_y(yVal), m_decayLength(decayLen), m_decayLengthXY(decayLenXY), m_decayLengthNormalised(decayLenNorm), m_decayLengthXYNormalised(decayLenXYNorm), m_cpa(cpaVal), m_impactParameterXY(impParXY) { - float xPVertex = 0.0f; - float yPVertex = 0.0f; - float zPVertex = 0.0f; - float xSecondaryVertex = 0.0f; - float ySecondaryVertex = 0.0f; - float zSecondaryVertex = 0.0f; - float px = 0.0f; - float py = 0.0f; - float pz = 0.0f; - float e = 0.0f; - float m = 0.0f; - float chi2PCA = 0.0f; - float errorDecayLength = 0.0f; - float errorDecayLengthXY = 0.0f; - - float rSecondaryVertex = 0.0f; - float pt = 0.0f; - float p = 0.0f; - std::array pVector = {0.0f, 0.0f, 0.0f}; - float eta = 0.0f; - float phi = 0.0f; - float y = 0.0f; - float decayLength = 0.0f; - float decayLengthXY = 0.0f; - float decayLengthNormalised = 0.0f; - float decayLengthXYNormalised = 0.0f; - float cpa = 0.0f; - float impactParameterXY = 0.0f; - - float maxSxy = -1.0f; - - for (const auto& prong : jet.template secondaryVertices_as()) { - float Sxy = -1.; - if (!doXYZ) { - Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); - } else { - Sxy = prong.decayLength() / prong.errorDecayLength(); - } - //if (!prongAcceptance(prong, doXYZ)) continue; - - if (maxSxy < Sxy) { - maxSxy = Sxy; - - xPVertex = prong.xPVertex(); - yPVertex = prong.yPVertex(); - zPVertex = prong.zPVertex(); - xSecondaryVertex = prong.xSecondaryVertex(); - ySecondaryVertex = prong.ySecondaryVertex(); - zSecondaryVertex = prong.zSecondaryVertex(); - px = prong.px(); - py = prong.py(); - pz = prong.pz(); - e = prong.e(); - m = prong.m(); - chi2PCA = prong.chi2PCA(); - errorDecayLength = prong.errorDecayLength(); - errorDecayLengthXY = prong.errorDecayLengthXY(); - rSecondaryVertex = prong.rSecondaryVertex(); - pt = prong.pt(); - p = prong.p(); - pVector = prong.pVector(); - eta = prong.eta(); - phi = prong.phi(); - y = prong.y(); - decayLength = prong.decayLength(); - decayLengthXY = prong.decayLengthXY(); - decayLengthNormalised = prong.decayLengthNormalised(); - decayLengthXYNormalised = prong.decayLengthXYNormalised(); - cpa = prong.cpa(); - impactParameterXY = prong.impactParameterXY(); - } - } - - return bjetCandSV( - xPVertex, yPVertex, zPVertex, - xSecondaryVertex, ySecondaryVertex, zSecondaryVertex, - px, py, pz, e, m, chi2PCA, - errorDecayLength, errorDecayLengthXY, - rSecondaryVertex, pt, p, - pVector, eta, phi, - y, decayLength, decayLengthXY, - decayLengthNormalised, decayLengthXYNormalised, - cpa, impactParameterXY - ); } - template - bool isTaggedJetSV (T const jet, U const& /*prongs*/, float const& doXYZ=false, float const& tagPointForSV = 15.) - { - auto bjetCand = jetFromProngMaxDecayLength(jet); + float xPVertex() const { return m_xPVertex; } + float yPVertex() const { return m_yPVertex; } + float zPVertex() const { return m_zPVertex; } + + float xSecondaryVertex() const { return m_xSecondaryVertex; } + float ySecondaryVertex() const { return m_ySecondaryVertex; } + float zSecondaryVertex() const { return m_zSecondaryVertex; } + + float px() const { return m_px; } + float py() const { return m_py; } + float pz() const { return m_pz; } + float e() const { return m_e; } + float m() const { return m_m; } + float chi2PCA() const { return m_chi2PCA; } + + float errorDecayLength() const { return m_errorDecayLength; } + float errorDecayLengthXY() const { return m_errorDecayLengthXY; } + + float rSecondaryVertex() const { return m_rSecondaryVertex; } + float pt() const { return m_pt; } + float p() const { return m_p; } + + std::array pVector() const { return m_pVector; } + + float eta() const { return m_eta; } + float phi() const { return m_phi; } + float y() const { return m_y; } + + float decayLength() const { return m_decayLength; } + float decayLengthXY() const { return m_decayLengthXY; } + float decayLengthNormalised() const { return m_decayLengthNormalised; } + float decayLengthXYNormalised() const { return m_decayLengthXYNormalised; } + + float cpa() const { return m_cpa; } + float impactParameterXY() const { return m_impactParameterXY; } + + private: + float m_xPVertex, m_yPVertex, m_zPVertex; + float m_xSecondaryVertex, m_ySecondaryVertex, m_zSecondaryVertex; + float m_px, m_py, m_pz, m_e, m_m, m_chi2PCA; + float m_errorDecayLength, m_errorDecayLengthXY; + float m_rSecondaryVertex, m_pt, m_p; + std::array m_pVector; + float m_eta, m_phi, m_y; + float m_decayLength, m_decayLengthXY, m_decayLengthNormalised, m_decayLengthXYNormalised; + float m_cpa, m_impactParameterXY; +}; + +template +bjetCandSV jetFromProngMaxDecayLength(const JetType& jet, bool doXYZ = false) +{ + float xPVertex = 0.0f; + float yPVertex = 0.0f; + float zPVertex = 0.0f; + float xSecondaryVertex = 0.0f; + float ySecondaryVertex = 0.0f; + float zSecondaryVertex = 0.0f; + float px = 0.0f; + float py = 0.0f; + float pz = 0.0f; + float e = 0.0f; + float m = 0.0f; + float chi2PCA = 0.0f; + float errorDecayLength = 0.0f; + float errorDecayLengthXY = 0.0f; + + float rSecondaryVertex = 0.0f; + float pt = 0.0f; + float p = 0.0f; + std::array pVector = {0.0f, 0.0f, 0.0f}; + float eta = 0.0f; + float phi = 0.0f; + float y = 0.0f; + float decayLength = 0.0f; + float decayLengthXY = 0.0f; + float decayLengthNormalised = 0.0f; + float decayLengthXYNormalised = 0.0f; + float cpa = 0.0f; + float impactParameterXY = 0.0f; + + float maxSxy = -1.0f; + + for (const auto& prong : jet.template secondaryVertices_as()) { + float Sxy = -1.; if (!doXYZ) { - auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); - if (maxSxy < tagPointForSV) return false; + Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); } else { - auto maxSxyz = bjetCand.decayLength() / bjetCand.errorDecayLength(); - if (maxSxyz < tagPointForSV) return false; + Sxy = prong.decayLength() / prong.errorDecayLength(); + } + // if (!prongAcceptance(prong, doXYZ)) continue; + + if (maxSxy < Sxy) { + maxSxy = Sxy; + + xPVertex = prong.xPVertex(); + yPVertex = prong.yPVertex(); + zPVertex = prong.zPVertex(); + xSecondaryVertex = prong.xSecondaryVertex(); + ySecondaryVertex = prong.ySecondaryVertex(); + zSecondaryVertex = prong.zSecondaryVertex(); + px = prong.px(); + py = prong.py(); + pz = prong.pz(); + e = prong.e(); + m = prong.m(); + chi2PCA = prong.chi2PCA(); + errorDecayLength = prong.errorDecayLength(); + errorDecayLengthXY = prong.errorDecayLengthXY(); + rSecondaryVertex = prong.rSecondaryVertex(); + pt = prong.pt(); + p = prong.p(); + pVector = prong.pVector(); + eta = prong.eta(); + phi = prong.phi(); + y = prong.y(); + decayLength = prong.decayLength(); + decayLengthXY = prong.decayLengthXY(); + decayLengthNormalised = prong.decayLengthNormalised(); + decayLengthXYNormalised = prong.decayLengthXYNormalised(); + cpa = prong.cpa(); + impactParameterXY = prong.impactParameterXY(); } - return true; } + return bjetCandSV( + xPVertex, yPVertex, zPVertex, + xSecondaryVertex, ySecondaryVertex, zSecondaryVertex, + px, py, pz, e, m, chi2PCA, + errorDecayLength, errorDecayLengthXY, + rSecondaryVertex, pt, p, + pVector, eta, phi, + y, decayLength, decayLengthXY, + decayLengthNormalised, decayLengthXYNormalised, + cpa, impactParameterXY); +} + +template +bool isTaggedJetSV(T const jet, U const& /*prongs*/, float const& doXYZ = false, float const& tagPointForSV = 15.) +{ + auto bjetCand = jetFromProngMaxDecayLength(jet); + if (!doXYZ) { + auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); + if (maxSxy < tagPointForSV) + return false; + } else { + auto maxSxyz = bjetCand.decayLength() / bjetCand.errorDecayLength(); + if (maxSxyz < tagPointForSV) + return false; + } + return true; +} }; // namespace jettaggingutilities diff --git a/PWGJE/TableProducer/jettaggerhf.cxx b/PWGJE/TableProducer/jettaggerhf.cxx index 7f281dc639c..ba5ea8b34b0 100644 --- a/PWGJE/TableProducer/jettaggerhf.cxx +++ b/PWGJE/TableProducer/jettaggerhf.cxx @@ -93,9 +93,9 @@ struct JetTaggerHFTask { std::unique_ptr fSignImpXYSigBeautyJetMC = nullptr; std::unique_ptr fSignImpXYSigLfJetMC = nullptr; - template - void calculateJetProbabilityMCD(int origin, T const& collision, U const& mcdjet, JetTagTracksMCD const& jtracks, OriTracksMCD const& tracks, std::vector& jetProb) { + void calculateJetProbabilityMCD(int origin, T const& collision, U const& mcdjet, JetTagTracksMCD const& jtracks, OriTracksMCD const& tracks, std::vector& jetProb) + { jetProb.clear(); jetProb.reserve(maxOrder); for (int order = 0; order < maxOrder; order++) { @@ -116,11 +116,14 @@ struct JetTaggerHFTask { } template - void evaluateTrackProbQA(int origin, T const& collision, U const& mcdjet, JetTagTracksMCD const& /*jtracks*/, OriTracksMCD const& /*tracks*/) { + void evaluateTrackProbQA(int origin, T const& collision, U const& mcdjet, JetTagTracksMCD const& /*jtracks*/, OriTracksMCD const& /*tracks*/) + { for (auto& jtrack : mcdjet.template tracks_as()) { auto track = jtrack.template track_as(); - if (!track.has_mcParticle()) continue; - if (!jettaggingutilities::trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) continue; + if (!track.has_mcParticle()) + continue; + if (!jettaggingutilities::trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) + continue; auto geoSign = jettaggingutilities::getGeoSign(collision, mcdjet, track); float probTrack = -1; if (useResoFuncFromIncJet) { @@ -144,7 +147,6 @@ struct JetTaggerHFTask { } } - HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; void init(InitContext const&) { diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index 7d501a73e21..d1565ffb602 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -325,7 +325,8 @@ struct JetTaggerHFQA { auto track = jtrack.template track_as(); if (!trackAcceptance(track)) continue; - if (!jettaggingutilities::trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) continue; + if (!jettaggingutilities::trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) + continue; // General parameters registry.fill(HIST("h3_jet_pt_track_pt_track_eta"), jet.pt(), track.pt(), track.eta()); registry.fill(HIST("h3_jet_pt_track_pt_track_phi"), jet.pt(), track.pt(), track.phi()); @@ -423,7 +424,8 @@ struct JetTaggerHFQA { auto track = jtrack.template track_as(); if (!trackAcceptance(track)) continue; - if (!jettaggingutilities::trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) continue; + if (!jettaggingutilities::trackAcceptanceWithDca(track, trackDcaXYMax, trackDcaZMax)) + continue; // General parameters registry.fill(HIST("h3_jet_pt_track_pt_flavour"), mcdjet.pt(), track.pt(), jetflavour); @@ -577,7 +579,7 @@ struct JetTaggerHFQA { float eventWeight = mcdjet.eventWeight(); int jetflavour = mcdjet.origin(); int jetflavourRun2Def = -1; - //if (!mcdjet.has_matchedJetGeo()) continue; + // if (!mcdjet.has_matchedJetGeo()) continue; for (auto& mcpjet : mcdjet.template matchedJetGeo_as()) { jetflavourRun2Def = jettaggingutilities::getJetFlavor(mcpjet, particlesPerColl); registry.fill(HIST("h3_response_matrix_jet_pt_jet_pt_part_flavour"), mcdjet.pt(), mcpjet.pt(), jetflavour, eventWeight); @@ -689,7 +691,8 @@ struct JetTaggerHFQA { } auto origin = mcdjet.origin(); registry.fill(HIST("h2_2prong_nprongs_flavour"), mcdjet.template secondaryVertices_as().size(), origin); - if (mcdjet.template secondaryVertices_as().size() < 1) continue; + if (mcdjet.template secondaryVertices_as().size() < 1) + continue; for (const auto& prong : mcdjet.template secondaryVertices_as()) { auto Lxy = prong.decayLengthXY(); auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); @@ -722,7 +725,8 @@ struct JetTaggerHFQA { } auto origin = mcdjet.origin(); registry.fill(HIST("h2_3prong_nprongs_flavour"), mcdjet.template secondaryVertices_as().size(), origin); - if (mcdjet.template secondaryVertices_as().size() < 1) continue; + if (mcdjet.template secondaryVertices_as().size() < 1) + continue; for (const auto& prong : mcdjet.template secondaryVertices_as()) { auto Lxy = prong.decayLengthXY(); auto Sxy = prong.decayLengthXY() / prong.errorDecayLengthXY(); From 73f65a0e6afff5846b982dc7a67c9a857f13ddd5 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Fri, 23 Aug 2024 15:39:39 +0200 Subject: [PATCH 20/21] fix clang of datamodel --- PWGJE/DataModel/JetTagging.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGJE/DataModel/JetTagging.h b/PWGJE/DataModel/JetTagging.h index 860572af981..7229fabd6ce 100644 --- a/PWGJE/DataModel/JetTagging.h +++ b/PWGJE/DataModel/JetTagging.h @@ -123,8 +123,8 @@ JETSV_TABLES_DEF(Charged, SecondaryVertex2Prong, "2PRONG"); { \ DECLARE_SOA_COLUMN(Origin, origin, int); \ DECLARE_SOA_COLUMN(JetProb, jetProb, std::vector); \ - DECLARE_SOA_COLUMN(FlagtaggedjetIP, flagtaggedjetIP, bool); \ - DECLARE_SOA_COLUMN(FlagtaggedjetSV, flagtaggedjetSV, bool); \ + DECLARE_SOA_COLUMN(FlagtaggedjetIP, flagtaggedjetIP, bool); \ + DECLARE_SOA_COLUMN(FlagtaggedjetSV, flagtaggedjetSV, bool); \ } \ DECLARE_SOA_TABLE(_jet_type_##Tags, "AOD", _description_ "Tags", _name_##tagging::Origin, _name_##tagging::JetProb, _name_##tagging::FlagtaggedjetIP, _name_##tagging::FlagtaggedjetSV); From f066d186d8bb79cede1e4a22e2aca722e55b96a7 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Fri, 23 Aug 2024 16:09:53 +0200 Subject: [PATCH 21/21] Add prong acceptance --- PWGJE/Core/JetTaggingUtilities.h | 20 +++++++++++++------- PWGJE/TableProducer/jettaggerhf.cxx | 12 ++++++++++-- PWGJE/Tasks/jettaggerhfQA.cxx | 16 ++++++++-------- 3 files changed, 31 insertions(+), 17 deletions(-) diff --git a/PWGJE/Core/JetTaggingUtilities.h b/PWGJE/Core/JetTaggingUtilities.h index 56c9fade744..590b495f570 100644 --- a/PWGJE/Core/JetTaggingUtilities.h +++ b/PWGJE/Core/JetTaggingUtilities.h @@ -340,12 +340,17 @@ bool trackAcceptanceWithDca(T const& track, float trackDcaXYMax, float trackDcaZ * retrun acceptance of prong about chi2 and error of decay length due to cut for high quality secondary vertex */ template -bool prongAcceptance(T const& prong, float prongChi2PCAMax, float prongsigmaLxyMax) +bool prongAcceptance(T const& prong, float prongChi2PCAMax, float prongsigmaLxyMax, bool doXYZ) { if (prong.chi2PCA() > prongChi2PCAMax) return false; - if (prong.errorDecayLengthXY() > prongsigmaLxyMax) - return false; + if (!doXYZ) { + if (prong.errorDecayLengthXY() > prongsigmaLxyMax) + return false; + } else { + if (prong.errorDecayLength() > prongsigmaLxyMax) + return false; + } return true; } @@ -578,7 +583,7 @@ class bjetCandSV }; template -bjetCandSV jetFromProngMaxDecayLength(const JetType& jet, bool doXYZ = false) +bjetCandSV jetFromProngMaxDecayLength(const JetType& jet, float const& prongChi2PCAMax, float const& prongsigmaLxyMax, const bool& doXYZ = false) { float xPVertex = 0.0f; float yPVertex = 0.0f; @@ -618,7 +623,8 @@ bjetCandSV jetFromProngMaxDecayLength(const JetType& jet, bool doXYZ = false) } else { Sxy = prong.decayLength() / prong.errorDecayLength(); } - // if (!prongAcceptance(prong, doXYZ)) continue; + if (!prongAcceptance(prong, prongChi2PCAMax, prongsigmaLxyMax, doXYZ)) + continue; if (maxSxy < Sxy) { maxSxy = Sxy; @@ -666,9 +672,9 @@ bjetCandSV jetFromProngMaxDecayLength(const JetType& jet, bool doXYZ = false) } template -bool isTaggedJetSV(T const jet, U const& /*prongs*/, float const& doXYZ = false, float const& tagPointForSV = 15.) +bool isTaggedJetSV(T const jet, U const& /*prongs*/, float const& prongChi2PCAMax, float const& prongsigmaLxyMax, float const& doXYZ = false, float const& tagPointForSV = 15.) { - auto bjetCand = jetFromProngMaxDecayLength(jet); + auto bjetCand = jetFromProngMaxDecayLength(jet, prongChi2PCAMax, prongsigmaLxyMax, doXYZ); if (!doXYZ) { auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); if (maxSxy < tagPointForSV) diff --git a/PWGJE/TableProducer/jettaggerhf.cxx b/PWGJE/TableProducer/jettaggerhf.cxx index ba5ea8b34b0..28a08a5aada 100644 --- a/PWGJE/TableProducer/jettaggerhf.cxx +++ b/PWGJE/TableProducer/jettaggerhf.cxx @@ -252,7 +252,11 @@ struct JetTaggerHFTask { } if (jettaggingutilities::isGreaterThanTaggingPoint(collision, jet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount)) flagtaggedjetIP = true; - flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(jet, prongs, useXYZForTagging, tagPointForSV); + if (!useXYZForTagging) { + flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(jet, prongs, prongChi2PCAMax, prongsigmaLxyMax, useXYZForTagging, tagPointForSV); + } else { + flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(jet, prongs, prongChi2PCAMax, prongsigmaLxyzMax, useXYZForTagging, tagPointForSV); + } taggingTableData(0, jetProb, flagtaggedjetIP, flagtaggedjetSV); } } @@ -301,7 +305,11 @@ struct JetTaggerHFTask { } if (jettaggingutilities::isGreaterThanTaggingPoint(collision, mcdjet, jtracks, tracks, trackDcaXYMax, trackDcaZMax, tagPointForIP, minIPCount)) flagtaggedjetIP = true; - flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(mcdjet, prongs, useXYZForTagging, tagPointForSV); + if (!useXYZForTagging) { + flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(mcdjet, prongs, prongChi2PCAMax, prongsigmaLxyMax, useXYZForTagging, tagPointForSV); + } else { + flagtaggedjetSV = jettaggingutilities::isTaggedJetSV(mcdjet, prongs, prongChi2PCAMax, prongsigmaLxyzMax, useXYZForTagging, tagPointForSV); + } taggingTableMCD(origin, jetProb, flagtaggedjetIP, flagtaggedjetSV); } } diff --git a/PWGJE/Tasks/jettaggerhfQA.cxx b/PWGJE/Tasks/jettaggerhfQA.cxx index d1565ffb602..255a09ad3b6 100644 --- a/PWGJE/Tasks/jettaggerhfQA.cxx +++ b/PWGJE/Tasks/jettaggerhfQA.cxx @@ -644,9 +644,9 @@ struct JetTaggerHFQA { registry.fill(HIST("h2_jet_pt_2prong_sigmaLxy"), jet.pt(), prong.errorDecayLengthXY()); registry.fill(HIST("h2_jet_pt_2prong_sigmaLxyz"), jet.pt(), prong.errorDecayLength()); } - auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(jet); + auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(jet, prongChi2PCAMax, prongsigmaLxyMax); auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); - auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(jet, true); + auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(jet, prongChi2PCAMax, prongsigmaLxyzMax, true); auto maxSxyz = bjetCandForXYZ.decayLength() / bjetCandForXYZ.errorDecayLength(); registry.fill(HIST("h2_jet_pt_2prong_Sxy_N1"), jet.pt(), maxSxy); registry.fill(HIST("h2_jet_pt_2prong_Sxyz_N1"), jet.pt(), maxSxyz); @@ -673,9 +673,9 @@ struct JetTaggerHFQA { registry.fill(HIST("h2_jet_pt_3prong_sigmaLxy"), jet.pt(), prong.errorDecayLengthXY()); registry.fill(HIST("h2_jet_pt_3prong_sigmaLxyz"), jet.pt(), prong.errorDecayLength()); } - auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(jet); + auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(jet, prongChi2PCAMax, prongsigmaLxyMax); auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); - auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(jet, true); + auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(jet, prongChi2PCAMax, prongsigmaLxyzMax, true); auto maxSxyz = bjetCandForXYZ.decayLength() / bjetCandForXYZ.errorDecayLength(); registry.fill(HIST("h2_jet_pt_3prong_Sxy_N1"), jet.pt(), maxSxy); registry.fill(HIST("h2_jet_pt_3prong_Sxyz_N1"), jet.pt(), maxSxyz); @@ -705,10 +705,10 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_jet_pt_2prong_sigmaLxy_flavour"), mcdjet.pt(), prong.errorDecayLengthXY(), origin); registry.fill(HIST("h3_jet_pt_2prong_sigmaLxyz_flavour"), mcdjet.pt(), prong.errorDecayLength(), origin); } - auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet); + auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, prongChi2PCAMax, prongsigmaLxyMax); auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); auto massSV = bjetCand.m(); - auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, true); + auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, prongChi2PCAMax, prongsigmaLxyzMax, true); auto maxSxyz = bjetCandForXYZ.decayLength() / bjetCandForXYZ.errorDecayLength(); registry.fill(HIST("h3_jet_pt_2prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); registry.fill(HIST("h3_jet_pt_2prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin); @@ -739,10 +739,10 @@ struct JetTaggerHFQA { registry.fill(HIST("h3_jet_pt_3prong_sigmaLxy_flavour"), mcdjet.pt(), prong.errorDecayLengthXY(), origin); registry.fill(HIST("h3_jet_pt_3prong_sigmaLxyz_flavour"), mcdjet.pt(), prong.errorDecayLength(), origin); } - auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet); + auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, prongChi2PCAMax, prongsigmaLxyMax); auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); auto massSV = bjetCand.m(); - auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, true); + auto bjetCandForXYZ = jettaggingutilities::jetFromProngMaxDecayLength(mcdjet, prongChi2PCAMax, prongsigmaLxyzMax, true); auto maxSxyz = bjetCandForXYZ.decayLength() / bjetCandForXYZ.errorDecayLength(); registry.fill(HIST("h3_jet_pt_3prong_Sxy_N1_flavour"), mcdjet.pt(), maxSxy, origin); registry.fill(HIST("h3_jet_pt_3prong_Sxyz_N1_flavour"), mcdjet.pt(), maxSxyz, origin);