diff --git a/PWGUD/Core/UPCHelpers.h b/PWGUD/Core/UPCHelpers.h index a7ca4b63af4..9dfa5fe3776 100644 --- a/PWGUD/Core/UPCHelpers.h +++ b/PWGUD/Core/UPCHelpers.h @@ -22,8 +22,6 @@ #include "CommonConstants/LHCConstants.h" #include "Framework/AnalysisDataModel.h" -#include "TLorentzVector.h" - using BCsWithBcSels = o2::soa::Join; using ForwardTracks = o2::soa::Join; @@ -96,14 +94,16 @@ struct FITInfo { int32_t distClosestBcT0A = 999; }; +constexpr double AbsorberMid = 26.5; + template void applyFwdCuts(UPCCutparHolder& upcCuts, const T& track, TSelectorsArray& fwdSelectors) { - fwdSelectors[kFwdSelPt] = track.pt() > upcCuts.getFwdPtLow() && track.pt() < upcCuts.getFwdPtHigh(); // check pt - fwdSelectors[kFwdSelEta] = track.eta() > upcCuts.getFwdEtaLow() && track.eta() < upcCuts.getFwdEtaHigh(); // check pseudorapidity - fwdSelectors[kFwdSelRabs] = track.rAtAbsorberEnd() > upcCuts.getMuonRAtAbsorberEndLow() && track.rAtAbsorberEnd() < upcCuts.getMuonRAtAbsorberEndHigh(); // check muon R - fwdSelectors[kFwdSelpDCA] = track.pDca() < 26.5 ? track.pDca() < upcCuts.getMuonPDcaHighFirst() : track.pDca() < upcCuts.getMuonPDcaHighSecond(); // check pDCA - fwdSelectors[kFwdSelChi2] = track.chi2() > upcCuts.getFwdChi2Low() && track.chi2() < upcCuts.getFwdChi2High(); // check chi2 + fwdSelectors[kFwdSelPt] = track.pt() > upcCuts.getFwdPtLow() && track.pt() < upcCuts.getFwdPtHigh(); // check pt + fwdSelectors[kFwdSelEta] = track.eta() > upcCuts.getFwdEtaLow() && track.eta() < upcCuts.getFwdEtaHigh(); // check pseudorapidity + fwdSelectors[kFwdSelRabs] = track.rAtAbsorberEnd() > upcCuts.getMuonRAtAbsorberEndLow() && track.rAtAbsorberEnd() < upcCuts.getMuonRAtAbsorberEndHigh(); // check muon R + fwdSelectors[kFwdSelpDCA] = track.rAtAbsorberEnd() < AbsorberMid ? track.pDca() < upcCuts.getMuonPDcaHighFirst() : track.pDca() < upcCuts.getMuonPDcaHighSecond(); // check pDCA + fwdSelectors[kFwdSelChi2] = track.chi2() > upcCuts.getFwdChi2Low() && track.chi2() < upcCuts.getFwdChi2High(); // check chi2 } template @@ -131,7 +131,7 @@ void applyBarrelCuts(UPCCutparHolder& upcCuts, const T& track, TSelectorsArray& if (upcCuts.getCheckMaxDcaXY()) { float dca = track.dcaXY(); - float maxDCA = 0.0105f + 0.0350f / pow(track.pt(), 1.1f); + float maxDCA = 0.0105f + 0.0350f / std::pow(track.pt(), 1.1f); barrelSelectors[kBarrelSelDCAXY] = dca < maxDCA; } else { barrelSelectors[kBarrelSelDCAXY] = true; diff --git a/PWGUD/TableProducer/CMakeLists.txt b/PWGUD/TableProducer/CMakeLists.txt index 651a51bddc6..c2fa718ddd0 100644 --- a/PWGUD/TableProducer/CMakeLists.txt +++ b/PWGUD/TableProducer/CMakeLists.txt @@ -41,7 +41,7 @@ o2physics_add_dpl_workflow(two-tracks-event-table-producer PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(fwdtrack-propagation +o2physics_add_dpl_workflow(fwd-track-propagation SOURCES fwdTrackPropagation.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::GlobalTracking COMPONENT_NAME Analysis) @@ -60,3 +60,8 @@ o2physics_add_dpl_workflow(tau-three-prong-event-table-producer SOURCES tauThreeProngEventTableProducer.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(upc-cand-producer-muon + SOURCES upcCandProducerMuon.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::UPCCutparHolder O2::GlobalTracking + COMPONENT_NAME Analysis) diff --git a/PWGUD/TableProducer/fwdTrackPropagation.cxx b/PWGUD/TableProducer/fwdTrackPropagation.cxx index 52fa510b176..64099fac4ed 100644 --- a/PWGUD/TableProducer/fwdTrackPropagation.cxx +++ b/PWGUD/TableProducer/fwdTrackPropagation.cxx @@ -8,29 +8,37 @@ // 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 fwdTrackPropagation.cxx +/// \brief Forward track propagation task for UD tasks /// \author Nazar Burmasov, nazar.burmasov@cern.ch /// \author Diana Krupova, diana.krupova@cern.ch /// \since 04.06.2024 -#include "Framework/AnalysisDataModel.h" -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" +#include "PWGUD/DataModel/UDTables.h" #include "CCDB/BasicCCDBManager.h" #include "DataFormatsParameters/GRPMagField.h" -#include "TGeoGlobalMagField.h" -#include "Field/MagneticField.h" #include "DetectorsBase/Propagator.h" +#include "Field/MagneticField.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/DataTypes.h" +#include "Framework/runDataProcessing.h" #include "GlobalTracking/MatchGlobalFwd.h" #include "MCHTracking/TrackExtrap.h" -#include "MCHTracking/TrackParam.h" -#include "Math/SMatrix.h" #include "ReconstructionDataFormats/TrackFwd.h" -#include "PWGUD/DataModel/UDTables.h" +#include "Math/SMatrix.h" +#include "TGeoGlobalMagField.h" + +#include +#include +#include using namespace o2::framework; using namespace o2::framework::expressions; +using o2::aod::fwdtrack::ForwardTrackTypeEnum; struct FwdTrackPropagation { using ForwardTracks = o2::soa::Join; @@ -65,13 +73,12 @@ struct FwdTrackPropagation { 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::track::TrackParCovFwd fwdtrack{muon.z(), tpars, tcovs, chi2}; o2::dataformats::GlobalFwdTrack propmuon; - if (static_cast(muon.trackType()) > 2) { + if (muon.trackType() > ForwardTrackTypeEnum::GlobalForwardTrack) { // tracks without MFT o2::dataformats::GlobalFwdTrack track; track.setParameters(tpars); - track.setZ(fwdtrack.getZ()); + track.setZ(muon.z()); track.setCovariances(tcovs); auto mchTrack = fMatching.FwdtoMCH(track); o2::mch::TrackExtrap::extrapToVertex(mchTrack, vtx[0], vtx[1], vtx[2], vtxCov[0], vtxCov[1]); @@ -79,13 +86,14 @@ struct FwdTrackPropagation { propmuon.setParameters(proptrack.getParameters()); propmuon.setZ(proptrack.getZ()); propmuon.setCovariances(proptrack.getCovariances()); - } else if (static_cast(muon.trackType()) < 2) { + } else if (muon.trackType() <= ForwardTrackTypeEnum::GlobalForwardTrack) { // tracks with MFT double centerMFT[3] = {0, 0, -61.4}; - o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); - auto Bz = field->getBz(centerMFT); // Get field at centre of MFT + o2::track::TrackParCovFwd fwdtrack{muon.z(), tpars, tcovs, chi2}; + auto* field = dynamic_cast(TGeoGlobalMagField::Instance()->GetField()); + auto bz = field->getBz(centerMFT); // Get field at centre of MFT auto geoMan = o2::base::GeometryManager::meanMaterialBudget(muon.x(), muon.y(), muon.z(), vtx[0], vtx[1], vtx[2]); auto x2x0 = static_cast(geoMan.meanX2X0); - fwdtrack.propagateToVtxhelixWithMCS(vtx[2], {vtx[0], vtx[1]}, {vtxCov[0], vtxCov[1]}, Bz, x2x0); + fwdtrack.propagateToVtxhelixWithMCS(vtx[2], {vtx[0], vtx[1]}, {vtxCov[0], vtxCov[1]}, bz, x2x0); propmuon.setParameters(fwdtrack.getParameters()); propmuon.setZ(fwdtrack.getZ()); propmuon.setCovariances(fwdtrack.getCovariances()); @@ -112,11 +120,12 @@ struct FwdTrackPropagation { propFwdTracks.reserve(fwdTracks.size()); propFwdTracksCov.reserve(fwdTracks.size()); + constexpr float ZAbsorberFront = -90.f; for (const auto& t : fwdTracks) { - if (t.z() < -90.f) { + if (t.z() < ZAbsorberFront) { std::array vtx = {0.f, 0.f, 0.f}; std::array vtxCov = {0.f, 0.f}; - if (t.has_collision()) { + if (t.has_collision() && t.trackType() <= ForwardTrackTypeEnum::GlobalForwardTrack) { // propagate only global muon tracks to collision vtx auto col = cols.iteratorAt(t.collisionId()); vtx[0] = col.posX(); vtx[1] = col.posY(); @@ -129,7 +138,7 @@ struct FwdTrackPropagation { pft.getX(), pft.getY(), pft.getZ(), pft.getPhi(), pft.getTgl(), pft.getInvQPt(), pft.getEta(), pft.getPt(), pft.getP(), t.nClusters(), t.pDca(), t.rAtAbsorberEnd(), - pft.getTrackChi2(), t.chi2MatchMCHMID(), t.chi2MatchMCHMFT(), + t.chi2(), t.chi2MatchMCHMID(), t.chi2MatchMCHMFT(), t.matchScoreMCHMFT(), t.matchMFTTrackId(), t.matchMCHTrackId(), t.mchBitMap(), t.midBoards(), t.midBitMap(), t.trackTime(), t.trackTimeRes()); @@ -137,11 +146,11 @@ struct FwdTrackPropagation { // LOGP(info, "track {}, before: {} {} {} {} {} {}", t.globalIndex(), t.x(), t.y(), t.z(), t.phi(), t.tgl(), t.signed1Pt()); // LOGP(info, "track {}, after: {} {} {} {} {} {}", t.globalIndex(), pft.getX(), pft.getY(), pft.getZ(), pft.getPhi(), pft.getTgl(), pft.getInvQPt()); SMatrix55 cov = pft.getCovariances(); - float sigX = TMath::Sqrt(cov(0, 0)); - float sigY = TMath::Sqrt(cov(1, 1)); - float sigPhi = TMath::Sqrt(cov(2, 2)); - float sigTgl = TMath::Sqrt(cov(3, 3)); - float sig1Pt = TMath::Sqrt(cov(4, 4)); + float sigX = std::sqrt(cov(0, 0)); + float sigY = std::sqrt(cov(1, 1)); + float sigPhi = std::sqrt(cov(2, 2)); + float sigTgl = std::sqrt(cov(3, 3)); + float sig1Pt = std::sqrt(cov(4, 4)); auto rhoXY = static_cast(128. * cov(0, 1) / (sigX * sigY)); auto rhoPhiX = static_cast(128. * cov(0, 2) / (sigPhi * sigX)); auto rhoPhiY = static_cast(128. * cov(1, 2) / (sigPhi * sigY)); diff --git a/PWGUD/TableProducer/upcCandProducerMuon.cxx b/PWGUD/TableProducer/upcCandProducerMuon.cxx new file mode 100644 index 00000000000..329a83b432c --- /dev/null +++ b/PWGUD/TableProducer/upcCandProducerMuon.cxx @@ -0,0 +1,536 @@ +// Copyright 2019-2020 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 upcCandProducerMuon.cxx +/// \brief UPC candidate producer for forward muons without MFT +/// \author Nazar Burmasov, nazar.burmasov@cern.ch +/// \since 19.12.2025 + +#include "PWGUD/Core/UPCCutparHolder.h" +#include "PWGUD/Core/UPCHelpers.h" + +#include "CCDB/BasicCCDBManager.h" +#include "CommonConstants/LHCConstants.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DetectorsBase/Propagator.h" +#include "Field/MagneticField.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "GlobalTracking/MatchGlobalFwd.h" +#include "MCHTracking/TrackExtrap.h" +#include "ReconstructionDataFormats/TrackFwd.h" + +#include "Math/SMatrix.h" +#include "TGeoGlobalMagField.h" + +#include +#include +#include +#include + +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct UpcCandProducerMuon { + bool fDoMC{false}; + + std::map fNewPartIDs; + + Produces udMCCollisions; + Produces udMCParticles; + Produces udFwdTrackLabels; + + Produces udFwdTracks; + Produces udFwdTracksExtra; + Produces udFwdIndices; + 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 fBcWindowMCH{"fBcWindowMCH", 20, "Time window for MCH-MID to MCH-only matching for Muon candidates"}; + Configurable fMaxFV0Amp{"fMaxFV0Amp", 100.f, "Max FV0 amplitude in the same BC"}; + + using ForwardTracks = 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"); + } + + bool cut(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; + } + + 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(); + // loop over MC particles to select only the ones from signal events + // and calculate new MC table IDs + 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{}; + + // storing MC particles + for (const auto& item : fNewPartIDs) { + int32_t mcPartID = item.first; + const auto& mcPart = mcParticles.iteratorAt(mcPartID); + int32_t mcEventID = mcPart.mcCollisionId(); + int32_t newEventID = newEventIDs[mcEventID]; + // collecting new mother IDs + 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); + } + } + } + // collecting new daughter IDs + 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(newEventID, mcPart.pdgCode(), mcPart.getHepMCStatusCode(), mcPart.flags(), newMotherIDs, newDaughterIDs, + mcPart.weight(), mcPart.px(), mcPart.py(), mcPart.pz(), mcPart.e()); + newMotherIDs.clear(); + } + + // storing MC events + 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; + } + + // find starting point for scrolling in some BC map -- closest to the input gbc + template + T::iterator getStartForScroll(uint64_t inGbc, T& gbcMap) + { + auto it1 = gbcMap.lower_bound(inGbc); + typename T::iterator it; + if (it1 != gbcMap.end()) { // found lower bound + 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 { // ended up in the end + it = it1; + --it; + } + return it; + } + + // scroll over gbcMap and do some operation on container + 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); + + // start scrolling backward + 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); // move back to the starting point + 1 + + if (it == gbcMap.end()) // ended up in the end of map + return; + + gbc = it->first; + dbc = bcDiff(inGbc, gbc); + + // start scrolling forward + while (dbc <= maxDbc) { + func(it, gbc); + ++it; + if (it == gbcMap.end()) + break; + gbc = it->first; + dbc = bcDiff(inGbc, gbc); + } + } + + void getMchTrackIds(uint64_t inGbc, std::map>& mchTracksPerBC, + uint64_t maxDbc, std::map& outMchTrkIds) + { + auto fillMchIds = [&outMchTrkIds](std::map>::iterator& inIt, uint64_t gbc) { + std::vector& ids = inIt->second; + for (const auto& id : ids) + outMchTrkIds[id] = gbc; + }; + scrollBackForth(inGbc, maxDbc, mchTracksPerBC, fillMchIds); + } + + 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) + { + const auto& track = fwdTracks.iteratorAt(trackId); + float px, py, pz; + int sign; + if (fUpcCuts.getUseFwdCuts()) { + auto pft = propagateToZero(track); + bool pass = cut(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(), -1., 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; + } + + void createCandidates(ForwardTracks const& fwdTracks, + o2::aod::AmbiguousFwdTracks const& ambFwdTracks, + 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) + { + using o2::aod::fwdtrack::ForwardTrackTypeEnum::MCHStandaloneTrack; + 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]; + } + + std::map> mapGlobalBcsWithMCHMIDTrackIds; + std::map> mapGlobalBcsWithMCHTrackIds; + for (const auto& fwdTrack : fwdTracks) { + if (fwdTrack.trackType() != MCHStandaloneTrack && 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.); + if (fwdTrack.trackType() == MuonStandaloneTrack) // MCH-MID + mapGlobalBcsWithMCHMIDTrackIds[globalBC].push_back(trackId); + if (fwdTrack.trackType() == MCHStandaloneTrack) // MCH-only + mapGlobalBcsWithMCHTrackIds[globalBC].push_back(trackId); + } + + int32_t candId = 0; + for (const auto& gbc_muids : mapGlobalBcsWithMCHMIDTrackIds) { + uint64_t globalBcMid = gbc_muids.first; + 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 = gbc_muids.second; + // writing MCH-MID tracks + for (const auto& imuon : vMuonIds) { + if (!addToFwdTable(candId, imuon, globalBcMid, 0., fwdTracks, mcFwdTrackLabels)) + continue; + numContrib++; + } + if (numContrib < 1) // didn't save any MCH-MID tracks + continue; + std::map mapMchIdBc{}; + 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)) + 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(); + mapGlobalBcsWithMCHMIDTrackIds.clear(); + mapGlobalBcsWithMCHTrackIds.clear(); + } + + void processFwd(ForwardTracks const& fwdTracks, + o2::aod::AmbiguousFwdTracks const& ambFwdTracks, + 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, bcs, collisions, fv0s, zdcs, (o2::aod::McFwdTrackLabels*)nullptr); + } + + void processFwdMC(ForwardTracks const& fwdTracks, + o2::aod::AmbiguousFwdTracks const& ambFwdTracks, + 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) + { + fDoMC = true; + skimMCInfo(mcCollisions, mcParticles); + createCandidates(fwdTracks, ambFwdTracks, bcs, collisions, fv0s, zdcs, &mcFwdTrackLabels); + fNewPartIDs.clear(); + } + + PROCESS_SWITCH(UpcCandProducerMuon, processFwd, "Produce candidates in forward region", true); + PROCESS_SWITCH(UpcCandProducerMuon, processFwdMC, "Produce candidates in forward region with MC information", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}