diff --git a/PWGUD/TableProducer/upcCandProducerGlobalMuon.cxx b/PWGUD/TableProducer/upcCandProducerGlobalMuon.cxx index 744b4e8523a..42aa8108dd1 100644 --- a/PWGUD/TableProducer/upcCandProducerGlobalMuon.cxx +++ b/PWGUD/TableProducer/upcCandProducerGlobalMuon.cxx @@ -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; @@ -132,30 +138,28 @@ struct UpcCandProducerGlobalMuon { histRegistry.get(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(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(1, "MuonStandalone"); - histRegistry.get(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(2, "MCHStandalone"); - histRegistry.get(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(3, "GlobalMuon"); - histRegistry.get(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(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(1, "MuonStandalone"); + histRegistry.get(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(2, "MCHStandalone"); + histRegistry.get(HIST("hTrackTypes"))->GetXaxis()->SetBinLabel(3, "GlobalMuon"); + histRegistry.get(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) @@ -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()) { @@ -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()); } @@ -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 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, @@ -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); @@ -631,11 +650,9 @@ struct UpcCandProducerGlobalMuon { // Map global BC to collisions for DCA-based vertex assignment std::map> 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 selTrackIds{}; // NEW: For cluster saving @@ -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); @@ -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; @@ -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(); @@ -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(fMaxDCAxy)) @@ -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.); } @@ -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(fMaxDCAxy)) @@ -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); @@ -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 mapMchIdBc{}; getMchTrackIds(globalBcMid, mapGlobalBcsWithMCHTrackIds, fBcWindowMCH, mapMchIdBc);