diff --git a/PWGUD/TableProducer/CMakeLists.txt b/PWGUD/TableProducer/CMakeLists.txt index 284b1af9fdc..a69671ab068 100644 --- a/PWGUD/TableProducer/CMakeLists.txt +++ b/PWGUD/TableProducer/CMakeLists.txt @@ -75,3 +75,8 @@ o2physics_add_dpl_workflow(upc-cand-producer-global-muon SOURCES upcCandProducerGlobalMuon.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::UPCCutparHolder O2::GlobalTracking COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(upc-cand-producer-semi-fwd + SOURCES upcCandProducerSemiFwd.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::UPCCutparHolder O2::GlobalTracking + COMPONENT_NAME Analysis) diff --git a/PWGUD/TableProducer/upcCandProducerSemiFwd.cxx b/PWGUD/TableProducer/upcCandProducerSemiFwd.cxx new file mode 100644 index 00000000000..958c2e8e796 --- /dev/null +++ b/PWGUD/TableProducer/upcCandProducerSemiFwd.cxx @@ -0,0 +1,681 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file upcCandProducerSemiFwd.cxx +/// \brief UPC candidate producer combining forward MCH-MID muons with central-barrel tracks. +/// Anchor: MCH-MID forward track (precise MID timing). Companions: barrel tracks +/// whose BC falls within fBcWindowBarrel of the anchor BC. One UD candidate per anchor BC. +/// \author Roman Lavička, lavicka.roman@gmail.com +/// \since 10.05.2026 + +#include "PWGUD/Core/UPCCutparHolder.h" +#include "PWGUD/Core/UPCHelpers.h" +#include "PWGUD/DataModel/UDTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct UpcCandProducerSemiFwd { + bool fDoMC{false}; + + std::map fNewPartIDs; + + Produces udMCCollisions; + Produces udMCParticles; + Produces udFwdTrackLabels; + Produces udTrackLabels; + + Produces udFwdTracks; + Produces udFwdTracksExtra; + Produces udFwdIndices; + + Produces udTracks; + Produces udTracksExtra; + Produces udTracksDCA; + Produces udTracksFlags; + Produces udTracksPID; + + Produces eventCandidates; + Produces eventCandidatesSelsFwd; + Produces udZdcsReduced; + + Configurable fSignalGenID{"fSignalGenID", 1, "Signal generator ID"}; + + UPCCutparHolder fUpcCuts = UPCCutparHolder(); + MutableConfigurable fUpcCutsConf{"fUpcCutsConf", {}, "UPC event cuts"}; + Configurable fBcWindowFITAmps{"fBcWindowFITAmps", 20, "BC range for T0A/V0A amplitudes array [-range, +(range-1)]"}; + Configurable fBcWindowBarrel{"fBcWindowBarrel", 20, "Time window for barrel-track to MCH-MID anchor matching"}; + Configurable fMaxFV0Amp{"fMaxFV0Amp", 100.f, "Max FV0 amplitude in the same BC"}; + + using ForwardTracks = o2::soa::Join; + + using BarrelTracks = o2::soa::Join; + + HistogramRegistry histRegistry{"HistRegistry", {}, OutputObjHandlingPolicy::AnalysisObject}; + + int fRun{0}; + Service fCCDB; + o2::ccdb::CcdbApi fCCDBApi; + o2::globaltracking::MatchGlobalFwd fMatching; + + void init(InitContext&) + { + fUpcCuts = (UPCCutparHolder)fUpcCutsConf; + if (fUpcCuts.getUseFwdCuts()) { + fCCDB->setURL("http://alice-ccdb.cern.ch"); + fCCDB->setCaching(true); + fCCDB->setLocalObjectValidityChecking(); + fCCDBApi.init("http://alice-ccdb.cern.ch"); + } + const AxisSpec axisSelFwd{upchelpers::kNFwdSels, 0., static_cast(upchelpers::kNFwdSels), ""}; + histRegistry.add("MuonsSelCounter", "", kTH1F, {axisSelFwd}); + histRegistry.get(HIST("MuonsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kFwdSelAll + 1, "All"); + histRegistry.get(HIST("MuonsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kFwdSelPt + 1, "Pt"); + histRegistry.get(HIST("MuonsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kFwdSelEta + 1, "Eta"); + histRegistry.get(HIST("MuonsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kFwdSelRabs + 1, "Rabs"); + histRegistry.get(HIST("MuonsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kFwdSelpDCA + 1, "pDCA"); + histRegistry.get(HIST("MuonsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kFwdSelChi2 + 1, "Chi2"); + + const AxisSpec axisSelBar{upchelpers::kNBarrelSels, 0., static_cast(upchelpers::kNBarrelSels), ""}; + histRegistry.add("BarrelsSelCounter", "", kTH1F, {axisSelBar}); + histRegistry.get(HIST("BarrelsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kBarrelSelAll + 1, "All"); + histRegistry.get(HIST("BarrelsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kBarrelSelHasTOF + 1, "HasTOF"); + histRegistry.get(HIST("BarrelsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kBarrelSelPt + 1, "Pt"); + histRegistry.get(HIST("BarrelsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kBarrelSelEta + 1, "Eta"); + histRegistry.get(HIST("BarrelsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kBarrelSelITSNCls + 1, "ITSNCls"); + histRegistry.get(HIST("BarrelsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kBarrelSelITSChi2 + 1, "ITSChi2"); + histRegistry.get(HIST("BarrelsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kBarrelSelTPCNCls + 1, "TPCNCls"); + histRegistry.get(HIST("BarrelsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kBarrelSelTPCChi2 + 1, "TPCChi2"); + histRegistry.get(HIST("BarrelsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kBarrelSelDCAXY + 1, "DCAXY"); + histRegistry.get(HIST("BarrelsSelCounter"))->GetXaxis()->SetBinLabel(upchelpers::kBarrelSelDCAZ + 1, "DCAZ"); + } + + bool fwdCut(const o2::dataformats::GlobalFwdTrack& pft, const ForwardTracks::iterator& fwdTrack) + { + histRegistry.fill(HIST("MuonsSelCounter"), upchelpers::kFwdSelAll, 1); + auto pt = pft.getPt(); + auto eta = pft.getEta(); + auto pdca = fwdTrack.pDca(); + auto rabs = fwdTrack.rAtAbsorberEnd(); + auto chi2 = fwdTrack.chi2(); + bool passPt = pt > fUpcCuts.getFwdPtLow() && pt < fUpcCuts.getFwdPtHigh(); + bool passEta = eta > fUpcCuts.getFwdEtaLow() && eta < fUpcCuts.getFwdEtaHigh(); + bool passRabs = rabs > fUpcCuts.getMuonRAtAbsorberEndLow() && rabs < fUpcCuts.getMuonRAtAbsorberEndHigh(); + bool passPDca = rabs < upchelpers::AbsorberMid ? pdca < fUpcCuts.getMuonPDcaHighFirst() : pdca < fUpcCuts.getMuonPDcaHighSecond(); + bool passChi2 = chi2 > fUpcCuts.getFwdChi2Low() && chi2 < fUpcCuts.getFwdChi2High(); + if (passPt) + histRegistry.fill(HIST("MuonsSelCounter"), upchelpers::kFwdSelPt, 1); + if (passEta) + histRegistry.fill(HIST("MuonsSelCounter"), upchelpers::kFwdSelEta, 1); + if (passRabs) + histRegistry.fill(HIST("MuonsSelCounter"), upchelpers::kFwdSelRabs, 1); + if (passPDca) + histRegistry.fill(HIST("MuonsSelCounter"), upchelpers::kFwdSelpDCA, 1); + if (passChi2) + histRegistry.fill(HIST("MuonsSelCounter"), upchelpers::kFwdSelChi2, 1); + return passPt && passEta && passRabs && passPDca && passChi2; + } + + bool barCut(const BarrelTracks::iterator& track) + { + histRegistry.fill(HIST("BarrelsSelCounter"), upchelpers::kBarrelSelAll, 1); + if (!fUpcCuts.getUseBarCuts()) + return true; + if (fUpcCuts.getAmbigSwitch() == 1 && !track.isPVContributor()) + return false; + if (fUpcCuts.getRequireTOF() && !track.hasTOF()) + return false; + bool passHasTOF = !fUpcCuts.getRequireTOF() || track.hasTOF(); + bool passPt = track.pt() > fUpcCuts.getBarPtLow() && track.pt() < fUpcCuts.getBarPtHigh(); + bool passEta = track.eta() > fUpcCuts.getBarEtaLow() && track.eta() < fUpcCuts.getBarEtaHigh(); + bool passITSNCls = track.itsNCls() >= static_cast(fUpcCuts.getITSNClusLow()) && + track.itsNCls() <= static_cast(fUpcCuts.getITSNClusHigh()); + bool passITSChi2 = track.itsChi2NCl() > fUpcCuts.getITSChi2Low() && track.itsChi2NCl() < fUpcCuts.getITSChi2High(); + bool passTPCNCls = track.tpcNClsFound() >= static_cast(fUpcCuts.getTPCNClsLow()) && + track.tpcNClsFound() <= static_cast(fUpcCuts.getTPCNClsHigh()); + bool passTPCChi2 = track.tpcChi2NCl() > fUpcCuts.getTPCChi2Low() && track.tpcChi2NCl() < fUpcCuts.getTPCChi2High(); + bool passDcaZ = track.dcaZ() > fUpcCuts.getDcaZLow() && track.dcaZ() < fUpcCuts.getDcaZHigh(); + bool passDcaXY = true; + if (fUpcCuts.getCheckMaxDcaXY()) { + constexpr float DcaXYConst = 0.0105f; + constexpr float DcaXYSlope = 0.0350f; + constexpr float DcaXYExp = 1.1f; + float maxDCA = DcaXYConst + DcaXYSlope / std::pow(track.pt(), DcaXYExp); + passDcaXY = track.dcaXY() < maxDCA; + } + if (passHasTOF) + histRegistry.fill(HIST("BarrelsSelCounter"), upchelpers::kBarrelSelHasTOF, 1); + if (passPt) + histRegistry.fill(HIST("BarrelsSelCounter"), upchelpers::kBarrelSelPt, 1); + if (passEta) + histRegistry.fill(HIST("BarrelsSelCounter"), upchelpers::kBarrelSelEta, 1); + if (passITSNCls) + histRegistry.fill(HIST("BarrelsSelCounter"), upchelpers::kBarrelSelITSNCls, 1); + if (passITSChi2) + histRegistry.fill(HIST("BarrelsSelCounter"), upchelpers::kBarrelSelITSChi2, 1); + if (passTPCNCls) + histRegistry.fill(HIST("BarrelsSelCounter"), upchelpers::kBarrelSelTPCNCls, 1); + if (passTPCChi2) + histRegistry.fill(HIST("BarrelsSelCounter"), upchelpers::kBarrelSelTPCChi2, 1); + if (passDcaXY) + histRegistry.fill(HIST("BarrelsSelCounter"), upchelpers::kBarrelSelDCAXY, 1); + if (passDcaZ) + histRegistry.fill(HIST("BarrelsSelCounter"), upchelpers::kBarrelSelDCAZ, 1); + return passPt && passEta && passITSNCls && passITSChi2 && passTPCNCls && passTPCChi2 && passDcaXY && passDcaZ; + } + + void skimMCInfo(o2::aod::McCollisions const& mcCollisions, o2::aod::McParticles const& mcParticles) + { + std::vector newEventIDs(mcCollisions.size(), -1); + + int32_t newPartID = 0; + int32_t newEventID = 0; + int32_t nMCParticles = mcParticles.size(); + for (int32_t mcPartID = 0; mcPartID < nMCParticles; mcPartID++) { + const auto& mcPart = mcParticles.iteratorAt(mcPartID); + if (!mcPart.has_mcCollision()) + continue; + int32_t mcEventID = mcPart.mcCollisionId(); + const auto& mcEvent = mcCollisions.iteratorAt(mcEventID); + bool isSignal = mcEvent.getGeneratorId() == fSignalGenID; + if (!isSignal) { + continue; + } + fNewPartIDs[mcPartID] = newPartID; + newPartID++; + if (newEventIDs[mcEventID] == -1) { + newEventIDs[mcEventID] = newEventID; + newEventID++; + } + } + + std::vector newMotherIDs{}; + + for (const auto& item : fNewPartIDs) { + int32_t mcPartID = item.first; + const auto& mcPart = mcParticles.iteratorAt(mcPartID); + int32_t mcEventID = mcPart.mcCollisionId(); + int32_t mappedEventID = newEventIDs[mcEventID]; + if (mcPart.has_mothers()) { + const auto& motherIDs = mcPart.mothersIds(); + for (const auto& motherID : motherIDs) { + if (motherID >= nMCParticles) { + continue; + } + auto it = fNewPartIDs.find(motherID); + if (it != fNewPartIDs.end()) { + newMotherIDs.push_back(it->second); + } + } + } + int32_t newDaughterIDs[2] = {-1, -1}; + if (mcPart.has_daughters()) { + const auto& daughterIDs = mcPart.daughtersIds(); + int32_t firstDaughter = daughterIDs.front(); + int32_t lastDaughter = daughterIDs.back(); + if (firstDaughter >= nMCParticles || lastDaughter >= nMCParticles) { + continue; + } + auto itFirst = fNewPartIDs.find(firstDaughter); + auto itLast = fNewPartIDs.find(lastDaughter); + if (itFirst != fNewPartIDs.end() && itLast != fNewPartIDs.end()) { + newDaughterIDs[0] = fNewPartIDs.at(daughterIDs.front()); + newDaughterIDs[1] = fNewPartIDs.at(daughterIDs.back()); + } + } + udMCParticles(mappedEventID, mcPart.pdgCode(), mcPart.getHepMCStatusCode(), mcPart.flags(), newMotherIDs, newDaughterIDs, + mcPart.weight(), mcPart.px(), mcPart.py(), mcPart.pz(), mcPart.e()); + newMotherIDs.clear(); + } + + for (int32_t i = 0; i < mcCollisions.size(); i++) { + if (newEventIDs[i] == -1) { + continue; + } + const auto& mcEvent = mcCollisions.iteratorAt(i); + const auto& bc = mcEvent.bc(); + udMCCollisions(bc.globalBC(), mcEvent.generatorsID(), mcEvent.posX(), mcEvent.posY(), mcEvent.posZ(), + mcEvent.t(), mcEvent.weight(), mcEvent.impactParameter()); + } + + newEventIDs.clear(); + } + + int64_t bcDiff(uint64_t bc1, uint64_t bc2) + { + return bc1 > bc2 ? bc1 - bc2 : bc2 - bc1; + } + + template + T::iterator getStartForScroll(uint64_t inGbc, T& gbcMap) + { + auto it1 = gbcMap.lower_bound(inGbc); + typename T::iterator it; + if (it1 != gbcMap.end()) { + auto it2 = it1; + uint64_t bc1 = it1->first; + if (it2 != gbcMap.begin()) + --it2; + uint64_t bc2 = it2->first; + uint64_t dbc1 = bcDiff(bc1, inGbc); + uint64_t dbc2 = bcDiff(bc2, inGbc); + it = (dbc1 <= dbc2) ? it1 : it2; + } else { + it = it1; + --it; + } + return it; + } + + template + void scrollBackForth(uint64_t inGbc, uint64_t maxDbc, T& gbcMap, F&& func) + { + auto it = getStartForScroll(inGbc, gbcMap); + uint64_t gbc = it->first; + uint64_t dbc = bcDiff(inGbc, gbc); + + int count = 0; + while (dbc <= maxDbc) { + func(it, gbc); + count++; + if (it == gbcMap.begin()) + break; + --it; + gbc = it->first; + dbc = bcDiff(inGbc, gbc); + } + + std::advance(it, count + 1); + + if (it == gbcMap.end()) + return; + + gbc = it->first; + dbc = bcDiff(inGbc, gbc); + + while (dbc <= maxDbc) { + func(it, gbc); + ++it; + if (it == gbcMap.end()) + break; + gbc = it->first; + dbc = bcDiff(inGbc, gbc); + } + } + + void getBarrelTrackIds(uint64_t inGbc, std::map>& tracksPerBC, + uint64_t maxDbc, std::map& outTrkIds) + { + auto fillIds = [&outTrkIds](std::map>::iterator& inIt, uint64_t gbc) { + std::vector& ids = inIt->second; + for (const auto& id : ids) + outTrkIds[id] = gbc; + }; + scrollBackForth(inGbc, maxDbc, tracksPerBC, fillIds); + } + + void getFV0Amplitudes(uint64_t inGbc, o2::aod::FV0As const& fv0s, uint64_t maxDbc, + std::map& mapBcs, std::vector& amps, std::vector& relBcs) + { + auto fillAmps = [this, &fv0s, &s, &relBcs, inGbc](std::map::iterator& inIt, uint64_t gbc) { + int64_t fv0Id = inIt->second; + const auto& fv0 = fv0s.iteratorAt(fv0Id); + const auto& amplitudes = fv0.amplitude(); + float totalAmp = std::accumulate(amplitudes.begin(), amplitudes.end(), 0.f); + if (totalAmp > 0.f) { + amps.push_back(totalAmp); + auto relBc = static_cast(bcDiff(gbc, inGbc)); + if (gbc < inGbc) + relBc *= -1; + relBcs.push_back(relBc); + } + }; + scrollBackForth(inGbc, maxDbc, mapBcs, fillAmps); + } + + auto propagateToZero(ForwardTracks::iterator const& muon) + { + using SMatrix55 = ROOT::Math::SMatrix>; + using SMatrix5 = ROOT::Math::SVector; + SMatrix5 tpars(muon.x(), muon.y(), muon.phi(), muon.tgl(), muon.signed1Pt()); + std::vector v1{muon.cXX(), muon.cXY(), muon.cYY(), muon.cPhiX(), muon.cPhiY(), + muon.cPhiPhi(), muon.cTglX(), muon.cTglY(), muon.cTglPhi(), muon.cTglTgl(), + muon.c1PtX(), muon.c1PtY(), muon.c1PtPhi(), muon.c1PtTgl(), muon.c1Pt21Pt2()}; + SMatrix55 tcovs(v1.begin(), v1.end()); + o2::dataformats::GlobalFwdTrack propmuon; + o2::dataformats::GlobalFwdTrack track; + track.setParameters(tpars); + track.setZ(muon.z()); + track.setCovariances(tcovs); + auto mchTrack = fMatching.FwdtoMCH(track); + o2::mch::TrackExtrap::extrapToVertex(mchTrack, 0., 0., 0., 0., 0.); + auto proptrack = fMatching.MCHtoFwd(mchTrack); + propmuon.setParameters(proptrack.getParameters()); + propmuon.setZ(proptrack.getZ()); + propmuon.setCovariances(proptrack.getCovariances()); + return propmuon; + } + + bool addToFwdTable(int64_t candId, int64_t trackId, uint64_t gbc, float trackTime, ForwardTracks const& fwdTracks, + const o2::aod::McFwdTrackLabels* mcFwdTrackLabels) + { + constexpr float NoMchMftMatch = -1.f; + const auto& track = fwdTracks.iteratorAt(trackId); + float px; + float py; + float pz; + int sign; + if (fUpcCuts.getUseFwdCuts()) { + auto pft = propagateToZero(track); + bool pass = fwdCut(pft, track); + if (!pass) + return false; + px = pft.getPx(); + py = pft.getPy(); + pz = pft.getPz(); + sign = (pft.getInvQPt() > 0) ? 1 : -1; + } else { + px = track.px(); + py = track.py(); + pz = track.pz(); + sign = track.sign(); + } + udFwdTracks(candId, px, py, pz, sign, gbc, trackTime, track.trackTimeRes()); + udFwdTracksExtra(track.trackType(), track.nClusters(), track.pDca(), track.rAtAbsorberEnd(), track.chi2(), + track.chi2MatchMCHMID(), NoMchMftMatch, track.mchBitMap(), track.midBitMap(), track.midBoards()); + udFwdIndices(candId, trackId, track.matchMCHTrackId(), -1); + if (fDoMC) { + const auto& label = mcFwdTrackLabels->iteratorAt(trackId); + uint16_t mcMask = label.mcMask(); + auto it = fNewPartIDs.find(label.mcParticleId()); + int32_t newPartID = it != fNewPartIDs.end() ? it->second : -1; + udFwdTrackLabels(newPartID, mcMask); + } + return true; + } + + bool addToBarrelTable(int64_t candId, int64_t trackId, uint64_t gbc, float trackTime, BarrelTracks const& barrelTracks, + const o2::aod::McTrackLabels* mcTrackLabels) + { + const auto& track = barrelTracks.iteratorAt(trackId); + if (!barCut(track)) + return false; + int32_t colId = track.collisionId() >= 0 ? track.collisionId() : -1; + udTracks(candId, track.px(), track.py(), track.pz(), track.sign(), gbc, trackTime, track.trackTimeRes()); + udTracksExtra(track.tpcInnerParam(), track.itsClusterSizes(), track.tpcNClsFindable(), track.tpcNClsFindableMinusFound(), + track.tpcNClsFindableMinusCrossedRows(), track.tpcNClsShared(), track.trdPattern(), track.itsChi2NCl(), + track.tpcChi2NCl(), track.trdChi2(), track.tofChi2(), track.tpcSignal(), track.tofSignal(), track.trdSignal(), + track.length(), track.tofExpMom(), track.detectorMap()); + udTracksDCA(track.dcaZ(), track.dcaXY()); + udTracksFlags(colId, track.isPVContributor()); + udTracksPID(track.tpcNSigmaEl(), track.tpcNSigmaMu(), track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr(), + track.beta(), track.betaerror(), + track.tofNSigmaEl(), track.tofNSigmaMu(), track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()); + if (fDoMC) { + const auto& label = mcTrackLabels->iteratorAt(trackId); + uint16_t mcMask = label.mcMask(); + auto it = fNewPartIDs.find(label.mcParticleId()); + int32_t newPartID = it != fNewPartIDs.end() ? it->second : -1; + udTrackLabels(newPartID, mcMask); + } + return true; + } + + void createCandidates(ForwardTracks const& fwdTracks, + o2::aod::AmbiguousFwdTracks const& ambFwdTracks, + BarrelTracks const& barrelTracks, + o2::aod::AmbiguousTracks const& ambBarrelTracks, + o2::aod::BCs const& bcs, + o2::aod::Collisions const& collisions, + o2::aod::FV0As const& fv0s, + o2::aod::Zdcs const& zdcs, + const o2::aod::McFwdTrackLabels* mcFwdTrackLabels, + const o2::aod::McTrackLabels* mcBarrelTrackLabels) + { + using o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack; + + int32_t runNumber = bcs.iteratorAt(0).runNumber(); + if (fUpcCuts.getUseFwdCuts()) { + if (runNumber != fRun) { + fRun = runNumber; + std::map metadata; + auto soreor = o2::ccdb::BasicCCDBManager::getRunDuration(fCCDBApi, fRun); + auto ts = soreor.second; + auto grpmag = fCCDBApi.retrieveFromTFileAny("GLO/Config/GRPMagField", metadata, ts); + o2::base::Propagator::initFieldFromGRP(grpmag); + if (!o2::base::GeometryManager::isGeometryLoaded()) + fCCDB->get("GLO/Config/GeometryAligned"); + o2::mch::TrackExtrap::setField(); + } + } + + auto nBcs = bcs.size(); + std::vector vGlobalBCs(nBcs, 0); + for (const auto& bc : bcs) { + vGlobalBCs[bc.globalIndex()] = bc.globalBC(); + } + + auto nCols = collisions.size(); + std::vector vColIndexBCs(nCols, -1); + for (const auto& col : collisions) { + vColIndexBCs[col.globalIndex()] = col.bcId(); + } + + std::map mapGlobalBcWithV0A{}; + constexpr float FV0ValidTime = 15.f; + for (const auto& fv0 : fv0s) { + if (std::abs(fv0.time()) > FV0ValidTime) + continue; + uint64_t globalBC = vGlobalBCs[fv0.bcId()]; + mapGlobalBcWithV0A[globalBC] = fv0.globalIndex(); + } + auto nFV0s = mapGlobalBcWithV0A.size(); + + std::map mapGlobalBcWithZdc{}; + constexpr float ZDCValidTime = 2.f; + for (const auto& zdc : zdcs) { + if (std::abs(zdc.timeZNA()) > ZDCValidTime && std::abs(zdc.timeZNC()) > ZDCValidTime) + continue; + uint64_t globalBC = vGlobalBCs[zdc.bcId()]; + mapGlobalBcWithZdc[globalBC] = zdc.globalIndex(); + } + auto nZdcs = mapGlobalBcWithZdc.size(); + + auto nFwdTracks = fwdTracks.size(); + auto nAmbFwdTracks = ambFwdTracks.size(); + std::vector vAmbFwdTrackIndex(nFwdTracks, -1); + std::vector vAmbFwdTrackIndexBCs(nAmbFwdTracks, -1); + for (const auto& ambTr : ambFwdTracks) { + vAmbFwdTrackIndex[ambTr.fwdtrackId()] = ambTr.globalIndex(); + vAmbFwdTrackIndexBCs[ambTr.globalIndex()] = ambTr.bcIds()[0]; + } + + auto nBarrelTracks = barrelTracks.size(); + auto nAmbBarrelTracks = ambBarrelTracks.size(); + std::vector vAmbBarTrackIndex(nBarrelTracks, -1); + std::vector vAmbBarTrackIndexBCs(nAmbBarrelTracks, -1); + for (const auto& ambTr : ambBarrelTracks) { + vAmbBarTrackIndex[ambTr.trackId()] = ambTr.globalIndex(); + vAmbBarTrackIndexBCs[ambTr.globalIndex()] = ambTr.bcIds()[0]; + } + + // anchor map: MCH-MID forward muons keyed by their global BC + std::map> mapGlobalBcsWithMCHMIDTrackIds; + for (const auto& fwdTrack : fwdTracks) { + if (fwdTrack.trackType() != MuonStandaloneTrack) + 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 + 1.); + mapGlobalBcsWithMCHMIDTrackIds[globalBC].push_back(trackId); + } + + // companion map: barrel tracks keyed by their global BC + std::map> mapGlobalBcsWithBarrelTrackIds; + for (const auto& barTrack : barrelTracks) { + auto trackId = barTrack.globalIndex(); + int64_t bcId = vAmbBarTrackIndex[trackId] < 0 + ? (barTrack.has_collision() ? vColIndexBCs[barTrack.collisionId()] : -1) + : vAmbBarTrackIndexBCs[vAmbBarTrackIndex[trackId]]; + if (bcId < 0) + continue; + auto globalBC = vGlobalBCs[bcId] + TMath::FloorNint(barTrack.trackTime() / o2::constants::lhc::LHCBunchSpacingNS + 1.); + mapGlobalBcsWithBarrelTrackIds[globalBC].push_back(trackId); + } + + int32_t candId = 0; + for (const auto& gbcMuids : mapGlobalBcsWithMCHMIDTrackIds) { + uint64_t globalBcMid = gbcMuids.first; + + // veto on FV0 in the same BC as the anchor + auto itFv0Id = mapGlobalBcWithV0A.find(globalBcMid); + if (itFv0Id != mapGlobalBcWithV0A.end()) { + auto fv0Id = itFv0Id->second; + const auto& fv0 = fv0s.iteratorAt(fv0Id); + float fv0Amp = 0.f; + for (const auto& amp : fv0.amplitude()) + fv0Amp += amp; + if (fv0Amp > fMaxFV0Amp) + continue; + } + + uint16_t numContrib = 0; + auto& vMuonIds = gbcMuids.second; + for (const auto& imuon : vMuonIds) { + if (!addToFwdTable(candId, imuon, globalBcMid, 0., fwdTracks, mcFwdTrackLabels)) + continue; + numContrib++; + } + if (numContrib < 1) // anchor produced no surviving fwd tracks + continue; + + // gather barrel tracks within ±fBcWindowBarrel of the anchor + std::map mapBarIdBc{}; + getBarrelTrackIds(globalBcMid, mapGlobalBcsWithBarrelTrackIds, fBcWindowBarrel, mapBarIdBc); + for (const auto& [ibar, gbc] : mapBarIdBc) { + float trackTime = (static_cast(gbc) - static_cast(globalBcMid)) * o2::constants::lhc::LHCBunchSpacingNS; + if (!addToBarrelTable(candId, ibar, gbc, trackTime, barrelTracks, mcBarrelTrackLabels)) + continue; + numContrib++; + } + + eventCandidates(globalBcMid, runNumber, 0., 0., 0., 0, numContrib, 0, 0); + std::vector amplitudesV0A{}; + std::vector relBCsV0A{}; + std::vector amplitudesT0A{}; + std::vector relBCsT0A{}; + if (nFV0s > 0) { + getFV0Amplitudes(globalBcMid, fv0s, fBcWindowFITAmps, mapGlobalBcWithV0A, amplitudesV0A, relBCsV0A); + } + eventCandidatesSelsFwd(0., 0., amplitudesT0A, relBCsT0A, amplitudesV0A, relBCsV0A); + if (nZdcs > 0) { + auto itZDC = mapGlobalBcWithZdc.find(globalBcMid); + if (itZDC != mapGlobalBcWithZdc.end()) { + const auto& zdc = zdcs.iteratorAt(itZDC->second); + float timeZNA = zdc.timeZNA(); + float timeZNC = zdc.timeZNC(); + float eComZNA = zdc.energyCommonZNA(); + float eComZNC = zdc.energyCommonZNC(); + udZdcsReduced(candId, timeZNA, timeZNC, eComZNA, eComZNC); + } + } + candId++; + } + + vGlobalBCs.clear(); + vColIndexBCs.clear(); + mapGlobalBcWithV0A.clear(); + mapGlobalBcWithZdc.clear(); + vAmbFwdTrackIndex.clear(); + vAmbFwdTrackIndexBCs.clear(); + vAmbBarTrackIndex.clear(); + vAmbBarTrackIndexBCs.clear(); + mapGlobalBcsWithMCHMIDTrackIds.clear(); + mapGlobalBcsWithBarrelTrackIds.clear(); + } + + void processSemiFwd(ForwardTracks const& fwdTracks, + o2::aod::AmbiguousFwdTracks const& ambFwdTracks, + BarrelTracks const& barrelTracks, + o2::aod::AmbiguousTracks const& ambBarrelTracks, + o2::aod::BCs const& bcs, + o2::aod::Collisions const& collisions, + o2::aod::FV0As const& fv0s, + o2::aod::Zdcs const& zdcs) + { + fDoMC = false; + createCandidates(fwdTracks, ambFwdTracks, barrelTracks, ambBarrelTracks, bcs, collisions, fv0s, zdcs, + (o2::aod::McFwdTrackLabels*)nullptr, (o2::aod::McTrackLabels*)nullptr); + } + + void processSemiFwdMC(ForwardTracks const& fwdTracks, + o2::aod::AmbiguousFwdTracks const& ambFwdTracks, + BarrelTracks const& barrelTracks, + o2::aod::AmbiguousTracks const& ambBarrelTracks, + o2::aod::BCs const& bcs, + o2::aod::Collisions const& collisions, + o2::aod::FV0As const& fv0s, + o2::aod::Zdcs const& zdcs, + o2::aod::McCollisions const& mcCollisions, + o2::aod::McParticles const& mcParticles, + o2::aod::McFwdTrackLabels const& mcFwdTrackLabels, + o2::aod::McTrackLabels const& mcBarrelTrackLabels) + { + fDoMC = true; + skimMCInfo(mcCollisions, mcParticles); + createCandidates(fwdTracks, ambFwdTracks, barrelTracks, ambBarrelTracks, bcs, collisions, fv0s, zdcs, + &mcFwdTrackLabels, &mcBarrelTrackLabels); + fNewPartIDs.clear(); + } + + PROCESS_SWITCH(UpcCandProducerSemiFwd, processSemiFwd, "Produce candidates combining forward MCH-MID and central-barrel tracks", true); + PROCESS_SWITCH(UpcCandProducerSemiFwd, processSemiFwdMC, "Produce semi-forward candidates with MC information", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGUD/Tasks/CMakeLists.txt b/PWGUD/Tasks/CMakeLists.txt index 5fe01ad5980..b30af171280 100644 --- a/PWGUD/Tasks/CMakeLists.txt +++ b/PWGUD/Tasks/CMakeLists.txt @@ -209,6 +209,11 @@ o2physics_add_dpl_workflow(upc-fwd-jpsi-rl PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(upc-semi-fwd-jpsi-rl + SOURCES upcSemiFwdJpsiRl.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(upc-event-itsrof-counter SOURCES upcEventITSROFcounter.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB diff --git a/PWGUD/Tasks/upcSemiFwdJpsiRl.cxx b/PWGUD/Tasks/upcSemiFwdJpsiRl.cxx new file mode 100644 index 00000000000..eb77adce860 --- /dev/null +++ b/PWGUD/Tasks/upcSemiFwdJpsiRl.cxx @@ -0,0 +1,532 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file upcSemiFwdJpsiRl.cxx +/// \brief UPC semi-forward J/psi -> mu+mu- analysis pairing one forward MCH-MID muon +/// with one central-barrel track. Consumes UD tables produced by +/// upcCandProducerSemiFwd. +/// \author Roman Lavicka, roman.lavicka@cern.ch +/// \since 10.05.2026 +/// +/// executable name o2-analysis-ud-upc-semi-fwd-jpsi-rl + +#include "PWGUD/DataModel/UDTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +// table for saving tree with semi-forward dimuon info +namespace jpsisfw +{ +// dimuon +DECLARE_SOA_COLUMN(RunNumber, runNumber, int); +DECLARE_SOA_COLUMN(M, m, float); +DECLARE_SOA_COLUMN(Energy, energy, float); +DECLARE_SOA_COLUMN(Px, px, float); +DECLARE_SOA_COLUMN(Py, py, float); +DECLARE_SOA_COLUMN(Pz, pz, float); +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(Rap, rap, float); +DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(PhiAv, phiAv, float); +DECLARE_SOA_COLUMN(PhiCh, phiCh, float); +// positive (p) and negative (n) tracks +DECLARE_SOA_COLUMN(EnergyP, energyP, float); +DECLARE_SOA_COLUMN(Pxp, pxp, float); +DECLARE_SOA_COLUMN(Pyp, pyp, float); +DECLARE_SOA_COLUMN(Pzp, pzp, float); +DECLARE_SOA_COLUMN(Ptp, ptp, float); +DECLARE_SOA_COLUMN(Etap, etap, float); +DECLARE_SOA_COLUMN(Phip, phip, float); +DECLARE_SOA_COLUMN(IsFwdP, isFwdP, int); // 1 if the positive track is forward, 0 if barrel +DECLARE_SOA_COLUMN(EnergyN, energyN, float); +DECLARE_SOA_COLUMN(Pxn, pxn, float); +DECLARE_SOA_COLUMN(Pyn, pyn, float); +DECLARE_SOA_COLUMN(Pzn, pzn, float); +DECLARE_SOA_COLUMN(Ptn, ptn, float); +DECLARE_SOA_COLUMN(Etan, etan, float); +DECLARE_SOA_COLUMN(Phin, phin, float); +DECLARE_SOA_COLUMN(IsFwdN, isFwdN, int); +// zn +DECLARE_SOA_COLUMN(Tzna, tzna, float); +DECLARE_SOA_COLUMN(Ezna, ezna, float); +DECLARE_SOA_COLUMN(Tznc, tznc, float); +DECLARE_SOA_COLUMN(Eznc, eznc, float); +DECLARE_SOA_COLUMN(Nclass, nclass, int); +} // namespace jpsisfw + +namespace o2::aod +{ +DECLARE_SOA_TABLE(JpsiSemiFwdRL, "AOD", "JPSISFW", + jpsisfw::RunNumber, + jpsisfw::M, jpsisfw::Energy, jpsisfw::Px, jpsisfw::Py, jpsisfw::Pz, jpsisfw::Pt, jpsisfw::Rap, jpsisfw::Phi, + jpsisfw::PhiAv, jpsisfw::PhiCh, + jpsisfw::EnergyP, jpsisfw::Pxp, jpsisfw::Pyp, jpsisfw::Pzp, jpsisfw::Ptp, jpsisfw::Etap, jpsisfw::Phip, jpsisfw::IsFwdP, + jpsisfw::EnergyN, jpsisfw::Pxn, jpsisfw::Pyn, jpsisfw::Pzn, jpsisfw::Ptn, jpsisfw::Etan, jpsisfw::Phin, jpsisfw::IsFwdN, + jpsisfw::Tzna, jpsisfw::Ezna, jpsisfw::Tznc, jpsisfw::Eznc, jpsisfw::Nclass); +} // namespace o2::aod + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +// constants used in the forward-muon selection (mirror upcFwdJpsiRl) +const float kRAbsMin = 17.6; +const float kRAbsMid = 26.5; +const float kRAbsMax = 89.5; +const float kPDca1 = 200.; +const float kPDca2 = 200.; +const float kFwdEtaMin = -4.0; +const float kFwdEtaMax = -2.5; +const float kFwdPtMin = 0.; + +const float kMaxAmpV0A = 100.; +const float kMaxZDCTime = 2.; +const float kMaxZDCTimeHisto = 10.; +const float kInvalidFloat = -999.; +const int kMaxRelBCsV0A = 1; + +struct UpcSemiFwdJpsiRl { + + using CandidatesSemiFwd = soa::Join; + using ForwardTracks = soa::Join; + using BarrelTracks = soa::Join; + + Produces dimuSel; + + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry reg0n0n{"reg0n0n", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry regXn0n{"regXn0n", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry regXnXn{"regXnXn", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + + // CONFIGURABLES + static constexpr double Pi = o2::constants::math::PI; + // pair kinematics + Configurable nBinsPt{"nBinsPt", 250, "N bins in pT histo"}; + Configurable lowPt{"lowPt", 0., "lower limit in pT histo"}; + Configurable highPt{"highPt", 2., "upper limit in pT histo"}; + Configurable nBinsMass{"nBinsMass", 500, "N bins in mass histo"}; + Configurable lowMass{"lowMass", 0., "lower limit in mass histo"}; + Configurable highMass{"highMass", 10., "upper limit in mass histo"}; + Configurable nBinsEta{"nBinsEta", 600, "N bins in eta histo"}; + Configurable lowEta{"lowEta", -10., "lower limit in eta histo"}; + Configurable highEta{"highEta", 5., "upper limit in eta histo"}; + Configurable nBinsRapidity{"nBinsRapidity", 500, "N bins in rapidity histo"}; + Configurable lowRapidity{"lowRapidity", -4.5, "lower limit in rapidity histo"}; + Configurable highRapidity{"highRapidity", 0., "upper limit in rapidity histo"}; + Configurable nBinsPhi{"nBinsPhi", 600, "N bins in phi histo"}; + Configurable lowPhi{"lowPhi", -Pi, "lower limit in phi histo"}; + Configurable highPhi{"highPhi", Pi, "upper limit in phi histo"}; + // single-track histos + Configurable nBinsPtSingle{"nBinsPtSingle", 500, "N bins in pT histo single track"}; + Configurable lowPtSingle{"lowPtSingle", 0., "lower limit in pT histo single track"}; + Configurable highPtSingle{"highPtSingle", 5., "upper limit in pT histo single track"}; + Configurable nBinsEtaSingle{"nBinsEtaSingle", 250, "N bins in eta histo single track"}; + Configurable lowEtaSingle{"lowEtaSingle", -4.5, "lower limit in eta histo single track"}; + Configurable highEtaSingle{"highEtaSingle", 1.5, "upper limit in eta histo single track"}; + Configurable nBinsPhiSingle{"nBinsPhiSingle", 600, "N bins in phi histo single track"}; + Configurable lowPhiSingle{"lowPhiSingle", -Pi, "lower limit in phi histo single track"}; + Configurable highPhiSingle{"highPhiSingle", Pi, "upper limit in phi histo single track"}; + // ZDC + Configurable nBinsZDCen{"nBinsZDCen", 200, "N bins in ZN energy"}; + Configurable lowEnZN{"lowEnZN", -50., "lower limit in ZN energy histo"}; + Configurable highEnZN{"highEnZN", 250., "upper limit in ZN energy histo"}; + + // central-barrel-muon selection + Configurable barEtaMin{"barEtaMin", -0.9, "barrel min eta"}; + Configurable barEtaMax{"barEtaMax", 0.9, "barrel max eta"}; + Configurable barPtMin{"barPtMin", 0.7, "barrel min pT (GeV/c)"}; + Configurable barMinITSNCls{"barMinITSNCls", 4, "barrel min number of ITS clusters"}; + Configurable barMinTPCNClsCR{"barMinTPCNClsCR", 70, "barrel min number of TPC crossed rows"}; + Configurable barMaxTpcNSigmaMu{"barMaxTpcNSigmaMu", 4., "barrel max |TPC nSigma muon|"}; + Configurable barRequireTOF{"barRequireTOF", false, "require TOF on the barrel track"}; + Configurable barMaxTofNSigmaMu{"barMaxTofNSigmaMu", 3., "barrel max |TOF nSigma muon| when TOF is present"}; + + void init(InitContext&) + { + std::vector ptFitBinning = { + 0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, + 0.11, 0.12, 0.13, 0.14, 0.15, 0.175, 0.20, 0.25, 0.30, 0.40, 0.50, + 0.60, 0.70, 0.80, 0.90, 1.00, 1.20, 1.40, 1.60, 1.80, 2.00, 2.50, + 3.00, 3.50}; + + const AxisSpec axisPt{nBinsPt, lowPt, highPt, "#it{p}_{T} GeV/#it{c}"}; + const AxisSpec axisPtFit = {ptFitBinning, "#it{p}_{T} (GeV/c)"}; + const AxisSpec axisMass{nBinsMass, lowMass, highMass, "m_{#mu#mu} GeV/#it{c}^{2}"}; + const AxisSpec axisEta{nBinsEta, lowEta, highEta, "#eta"}; + const AxisSpec axisRapidity{nBinsRapidity, lowRapidity, highRapidity, "Rapidity"}; + const AxisSpec axisPhi{nBinsPhi, lowPhi, highPhi, "#varphi"}; + const AxisSpec axisPtSingle{nBinsPtSingle, lowPtSingle, highPtSingle, "#it{p}_{T}_{ trk} GeV/#it{c}"}; + const AxisSpec axisTimeZN{200, -10, 10, "ZDC time (ns)"}; + const AxisSpec axisEnergyZNA{nBinsZDCen, lowEnZN, highEnZN, "ZNA energy (TeV)"}; + const AxisSpec axisEnergyZNC{nBinsZDCen, lowEnZN, highEnZN, "ZNC energy (TeV)"}; + const AxisSpec axisEtaSingle{nBinsEtaSingle, lowEtaSingle, highEtaSingle, "#eta_{trk}"}; + const AxisSpec axisPhiSingle{nBinsPhiSingle, lowPhiSingle, highPhiSingle, "#varphi_{trk}"}; + + // pair histos + registry.add("hMass", "Invariant mass of mu pairs;;#counts", kTH1D, {axisMass}); + registry.add("hPt", "Transverse momentum of mu pairs;;#counts", kTH1D, {axisPt}); + registry.add("hPtFit", "Transverse momentum of mu pairs;;#counts", kTH1D, {axisPtFit}); + registry.add("hEta", "Pseudorapidity of mu pairs;;#counts", kTH1D, {axisEta}); + registry.add("hRapidity", "Rapidity of mu pairs;;#counts", kTH1D, {axisRapidity}); + registry.add("hPhi", "#varphi of mu pairs;;#counts", kTH1D, {axisPhi}); + registry.add("hCharge", "Charge;;;#counts", kTH1D, {{5, -2.5, 2.5}}); + registry.add("hContrib", "hContrib;;#counts", kTH1D, {{20, -0.5, 19.5}}); + registry.add("hEvSign", "Sum of the charges of all the tracks in each event;;#counts", kTH1D, {{5, -2.5, 2.5}}); + registry.add("hSameSign", "hSameSign;;#counts", kTH1D, {{20, -0.5, 19.5}}); + registry.add("hPhiCharge", "#phi #it{charge}", kTH1D, {axisPhi}); + registry.add("hPhiAverage", "#phi #it{average}", kTH1D, {axisPhi}); + + // forward / barrel single-track histos + registry.add("hPtTrkFwd", "Pt of forward muons;;#counts", kTH1D, {axisPtSingle}); + registry.add("hPtTrkBar", "Pt of barrel muons;;#counts", kTH1D, {axisPtSingle}); + registry.add("hEtaTrkFwd", "#eta of forward muons;;#counts", kTH1D, {axisEtaSingle}); + registry.add("hEtaTrkBar", "#eta of barrel muons;;#counts", kTH1D, {axisEtaSingle}); + registry.add("hPhiTrkFwd", "#varphi of forward muons;;#counts", kTH1D, {axisPhiSingle}); + registry.add("hPhiTrkBar", "#varphi of barrel muons;;#counts", kTH1D, {axisPhiSingle}); + + // ZDC + registry.add("hTimeZNA", "ZNA Times;;#counts", kTH1D, {axisTimeZN}); + registry.add("hTimeZNC", "ZNC Times;;#counts", kTH1D, {axisTimeZN}); + registry.add("hEnergyZN", "ZNA vs ZNC energy", kTH2D, {axisEnergyZNA, axisEnergyZNC}); + + // neutron classes + reg0n0n.add("hMass", "Invariant mass of mu pairs - 0n0n;;#counts", kTH1D, {axisMass}); + reg0n0n.add("hPt", "Transverse momentum of mu pairs - 0n0n;;#counts", kTH1D, {axisPt}); + reg0n0n.add("hEta", "Pseudorapidity of mu pairs - 0n0n;;#counts", kTH1D, {axisEta}); + reg0n0n.add("hRapidity", "Rapidity of mu pairs - 0n0n;;#counts", kTH1D, {axisRapidity}); + reg0n0n.add("hPtFit", "Transverse momentum of mu pairs - 0n0n;;#counts", kTH1D, {axisPtFit}); + + regXn0n.add("hMass", "Invariant mass of mu pairs - Xn0n;;#counts", kTH1D, {axisMass}); + regXn0n.add("hPt", "Transverse momentum of mu pairs - Xn0n;;#counts", kTH1D, {axisPt}); + regXn0n.add("hEta", "Pseudorapidity of mu pairs - Xn0n;;#counts", kTH1D, {axisEta}); + regXn0n.add("hRapidity", "Rapidity of mu pairs - Xn0n;;#counts", kTH1D, {axisRapidity}); + regXn0n.add("hPtFit", "Transverse momentum of mu pairs - Xn0n;;#counts", kTH1D, {axisPtFit}); + + regXnXn.add("hMass", "Invariant mass of mu pairs - XnXn;;#counts", kTH1D, {axisMass}); + regXnXn.add("hPt", "Transverse momentum of mu pairs - XnXn;;#counts", kTH1D, {axisPt}); + regXnXn.add("hEta", "Pseudorapidity of mu pairs - XnXn;;#counts", kTH1D, {axisEta}); + regXnXn.add("hRapidity", "Rapidity of mu pairs - XnXn;;#counts", kTH1D, {axisRapidity}); + regXnXn.add("hPtFit", "Transverse momentum of mu pairs - XnXn;;#counts", kTH1D, {axisPtFit}); + } + + // FUNCTIONS + using LorentzVec = ROOT::Math::PxPyPzMVector; + + // collect tracks per candidate into a map + template + void collectCandIDs(std::unordered_map>& tracksPerCand, TTracks& tracks) + { + for (const auto& tr : tracks) { + int32_t candId = tr.udCollisionId(); + if (candId < 0) { + continue; + } + tracksPerCand[candId].push_back(tr.globalIndex()); + } + } + + struct ZDCinfo { + float timeA; + float timeC; + float enA; + float enC; + }; + + void collectCandZDCInfo(std::unordered_map& zdcPerCand, o2::aod::UDZdcsReduced const& ZDCs) + { + for (const auto& zdc : ZDCs) { + int32_t candId = zdc.udCollisionId(); + if (candId < 0) { + continue; + } + zdcPerCand[candId].timeA = zdc.timeZNA(); + zdcPerCand[candId].timeC = zdc.timeZNC(); + zdcPerCand[candId].enA = zdc.energyCommonZNA(); + zdcPerCand[candId].enC = zdc.energyCommonZNC(); + + if (std::isinf(zdcPerCand[candId].timeA)) + zdcPerCand[candId].timeA = kInvalidFloat; + if (std::isinf(zdcPerCand[candId].timeC)) + zdcPerCand[candId].timeC = kInvalidFloat; + if (std::isinf(zdcPerCand[candId].enA)) + zdcPerCand[candId].enA = kInvalidFloat; + if (std::isinf(zdcPerCand[candId].enC)) + zdcPerCand[candId].enC = kInvalidFloat; + } + } + + // forward-muon selection (MCH-MID muon, like upcFwdJpsiRl) + template + bool isFwdMuonSelected(const TTrack& fwdTrack) + { + float rAbs = fwdTrack.rAtAbsorberEnd(); + float pDca = fwdTrack.pDca(); + auto mMu = o2::constants::physics::MassMuon; + LorentzVec p(fwdTrack.px(), fwdTrack.py(), fwdTrack.pz(), mMu); + float eta = p.Eta(); + float pt = p.Pt(); + float pDcaMax = rAbs < kRAbsMid ? kPDca1 : kPDca2; + + if (eta < kFwdEtaMin || eta > kFwdEtaMax) + return false; + if (pt < kFwdPtMin) + return false; + if (rAbs < kRAbsMin || rAbs > kRAbsMax) + return false; + if (pDca > pDcaMax) + return false; + if (fwdTrack.chi2MatchMCHMID() <= 0) // require MID match + return false; + return true; + } + + // barrel-muon selection + template + bool isBarrelMuonSelected(const TTrack& barTrack) + { + auto mMu = o2::constants::physics::MassMuon; + LorentzVec p(barTrack.px(), barTrack.py(), barTrack.pz(), mMu); + float eta = p.Eta(); + float pt = p.Pt(); + + if (pt < barPtMin) + return false; + if (eta < barEtaMin || eta > barEtaMax) + return false; + if (barTrack.itsNCls() < static_cast(barMinITSNCls)) + return false; + if (barTrack.tpcNClsCrossedRows() < static_cast(barMinTPCNClsCR)) + return false; + if (std::abs(barTrack.tpcNSigmaMu()) > barMaxTpcNSigmaMu) + return false; + if (barRequireTOF && !barTrack.hasTOF()) + return false; + if (barTrack.hasTOF() && std::abs(barTrack.tofNSigmaMu()) > barMaxTofNSigmaMu) + return false; + return true; + } + + // azimuth anisotropy phi + void computePhiAnis(LorentzVec p1, LorentzVec p2, int sign1, float& phiAverage, float& phiCharge) + { + auto tSum = p1 + p2; + float halfUnity = 0.5; + decltype(tSum) tDiffCh, tDiffAv; + if (sign1 > 0) { + tDiffCh = p1 - p2; + if (gRandom->Rndm() > halfUnity) + tDiffAv = p1 - p2; + else + tDiffAv = p2 - p1; + } else { + tDiffCh = p2 - p1; + if (gRandom->Rndm() > halfUnity) + tDiffAv = p2 - p1; + else + tDiffAv = p1 - p2; + } + phiAverage = ROOT::Math::VectorUtil::DeltaPhi(tSum, tDiffAv); + phiCharge = ROOT::Math::VectorUtil::DeltaPhi(tSum, tDiffCh); + } + + // process a single (forward, barrel) pair + void processPair(CandidatesSemiFwd::iterator const& cand, + ForwardTracks::iterator const& fwd, BarrelTracks::iterator const& bar, + ZDCinfo const& zdc) + { + // V0 selection + const auto& ampsV0A = cand.amplitudesV0A(); + const auto& ampsRelBCsV0A = cand.ampRelBCsV0A(); + for (unsigned int i = 0; i < ampsV0A.size(); ++i) { + if (std::abs(ampsRelBCsV0A[i]) <= kMaxRelBCsV0A) { + if (ampsV0A[i] > kMaxAmpV0A) + return; + } + } + + // opposite charge only + int pairCharge = fwd.sign() + bar.sign(); + if (pairCharge != 0) { + registry.fill(HIST("hSameSign"), cand.numContrib()); + return; + } + + // track selection + if (!isFwdMuonSelected(fwd)) + return; + if (!isBarrelMuonSelected(bar)) + return; + + // form Lorentz vectors + auto mMu = o2::constants::physics::MassMuon; + LorentzVec pFwd(fwd.px(), fwd.py(), fwd.pz(), mMu); + LorentzVec pBar(bar.px(), bar.py(), bar.pz(), mMu); + auto p = pFwd + pBar; + + // pair-level cuts + if (p.M() < lowMass || p.M() > highMass) + return; + if (p.Pt() < lowPt || p.Pt() > highPt) + return; + if (p.Rapidity() < lowRapidity || p.Rapidity() > highRapidity) + return; + + // azimuth anisotropy + float phiAverage = 0; + float phiCharge = 0; + computePhiAnis(pFwd, pBar, fwd.sign(), phiAverage, phiCharge); + + // ZDC info histograms + if (std::abs(zdc.timeA) < kMaxZDCTimeHisto) + registry.fill(HIST("hTimeZNA"), zdc.timeA); + if (std::abs(zdc.timeC) < kMaxZDCTimeHisto) + registry.fill(HIST("hTimeZNC"), zdc.timeC); + registry.fill(HIST("hEnergyZN"), zdc.enA, zdc.enC); + + // neutron classes + bool neutronA = std::abs(zdc.timeA) < kMaxZDCTime && !std::isinf(zdc.timeA); + bool neutronC = std::abs(zdc.timeC) < kMaxZDCTime && !std::isinf(zdc.timeC); + int znClass = -1; + + if (!neutronC && !neutronA) { + znClass = 1; + reg0n0n.fill(HIST("hMass"), p.M()); + reg0n0n.fill(HIST("hPt"), p.Pt()); + reg0n0n.fill(HIST("hPtFit"), p.Pt()); + reg0n0n.fill(HIST("hEta"), p.Eta()); + reg0n0n.fill(HIST("hRapidity"), p.Rapidity()); + } else if (neutronA ^ neutronC) { + znClass = neutronA ? 2 : 3; + regXn0n.fill(HIST("hMass"), p.M()); + regXn0n.fill(HIST("hPt"), p.Pt()); + regXn0n.fill(HIST("hPtFit"), p.Pt()); + regXn0n.fill(HIST("hEta"), p.Eta()); + regXn0n.fill(HIST("hRapidity"), p.Rapidity()); + } else if (neutronA && neutronC) { + znClass = 4; + regXnXn.fill(HIST("hMass"), p.M()); + regXnXn.fill(HIST("hPt"), p.Pt()); + regXnXn.fill(HIST("hPtFit"), p.Pt()); + regXnXn.fill(HIST("hEta"), p.Eta()); + regXnXn.fill(HIST("hRapidity"), p.Rapidity()); + } + + // single-track histos by side + registry.fill(HIST("hPtTrkFwd"), pFwd.Pt()); + registry.fill(HIST("hPtTrkBar"), pBar.Pt()); + registry.fill(HIST("hEtaTrkFwd"), pFwd.Eta()); + registry.fill(HIST("hEtaTrkBar"), pBar.Eta()); + registry.fill(HIST("hPhiTrkFwd"), pFwd.Phi()); + registry.fill(HIST("hPhiTrkBar"), pBar.Phi()); + + // pair histos + registry.fill(HIST("hContrib"), cand.numContrib()); + registry.fill(HIST("hEvSign"), cand.netCharge()); + registry.fill(HIST("hMass"), p.M()); + registry.fill(HIST("hPt"), p.Pt()); + registry.fill(HIST("hPtFit"), p.Pt()); + registry.fill(HIST("hEta"), p.Eta()); + registry.fill(HIST("hRapidity"), p.Rapidity()); + registry.fill(HIST("hPhi"), p.Phi()); + registry.fill(HIST("hCharge"), fwd.sign()); + registry.fill(HIST("hCharge"), bar.sign()); + registry.fill(HIST("hPhiAverage"), phiAverage); + registry.fill(HIST("hPhiCharge"), phiCharge); + + // assign positive/negative leg + bool fwdIsPos = fwd.sign() > 0; + const LorentzVec& pPos = fwdIsPos ? pFwd : pBar; + const LorentzVec& pNeg = fwdIsPos ? pBar : pFwd; + int isFwdPos = fwdIsPos ? 1 : 0; + int isFwdNeg = fwdIsPos ? 0 : 1; + + dimuSel(cand.runNumber(), + p.M(), p.E(), p.Px(), p.Py(), p.Pz(), p.Pt(), p.Rapidity(), p.Phi(), + phiAverage, phiCharge, + pPos.E(), pPos.Px(), pPos.Py(), pPos.Pz(), pPos.Pt(), pPos.Eta(), pPos.Phi(), isFwdPos, + pNeg.E(), pNeg.Px(), pNeg.Py(), pNeg.Pz(), pNeg.Pt(), pNeg.Eta(), pNeg.Phi(), isFwdNeg, + zdc.timeA, zdc.enA, zdc.timeC, zdc.enC, znClass); + } + + // PROCESS FUNCTION + void processData(CandidatesSemiFwd const& eventCandidates, + o2::aod::UDZdcsReduced const& ZDCs, + ForwardTracks const& fwdTracks, + BarrelTracks const& barrelTracks) + { + std::unordered_map> fwdTracksPerCand; + collectCandIDs(fwdTracksPerCand, fwdTracks); + + std::unordered_map> barTracksPerCand; + collectCandIDs(barTracksPerCand, barrelTracks); + + std::unordered_map zdcPerCand; + collectCandZDCInfo(zdcPerCand, ZDCs); + + // loop over candidates that have at least one forward track + for (const auto& [candID, fwdIds] : fwdTracksPerCand) { + auto itBar = barTracksPerCand.find(candID); + if (itBar == barTracksPerCand.end()) + continue; // need at least one barrel companion + const auto& barIds = itBar->second; + auto cand = eventCandidates.iteratorAt(candID); + + ZDCinfo zdc; + auto itZdc = zdcPerCand.find(candID); + if (itZdc != zdcPerCand.end()) { + zdc = itZdc->second; + } else { + zdc.timeA = kInvalidFloat; + zdc.timeC = kInvalidFloat; + zdc.enA = kInvalidFloat; + zdc.enC = kInvalidFloat; + } + + for (const auto& fwdIdx : fwdIds) { + for (const auto& barIdx : barIds) { + auto fwd = fwdTracks.iteratorAt(fwdIdx); + auto bar = barrelTracks.iteratorAt(barIdx); + processPair(cand, fwd, bar, zdc); + } + } + } + } + + PROCESS_SWITCH(UpcSemiFwdJpsiRl, processData, "Process semi-forward J/psi candidates", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + }; +}