Skip to content

[PWGUD] Modification to global muon producer #15697

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 4 commits into from
Apr 9, 2026
+69 −52
Merged

[PWGUD] Modification to global muon producer #15697

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a2ddb40
all histos for fwd tracks only too
rolavick Apr 9, 2026
e7f1c5c
fixing performance issues
rolavick Apr 9, 2026
d3121e9
Please consider the following formatting changes
alibuild Apr 9, 2026
398b0b1
Merge pull request #75 from alibuild/alibot-cleanup-15697
rolavick Apr 9, 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
121 changes: 69 additions & 52 deletions PWGUD/TableProducer/upcCandProducerGlobalMuon.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,12 @@ struct UpcCandProducerGlobalMuon {
static constexpr double fCcenterMFT[3] = {0, 0, -61.4}; // Field evaluation point at center of MFT
float fZShift{0}; // z-vertex shift for forward track propagation

// Named constants (avoid magic numbers in expressions)
static constexpr double kInvalidDCA = 999.; // Sentinel for "no valid DCA computed yet"
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

void init(InitContext&)
{
fUpcCuts = (UPCCutparHolder)fUpcCutsConf;
Expand All @@ -132,30 +138,28 @@ struct UpcCandProducerGlobalMuon {
histRegistry.get<TH1>(HIST("MuonsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kFwdSelChi2 + 1, "Chi2");

// NEW: Add histograms for global track monitoring
if (fEnableMFT) {
const AxisSpec axisTrackType{5, -0.5, 4.5, "Track Type"};
histRegistry.add("hTrackTypes", "Track type distribution", kTH1F, {axisTrackType});
histRegistry.get<TH1>(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(1, "MuonStandalone");
histRegistry.get<TH1>(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(2, "MCHStandalone");
histRegistry.get<TH1>(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(3, "GlobalMuon");
histRegistry.get<TH1>(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(4, "GlobalFwd");

const AxisSpec axisEta{100, -4.0, -2.0, "#eta"};
histRegistry.add("hEtaGlobal", "Global track eta", kTH1F, {axisEta});

const AxisSpec axisDCAxy{200, 0., 10., "DCA_{xy} (cm)"};
const AxisSpec axisDCAz{200, -10., 10., "DCA_{z} (cm)"};
histRegistry.add("hDCAxyGlobal", "DCAxy of global tracks to best collision", kTH1F, {axisDCAxy});
histRegistry.add("hDCAzGlobal", "DCAz of global tracks to best collision", kTH1F, {axisDCAz});
histRegistry.add("hNCompatColls", "Number of compatible collisions per global track", kTH1F, {{21, -0.5, 20.5}});

const AxisSpec axisChi2Match{200, 0., 100., "#chi^{2}_{MCH-MFT}"};
histRegistry.add("hChi2MatchMCHMFT", "Chi2 of MCH-MFT matching (before cut)", kTH1F, {axisChi2Match});

const AxisSpec axisMass{500, 0., 10., "m_{inv} (GeV/c^{2})"};
histRegistry.add("hMassGlobalMuon", "Invariant mass from MCH-MID-MFT tracks only", kTH1F, {axisMass});
histRegistry.add("hMassGlobalMuonWithMCHMFT", "Invariant mass from MCH-MID-MFT + MCH-MFT tracks", kTH1F, {axisMass});
}
const AxisSpec axisTrackType{5, -0.5, 4.5, "Track Type"};
histRegistry.add("hTrackTypes", "Track type distribution", kTH1F, {axisTrackType});
histRegistry.get<TH1>(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(1, "MuonStandalone");
histRegistry.get<TH1>(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(2, "MCHStandalone");
histRegistry.get<TH1>(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(3, "GlobalMuon");
histRegistry.get<TH1>(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(4, "GlobalFwd");

const AxisSpec axisEta{100, -4.0, -2.0, "#eta"};
histRegistry.add("hEtaGlobal", "Global track eta", kTH1F, {axisEta});

const AxisSpec axisDCAxy{200, 0., 10., "DCA_{xy} (cm)"};
const AxisSpec axisDCAz{200, -10., 10., "DCA_{z} (cm)"};
histRegistry.add("hDCAxyGlobal", "DCAxy of global tracks to best collision", kTH1F, {axisDCAxy});
histRegistry.add("hDCAzGlobal", "DCAz of global tracks to best collision", kTH1F, {axisDCAz});
histRegistry.add("hNCompatColls", "Number of compatible collisions per global track", kTH1F, {{21, -0.5, 20.5}});

const AxisSpec axisChi2Match{200, 0., 100., "#chi^{2}_{MCH-MFT}"};
histRegistry.add("hChi2MatchMCHMFT", "Chi2 of MCH-MFT matching (before cut)", kTH1F, {axisChi2Match});

const AxisSpec axisMass{500, 0., 10., "m_{inv} (GeV/c^{2})"};
histRegistry.add("hMassGlobalMuon", "Invariant mass from MCH-MID-MFT tracks only", kTH1F, {axisMass});
histRegistry.add("hMassGlobalMuonWithMCHMFT", "Invariant mass from MCH-MID-MFT + MCH-MFT tracks", kTH1F, {axisMass});
}

bool cut(const o2::dataformats::GlobalFwdTrack& pft, const ForwardTracks::iterator& fwdTrack)
Expand Down Expand Up @@ -393,13 +397,11 @@ struct UpcCandProducerGlobalMuon {
float px, py, pz;
int sign;

// NEW: Fill track type histogram if MFT enabled
if (fEnableMFT) {
histRegistry.fill(HIST("hTrackTypes"), track.trackType());
if (track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack ||
track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalForwardTrack) {
histRegistry.fill(HIST("hEtaGlobal"), track.eta());
}
// NEW: 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) {
histRegistry.fill(HIST("hEtaGlobal"), track.eta());
}

if (fUpcCuts.getUseFwdCuts()) {
Expand Down Expand Up @@ -476,11 +478,11 @@ struct UpcCandProducerGlobalMuon {
}

int newId = 0;
for (auto trackId : trackIds) {
for (const auto& trackId : trackIds) {
auto it = clustersPerTrack.find(trackId);
if (it != clustersPerTrack.end()) {
const auto& clusters = it->second;
for (auto clsId : clusters) {
for (const auto& clsId : clusters) {
const auto& clsInfo = fwdTrkCls.iteratorAt(clsId);
udFwdTrkClusters(newId, clsInfo.x(), clsInfo.y(), clsInfo.z(), clsInfo.clInfo());
}
Expand All @@ -502,6 +504,23 @@ struct UpcCandProducerGlobalMuon {
return {dcaXY, dcaOrig[2], dcaOrig[0], dcaOrig[1]};
}

// Dispatch DCA computation: use MFT helix propagation when fEnableMFT is true,
// otherwise fall back to MCH extrapolation to (0,0,0) via propagateToZero.
// Returns {DCAxy, DCAz, DCAx, DCAy}
std::array<double, 4> propagateFwdToDCA(ForwardTracks::iterator const& track,
double colX, double colY, double colZ)
{
if (fEnableMFT) {
return propagateGlobalToDCA(track, colX, colY, colZ);
}
auto pft = propagateToZero(track);
double dcaX = pft.getX() - colX;
double dcaY = pft.getY() - colY;
double dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY);
double dcaZ = pft.getZ() - colZ;
return {dcaXY, dcaZ, dcaX, dcaY};
}

void createCandidates(ForwardTracks const& fwdTracks,
o2::aod::FwdTrkCls const& fwdTrkCls,
o2::aod::AmbiguousFwdTracks const& ambFwdTracks,
Expand Down Expand Up @@ -610,18 +629,18 @@ struct UpcCandProducerGlobalMuon {

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 + 1.);
auto globalBC = vGlobalBCs[indexBC] + TMath::FloorNint(fwdTrack.trackTime() / o2::constants::lhc::LHCBunchSpacingNS + kBcTimeRoundingOffset);

if (trackType == MuonStandaloneTrack) { // MCH-MID
mapGlobalBcsWithMCHMIDTrackIds[globalBC].push_back(trackId);
} else if (trackType == MCHStandaloneTrack) { // MCH-only
mapGlobalBcsWithMCHTrackIds[globalBC].push_back(trackId);
} else if (fEnableMFT && trackType == GlobalMuonTrack) { // MCH-MID-MFT: good timing, used as anchor
} else if (trackType == GlobalMuonTrack) { // MCH-MID-MFT: good timing, used as anchor
histRegistry.fill(HIST("hChi2MatchMCHMFT"), fwdTrack.chi2MatchMCHMFT());
if (fwdTrack.chi2MatchMCHMFT() > 0 && fwdTrack.chi2MatchMCHMFT() < fMaxChi2MatchMCHMFT) {
mapGlobalBcsWithGlobalMuonTrackIds[globalBC].push_back(trackId);
}
} else if (fEnableMFT && trackType == GlobalForwardTrack) { // MCH-MFT: poor timing, matched to anchors
} else if (trackType == GlobalForwardTrack) { // MCH-MFT: poor timing, matched to anchors
histRegistry.fill(HIST("hChi2MatchMCHMFT"), fwdTrack.chi2MatchMCHMFT());
if (fwdTrack.chi2MatchMCHMFT() > 0 && fwdTrack.chi2MatchMCHMFT() < fMaxChi2MatchMCHMFT) {
mapGlobalBcsWithMCHMFTTrackIds[globalBC].push_back(trackId);
Expand All @@ -631,11 +650,9 @@ struct UpcCandProducerGlobalMuon {

// Map global BC to collisions for DCA-based vertex assignment
std::map<uint64_t, std::vector<int64_t>> mapGlobalBCtoCollisions;
if (fEnableMFT) {
for (const auto& col : collisions) {
uint64_t gbc = vGlobalBCs[col.bcId()];
mapGlobalBCtoCollisions[gbc].push_back(col.globalIndex());
}
for (const auto& col : collisions) {
uint64_t gbc = vGlobalBCs[col.bcId()];
mapGlobalBCtoCollisions[gbc].push_back(col.globalIndex());
}

std::vector<int64_t> selTrackIds{}; // NEW: For cluster saving
Expand All @@ -645,7 +662,7 @@ struct UpcCandProducerGlobalMuon {
// Process global tracks: MCH-MID-MFT anchors + MCH-MFT in BC window
// MCH-MID-MFT tracks have good timing (from MID) and serve as anchors.
// MCH-MFT tracks have poor timing and are searched in a BC window around anchors.
if (fEnableMFT && !mapGlobalBcsWithGlobalMuonTrackIds.empty()) {
if (!mapGlobalBcsWithGlobalMuonTrackIds.empty()) {
for (const auto& gbc_anchorids : mapGlobalBcsWithGlobalMuonTrackIds) {
uint64_t globalBcAnchor = gbc_anchorids.first;
auto itFv0Id = mapGlobalBcWithV0A.find(globalBcAnchor);
Expand Down Expand Up @@ -675,7 +692,7 @@ struct UpcCandProducerGlobalMuon {

// Step 1: Find best collision vertex using DCA-based propagation
float bestVtxX = 0., bestVtxY = 0., bestVtxZ = 0.;
double bestAvgDCA = 999.;
double bestAvgDCA = kInvalidDCA;
bool hasVertex = false;
int nCompatColls = 0;

Expand All @@ -684,18 +701,18 @@ struct UpcCandProducerGlobalMuon {
auto itCol = mapGlobalBCtoCollisions.find(searchBC);
if (itCol == mapGlobalBCtoCollisions.end())
continue;
for (auto colIdx : itCol->second) {
for (const auto& colIdx : itCol->second) {
nCompatColls++;
const auto& col = collisions.iteratorAt(colIdx);
double sumDCAxy = 0.;
int nTracks = 0;
for (const auto& iglobal : allTrackIds) {
const auto& trk = fwdTracks.iteratorAt(iglobal);
auto dca = propagateGlobalToDCA(trk, col.posX(), col.posY(), col.posZ());
auto dca = propagateFwdToDCA(trk, col.posX(), col.posY(), col.posZ());
sumDCAxy += dca[0];
nTracks++;
}
double avgDCA = nTracks > 0 ? sumDCAxy / nTracks : 999.;
double avgDCA = nTracks > 0 ? sumDCAxy / nTracks : kInvalidDCA;
if (!hasVertex || avgDCA < bestAvgDCA) {
bestAvgDCA = avgDCA;
bestVtxX = col.posX();
Expand All @@ -715,7 +732,7 @@ struct UpcCandProducerGlobalMuon {
for (const auto& ianchor : vAnchorIds) {
if (hasVertex) {
const auto& trk = fwdTracks.iteratorAt(ianchor);
auto dca = propagateGlobalToDCA(trk, bestVtxX, bestVtxY, bestVtxZ);
auto dca = propagateFwdToDCA(trk, bestVtxX, bestVtxY, bestVtxZ);
histRegistry.fill(HIST("hDCAxyGlobal"), dca[0]);
histRegistry.fill(HIST("hDCAzGlobal"), dca[1]);
if (dca[0] > static_cast<double>(fMaxDCAxy))
Expand All @@ -735,7 +752,7 @@ struct UpcCandProducerGlobalMuon {

// Fill invariant mass from MCH-MID-MFT anchors only
uint16_t numContribAnch = numContrib;
if (numContribAnch >= 2) {
if (numContribAnch >= kMinTracksForPair) {
double mass2 = sumE * sumE - sumPx * sumPx - sumPy * sumPy - sumPz * sumPz;
histRegistry.fill(HIST("hMassGlobalMuon"), mass2 > 0. ? std::sqrt(mass2) : 0.);
}
Expand All @@ -744,7 +761,7 @@ struct UpcCandProducerGlobalMuon {
for (const auto& [imchMft, gbc] : mapMchMftIdBc) {
if (hasVertex) {
const auto& trk = fwdTracks.iteratorAt(imchMft);
auto dca = propagateGlobalToDCA(trk, bestVtxX, bestVtxY, bestVtxZ);
auto dca = propagateFwdToDCA(trk, bestVtxX, bestVtxY, bestVtxZ);
histRegistry.fill(HIST("hDCAxyGlobal"), dca[0]);
histRegistry.fill(HIST("hDCAzGlobal"), dca[1]);
if (dca[0] > static_cast<double>(fMaxDCAxy))
Expand All @@ -763,12 +780,12 @@ struct UpcCandProducerGlobalMuon {
}

// Fill invariant mass including MCH-MFT tracks (only if MCH-MFT tracks were added)
if (numContrib > numContribAnch && numContrib >= 2) {
if (numContrib > numContribAnch && numContrib >= kMinTracksForPair) {
double mass2 = sumE * sumE - sumPx * sumPx - sumPy * sumPy - sumPz * sumPz;
histRegistry.fill(HIST("hMassGlobalMuonWithMCHMFT"), mass2 > 0. ? std::sqrt(mass2) : 0.);
}

if (numContrib < 1)
if (numContrib < kMinTracksForCandidate)
continue;

eventCandidates(globalBcAnchor, runNumber, bestVtxX, bestVtxY, bestVtxZ, 0, numContrib, 0, 0);
Expand Down Expand Up @@ -817,7 +834,7 @@ struct UpcCandProducerGlobalMuon {
numContrib++;
selTrackIds.push_back(imuon);
}
if (numContrib < 1) // didn't save any MCH-MID tracks
if (numContrib < kMinTracksForCandidate) // didn't save any MCH-MID tracks
continue;
std::map<int64_t, uint64_t> mapMchIdBc{};
getMchTrackIds(globalBcMid, mapGlobalBcsWithMCHTrackIds, fBcWindowMCH, mapMchIdBc);
Expand Down
Loading