Skip to content

[PWGUD] GlobalMuon matching based on chi2matching #15938

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
rolavick merged 5 commits into from
Apr 23, 2026
+70 −17
Merged

[PWGUD] GlobalMuon matching based on chi2matching #15938

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5023837
GlobalMuon matching based on chi2matching
rolavick Apr 23, 2026
b4d397b
Please consider the following formatting changes
alibuild Apr 23, 2026
9d512de
Merge pull request #76 from alibuild/alibot-cleanup-15938
rolavick Apr 23, 2026
cd98f9e
Making MegaLinter and O2linter happy
rolavick Apr 23, 2026
28047cf
fix unused variable
rolavick Apr 23, 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
87 changes: 70 additions & 17 deletions PWGUD/TableProducer/upcCandProducerGlobalMuon.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@
#include "PWGUD/Core/UPCHelpers.h"
#include "PWGUD/DataModel/UDTables.h"

#include "Common/Core/fwdtrackUtilities.h"

#include <CCDB/BasicCCDBManager.h>
#include <CCDB/CcdbApi.h>
#include <CommonConstants/LHCConstants.h>
Expand Down Expand Up @@ -55,6 +57,8 @@
#include <map>
#include <numeric>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

using namespace o2::framework;
Expand All @@ -64,6 +68,7 @@ struct UpcCandProducerGlobalMuon {
bool fDoMC{false};

std::map<int32_t, int32_t> fNewPartIDs;
std::map<uint32_t, bool> fBestMuonMatch;

Produces<o2::aod::UDMcCollisions> udMCCollisions;
Produces<o2::aod::UDMcParticles> udMCParticles;
Expand Down Expand Up @@ -97,6 +102,7 @@ struct UpcCandProducerGlobalMuon {
Configurable<bool> fUseAbsDCA{"fUseAbsDCA", true, "Use absolute DCA in FwdDCAFitter"};
Configurable<float> fMaxChi2MatchMCHMFT{"fMaxChi2MatchMCHMFT", 4.f, "Maximum chi2 for MCH-MFT matching quality filter"};
Configurable<int> fBcWindowMCHMFT{"fBcWindowMCHMFT", 20, "BC window for searching MCH-MFT tracks around MCH-MID-MFT anchors"};
Configurable<bool> fKeepBestMuonMatch{"fKeepBestMuonMatch", true, "Keep only the best MCH-MFT match per MCH track (lowest chi2)"};

using ForwardTracks = o2::soa::Join<o2::aod::FwdTracks, o2::aod::FwdTracksCov>;

Expand All @@ -111,9 +117,10 @@ struct UpcCandProducerGlobalMuon {
o2::vertexing::FwdDCAFitterN<2> fFwdFitter;

// Named constants (avoid magic numbers in expressions)
static constexpr double kBcTimeRoundingOffset = 1.; // Offset used when rounding trackTime to BC units
static constexpr uint16_t kMinTracksForPair = 2; // Minimum tracks required to compute a pair invariant mass
static constexpr uint16_t kMinTracksForCandidate = 1; // Minimum contributors required to save a candidate
static constexpr double kBcTimeRoundingOffset = 1.; // Offset used when rounding trackTime to BC units
static constexpr uint16_t kMinTracksForPair = 2; // Minimum tracks required to compute a pair invariant mass
static constexpr uint16_t kMinTracksForCandidate = 1; // Minimum contributors required to save a candidate
static constexpr int kUpperBoundaryToTrackTypeEnum = 2; // Make sure you use MFT tracks

void init(InitContext&)
{
Expand Down Expand Up @@ -384,13 +391,15 @@ struct UpcCandProducerGlobalMuon {
return propmuon;
}

bool addToFwdTable(int64_t candId, int64_t trackId, uint64_t gbc, float trackTime, ForwardTracks const& fwdTracks, const o2::aod::McFwdTrackLabels* mcFwdTrackLabels)
bool addToFwdTable(int64_t candId, int64_t trackId, uint64_t gbc, float trackTime,
ForwardTracks const& fwdTracks, o2::aod::MFTTracks const& /*mftTracks*/,
const o2::aod::McFwdTrackLabels* mcFwdTrackLabels)
{
const auto& track = fwdTracks.iteratorAt(trackId);
float px, py, pz;
int sign;

// NEW: Fill track type histogram
// Fill track type histogram
histRegistry.fill(HIST("hTrackTypes"), track.trackType());
if (track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack ||
track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalForwardTrack) {
Expand All @@ -402,12 +411,16 @@ struct UpcCandProducerGlobalMuon {
track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalForwardTrack;
o2::dataformats::GlobalFwdTrack pft;
if (isGlobal && fEnableMFT) {
// For global tracks, use raw kinematics for cuts;
// the FwdDCAFitter handles proper vertex propagation
auto trackPar = fwdToTrackPar(track);
pft.setParameters(trackPar.getParameters());
pft.setZ(trackPar.getZ());
pft.setCovariances(trackPar.getCovariances());
// Refit global muon: propagate MCH component to vertex, combine with MFT spatial info
auto mchTrack = track.matchMCHTrack_as<ForwardTracks>();
auto propMuon = propagateToZero(mchTrack);
auto mfttrack = track.matchMFTTrack_as<o2::aod::MFTTracks>();
using SMatrix5 = ROOT::Math::SVector<double, 5>;
using SMatrix55 = ROOT::Math::SMatrix<double, 5, 5, ROOT::Math::MatRepSym<double, 5>>;
SMatrix5 tpars(mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt());
SMatrix55 tcovs{};
o2::track::TrackParCovFwd mft{mfttrack.z(), tpars, tcovs, mfttrack.chi2()};
pft = o2::aod::fwdtrackutils::refitGlobalMuonCov(propMuon, mft);
} else {
pft = propagateToZero(track);
}
Expand Down Expand Up @@ -495,7 +508,33 @@ struct UpcCandProducerGlobalMuon {
return o2::track::TrackParCovFwd{track.z(), tpars, tcovs, track.chi2()};
}

// Select the best MCH-MFT match per MCH track based on lowest chi2MatchMCHMFT.
// Multiple global tracks can share the same MCH track with different MFT matches;
// this function keeps only the best one to reduce combinatorial background.
void selectBestMuonMatches(ForwardTracks const& fwdTracks)
{
fBestMuonMatch.clear();
std::unordered_map<int, std::pair<float, int>> mCandidates;
for (const auto& muon : fwdTracks) {
if (static_cast<int>(muon.trackType()) < kUpperBoundaryToTrackTypeEnum) {
auto muonID = muon.matchMCHTrackId();
auto chi2 = muon.chi2MatchMCHMFT();
if (mCandidates.find(muonID) == mCandidates.end()) {
mCandidates[muonID] = {chi2, muon.globalIndex()};
} else {
if (chi2 < mCandidates[muonID].first) {
mCandidates[muonID] = {chi2, muon.globalIndex()};
}
}
}
}
for (const auto& pairCand : mCandidates) {
fBestMuonMatch[pairCand.second.second] = true;
}
}

void createCandidates(ForwardTracks const& fwdTracks,
o2::aod::MFTTracks const& mftTracks,
o2::aod::FwdTrkCls const& fwdTrkCls,
o2::aod::AmbiguousFwdTracks const& ambFwdTracks,
o2::aod::BCs const& bcs,
Expand Down Expand Up @@ -580,6 +619,11 @@ struct UpcCandProducerGlobalMuon {
vAmbFwdTrackIndexBCs[ambTr.globalIndex()] = ambTr.bcIds()[0];
}

// Select best MCH-MFT match per MCH track before sorting into BC maps
if (fKeepBestMuonMatch) {
selectBestMuonMatches(fwdTracks);
}

std::map<uint64_t, std::vector<int64_t>> mapGlobalBcsWithMCHMIDTrackIds;
std::map<uint64_t, std::vector<int64_t>> mapGlobalBcsWithMCHTrackIds;
std::map<uint64_t, std::vector<int64_t>> mapGlobalBcsWithGlobalMuonTrackIds; // MCH-MID-MFT (good timing from MID)
Expand All @@ -595,6 +639,13 @@ struct UpcCandProducerGlobalMuon {
trackType != GlobalForwardTrack)
continue;

// For global tracks, skip if not the best match for this MCH track
if (fKeepBestMuonMatch && static_cast<int>(trackType) < kUpperBoundaryToTrackTypeEnum) {
if (fBestMuonMatch.find(fwdTrack.globalIndex()) == fBestMuonMatch.end()) {
continue;
}
}

auto trackId = fwdTrack.globalIndex();
int64_t indexBC = vAmbFwdTrackIndex[trackId] < 0 ? vColIndexBCs[fwdTrack.collisionId()] : vAmbFwdTrackIndexBCs[vAmbFwdTrackIndex[trackId]];
auto globalBC = vGlobalBCs[indexBC] + TMath::FloorNint(fwdTrack.trackTime() / o2::constants::lhc::LHCBunchSpacingNS + kBcTimeRoundingOffset);
Expand Down Expand Up @@ -675,7 +726,7 @@ struct UpcCandProducerGlobalMuon {
uint16_t numContrib = 0;
double sumPx = 0., sumPy = 0., sumPz = 0., sumE = 0.;
for (const auto& ianchor : vAnchorIds) {
if (!addToFwdTable(candId, ianchor, globalBcAnchor, 0., fwdTracks, mcFwdTrackLabels))
if (!addToFwdTable(candId, ianchor, globalBcAnchor, 0., fwdTracks, mftTracks, mcFwdTrackLabels))
continue;
const auto& trk = fwdTracks.iteratorAt(ianchor);
double p2 = trk.px() * trk.px() + trk.py() * trk.py() + trk.pz() * trk.pz();
Expand All @@ -696,7 +747,7 @@ struct UpcCandProducerGlobalMuon {

// Step 3: Write matched MCH-MFT tracks with adjusted track time
for (const auto& [imchMft, gbc] : mapMchMftIdBc) {
if (!addToFwdTable(candId, imchMft, gbc, (gbc - globalBcAnchor) * o2::constants::lhc::LHCBunchSpacingNS, fwdTracks, mcFwdTrackLabels))
if (!addToFwdTable(candId, imchMft, gbc, (gbc - globalBcAnchor) * o2::constants::lhc::LHCBunchSpacingNS, fwdTracks, mftTracks, mcFwdTrackLabels))
continue;
const auto& trk = fwdTracks.iteratorAt(imchMft);
double p2 = trk.px() * trk.px() + trk.py() * trk.py() + trk.pz() * trk.pz();
Expand Down Expand Up @@ -758,7 +809,7 @@ struct UpcCandProducerGlobalMuon {
auto& vMuonIds = gbc_muids.second;
// writing MCH-MID tracks
for (const auto& imuon : vMuonIds) {
if (!addToFwdTable(candId, imuon, globalBcMid, 0., fwdTracks, mcFwdTrackLabels))
if (!addToFwdTable(candId, imuon, globalBcMid, 0., fwdTracks, mftTracks, mcFwdTrackLabels))
continue;
numContrib++;
selTrackIds.push_back(imuon);
Expand All @@ -769,7 +820,7 @@ struct UpcCandProducerGlobalMuon {
getMchTrackIds(globalBcMid, mapGlobalBcsWithMCHTrackIds, fBcWindowMCH, mapMchIdBc);
// writing MCH-only tracks
for (const auto& [imch, gbc] : mapMchIdBc) {
if (!addToFwdTable(candId, imch, gbc, (gbc - globalBcMid) * o2::constants::lhc::LHCBunchSpacingNS, fwdTracks, mcFwdTrackLabels))
if (!addToFwdTable(candId, imch, gbc, (gbc - globalBcMid) * o2::constants::lhc::LHCBunchSpacingNS, fwdTracks, mftTracks, mcFwdTrackLabels))
continue;
numContrib++;
selTrackIds.push_back(imch);
Expand Down Expand Up @@ -816,6 +867,7 @@ struct UpcCandProducerGlobalMuon {
}

void processFwd(ForwardTracks const& fwdTracks,
o2::aod::MFTTracks const& mftTracks,
o2::aod::FwdTrkCls const& fwdTrkCls,
o2::aod::AmbiguousFwdTracks const& ambFwdTracks,
o2::aod::BCs const& bcs,
Expand All @@ -824,10 +876,11 @@ struct UpcCandProducerGlobalMuon {
o2::aod::Zdcs const& zdcs)
{
fDoMC = false;
createCandidates(fwdTracks, fwdTrkCls, ambFwdTracks, bcs, collisions, fv0s, zdcs, (o2::aod::McFwdTrackLabels*)nullptr);
createCandidates(fwdTracks, mftTracks, fwdTrkCls, ambFwdTracks, bcs, collisions, fv0s, zdcs, (o2::aod::McFwdTrackLabels*)nullptr);
}

void processFwdMC(ForwardTracks const& fwdTracks,
o2::aod::MFTTracks const& mftTracks,
o2::aod::FwdTrkCls const& fwdTrkCls,
o2::aod::AmbiguousFwdTracks const& ambFwdTracks,
o2::aod::BCs const& bcs,
Expand All @@ -840,7 +893,7 @@ struct UpcCandProducerGlobalMuon {
{
fDoMC = true;
skimMCInfo(mcCollisions, mcParticles);
createCandidates(fwdTracks, fwdTrkCls, ambFwdTracks, bcs, collisions, fv0s, zdcs, &mcFwdTrackLabels);
createCandidates(fwdTracks, mftTracks, fwdTrkCls, ambFwdTracks, bcs, collisions, fv0s, zdcs, &mcFwdTrackLabels);
fNewPartIDs.clear();
}

Expand Down
Loading