From 9415150dda61acf60965ce219718967427866741 Mon Sep 17 00:00:00 2001 From: Francesco Mazzaschi Date: Tue, 15 Sep 2026 17:33:53 +0200 Subject: [PATCH] DoubleOmega:Improve findable study + kaon PID information --- PWGLF/DataModel/LFDoubleOmegaTables.h | 19 + .../Strangeness/doubleOmegaTreeCreator.cxx | 404 ++++++++++++------ 2 files changed, 291 insertions(+), 132 deletions(-) diff --git a/PWGLF/DataModel/LFDoubleOmegaTables.h b/PWGLF/DataModel/LFDoubleOmegaTables.h index 95a77ba2181..9cc38e08f43 100644 --- a/PWGLF/DataModel/LFDoubleOmegaTables.h +++ b/PWGLF/DataModel/LFDoubleOmegaTables.h @@ -14,6 +14,8 @@ #include +#include + namespace o2::aod { @@ -25,6 +27,9 @@ DECLARE_SOA_COLUMN(PhiDoubleOmega, phiDoubleOmega, float); DECLARE_SOA_COLUMN(DecayVtxX, decayVtxX, float); DECLARE_SOA_COLUMN(DecayVtxY, decayVtxY, float); DECLARE_SOA_COLUMN(DecayVtxZ, decayVtxZ, float); +DECLARE_SOA_COLUMN(PrimaryVtxX, primaryVtxX, float); +DECLARE_SOA_COLUMN(PrimaryVtxY, primaryVtxY, float); +DECLARE_SOA_COLUMN(PrimaryVtxZ, primaryVtxZ, float); DECLARE_SOA_COLUMN(CosPAOmega, cosPAOmega, float); DECLARE_SOA_COLUMN(CosPADirectLambda, cosPADirectLambda, float); DECLARE_SOA_COLUMN(CosPADoubleOmega, cosPADoubleOmega, float); @@ -34,6 +39,12 @@ DECLARE_SOA_COLUMN(DCAxyDirectLambdaToPV, dcaXYDirectLambdaToPV, float); DECLARE_SOA_COLUMN(DCAzDirectLambdaToPV, dcaZDirectLambdaToPV, float); DECLARE_SOA_COLUMN(DCAxyDirectKaonToPV, dcaXYDirectKaonToPV, float); DECLARE_SOA_COLUMN(DCAzDirectKaonToPV, dcaZDirectKaonToPV, float); +DECLARE_SOA_COLUMN(TPCNSigmaDirectKaon, tpcNSigmaDirectKaon, float); // -999 when TPC is unavailable or the candidate is not reconstructed. +DECLARE_SOA_COLUMN(TOFNSigmaDirectKaon, tofNSigmaDirectKaon, float); // -999 when TOF is unavailable or the candidate is not reconstructed. +// Three bits per daughter (ITS, TPC, TOF), ordered as: Omega kaon, Omega +// proton, Omega pion, direct proton, direct pion, direct kaon (also for c.c.). +// Bits 18-31 are zero; generated-only rows have a zero bitmap. +DECLARE_SOA_COLUMN(DaughterDetectorMap, daughterDetectorMap, uint32_t); DECLARE_SOA_COLUMN(MassDoubleOmega, massDoubleOmega, float); DECLARE_SOA_COLUMN(MassOmega, massOmega, float); DECLARE_SOA_COLUMN(MassXi, massXi, float); @@ -43,6 +54,7 @@ DECLARE_SOA_COLUMN(GenEta, genEta, float); DECLARE_SOA_COLUMN(GenPhi, genPhi, float); DECLARE_SOA_COLUMN(GenDecayLength, genDecayLength, float); DECLARE_SOA_COLUMN(PdgDoubleOmega, pdgDoubleOmega, int); +DECLARE_SOA_COLUMN(GenMotherId, genMotherId, int64_t); // MC particle index within the input dataframe. DECLARE_SOA_COLUMN(IsReco, isReco, bool); } // namespace DoubleOmegaTables @@ -53,6 +65,9 @@ DECLARE_SOA_COLUMN(IsReco, isReco, bool); DoubleOmegaTables::DecayVtxX, \ DoubleOmegaTables::DecayVtxY, \ DoubleOmegaTables::DecayVtxZ, \ + DoubleOmegaTables::PrimaryVtxX, \ + DoubleOmegaTables::PrimaryVtxY, \ + DoubleOmegaTables::PrimaryVtxZ, \ DoubleOmegaTables::CosPAOmega, \ DoubleOmegaTables::CosPADirectLambda, \ DoubleOmegaTables::CosPADoubleOmega, \ @@ -62,6 +77,9 @@ DECLARE_SOA_COLUMN(IsReco, isReco, bool); DoubleOmegaTables::DCAzDirectLambdaToPV, \ DoubleOmegaTables::DCAxyDirectKaonToPV, \ DoubleOmegaTables::DCAzDirectKaonToPV, \ + DoubleOmegaTables::TPCNSigmaDirectKaon, \ + DoubleOmegaTables::TOFNSigmaDirectKaon, \ + DoubleOmegaTables::DaughterDetectorMap, \ DoubleOmegaTables::MassDoubleOmega, \ DoubleOmegaTables::MassOmega, \ DoubleOmegaTables::MassXi @@ -76,6 +94,7 @@ DECLARE_SOA_TABLE(DoubleOmegaTableMC, "AOD", "DBLOMEGAMCTABLE", DoubleOmegaTables::GenPhi, DoubleOmegaTables::GenDecayLength, DoubleOmegaTables::PdgDoubleOmega, + DoubleOmegaTables::GenMotherId, DoubleOmegaTables::IsReco); #undef DOUBLE_OMEGA_RECO_COLUMNS diff --git a/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx b/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx index 021fee63251..0dc52f1beee 100644 --- a/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx +++ b/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx @@ -19,6 +19,7 @@ #include "Common/Core/trackUtilities.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -35,6 +36,7 @@ #include #include #include +#include #include #include #include @@ -49,6 +51,7 @@ #include #include #include +#include #include #include @@ -60,10 +63,9 @@ using CollisionsTable = soa::Join; using FullCascades = aod::CascDataExt; -using TracksFull = soa::Join; -using TracksFullIU = soa::Join; -using TracksFullIUMC = soa::Join; -using TracksFullMC = soa::Join; +using TracksFull = soa::Join; +using TracksFullIU = soa::Join; +using TracksFullIUMC = soa::Join; struct DoubleOmegaCandidate { float pt = -999.f; @@ -72,6 +74,9 @@ struct DoubleOmegaCandidate { float x = -999.f; float y = -999.f; float z = -999.f; + float pvX = -999.f; + float pvY = -999.f; + float pvZ = -999.f; float cosPAOmega = -999.f; float cosPADirectLambda = -999.f; float cosPADoubleOmega = -999.f; @@ -81,6 +86,9 @@ struct DoubleOmegaCandidate { float dcaZDirectLambdaToPV = -999.f; float dcaXYDirectKaonToPV = -999.f; float dcaZDirectKaonToPV = -999.f; + float tpcNSigmaDirectKaon = -999.f; + float tofNSigmaDirectKaon = -999.f; + uint32_t daughterDetectorMap = 0; float mass = -999.f; float massOmega = -999.f; float massXi = -999.f; @@ -93,6 +101,11 @@ struct DoubleOmegaMCInfo { float phi = -999.f; float decayLength = -999.f; int pdgCode = 0; + // Generated vertices: double-Omega production point (primary vertex), double-Omega + // decay point, and Omega decay point. + std::array genProductionVertex{-999.f, -999.f, -999.f}; + std::array genDecayVertex{-999.f, -999.f, -999.f}; + std::array genOmegaDecayVertex{-999.f, -999.f, -999.f}; }; struct LambdaCandidate { @@ -106,9 +119,9 @@ struct BuiltLambdaCandidate { LambdaCandidate candidate; o2::track::TrackParCov parentTrack; std::array decayVertex{}; - int64_t v0Id = -1; int64_t posTrackId = -1; int64_t negTrackId = -1; + uint32_t daughterDetectorMap = 0; // Proton, pion; three detector bits each. }; struct BuiltOmegaCandidate { @@ -121,11 +134,11 @@ struct BuiltOmegaCandidate { float massOmega = -999.f; float massXi = -999.f; int8_t sign = 0; - int64_t cascadeId = -1; int64_t posTrackId = -1; int64_t negTrackId = -1; int64_t bachelorId = -1; o2::track::TrackParCov parentTrack; + uint32_t daughterDetectorMap = 0; // Kaon, proton, pion; three detector bits each. }; struct doubleOmegaTreeCreator { @@ -198,6 +211,14 @@ struct doubleOmegaTreeCreator { HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + template + static uint32_t detectorMap(Track const& track) + { + return static_cast(track.hasITS()) | + (static_cast(track.hasTPC()) << 1) | + (static_cast(track.hasTOF()) << 2); + } + template bool selectTrack(T const& track) { @@ -209,6 +230,7 @@ struct doubleOmegaTreeCreator { track.tpcNClsCrossedRows() < minNCrossedRows || track.tpcNClsCrossedRows() < 0.8f * track.tpcNClsFindable() || track.tpcNClsShared() > maxNSharedTPCClus) { + LOG(debug) << "Track failed selection: itsNCls=" << static_cast(track.itsNCls()) << ", tpcNClsFound=" << track.tpcNClsFound() << ", tpcNClsCrossedRows=" << static_cast(track.tpcNClsCrossedRows()) << ", tpcNClsFindable=" << static_cast(track.tpcNClsFindable()) << ", tpcNClsShared=" << static_cast(track.tpcNClsShared()); return false; } return true; @@ -239,7 +261,6 @@ struct doubleOmegaTreeCreator { fitter.setMatCorrType(static_cast(cfgMaterialCorrection.value)); fitter3Body.setMatCorrType(static_cast(cfgMaterialCorrection.value)); - LOG(info) << "Retrieved GRP for timestamp " << bc.timestamp() << " with magnetic field " << mBz << " kG"; mRunNumber = bc.runNumber(); if (cfgSkimmedProcessing) { zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), "fDoubleOmega,fOmegaXi"); @@ -247,8 +268,68 @@ struct doubleOmegaTreeCreator { } } + void printConfig() const + { + LOG(info) << "doubleOmegaTreeCreator runtime configuration (boolean values: 0=false, 1=true):"; + const auto printValue = [](auto const& config) { + LOG(info) << config.name << " = " << config.value; + }; + const auto printAxis = [](ConfigurableAxis const& axis) { + std::ostringstream values; + values << '['; + for (size_t i = 0; i < axis.value.size(); ++i) { + if (i > 0) { + values << ", "; + } + values << axis.value[i]; + } + values << ']'; + LOG(info) << axis.name << " = " << values.str(); + }; + + printValue(cfgSkimmedProcessing); + printValue(cfgApplyEventSelection); + printValue(ccdburl); + printValue(cfgMaterialCorrection); + printAxis(zVtxAxis); + printAxis(massOmegaAxis); + printAxis(massXiAxis); + printAxis(massLambdaAxis); + printAxis(momAxis); + printValue(zVtxMax); + printValue(etaMax); + printValue(cascPtMin); + printValue(cascPtMax); + printValue(minNCrossedRows); + printValue(minNITSClus); + printValue(minNTPCClus); + printValue(maxNSharedTPCClus); + printValue(nSigmaTPCCut); + printValue(dcaBachToPV); + printValue(dcaKaonToPV); + printValue(dcaOmegaToPV); + printValue(dcaDirectLambdaToPV); + printValue(dcaV0DauToPV); + printValue(dcaV0Bach); + printValue(dcaLambdaDaughters); + printValue(mXiWindow); + printValue(mOmegaWindow); + printValue(mLambdaWindow); + printValue(maxDoubleOmegaMass); + printValue(minCosPAOmega); + printValue(minCosPADirectLambda); + printValue(minCosPADoubleOmega); + printValue(minDoubleOmegaDecayRadius); + printValue(doprocessData); + printValue(doprocessDataFromCascades); + printValue(doprocessMC); + printValue(doprocessFindableTracks); + } + void init(InitContext const&) { + printConfig(); + ccdb->setURL(ccdburl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -305,12 +386,11 @@ struct doubleOmegaTreeCreator { return std::sqrt(std::max(0.f, massSquared)); } - template - bool buildLambda(T const&, V0 const& v0, int8_t charge, BuiltLambdaCandidate& builtLambda) + template + bool buildLambda(Track const& posTrack, Track const& negTrack, int8_t charge, BuiltLambdaCandidate& builtLambda) { - auto posTrack = v0.template posTrack_as(); - auto negTrack = v0.template negTrack_as(); - if (!selectTrack(posTrack) || !selectTrack(negTrack)) { + if (posTrack.sign() <= 0 || negTrack.sign() >= 0 || + !selectTrack(posTrack) || !selectTrack(negTrack)) { return false; } @@ -319,6 +399,7 @@ struct doubleOmegaTreeCreator { const float pionNSigma = isMatter ? negTrack.tpcNSigmaPi() : posTrack.tpcNSigmaPi(); if (std::abs(protonNSigma) > nSigmaTPCCut || std::abs(pionNSigma) > nSigmaTPCCut || std::abs(posTrack.dcaXY()) < dcaV0DauToPV || std::abs(negTrack.dcaXY()) < dcaV0DauToPV) { + LOG(debug) << "Lambda dau don't pass selections: " << "nsigmaPR: " << protonNSigma << ", pionNSigma: " << pionNSigma << ", posTrack.dcaXY(): " << posTrack.dcaXY() << ", negTrack.dcaXY()" << negTrack.dcaXY(); return false; } @@ -328,7 +409,7 @@ struct doubleOmegaTreeCreator { try { nCandidates = fitter.process(posTrackParCov, negTrackParCov); } catch (...) { - LOG(error) << "Exception while fitting raw V0 " << v0.globalIndex(); + LOG(error) << "Exception while fitting Lambda daughter tracks " << posTrack.globalIndex() << ", " << negTrack.globalIndex(); return false; } if (nCandidates == 0) { @@ -365,87 +446,17 @@ struct doubleOmegaTreeCreator { .mass = mass}; builtLambda.parentTrack = fitter.createParentTrackParCov(0); builtLambda.decayVertex = decayVertex; - builtLambda.v0Id = v0.globalIndex(); - builtLambda.posTrackId = v0.posTrackId(); - builtLambda.negTrackId = v0.negTrackId(); - return true; - } - - template - bool buildOmega(T const& tracks, Casc const& cascade, BuiltOmegaCandidate& builtOmega) - { - auto bachelor = cascade.template bachelor_as(); - if (bachelor.sign() == 0 || - !selectTrack(bachelor) || - std::abs(bachelor.tpcNSigmaKa()) > nSigmaTPCCut || - std::abs(bachelor.dcaXY()) < dcaBachToPV) { - return false; - } - - auto v0 = cascade.template v0_as(); - BuiltLambdaCandidate lambda; - if (!buildLambda(tracks, v0, bachelor.sign(), lambda)) { - return false; - } - - auto v0TrackParCov = lambda.parentTrack; - auto bachelorTrackParCov = getTrackParCov(bachelor); - int nCandidates = 0; - try { - nCandidates = fitter.process(v0TrackParCov, bachelorTrackParCov); - } catch (...) { - LOG(error) << "Exception while fitting raw cascade " << cascade.globalIndex(); - return false; - } - if (nCandidates == 0) { - return false; - } - - std::array lambdaMomentum{}; - std::array bachelorMomentum{}; - fitter.getTrack(0).getPxPyPzGlo(lambdaMomentum); - fitter.getTrack(1).getPxPyPzGlo(bachelorMomentum); - const std::array omegaMomentum{ - lambdaMomentum[0] + bachelorMomentum[0], - lambdaMomentum[1] + bachelorMomentum[1], - lambdaMomentum[2] + bachelorMomentum[2]}; - const auto& fittedVertex = fitter.getPCACandidate(); - const std::array decayVertex{ - static_cast(fittedVertex[0]), static_cast(fittedVertex[1]), static_cast(fittedVertex[2])}; - const float massOmega = invariantMass2Body(lambdaMomentum, o2::constants::physics::MassLambda0, - bachelorMomentum, o2::constants::physics::MassKaonCharged); - const float massXi = invariantMass2Body(lambdaMomentum, o2::constants::physics::MassLambda0, - bachelorMomentum, o2::constants::physics::MassPionCharged); - const float dcaDaughters = std::sqrt(std::abs(fitter.getChi2AtPCACandidate())); - const float eta = etaFromMomentum(omegaMomentum[0], omegaMomentum[1], omegaMomentum[2]); - if (dcaDaughters > dcaV0Bach || - std::abs(eta) > etaMax || - std::abs(massOmega - o2::constants::physics::MassOmegaMinus) > mOmegaWindow) { - return false; - } - - builtOmega = { - omegaMomentum[0], - omegaMomentum[1], - omegaMomentum[2], - decayVertex[0], - decayVertex[1], - decayVertex[2], - massOmega, - massXi, - static_cast(bachelor.sign()), - cascade.globalIndex(), - v0.posTrackId(), - v0.negTrackId(), - cascade.bachelorId(), - fitter.createParentTrackParCov()}; + builtLambda.posTrackId = posTrack.globalIndex(); + builtLambda.negTrackId = negTrack.globalIndex(); + // Keep particle-role ordering for matter and antimatter alike. + builtLambda.daughterDetectorMap = isMatter ? detectorMap(posTrack) | (detectorMap(negTrack) << 3) + : detectorMap(negTrack) | (detectorMap(posTrack) << 3); return true; } - template - bool buildOmegaFromBuilder(T const& tracks, Casc const& cascade, BuiltOmegaCandidate& builtOmega) + template + bool buildOmega(Track const& posTrack, Track const& negTrack, Track const& bachelor, BuiltOmegaCandidate& builtOmega) { - auto bachelor = cascade.template bachelor_as(); if (bachelor.sign() == 0 || !selectTrack(bachelor) || std::abs(bachelor.tpcNSigmaKa()) > nSigmaTPCCut || @@ -454,7 +465,7 @@ struct doubleOmegaTreeCreator { } BuiltLambdaCandidate lambda; - if (!buildLambda(tracks, cascade, bachelor.sign(), lambda)) { + if (!buildLambda(posTrack, negTrack, bachelor.sign(), lambda)) { return false; } @@ -464,7 +475,7 @@ struct doubleOmegaTreeCreator { try { nCandidates = fitter.process(v0TrackParCov, bachelorTrackParCov); } catch (...) { - LOG(error) << "Exception while refitting cascade " << cascade.globalIndex(); + LOG(error) << "Exception while fitting Omega daughter tracks " << posTrack.globalIndex() << ", " << negTrack.globalIndex() << ", " << bachelor.globalIndex(); return false; } if (nCandidates == 0) { @@ -504,11 +515,11 @@ struct doubleOmegaTreeCreator { massOmega, massXi, static_cast(bachelor.sign()), - cascade.globalIndex(), - cascade.posTrackId(), - cascade.negTrackId(), - cascade.bachelorId(), - fitter.createParentTrackParCov()}; + posTrack.globalIndex(), + negTrack.globalIndex(), + bachelor.globalIndex(), + fitter.createParentTrackParCov(), + detectorMap(bachelor) | (lambda.daughterDetectorMap << 3)}; return true; } @@ -580,6 +591,9 @@ struct doubleOmegaTreeCreator { cand.x = decayVertex[0]; cand.y = decayVertex[1]; cand.z = decayVertex[2]; + cand.pvX = primaryVertex[0]; + cand.pvY = primaryVertex[1]; + cand.pvZ = primaryVertex[2]; cand.cosPAOmega = RecoDecay::cpa(decayVertex, omegaDecayVertex, omegaMomentum); cand.cosPADirectLambda = RecoDecay::cpa(decayVertex, directLambda.decayVertex, lambdaMomentum); cand.cosPADoubleOmega = RecoDecay::cpa(primaryVertex, decayVertex, totalMomentum); @@ -589,6 +603,11 @@ struct doubleOmegaTreeCreator { cand.dcaZDirectLambdaToPV = dcaDirectLambda[1]; cand.dcaXYDirectKaonToPV = dcaDirectKaon[0]; cand.dcaZDirectKaonToPV = dcaDirectKaon[1]; + cand.tpcNSigmaDirectKaon = kaon.hasTPC() ? kaon.tpcNSigmaKa() : -999.f; + cand.tofNSigmaDirectKaon = kaon.hasTOF() ? kaon.tofNSigmaKa() : -999.f; + cand.daughterDetectorMap = omega.daughterDetectorMap | + (directLambda.daughterDetectorMap << 9) | + (detectorMap(kaon) << 15); cand.mass = std::sqrt(std::max(0.f, massSquared)); if (cand.cosPAOmega < minCosPAOmega || cand.cosPADirectLambda < minCosPADirectLambda || @@ -726,11 +745,17 @@ struct doubleOmegaTreeCreator { bool foundLambda = false; bool foundKaon = false; std::array decayVertex{0.f, 0.f, 0.f}; + std::array omegaDecayVertex{-999.f, -999.f, -999.f}; for (const auto& daughter : particle.template daughters_as()) { if (daughter.pdgCode() == expectedOmegaPdg) { foundOmega = true; decayVertex = {daughter.vx(), daughter.vy(), daughter.vz()}; + // The Omega decay point is the production point of any of its daughters. + for (const auto& omegaDaughter : daughter.template daughters_as()) { + omegaDecayVertex = {omegaDaughter.vx(), omegaDaughter.vy(), omegaDaughter.vz()}; + break; + } } else if (daughter.pdgCode() == expectedLambdaPdg) { foundLambda = true; } else if (daughter.pdgCode() == expectedKaonPdg) { @@ -749,6 +774,9 @@ struct doubleOmegaTreeCreator { decayVertex[1] - particle.vy(), decayVertex[2] - particle.vz()); info.pdgCode = particle.pdgCode(); + info.genProductionVertex = {particle.vx(), particle.vy(), particle.vz()}; + info.genDecayVertex = decayVertex; + info.genOmegaDecayVertex = omegaDecayVertex; return true; } @@ -777,6 +805,9 @@ struct doubleOmegaTreeCreator { cand.x, cand.y, cand.z, + cand.pvX, + cand.pvY, + cand.pvZ, cand.cosPAOmega, cand.cosPADirectLambda, cand.cosPADoubleOmega, @@ -786,6 +817,9 @@ struct doubleOmegaTreeCreator { cand.dcaZDirectLambdaToPV, cand.dcaXYDirectKaonToPV, cand.dcaZDirectKaonToPV, + cand.tpcNSigmaDirectKaon, + cand.tofNSigmaDirectKaon, + cand.daughterDetectorMap, cand.mass, cand.massOmega, cand.massXi); @@ -799,6 +833,9 @@ struct doubleOmegaTreeCreator { cand.x, cand.y, cand.z, + cand.pvX, + cand.pvY, + cand.pvZ, cand.cosPAOmega, cand.cosPADirectLambda, cand.cosPADoubleOmega, @@ -808,6 +845,9 @@ struct doubleOmegaTreeCreator { cand.dcaZDirectLambdaToPV, cand.dcaXYDirectKaonToPV, cand.dcaZDirectKaonToPV, + cand.tpcNSigmaDirectKaon, + cand.tofNSigmaDirectKaon, + cand.daughterDetectorMap, cand.mass, cand.massOmega, cand.massXi, @@ -816,16 +856,42 @@ struct doubleOmegaTreeCreator { mcInfo.phi, mcInfo.decayLength, mcInfo.pdgCode, + mcInfo.motherId, isReco); } + template + void logVertexComparison(C const& collision, DoubleOmegaCandidate const& cand, + BuiltOmegaCandidate const& omega, DoubleOmegaMCInfo const& mcInfo) + { + auto logVertex = [](const char* label, std::array const& reco, std::array const& gen) { + LOG(debug) << label + << ": reco = (" << reco[0] << ", " << reco[1] << ", " << reco[2] << ") cm" + << " | gen = (" << gen[0] << ", " << gen[1] << ", " << gen[2] << ") cm" + << " | reco - gen = (" << reco[0] - gen[0] << ", " << reco[1] - gen[1] << ", " << reco[2] - gen[2] << ") cm"; + }; + LOG(debug) << "---------------------------------------- vertices, mother " << mcInfo.motherId + << ", pdgCode = " << mcInfo.pdgCode << ", gen pt = " << mcInfo.pt << " GeV/c"; + logVertex("PV ", {cand.pvX, cand.pvY, cand.pvZ}, mcInfo.genProductionVertex); + logVertex("doubleOmega ", {cand.x, cand.y, cand.z}, mcInfo.genDecayVertex); + logVertex("Omega ", {omega.x, omega.y, omega.z}, mcInfo.genOmegaDecayVertex); + LOG(debug) << "decay length: reco = " + << std::hypot(cand.x - cand.pvX, cand.y - cand.pvY, cand.z - cand.pvZ) + << " cm | gen = " << mcInfo.decayLength << " cm"; + LOG(debug) << "collision " << collision.globalIndex() << " (direct kaon), numContrib = " << collision.numContrib(); + } + template void fillFromRawTables(C const& collision, T const& tracks, Cascades const& cascades, V0s const& v0s, aod::McParticles const* mcParticles = nullptr) { for (const auto& rawCascade : cascades) { + const auto omegaV0 = rawCascade.template v0_as(); + const auto omegaPosTrack = omegaV0.template posTrack_as(); + const auto omegaNegTrack = omegaV0.template negTrack_as(); + const auto omegaBachelorTrack = rawCascade.template bachelor_as(); BuiltOmegaCandidate omega; - if (!buildOmega(tracks, rawCascade, omega)) { + if (!buildOmega(omegaPosTrack, omegaNegTrack, omegaBachelorTrack, omega)) { continue; } const float omegaPt = std::hypot(omega.px, omega.py); @@ -834,8 +900,10 @@ struct doubleOmegaTreeCreator { const std::array omegaTrackIds{omega.posTrackId, omega.negTrackId, omega.bachelorId}; for (const auto& rawV0 : v0s) { + const auto lambdaPosTrack = rawV0.template posTrack_as(); + const auto lambdaNegTrack = rawV0.template negTrack_as(); BuiltLambdaCandidate lambda; - if (!buildLambda(tracks, rawV0, omega.sign, lambda)) { + if (!buildLambda(lambdaPosTrack, lambdaNegTrack, omega.sign, lambda)) { continue; } if (std::find(omegaTrackIds.begin(), omegaTrackIds.end(), lambda.posTrackId) != omegaTrackIds.end() || @@ -859,12 +927,6 @@ struct doubleOmegaTreeCreator { if (mcParticles == nullptr) { continue; } - auto omegaV0 = rawCascade.template v0_as(); - auto omegaPosTrack = omegaV0.template posTrack_as(); - auto omegaNegTrack = omegaV0.template negTrack_as(); - auto omegaBachelorTrack = rawCascade.template bachelor_as(); - auto lambdaPosTrack = rawV0.template posTrack_as(); - auto lambdaNegTrack = rawV0.template negTrack_as(); DoubleOmegaMCInfo mcInfo; if (getMCInfo(mcInfo, omegaPosTrack, omegaNegTrack, omegaBachelorTrack, @@ -882,11 +944,11 @@ struct doubleOmegaTreeCreator { } template - void fillFromCascades(C const& collision, T const& tracks, Cascades const& cascades) + void fillFromCascades(C const& collision, T const&, Cascades const& cascades) { for (const auto& omegaRow : cascades) { BuiltOmegaCandidate omega; - if (!buildOmegaFromBuilder(tracks, omegaRow, omega)) { + if (!buildOmega(omegaRow.template posTrack_as(), omegaRow.template negTrack_as(), omegaRow.template bachelor_as(), omega)) { continue; } const float omegaPt = std::hypot(omega.px, omega.py); @@ -918,7 +980,7 @@ struct doubleOmegaTreeCreator { continue; } BuiltLambdaCandidate lambda; - if (!buildLambda(tracks, lambdaKaonSource, omega.sign, lambda)) { + if (!buildLambda(lambdaKaonSource.template posTrack_as(), lambdaKaonSource.template negTrack_as(), omega.sign, lambda)) { continue; } histos.fill(HIST("QA/massLambda"), std::hypot(lambda.candidate.px, lambda.candidate.py), lambda.candidate.mass); @@ -930,6 +992,63 @@ struct doubleOmegaTreeCreator { } } + // Always take the vertex the direct kaon is associated with. To be improved + template + int64_t selectFindableCollision(Kaon const& kaon) + { + return kaon.collisionId(); + } + + template + bool fillFromSelectedTracks(C const& collisions, T const& tracks, + std::array const& daughterTrackIds, + DoubleOmegaMCInfo const& mcInfo) + { + // Daughter slots: Omega kaon, Omega proton/pion, direct proton/pion, direct kaon. + // Swap proton and pion slots for antimatter to pass positive/negative tracks. + const bool isMatter = mcInfo.pdgCode > 0; + const int8_t charge = isMatter ? -1 : 1; + const auto bachelor = tracks.rawIteratorAt(daughterTrackIds[0]); + if (bachelor.sign() != charge) { + return false; + } + const auto omegaPosTrack = tracks.rawIteratorAt(daughterTrackIds[isMatter ? 1 : 2]); + const auto omegaNegTrack = tracks.rawIteratorAt(daughterTrackIds[isMatter ? 2 : 1]); + const auto lambdaPosTrack = tracks.rawIteratorAt(daughterTrackIds[isMatter ? 3 : 4]); + const auto lambdaNegTrack = tracks.rawIteratorAt(daughterTrackIds[isMatter ? 4 : 3]); + const auto kaon = tracks.rawIteratorAt(daughterTrackIds[5]); + if (!isSelectedKaon(kaon, charge)) { + return false; + } + const int64_t collisionId = selectFindableCollision(kaon); + if (collisionId < 0) { + return false; + } + const auto collision = collisions.rawIteratorAt(collisionId); + const auto bc = collision.template bc_as(); + initCCDB(bc); + + BuiltOmegaCandidate omega; + if (!buildOmega(omegaPosTrack, omegaNegTrack, bachelor, omega)) { + return false; + } + const float omegaPt = std::hypot(omega.px, omega.py); + histos.fill(HIST("QA/massXi"), omegaPt, omega.massXi); + histos.fill(HIST("QA/massOmega"), omegaPt, omega.massOmega); + BuiltLambdaCandidate lambda; + if (!buildLambda(lambdaPosTrack, lambdaNegTrack, charge, lambda)) { + return false; + } + histos.fill(HIST("QA/massLambda"), std::hypot(lambda.candidate.px, lambda.candidate.py), lambda.candidate.mass); + DoubleOmegaCandidate cand; + if (!buildDoubleOmega(collision, omega, lambda, kaon, cand)) { + return false; + } + logVertexComparison(collision, cand, omega, mcInfo); + writeMCCandidate(cand, mcInfo, true); + return true; + } + template bool acceptCollision(C const& collision, aod::BCsWithTimestamps const&) { @@ -996,7 +1115,6 @@ struct doubleOmegaTreeCreator { cascadesThisCollision.bindExternalIndices(&v0s); fillFromRawTables(collision, tracksThisCollision, cascadesThisCollision, v0sThisCollision, &mcParticles); } - LOG(info) << "Found " << reconstructedDoubleOmegaIds.size() << " reconstructed double-Omega candidates in MC"; for (const auto& mcParticle : mcParticles) { DoubleOmegaMCInfo mcInfo; if (!getGeneratedMCInfo(mcInfo, mcParticle)) { @@ -1011,10 +1129,13 @@ struct doubleOmegaTreeCreator { } PROCESS_SWITCH(doubleOmegaTreeCreator, processMC, "Reconstruct and truth-match Omega + Lambda + kaon in MC", false); - void processFindableTracks(TracksFullMC const& tracks, + void processFindableTracks(soa::Join const& collisions, + TracksFullIUMC const& tracks, aod::Cascades const& cascades, aod::V0s const& v0s, - aod::McParticles const& mcParticles) + aod::McParticles const& mcParticles, + aod::McCollisions const&, + aod::BCsWithTimestamps const&) { for (const auto& mcParticle : mcParticles) { if (std::abs(mcParticle.pdgCode()) != kDoubleOmegaPdg) { @@ -1080,8 +1201,11 @@ struct doubleOmegaTreeCreator { directKaonId}; std::array reconstructedDaughters{}; std::array selectedReconstructedDaughters{}; - std::array daughterTrackIds{}; + std::array daughterTrackIds{}; daughterTrackIds.fill(-1); + // Match by generated daughter, independently of reconstructed collision association. + std::array selectedTrackIds{}; + selectedTrackIds.fill(-1); for (const auto& track : tracks) { if (!track.has_mcParticle()) { continue; @@ -1095,9 +1219,34 @@ struct doubleOmegaTreeCreator { } if (selectTrack(track)) { selectedReconstructedDaughters[daughterIndex] = true; + if (selectedTrackIds[daughterIndex] < 0) { + selectedTrackIds[daughterIndex] = track.globalIndex(); + } } } } + auto printDaughterTracks = [&](const char* tag) { + LOG(debug) << "---------------------------------------- " << tag; + for (std::size_t iDaughter = 0; iDaughter < daughterTrackIds.size(); ++iDaughter) { + LOG(debug) << "+++++++"; + const auto daughterTrack = tracks.rawIteratorAt(daughterTrackIds[iDaughter]); + const auto daughterMCParticle = daughterTrack.mcParticle_as(); + LOG(debug) << "Dau" << iDaughter + 1 + << ": hasITS = " << daughterTrack.hasITS() + << ", hasTPC = " << daughterTrack.hasTPC() + << ", hasTOF = " << daughterTrack.hasTOF() + << ", hasTRD = " << daughterTrack.hasTRD() + << ", isITSAfterburner = " << daughterTrack.isITSAfterburner(); + LOG(debug) << "Dau" << iDaughter + 1 + << ": eta = " << daughterTrack.eta() + << ", pt = " << daughterTrack.pt() + << ", CollisionID = " << daughterTrack.collisionId() + << ", production point = (" << daughterMCParticle.vx() + << ", " << daughterMCParticle.vy() + << ", " << daughterMCParticle.vz() << ") cm"; + } + }; + bool allDaughtersReconstructed = true; bool allDaughtersSelected = true; for (std::size_t iDaughter = 0; iDaughter < reconstructedDaughters.size(); ++iDaughter) { @@ -1142,33 +1291,24 @@ struct doubleOmegaTreeCreator { if (omegaInCascade && directLambdaInV0) { histos.fill(HIST("MC/generatedAndFindable"), kFindableCascadeAndV0); + printDaughterTracks("findable cascade and V0"); + DoubleOmegaMCInfo mcInfo; + if (!getGeneratedMCInfo(mcInfo, mcParticle)) { + continue; + } + const bool isReco = fillFromSelectedTracks(collisions, tracks, daughterTrackIds, mcInfo); + if (!isReco) { + writeMCCandidate(DoubleOmegaCandidate{}, mcInfo, false); + } } } if (allDaughtersSelected) { histos.fill(HIST("MC/generatedAndFindable"), kFindableSelectedTracks); - LOG(debug) << "----------------------------------------"; - for (std::size_t iDaughter = 0; iDaughter < daughterTrackIds.size(); ++iDaughter) { - LOG(debug) << "+++++++"; - const auto daughterTrack = tracks.rawIteratorAt(daughterTrackIds[iDaughter]); - const auto daughterMCParticle = daughterTrack.mcParticle_as(); - LOG(debug) << "Dau" << iDaughter + 1 - << ": hasITS = " << daughterTrack.hasITS() - << ", hasTPC = " << daughterTrack.hasTPC() - << ", hasTOF = " << daughterTrack.hasTOF() - << ", hasTRD = " << daughterTrack.hasTRD() - << ", isITSAfterburner = " << daughterTrack.isITSAfterburner(); - LOG(debug) << "Dau" << iDaughter + 1 - << ": eta = " << daughterTrack.eta() - << ", pt = " << daughterTrack.pt() - << ", CollisionID = " << daughterTrack.collisionId() - << ", production point = (" << daughterMCParticle.vx() - << ", " << daughterMCParticle.vy() - << ", " << daughterMCParticle.vz() << ") cm"; - } + printDaughterTracks("findable selected tracks"); } } } - PROCESS_SWITCH(doubleOmegaTreeCreator, processFindableTracks, "Count generated and track-findable double-Omega decays", false); + PROCESS_SWITCH(doubleOmegaTreeCreator, processFindableTracks, "Count findable double-Omega decays and reconstruct selected daughter tracks", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)