From e961d30236fc50146cf67a4a9c30328d4e7469d1 Mon Sep 17 00:00:00 2001 From: nburmaso Date: Tue, 3 Nov 2026 17:47:01 +0300 Subject: [PATCH] Introduce compact MCH and MCH-MID efficiency table producer --- PWGUD/TableProducer/CMakeLists.txt | 5 + .../TableProducer/muonEfficiencyProducer.cxx | 259 ++++++++++++++++++ 2 files changed, 264 insertions(+) create mode 100644 PWGUD/TableProducer/muonEfficiencyProducer.cxx diff --git a/PWGUD/TableProducer/CMakeLists.txt b/PWGUD/TableProducer/CMakeLists.txt index a26219c2e6f..85923b12bc2 100644 --- a/PWGUD/TableProducer/CMakeLists.txt +++ b/PWGUD/TableProducer/CMakeLists.txt @@ -100,3 +100,8 @@ o2physics_add_dpl_workflow(upc-cand-producer-barrel SOURCES upcCandProducerBarrel.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::CCDB O2::DataFormatsParameters O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(muon-efficiency-producer + SOURCES muonEfficiencyProducer.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::GlobalTracking + COMPONENT_NAME Analysis) diff --git a/PWGUD/TableProducer/muonEfficiencyProducer.cxx b/PWGUD/TableProducer/muonEfficiencyProducer.cxx new file mode 100644 index 00000000000..9a7f0987bde --- /dev/null +++ b/PWGUD/TableProducer/muonEfficiencyProducer.cxx @@ -0,0 +1,259 @@ +// 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 muonEfficiencyProducer.cxx +/// \brief Compact forward-muon skim for offline MCH and MCH-MID efficiency studies +/// \author Nazar Burmasov (JINR) + +#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; +using namespace o2::framework; + +namespace o2::aod::mcheffskim +{ +DECLARE_SOA_COLUMN(RunNumber, runNumber, int32_t); +DECLARE_SOA_COLUMN(ReferenceBC, referenceBC, int64_t); +DECLARE_SOA_COLUMN(FT0AAmplitude, ft0aAmplitude, float); +DECLARE_SOA_COLUMN(FT0CAmplitude, ft0cAmplitude, float); +DECLARE_SOA_COLUMN(FV0AAmplitude, fv0aAmplitude, float); +DECLARE_SOA_COLUMN(CollisionId, collisionId, int32_t); // -1 for ambiguous tracks +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(Eta, eta, float); +DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(Sign, sign, int8_t); +DECLARE_SOA_COLUMN(RAbs, rAbs, float); +DECLARE_SOA_COLUMN(PDca, pDca, float); +DECLARE_SOA_COLUMN(Chi2Ndf, chi2Ndf, float); // -1 if propagation failed +DECLARE_SOA_COLUMN(NClusters, nClusters, int8_t); +DECLARE_SOA_COLUMN(MchBitMap, mchBitMap, uint16_t); +DECLARE_SOA_COLUMN(MidBitMap, midBitMap, uint8_t); +DECLARE_SOA_COLUMN(MidBoards, midBoards, uint32_t); +DECLARE_SOA_COLUMN(Chi2MatchMCHMID, chi2MatchMCHMID, float); // -1 for MCH-only tracks +DECLARE_SOA_COLUMN(TrackTime, trackTime, float); // ns relative to referenceBC +DECLARE_SOA_COLUMN(TrackTimeRes, trackTimeRes, float); +} // namespace o2::aod::mcheffskim + +namespace o2::aod +{ +DECLARE_SOA_TABLE(MchEffBCs, "AOD", "MCHEFFBC", + o2::soa::Index<>, + mcheffskim::RunNumber, mcheffskim::ReferenceBC, + mcheffskim::FT0AAmplitude, mcheffskim::FT0CAmplitude, + mcheffskim::FV0AAmplitude); +} // namespace o2::aod + +namespace o2::aod::mcheffskim +{ +DECLARE_SOA_INDEX_COLUMN(MchEffBC, mchEffBC); // links a muon to the compact BC table +} // namespace o2::aod::mcheffskim + +namespace o2::aod +{ +DECLARE_SOA_TABLE(MchEffSkims, "AOD", "MCHEFFSKIM", + o2::soa::Index<>, mcheffskim::MchEffBCId, + mcheffskim::CollisionId, + mcheffskim::Pt, mcheffskim::Eta, mcheffskim::Phi, mcheffskim::Sign, + mcheffskim::RAbs, mcheffskim::PDca, mcheffskim::Chi2Ndf, + mcheffskim::NClusters, mcheffskim::MchBitMap, + mcheffskim::MidBitMap, mcheffskim::MidBoards, + mcheffskim::Chi2MatchMCHMID, + mcheffskim::TrackTime, mcheffskim::TrackTimeRes); +} // namespace o2::aod + +struct MuonEfficiencyProducer { + static constexpr float MinEta = -4.f; + static constexpr float MaxEta = -2.5f; + static constexpr float NoAmplitude = -1.f; + static constexpr float MaxFITTime = 30.f; + + using Tracks = soa::Join; + using Matrix55 = ROOT::Math::SMatrix>; + using Vector5 = ROOT::Math::SVector; + + Produces outputBCs; + Produces output; + Service ccdb{}; + globaltracking::MatchGlobalFwd matcher; + int cachedFieldRun = -1; + + Configurable minPt{"minPt", 0.4f, "Minimum track pT (Gev/c)"}; + Configurable minRAbs{"minRAbs", 17.6f, "Minimum absorber-end radius (cm)"}; + Configurable maxRAbs{"maxRAbs", 89.5f, "Maximum absorber-end radius (cm)"}; + Configurable maxPDca{"maxPDca", 400.f, "Maximum stored pDCA (GeV cm/c)"}; + Configurable maxChi2{"maxChi2", 100.f, "Maximum MCH chi2/ndf"}; + + void init(InitContext&) + { + ccdb->setLocalObjectValidityChecking(); + } + + void setField(int run) + { + if (cachedFieldRun == run) { + return; + } + const auto runInfo = parameters::AggregatedRunInfo::buildAggregatedRunInfo(ccdb->instance(), run); + const auto ts = (runInfo.sor + runInfo.eor) / 2; + auto* grp = ccdb->getForTimeStamp("GLO/Config/GRPMagField", ts); + if (!grp) { + LOGF(fatal, "GRPMagField not found for run %d", run); + } + base::Propagator::initFieldFromGRP(grp); + if (!base::GeometryManager::isGeometryLoaded()) { + ccdb->get("GLO/Config/GeometryAligned"); + } + mch::TrackExtrap::setField(); + cachedFieldRun = run; + } + + struct Kinematics { + float pt, eta, phi; + }; + + template + bool propagateMCH(const Track& t, const std::array& vtx, const std::array& vtxCov, Kinematics& k) + { + const Vector5 params(t.x(), t.y(), t.phi(), t.tgl(), t.signed1Pt()); + const std::vector covElements{t.cXX(), t.cXY(), t.cYY(), t.cPhiX(), t.cPhiY(), t.cPhiPhi(), + t.cTglX(), t.cTglY(), t.cTglPhi(), t.cTglTgl(), + t.c1PtX(), t.c1PtY(), t.c1PtPhi(), t.c1PtTgl(), t.c1Pt21Pt2()}; + const Matrix55 cov(covElements.begin(), covElements.end()); + dataformats::GlobalFwdTrack source; + source.setParameters(params); + source.setZ(t.z()); + source.setCovariances(cov); + auto mchTrack = matcher.FwdtoMCH(source); + if (!mch::TrackExtrap::extrapToVertex(mchTrack, vtx[0], vtx[1], vtx[2], vtxCov[0], vtxCov[1])) { + return false; + } + const auto propagated = matcher.MCHtoFwd(mchTrack); + k = {.pt = static_cast(propagated.getPt()), .eta = static_cast(propagated.getEta()), .phi = static_cast(propagated.getPhi())}; + return true; + } + + void process(aod::BCs const& bcs, aod::FT0s const& ft0s, aod::FV0As const& fv0as, + Tracks const& tracks, aod::Collisions const& collisions, + aod::AmbiguousFwdTracks const& ambiguous) + { + if (bcs.size() == 0) { + return; + } + const int run = bcs.begin().runNumber(); + + std::vector ft0aByBC(bcs.size(), NoAmplitude); + std::vector ft0cByBC(bcs.size(), NoAmplitude); + std::vector fv0aByBC(bcs.size(), NoAmplitude); + std::vector skimBCIds(bcs.size(), -1); + for (const auto& ft0 : ft0s) { + if (ft0.timeA() < MaxFITTime) { + ft0aByBC[ft0.bcId()] = ft0.sumAmpA(); + } + if (ft0.timeC() < MaxFITTime) { + ft0cByBC[ft0.bcId()] = ft0.sumAmpC(); + } + } + for (const auto& fv0 : fv0as) { + if (fv0.time() >= MaxFITTime) { + continue; + } + float amplitude = 0.f; + for (const float a : fv0.amplitude()) { + if (a > 0.f) { + amplitude += a; + } + } + fv0aByBC[fv0.bcId()] = amplitude; + } + + std::vector ambiguousBCIds(tracks.size(), -1); + for (const auto& entry : ambiguous) { + const auto& ids = entry.bcIds(); + if (ids.empty()) { + continue; + } + ambiguousBCIds[entry.fwdtrackId()] = *ids.begin(); + } + + for (const auto& t : tracks) { + const auto type = t.trackType(); + if (type != aod::fwdtrack::MuonStandaloneTrack && type != aod::fwdtrack::MCHStandaloneTrack) { + continue; + } + const float rAbs = t.rAtAbsorberEnd(); + if (rAbs < minRAbs || rAbs > maxRAbs || t.pDca() > maxPDca || t.chi2() > maxChi2) { + continue; + } + Kinematics k{.pt = t.pt(), .eta = t.eta(), .phi = t.phi()}; + std::array vtx{0.f, 0.f, 0.f}; + std::array vtxCov{0.f, 0.f}; + int32_t referenceBCId = ambiguousBCIds[t.globalIndex()]; + if (t.has_collision()) { + const auto col = collisions.iteratorAt(t.collisionId()); + vtx = {col.posX(), col.posY(), col.posZ()}; + vtxCov = {col.covXX(), col.covYY()}; + referenceBCId = col.bcId(); + } + setField(run); + if (!propagateMCH(t, vtx, vtxCov, k)) { + continue; + } + if (k.eta < MinEta || k.eta > MaxEta) { + continue; + } + if (k.pt < minPt) { + continue; + } + int32_t skimBCId = -1; + if (referenceBCId >= 0) { + skimBCId = skimBCIds[referenceBCId]; + if (skimBCId < 0) { + outputBCs(run, bcs.iteratorAt(referenceBCId).globalBC(), + ft0aByBC[referenceBCId], ft0cByBC[referenceBCId], fv0aByBC[referenceBCId]); + skimBCId = outputBCs.lastIndex(); + skimBCIds[referenceBCId] = skimBCId; + } + } + output(skimBCId, t.collisionId(), k.pt, k.eta, k.phi, static_cast(t.signed1Pt() > 0 ? 1 : -1), + rAbs, t.pDca(), t.chi2(), t.nClusters(), t.mchBitMap(), t.midBitMap(), t.midBoards(), + type == aod::fwdtrack::MuonStandaloneTrack ? t.chi2MatchMCHMID() : -1.f, t.trackTime(), t.trackTimeRes()); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& context) +{ + return WorkflowSpec{adaptAnalysisTask(context)}; +}