Skip to content

[PWGJE] modification for systematic uncertainty of b-jets and multiplicity analysis #15830

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
nzardosh merged 12 commits into from
Apr 17, 2026
Merged

[PWGJE] modification for systematic uncertainty of b-jets and multiplicity analysis #15830

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
3254b2b
Being developing mult of b-jet
hanseopark Jan 22, 2026
9f9e67a
Add bjet multiplcity analysis
hanseopark Feb 2, 2026
5d87074
fix CMakeList
hanseopark Feb 2, 2026
b3167ef
update bkg sub b-jet mult
hanseopark Apr 16, 2026
f8dcc12
fix to complict rebase
hanseopark Apr 17, 2026
b7049a2
modificaiton validation of flavour and seconday contamition
hanseopark Apr 17, 2026
a670bf8
fix CMakeList
hanseopark Apr 17, 2026
4ac9922
update background subtraction of b-jets for multiplicity
hanseopark Apr 17, 2026
38b0955
fix o2linter
hanseopark Apr 17, 2026
e827389
fix o2linter
hanseopark Apr 17, 2026
2ec4bde
fix magic number
hanseopark Apr 17, 2026
b6a2eff
fix o2linter
hanseopark Apr 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 23 additions & 21 deletions PWGJE/TableProducer/secondaryVertexReconstruction.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,22 +22,22 @@
#include "Common/Core/trackUtilities.h"
#include "Common/DataModel/TrackSelectionTables.h"

#include "CommonConstants/PhysicsConstants.h"
#include "DCAFitter/DCAFitterN.h"
#include "Framework/ASoA.h"
#include "Framework/AnalysisDataModel.h"
#include "Framework/AnalysisTask.h"
#include "ReconstructionDataFormats/DCA.h"
Comment on lines +25 to +30
Copy link
Copy Markdown
Collaborator

@vkucera vkucera Apr 17, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is wrong. Why did you do this?

#include <CCDB/BasicCCDBManager.h>
#include <CommonConstants/PhysicsConstants.h>
#include <DCAFitter/DCAFitterN.h>
#include <DetectorsBase/MatLayerCylSet.h>
#include <DetectorsBase/Propagator.h>
#include <Framework/ASoA.h>
#include <Framework/AnalysisDataModel.h>
#include <Framework/AnalysisHelpers.h>
#include <Framework/AnalysisTask.h>
#include <Framework/Configurable.h>
#include <Framework/HistogramRegistry.h>
#include <Framework/HistogramSpec.h>
#include <Framework/InitContext.h>
#include <Framework/Logger.h>
#include <Framework/runDataProcessing.h>
#include <ReconstructionDataFormats/DCA.h>
#include <ReconstructionDataFormats/TrackParametrizationWithError.h>
#include <ReconstructionDataFormats/Vertex.h>

Expand Down Expand Up @@ -74,12 +74,12 @@ struct SecondaryVertexReconstruction {
Configurable<bool> propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"};
Configurable<bool> useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"};
Configurable<bool> useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"};
Configurable<double> maxR{"maxR", 200., "reject PCA's above this radius"};
Configurable<double> maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"};
Configurable<double> maxRsv{"maxRsv", 999., "max. radius of the reconstruced SV"};
Configurable<double> maxZsv{"maxZsv", 999., "max. Z coordinates of the reconstruced SV"};
Configurable<double> minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"};
Configurable<double> minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"};
Configurable<float> maxR{"maxR", 200., "reject PCA's above this radius"};
Configurable<float> maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"};
Configurable<float> maxRsv{"maxRsv", 999., "max. radius of the reconstruced SV"};
Configurable<float> maxZsv{"maxZsv", 999., "max. Z coordinates of the reconstruced SV"};
Configurable<float> minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"};
Configurable<float> minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"};
Configurable<float> ptMinTrack{"ptMinTrack", -1., "min. track pT"};
Configurable<float> etaMinTrack{"etaMinTrack", -99999., "min. pseudorapidity"};
Configurable<float> etaMaxTrack{"etaMaxTrack", 4., "max. pseudorapidity"};
Expand All @@ -94,13 +94,15 @@ struct SecondaryVertexReconstruction {
o2::vertexing::DCAFitterN<2> df2; // 2-prong vertex fitter
o2::vertexing::DCAFitterN<3> df3; // 3-prong vertex fitter
Service<o2::ccdb::BasicCCDBManager> ccdb;
o2::base::MatLayerCylSet* lut;
o2::base::MatLayerCylSet* lut = nullptr;

HistogramRegistry registry{"registry"};

int runNumber{0};
float toMicrometers = 10000.; // from cm to µm
double bz{0.};
static constexpr int TwoProngCount = 2;
static constexpr int ThreeProngCount = 3;

void init(InitContext const&)
{
Expand Down Expand Up @@ -216,7 +218,7 @@ struct SecondaryVertexReconstruction {
auto covMatrixPV = primaryVertex.getCov();

// Get track momenta and impact parameters
std::array<std::array<float, 3>, numProngs> arrayMomenta;
std::array<std::array<float, 3>, numProngs> arrayMomenta{};
std::array<o2::dataformats::DCA, numProngs> impactParameters;
for (unsigned int inum = 0; inum < numProngs; ++inum) {
trackParVars[inum].getPxPyPzGlo(arrayMomenta[inum]);
Expand All @@ -236,12 +238,12 @@ struct SecondaryVertexReconstruction {
auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.));

// calculate invariant mass
std::array<double, numProngs> massArray;
std::array<double, numProngs> massArray{};
std::fill(massArray.begin(), massArray.end(), o2::constants::physics::MassPiPlus);
double massSV = RecoDecay::m(std::move(arrayMomenta), massArray);
double massSV = RecoDecay::m(arrayMomenta, massArray);

// fill candidate table rows
if ((doprocessData3Prongs || doprocessData3ProngsExternalMagneticField) && numProngs == 3) {
if ((doprocessData3Prongs || doprocessData3ProngsExternalMagneticField) && numProngs == ThreeProngCount) {
sv3prongTableData(analysisJet.globalIndex(),
primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(),
secondaryVertex[0], secondaryVertex[1], secondaryVertex[2],
Expand All @@ -250,7 +252,7 @@ struct SecondaryVertexReconstruction {
arrayMomenta[0][2] + arrayMomenta[1][2] + arrayMomenta[2][2],
energySV, massSV, chi2PCA, dispersion, errorDecayLength, errorDecayLengthXY);
svIndices.push_back(sv3prongTableData.lastIndex());
} else if ((doprocessData2Prongs || doprocessData2ProngsExternalMagneticField) && numProngs == 2) {
} else if ((doprocessData2Prongs || doprocessData2ProngsExternalMagneticField) && numProngs == TwoProngCount) {
sv2prongTableData(analysisJet.globalIndex(),
primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(),
secondaryVertex[0], secondaryVertex[1], secondaryVertex[2],
Expand All @@ -259,7 +261,7 @@ struct SecondaryVertexReconstruction {
arrayMomenta[0][2] + arrayMomenta[1][2],
energySV, massSV, chi2PCA, dispersion, errorDecayLength, errorDecayLengthXY);
svIndices.push_back(sv2prongTableData.lastIndex());
} else if ((doprocessMCD3Prongs || doprocessMCD3ProngsExternalMagneticField) && numProngs == 3) {
} else if ((doprocessMCD3Prongs || doprocessMCD3ProngsExternalMagneticField) && numProngs == ThreeProngCount) {
sv3prongTableMCD(analysisJet.globalIndex(),
primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(),
secondaryVertex[0], secondaryVertex[1], secondaryVertex[2],
Expand All @@ -268,7 +270,7 @@ struct SecondaryVertexReconstruction {
arrayMomenta[0][2] + arrayMomenta[1][2] + arrayMomenta[2][2],
energySV, massSV, chi2PCA, dispersion, errorDecayLength, errorDecayLengthXY);
svIndices.push_back(sv3prongTableMCD.lastIndex());
} else if ((doprocessMCD2Prongs || doprocessMCD2ProngsExternalMagneticField) && numProngs == 2) {
} else if ((doprocessMCD2Prongs || doprocessMCD2ProngsExternalMagneticField) && numProngs == TwoProngCount) {
sv2prongTableMCD(analysisJet.globalIndex(),
primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(),
secondaryVertex[0], secondaryVertex[1], secondaryVertex[2],
Expand Down Expand Up @@ -399,5 +401,5 @@ struct SecondaryVertexReconstruction {

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
{
return WorkflowSpec{adaptAnalysisTask<SecondaryVertexReconstruction>(cfgc, TaskName{"jet-sv-reconstruction-charged"})}; // o2-linter: disable=name/o2-task
return WorkflowSpec{adaptAnalysisTask<SecondaryVertexReconstruction>(cfgc, TaskName{"jet-sv-reconstruction-charged"})}; // o2-linter: disable=name/o2-task (templated struct)
}
80 changes: 80 additions & 0 deletions PWGJE/Tasks/bjetCentMult.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,25 @@ struct BjetCentMultTask {
registry.add("hn_taggedjet_3prong_Sxyz_N1_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisSxyz}, {axisMass}, {axisCentrality}}});
}
}
if (doprocessRhoAreaSubSV3ProngData) {
registry.add("h_event_centrality", "", {HistType::kTH1F, {{axisCentrality}}});
registry.add("h2_jet_pt_rhoareasubtracted_centrality", "", {HistType::kTH2F, {{axisJetPt}, {axisCentrality}}});
registry.add("h2_jet_eta_rhoareasubtracted_centrality", "", {HistType::kTH2F, {{axisEta}, {axisCentrality}}});
registry.add("h2_jet_phi_rhoareasubtracted_centrality", "", {HistType::kTH2F, {{axisPhi}, {axisCentrality}}});
if (fillGeneralSVQA) {
registry.add("h2_3prong_nprongs_rhoareasubtracted_centrality", "", {HistType::kTH2F, {{axisNprongs}, {axisCentrality}}});
registry.add("hn_jet_3prong_Sxy_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisLxy}, {axisSigmaLxy}, {axisSxy}, {axisCentrality}}});
if (fillSVxyz) {
registry.add("hn_jet_3prong_Sxyz_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisLxyz}, {axisSigmaLxyz}, {axisSxyz}, {axisCentrality}}});
}
}
registry.add("hn_jet_3prong_Sxy_N1_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisSxy}, {axisMass}, {axisCentrality}}});
registry.add("hn_taggedjet_3prong_Sxy_N1_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisSxy}, {axisMass}, {axisCentrality}}});
if (fillSVxyz) {
registry.add("hn_jet_3prong_Sxyz_N1_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisSxyz}, {axisMass}, {axisCentrality}}});
registry.add("hn_taggedjet_3prong_Sxyz_N1_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisSxyz}, {axisMass}, {axisCentrality}}});
}
}
if (doprocessSV3ProngMCD || doprocessSV3ProngMCPMCDMatched) {
registry.add("h_event_centrality", "", {HistType::kTH1F, {{axisCentrality}}});
registry.add("h3_jet_pt_centrality_flavour", "", {HistType::kTH3F, {{axisJetPt}, {axisCentrality}, {axisJetFlavour}}});
Expand Down Expand Up @@ -312,6 +331,47 @@ struct BjetCentMultTask {
registry.fill(HIST("h2_jet_phi_part_flavour"), mcpjet.phi(), jetflavour, eventWeight);
}

template <typename T, typename U>
void fillRhoAreaSubtractedHistogramSV3ProngData(T const& jet, U const& /*prongs*/, float centrality, float rho)
{
if (jet.template secondaryVertices_as<U>().size() < 1)
return;
registry.fill(HIST("h2_jet_pt_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), centrality);
registry.fill(HIST("h2_jet_eta_rhoareasubtracted_centrality"), jet.eta(), centrality);
registry.fill(HIST("h2_jet_phi_rhoareasubtracted_centrality"), jet.phi(), centrality);
if (fillGeneralSVQA) {
registry.fill(HIST("h2_3prong_nprongs_rhoareasubtracted_centrality"), jet.template secondaryVertices_as<U>().size(), centrality);
for (const auto& prong : jet.template secondaryVertices_as<U>()) {
registry.fill(HIST("hn_jet_3prong_Sxy_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), prong.decayLengthXY(), prong.errorDecayLengthXY(), prong.decayLengthXY() / prong.errorDecayLengthXY(), centrality);
if (fillSVxyz) {
registry.fill(HIST("hn_jet_3prong_Sxyz_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), prong.decayLength(), prong.errorDecayLength(), prong.decayLength() / prong.errorDecayLength(), centrality);
}
}
}
bool checkSv = false;
auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength<U>(jet, prongCuts->at(0), prongCuts->at(1), prongCuts->at(2), prongCuts->at(4), prongCuts->at(5), false, &checkSv);
if (checkSv && jettaggingutilities::svAcceptance(bjetCand, svDispersionMax)) {
auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY();
auto massSV = bjetCand.m();
registry.fill(HIST("hn_jet_3prong_Sxy_N1_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), maxSxy, massSV, centrality);
if (jet.isTagged(BJetTaggingMethod::SV)) {
registry.fill(HIST("hn_taggedjet_3prong_Sxy_N1_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), maxSxy, massSV, centrality);
}
}
if (fillSVxyz) {
checkSv = false;
auto bjetCandXYZ = jettaggingutilities::jetFromProngMaxDecayLength<U>(jet, prongCuts->at(0), prongCuts->at(1), prongCuts->at(3), prongCuts->at(4), prongCuts->at(5), true, &checkSv);
if (checkSv && jettaggingutilities::svAcceptance(bjetCandXYZ, svDispersionMax)) {
auto maxSxyz = bjetCandXYZ.decayLength() / bjetCandXYZ.errorDecayLength();
auto massSV = bjetCandXYZ.m();
registry.fill(HIST("hn_jet_3prong_Sxyz_N1_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), maxSxyz, massSV, centrality);
if (jet.isTagged(BJetTaggingMethod::SV3D)) {
registry.fill(HIST("hn_taggedjet_3prong_Sxyz_N1_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), maxSxyz, massSV, centrality);
}
}
}
}

template <typename T, typename U>
void fillHistogramSV3ProngData(T const& jet, U const& /*prongs*/, float centrality)
{
Expand Down Expand Up @@ -513,6 +573,26 @@ struct BjetCentMultTask {
}
PROCESS_SWITCH(BjetCentMultTask, processSV3ProngData, "Fill 3prong imformation for data jets", false);

void processRhoAreaSubSV3ProngData(soa::Filtered<soa::Join<aod::JetCollisions, aod::BkgChargedRhos>>::iterator const& collision, soa::Join<JetTableData, TagTableData, aod::DataSecondaryVertex3ProngIndices> const& jets, aod::DataSecondaryVertex3Prongs const& prongs)
{
if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) {
return;
}
float centrality = collision.centFT0M();
float rho = collision.rho();
registry.fill(HIST("h_event_centrality"), centrality);
for (auto const& jet : jets) {
if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaCuts->at(0), jetEtaCuts->at(1), trackCuts->at(2), trackCuts->at(3))) {
continue;
}
if (!isAcceptedJet<aod::JetTracks>(jet)) {
continue;
}
fillRhoAreaSubtractedHistogramSV3ProngData(jet, prongs, centrality, rho);
}
}
PROCESS_SWITCH(BjetCentMultTask, processRhoAreaSubSV3ProngData, "Fill 3prong imformation for data jets with background subtraction", false);

void processSV3ProngMCD(soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::JCollisionOutliers>>::iterator const& collision, soa::Join<JetTableMCD, TagTableMCD, aod::MCDSecondaryVertex3ProngIndices> const& mcdjets, aod::MCDSecondaryVertex3Prongs const& prongs)
{
if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) {
Expand Down
9 changes: 4 additions & 5 deletions PWGJE/Tasks/jetTaggerHFQA.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,7 @@ struct JetTaggerHFQA {
}
}
if (doprocessTracksInJetsData) {
registry.add("h2_track_pt_impact_parameter_xy", "", {HistType::kTH2F, {{axisTrackPt}, {axisImpactParameterXY}}});
registry.add("h3_jet_pt_track_pt_impact_parameter_xy", "", {HistType::kTH3F, {{axisJetPt}, {axisTrackPt}, {axisImpactParameterXY}}});
}
if (doprocessSecondaryContaminationMCD) {
registry.add("hn_jet_pt_track_pt_impact_parameter_xy_physical_primary_flavour", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisTrackPt}, {axisImpactParameterXY}, {axisJetFlavour}}});
Expand Down Expand Up @@ -1158,7 +1158,7 @@ struct JetTaggerHFQA {
}
for (auto const& track : jet.template tracks_as<JetTagTracksData>()) {
float varImpXY = track.dcaXY() * jettaggingutilities::cmTomum;
registry.fill(HIST("h2_track_pt_impact_parameter_xy"), track.pt(), varImpXY);
registry.fill(HIST("h3_jet_pt_track_pt_impact_parameter_xy"), jet.pt(), track.pt(), varImpXY);
}
}
}
Expand All @@ -1179,7 +1179,7 @@ struct JetTaggerHFQA {
if (!isAcceptedJet<aod::JetTracks>(mcdjet)) {
continue;
}
auto eventWeight = collision.mcCollision_as<soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>>().weight();
auto eventWeight = collision.weight();
float pTHat = 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
if (mcdjet.pt() > pTHatMaxMCD * pTHat) {
continue;
Expand Down Expand Up @@ -1224,8 +1224,7 @@ struct JetTaggerHFQA {
if (!isAcceptedJet<aod::JetTracks>(mcdjet)) {
continue;
}
auto eventWeight = collision.mcCollision_as<soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>>().weight();
fillValidationFlavourDefMCD<soa::Join<JetTableMCP, JetTableMCPMCD>>(mcdjet, tracks, particles, particlesPerColl, eventWeight);
fillValidationFlavourDefMCD<soa::Join<JetTableMCP, JetTableMCPMCD>>(mcdjet, tracks, particles, particlesPerColl, collision.weight());
}
}
PROCESS_SWITCH(JetTaggerHFQA, processValFlavourDefMCD, "to check the validation of jet-flavour definition when compared to distance for mcd jets", false);
Expand Down
Loading