diff --git a/ALICE3/DataModel/OTFStrangeness.h b/ALICE3/DataModel/OTFStrangeness.h index 3bc6872eb4c..73bbf0fc39c 100644 --- a/ALICE3/DataModel/OTFStrangeness.h +++ b/ALICE3/DataModel/OTFStrangeness.h @@ -20,6 +20,8 @@ #define ALICE3_DATAMODEL_OTFSTRANGENESS_H_ // O2 includes +#include "Common/Core/RecoDecay.h" + #include "Framework/AnalysisDataModel.h" namespace o2::aod @@ -70,7 +72,7 @@ namespace otfv0 DECLARE_SOA_INDEX_COLUMN(Collision, collision); //! DECLARE_SOA_INDEX_COLUMN_FULL(PosTrack, posTrack, int, Tracks, "_Pos"); //! DECLARE_SOA_INDEX_COLUMN_FULL(NegTrack, negTrack, int, Tracks, "_Neg"); //! -DECLARE_SOA_INDEX_COLUMN(V0, v0); //! +DECLARE_SOA_INDEX_COLUMN(V0, v0); //! index of the mc particle corresponding to the V0 // topo vars DECLARE_SOA_COLUMN(DCAV0Daughters, dcaV0Daughters, float); @@ -86,6 +88,7 @@ DECLARE_SOA_COLUMN(Pt, pt, float); DECLARE_SOA_TABLE(UpgradeV0s, "AOD", "UPGRADEV0S", o2::soa::Index<>, otfv0::CollisionId, + otfv0::V0Id, otfv0::PosTrackId, otfv0::NegTrackId, otfv0::DCAV0Daughters, @@ -96,6 +99,152 @@ DECLARE_SOA_TABLE(UpgradeV0s, "AOD", "UPGRADEV0S", otfv0::Pt); using UpgradeV0 = UpgradeV0s::iterator; -} // namespace o2::aod +namespace candidatev0 +{ +DECLARE_SOA_INDEX_COLUMN(Collision, collision); //! +DECLARE_SOA_INDEX_COLUMN_FULL(PosTrack, posTrack, int, Tracks, "_Pos"); //! +DECLARE_SOA_INDEX_COLUMN_FULL(NegTrack, negTrack, int, Tracks, "_Neg"); //! +DECLARE_SOA_INDEX_COLUMN(V0, v0); //! + +// General V0 properties: position, momentum +DECLARE_SOA_COLUMN(PosX, posX, float); //! positive track X at min +DECLARE_SOA_COLUMN(NegX, negX, float); //! negative track X at min +DECLARE_SOA_COLUMN(PxPos, pxpos, float); //! positive track px at min +DECLARE_SOA_COLUMN(PyPos, pypos, float); //! positive track py at min +DECLARE_SOA_COLUMN(PzPos, pzpos, float); //! positive track pz at min +DECLARE_SOA_COLUMN(PxNeg, pxneg, float); //! negative track px at min +DECLARE_SOA_COLUMN(PyNeg, pyneg, float); //! negative track py at min +DECLARE_SOA_COLUMN(PzNeg, pzneg, float); //! negative track pz at min +DECLARE_SOA_COLUMN(X, x, float); //! decay position X +DECLARE_SOA_COLUMN(Y, y, float); //! decay position Y +DECLARE_SOA_COLUMN(Z, z, float); //! decay position Z + +// topo vars +DECLARE_SOA_COLUMN(DCAV0Daughters, dcaV0Daughters, float); +DECLARE_SOA_COLUMN(CosPA, cosPA, float); +DECLARE_SOA_COLUMN(DCAPosToPV, dcaPosToPV, float); +DECLARE_SOA_COLUMN(DCANegToPV, dcaNegToPV, float); +DECLARE_SOA_COLUMN(DCAV0ToPV, dcaV0ToPV, float); + +//______________________________________________________ +// DYNAMIC COLUMNS + +DECLARE_SOA_DYNAMIC_COLUMN(Px, px, //! V0 px + [](float pxPos, float pxNeg) -> float { return pxPos + pxNeg; }); +DECLARE_SOA_DYNAMIC_COLUMN(Py, py, //! V0 py + [](float pyPos, float pyNeg) -> float { return pyPos + pyNeg; }); +DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, //! V0 pz + [](float pzPos, float pzNeg) -> float { return pzPos + pzNeg; }); +DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, //! Transverse momentum in GeV/c + [](float pxPos, float pyPos, float pxNeg, float pyNeg) -> float { + return RecoDecay::sqrtSumOfSquares(pxPos + pxNeg, pyPos + pyNeg); + }); +DECLARE_SOA_DYNAMIC_COLUMN(P, p, //! Total momentum in GeV/c + [](float pxPos, float pyPos, float pzPos, float pxNeg, float pyNeg, float pzNeg) -> float { + return RecoDecay::sqrtSumOfSquares(pxPos + pxNeg, pyPos + pyNeg, pzPos + pzNeg); + }); +DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, //! Phi in the range [0, 2pi) + [](float pxPos, float pyPos, float pxNeg, float pyNeg) -> float { return RecoDecay::phi(pxPos + pxNeg, pyPos + pyNeg); }); +DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, //! Pseudorapidity, conditionally defined to avoid FPEs + [](float pxPos, float pyPos, float pzPos, float pxNeg, float pyNeg, float pzNeg) -> float { + return RecoDecay::eta(std::array{pxPos + pxNeg, pyPos + pyNeg, pzPos + pzNeg}); + }); +// Length quantities +DECLARE_SOA_DYNAMIC_COLUMN(V0Radius, v0radius, //! V0 decay radius (2D, centered at zero) + [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); + +// Distance Over To Mom +DECLARE_SOA_DYNAMIC_COLUMN(DistOverTotMom, distovertotmom, //! PV to V0decay distance over total momentum + [](float X, float Y, float Z, float pxPos, float pyPos, float pzPos, float pxNeg, float pyNeg, float pzNeg, float pvX, float pvY, float pvZ) { + float P = RecoDecay::sqrtSumOfSquares(pxPos + pxNeg, pyPos + pyNeg, pzPos + pzNeg); + return std::sqrt(std::pow(X - pvX, 2) + std::pow(Y - pvY, 2) + std::pow(Z - pvZ, 2)) / (P + 1E-10); + }); + +// Armenteros-Podolanski variables +DECLARE_SOA_DYNAMIC_COLUMN(Alpha, alpha, //! Armenteros Alpha + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) { + float momTot = RecoDecay::p(pxpos + pxneg, pypos + pyneg, pzpos + pzneg); + float lQlNeg = RecoDecay::dotProd(std::array{pxneg, pyneg, pzneg}, std::array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}) / momTot; + float lQlPos = RecoDecay::dotProd(std::array{pxpos, pypos, pzpos}, std::array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}) / momTot; + return (lQlPos - lQlNeg) / (lQlPos + lQlNeg); // alphav0 + }); + +DECLARE_SOA_DYNAMIC_COLUMN(QtArm, qtarm, //! Armenteros Qt + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) { + float momTot = RecoDecay::p2(pxpos + pxneg, pypos + pyneg, pzpos + pzneg); + float dp = RecoDecay::dotProd(std::array{pxneg, pyneg, pzneg}, std::array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}); + return std::sqrt(RecoDecay::p2(pxneg, pyneg, pzneg) - dp * dp / momTot); // qtarm + }); +// Mass assumption +DECLARE_SOA_DYNAMIC_COLUMN(MLambda, mLambda, //! mass under lambda hypothesis + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) -> float { + return RecoDecay::m(std::array{std::array{pxpos, pypos, pzpos}, std::array{pxneg, pyneg, pzneg}}, std::array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}); + }); +DECLARE_SOA_DYNAMIC_COLUMN(MAntiLambda, mAntiLambda, //! mass under antilambda hypothesis + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) -> float { + return RecoDecay::m(std::array{std::array{pxpos, pypos, pzpos}, std::array{pxneg, pyneg, pzneg}}, std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton}); + }); +DECLARE_SOA_DYNAMIC_COLUMN(MK0Short, mK0Short, //! mass under K0short hypothesis + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) -> float { + return RecoDecay::m(std::array{std::array{pxpos, pypos, pzpos}, std::array{pxneg, pyneg, pzneg}}, std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassPionCharged}); + }); +// Rapidity +DECLARE_SOA_DYNAMIC_COLUMN(YK0Short, yK0Short, //! V0 y with K0short hypothesis + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) -> float { + return RecoDecay::y(std::array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}, o2::constants::physics::MassKaonNeutral); + }); +DECLARE_SOA_DYNAMIC_COLUMN(YLambda, yLambda, //! V0 y with lambda or antilambda hypothesis + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) -> float { + return RecoDecay::y(std::array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg}, o2::constants::physics::MassLambda); + }); +// Daughter track momenta +DECLARE_SOA_DYNAMIC_COLUMN(NegativePt, negativept, //! negative daughter pT + [](float pxneg, float pyneg) -> float { return RecoDecay::sqrtSumOfSquares(pxneg, pyneg); }); +DECLARE_SOA_DYNAMIC_COLUMN(PositivePt, positivept, //! positive daughter pT + [](float pxpos, float pypos) -> float { return RecoDecay::sqrtSumOfSquares(pxpos, pypos); }); +DECLARE_SOA_DYNAMIC_COLUMN(NegativeEta, negativeeta, //! negative daughter eta + [](float PxNeg, float PyNeg, float PzNeg) -> float { return RecoDecay::eta(std::array{PxNeg, PyNeg, PzNeg}); }); +DECLARE_SOA_DYNAMIC_COLUMN(NegativePhi, negativephi, //! negative daughter phi + [](float PxNeg, float PyNeg) -> float { return RecoDecay::phi(PxNeg, PyNeg); }); +DECLARE_SOA_DYNAMIC_COLUMN(PositiveEta, positiveeta, //! positive daughter eta + [](float PxPos, float PyPos, float PzPos) -> float { return RecoDecay::eta(std::array{PxPos, PyPos, PzPos}); }); +DECLARE_SOA_DYNAMIC_COLUMN(PositivePhi, positivephi, //! positive daughter phi + [](float PxPos, float PyPos) -> float { return RecoDecay::phi(PxPos, PyPos); }); +} // namespace candidatev0 +DECLARE_SOA_TABLE(V0CandidateIndices, "AOD", "V0CANDIDATEINDEX", //! index table + o2::soa::Index<>, candidatev0::CollisionId, candidatev0::PosTrackId, candidatev0::NegTrackId); + +DECLARE_SOA_TABLE(V0CandidateCores, "AOD", "V0CANDIDATECORE", + o2::soa::Index<>, + candidatev0::X, candidatev0::Y, candidatev0::Z, + candidatev0::PxPos, candidatev0::PyPos, candidatev0::PzPos, + candidatev0::PxNeg, candidatev0::PyNeg, candidatev0::PzNeg, + candidatev0::DCAV0Daughters, candidatev0::DCAPosToPV, candidatev0::DCANegToPV, + candidatev0::CosPA, candidatev0::DCAV0ToPV, + candidatev0::Px, + candidatev0::Py, + candidatev0::Pz, + candidatev0::Pt, + candidatev0::P, + candidatev0::Phi, + candidatev0::Eta, + candidatev0::V0Radius, + candidatev0::DistOverTotMom, + candidatev0::Alpha, + candidatev0::QtArm, + candidatev0::MLambda, + candidatev0::MAntiLambda, + candidatev0::MK0Short, + candidatev0::YK0Short, + candidatev0::YLambda, + candidatev0::NegativePt, + candidatev0::PositivePt, + candidatev0::NegativeEta, + candidatev0::NegativePhi, + candidatev0::PositiveEta, + candidatev0::PositivePhi); + +using V0CandidateCore = V0CandidateCores::iterator; +} // namespace o2::aod #endif // ALICE3_DATAMODEL_OTFSTRANGENESS_H_ diff --git a/ALICE3/TableProducer/CMakeLists.txt b/ALICE3/TableProducer/CMakeLists.txt index 61f46fba7cd..a3ec83bd851 100644 --- a/ALICE3/TableProducer/CMakeLists.txt +++ b/ALICE3/TableProducer/CMakeLists.txt @@ -60,3 +60,8 @@ o2physics_add_dpl_workflow(alice3-tracking-translator SOURCES alice3TrackingTranslator.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(alice3-strangeness-finder + SOURCES alice3-strangenessFinder.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 0239be242d3..4e10b6a5779 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -248,6 +248,7 @@ struct OnTheFlyTracker { struct v0candidate { int positiveId; // track index in the Tracks table int negativeId; // track index in the Tracks table + int mcParticleId; // mc particle index float pt; @@ -1193,6 +1194,7 @@ struct OnTheFlyTracker { // n-2: negative Track from V0 thisV0.positiveId = lastTrackIndex + tracksAlice3.size() - 1; thisV0.negativeId = lastTrackIndex + tracksAlice3.size() - 2; + thisV0.mcParticleId = mcParticle.globalIndex(); // use DCA fitters int nCand = 0; bool dcaFitterOK_V0 = true; @@ -1535,6 +1537,7 @@ struct OnTheFlyTracker { // populate V0s for (const auto& v0 : v0sAlice3) { tableUpgradeV0s(tableCollisions.lastIndex(), // now we know the collision index -> populate table + v0.mcParticleId, v0.positiveId, v0.negativeId, v0.dcaV0dau, diff --git a/ALICE3/TableProducer/alice3-strangenessFinder.cxx b/ALICE3/TableProducer/alice3-strangenessFinder.cxx new file mode 100644 index 00000000000..29e05d5e7e7 --- /dev/null +++ b/ALICE3/TableProducer/alice3-strangenessFinder.cxx @@ -0,0 +1,197 @@ +// 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 alice3-strangenessFinder.cxx +/// +/// \brief finding of V0 and cascade candidates for ALICE 3 +/// +/// This task finds and build condidates for strange hadrons (K0s, Lambda, AntiLambda, Xi-, Xi+, Omega-, Omega+) +/// using the output of the on-the-fly tracker. +/// +/// \author Lucia Anna Tarasovičová, Pavol Jozef Šafárik University (SK) +/// + +#include "ALICE3/DataModel/OTFPIDTrk.h" +#include "ALICE3/DataModel/OTFRICH.h" +#include "ALICE3/DataModel/OTFStrangeness.h" +#include "ALICE3/DataModel/OTFTOF.h" +#include "ALICE3/DataModel/tracksAlice3.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "DCAFitter/DCAFitterN.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" +#include + +#include + +using namespace o2; +// using namespace o2::analysis; +using namespace o2::framework; +using namespace o2::constants::physics; + +using Alice3TracksWPid = soa::Join; +using Alice3Tracks = soa::Join; + +struct alice3strangenessFinder { + SliceCache cache; + + Produces v0CandidateIndices; // contains V0 candidate indices + Produces v0CandidateCores; // contains V0 candidate core information + + Configurable mcSameMotherCheck{"mcSameMotherCheck", true, "check if tracks come from the same MC mother"}; + + Configurable nSigmaTOF{"nSigmaTOF", 5.0f, "Nsigma for TOF PID (if enabled)"}; + Configurable dcaXYconstant{"dcaXYconstant", -1.0f, "[0] in |DCAxy| > [0]+[1]/pT"}; + Configurable dcaXYpTdep{"dcaXYpTdep", 0.0, "[1] in |DCAxy| > [0]+[1]/pT"}; + + o2::vertexing::DCAFitterN<2> fitter; + o2::vertexing::DCAFitterN<3> fitter3; + + // partitions for D mesons + Partition positiveSecondaryTracks = + aod::track::signed1Pt > 0.0f && nabs(aod::track::dcaXY) > dcaXYconstant + dcaXYpTdep* nabs(aod::track::signed1Pt); + Partition negativeSecondaryTracks = + aod::track::signed1Pt < 0.0f && nabs(aod::track::dcaXY) > dcaXYconstant + dcaXYpTdep* nabs(aod::track::signed1Pt); + // Partition negativeSecondaryPions = nabs(aod::upgrade_tof::nSigmaPionInnerTOF) < nSigmaTOF && nabs(aod::upgrade_tof::nSigmaPionOuterTOF) < nSigmaTOF && aod::track::signed1Pt < 0.0f && nabs(aod::track::dcaXY) > dcaXYconstant + dcaXYpTdep* nabs(aod::track::signed1Pt); + // Partition positiveSecondaryPions = nabs(aod::upgrade_tof::nSigmaPionInnerTOF) < nSigmaTOF && nabs(aod::upgrade_tof::nSigmaPionOuterTOF) < nSigmaTOF && aod::track::signed1Pt > 0.0f && nabs(aod::track::dcaXY) > dcaXYconstant + dcaXYpTdep* nabs(aod::track::signed1Pt); + // Partition secondaryProtons = nabs(aod::upgrade_tof::nSigmaProtonInnerTOF) < nSigmaTOF && nabs(aod::upgrade_tof::nSigmaProtonOuterTOF) < nSigmaTOF && aod::track::signed1Pt > 0.0f && nabs(aod::track::dcaXY) > dcaXYconstant + dcaXYpTdep* nabs(aod::track::signed1Pt); + // Partition secondaryAntiProtons = nabs(aod::upgrade_tof::nSigmaProtonInnerTOF) < nSigmaTOF && nabs(aod::upgrade_tof::nSigmaProtonOuterTOF) < nSigmaTOF && aod::track::signed1Pt < 0.0f && nabs(aod::track::dcaXY) > dcaXYconstant + dcaXYpTdep* nabs(aod::track::signed1Pt); + + struct { + float dcaDau; + std::array posSV; + std::array P; + std::array Pdaug; // positive track + std::array Ndaug; // negative track + float cosPA; + float dcaToPV; + } v0cand; + + void init(InitContext&) + { + // Initialization code here + } + /// function to check if tracks have the same mother in MC + template + bool checkSameMother(TTrackType const& track1, TTrackType const& track2) + { + bool sameMother = false; + if (!track1.has_mcParticle() || !track2.has_mcParticle()) + return sameMother; + auto mcParticle1 = track1.template mcParticle_as(); + auto mcParticle2 = track2.template mcParticle_as(); + if (!mcParticle1.has_mothers() || !mcParticle2.has_mothers()) + return sameMother; + for (auto& mcParticleMother1 : mcParticle1.template mothers_as()) { + for (auto& mcParticleMother2 : mcParticle2.template mothers_as()) { + if (mcParticleMother1.globalIndex() == mcParticleMother2.globalIndex()) { + sameMother = true; + } + } + } + return sameMother; + } + + template + bool buildDecayCandidateTwoBody(TTrackType const& posTrack, TTrackType const& negTrack) + { + o2::track::TrackParCov posTrackCov = getTrackParCov(posTrack); + o2::track::TrackParCov negTrackCov = getTrackParCov(negTrack); + + //}-{}-{}-{}-{}-{}-{}-{}-{}-{} + // Move close to minima + int nCand = 0; + try { + nCand = fitter.process(posTrackCov, negTrackCov); + } catch (...) { + return false; + } + if (nCand == 0) { + return false; + } + //}-{}-{}-{}-{}-{}-{}-{}-{}-{} + + posTrackCov = fitter.getTrack(0); + negTrackCov = fitter.getTrack(1); + std::array posP; + std::array negP; + posTrackCov.getPxPyPzGlo(posP); + negTrackCov.getPxPyPzGlo(negP); + v0cand.dcaDau = TMath::Sqrt(fitter.getChi2AtPCACandidate()); + v0cand.Pdaug[0] = posP[0]; + v0cand.Pdaug[1] = posP[1]; + v0cand.Pdaug[2] = posP[2]; + v0cand.Ndaug[0] = negP[0]; + v0cand.Ndaug[1] = negP[1]; + v0cand.Ndaug[2] = negP[2]; + v0cand.P[0] = posP[0] + negP[0]; + v0cand.P[1] = posP[1] + negP[1]; + v0cand.P[2] = posP[2] + negP[2]; + const auto posSV = fitter.getPCACandidate(); + v0cand.posSV[0] = posSV[0]; + v0cand.posSV[1] = posSV[1]; + v0cand.posSV[2] = posSV[2]; + + return true; + } + float CalculateDCAStraightToPV(float X, float Y, float Z, float Px, float Py, float Pz, float pvX, float pvY, float pvZ) + { + return std::sqrt((std::pow((pvY - Y) * Pz - (pvZ - Z) * Py, 2) + std::pow((pvX - X) * Pz - (pvZ - Z) * Px, 2) + std::pow((pvX - X) * Py - (pvY - Y) * Px, 2)) / (Px * Px + Py * Py + Pz * Pz)); + } + void processFindV0CandidateNoPid(aod::Collision const& collision, Alice3Tracks const&, aod::McParticles const&) + { + auto negativeSecondaryTracksGrouped = negativeSecondaryTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + auto positiveSecondaryTracksGrouped = positiveSecondaryTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + + for (auto const& posTrack : positiveSecondaryTracksGrouped) { + for (auto const& negTrack : negativeSecondaryTracksGrouped) { + + if (mcSameMotherCheck && !checkSameMother(posTrack, negTrack)) + continue; + if (!buildDecayCandidateTwoBody(posTrack, negTrack)) + continue; + v0cand.cosPA = RecoDecay::cpa(std::array{collision.posX(), collision.posY(), collision.posZ()}, std::array{v0cand.posSV[0], v0cand.posSV[1], v0cand.posSV[2]}, std::array{v0cand.P[0], v0cand.P[1], v0cand.P[2]}); + v0cand.dcaToPV = CalculateDCAStraightToPV( + v0cand.posSV[0], v0cand.posSV[1], v0cand.posSV[2], + v0cand.P[0], v0cand.P[1], v0cand.P[2], + collision.posX(), collision.posY(), collision.posZ()); + v0CandidateIndices(collision.globalIndex(), + posTrack.globalIndex(), + negTrack.globalIndex()); + v0CandidateCores( + v0cand.posSV[0], v0cand.posSV[1], v0cand.posSV[2], + v0cand.Pdaug[0], v0cand.Pdaug[1], v0cand.Pdaug[2], + v0cand.Ndaug[0], v0cand.Ndaug[1], v0cand.Ndaug[2], + v0cand.dcaDau, posTrack.dcaXY(), negTrack.dcaXY(), + v0cand.cosPA, v0cand.dcaToPV); + } + } + } + // void processFindV0CandidateWithPid(aod::Collision const& collision, aod::McParticles const& mcParticles, Alice3TracksWPid const&) + // { + // auto negativeSecondaryPionsGrouped = negativeSecondaryPions->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + // auto positiveSecondaryPionsGrouped = positiveSecondaryPions->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + // auto secondaryProtonsGrouped = secondaryProtons->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + // auto secondaryAntiProtonsGrouped = secondaryAntiProtons->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + // } + PROCESS_SWITCH(alice3strangenessFinder, processFindV0CandidateNoPid, "find V0 without PID", true); + // PROCESS_SWITCH(alice3strangenessFinder, processFindV0CandidateWithPid, "find V0 with PID", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/ALICE3/Tasks/alice3-strangeness.cxx b/ALICE3/Tasks/alice3-strangeness.cxx index c29e764e9c0..107356fc2ac 100644 --- a/ALICE3/Tasks/alice3-strangeness.cxx +++ b/ALICE3/Tasks/alice3-strangeness.cxx @@ -18,6 +18,7 @@ /// #include "ALICE3/DataModel/OTFStrangeness.h" +#include "ALICE3/DataModel/OTFTracks.h" #include "ALICE3/DataModel/tracksAlice3.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -32,10 +33,13 @@ #include #include #include +#include #include #include +#include #include #include +#include #include #include #include @@ -46,31 +50,124 @@ #include #include +#include + #include #include using namespace o2; using namespace o2::framework; +using namespace o2::constants::math; using alice3tracks = soa::Join; +using fullV0Candidates = soa::Join; +using fullCollisions = soa::Join; struct alice3strangeness { HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "pt axis for QA histograms"}; - ConfigurableAxis axisK0Mass{"axisK0Mass", {200, 0.4f, 0.6f}, ""}; - ConfigurableAxis axisVertexZ{"axisVertexZ", {40, -20, 20}, "vertex Z (cm)"}; + + Configurable idGeometry{"idGeometry", 0, "geometry ID used for propagation"}; + struct ConfigurableGroup { + ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "pt axis for QA histograms"}; + ConfigurableAxis axisK0Mass{"axisK0Mass", {200, 0.4f, 0.6f}, ""}; + ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.08f, 1.2f}, ""}; + ConfigurableAxis axisVertexZ{"axisVertexZ", {40, -20, 20}, "vertex Z (cm)"}; + ConfigurableAxis axisDCA{"axisDCA", {200, 0, 5}, "DCA (cm)"}; + ConfigurableAxis axisV0Radius{"axisV0Radius", {50, 0.0, 100}, "V0 radius (cm)"}; + ConfigurableAxis axisDCAV0Daughters{"axisDCAV0Daughters", {20, 0, 5}, "DCA V0 daughters"}; + ConfigurableAxis axisPointingAngle{"axisPointingAngle", {40, 0.0f, 0.4f}, "pointing angle "}; + ConfigurableAxis axisProperLifeTime{"axisProperLifeTime", {100, 0.0f, 100.0f}, "proper lifetime (cm)"}; + + } histAxes; + + struct : ConfigurableGroup { + Configurable applyRapiditySelection{"applyRapiditySelection", true, "apply rapidity selection"}; + Configurable applyDCAdaughterSelection{"applyDCADaughterSelection", true, "apply DCA daughter selection"}; + Configurable applyCosOfPAngleSelection{"applyCosOfPAngleSelection", true, "apply cosine of pointing angle selection"}; + Configurable applyDCAdaughtersToPVSelection{"applyDCADaughtersToPVSelection", true, "apply DCA daughters to primary vertex selection"}; + Configurable applyV0RadiusSelection{"applyV0RadiusSelection", true, "apply V0 radius selection"}; + Configurable applyArmenterosSelection{"applyArmenterosSelection", true, "apply Armenteros-Podolanski selection"}; + Configurable applyCompetingMassRejection{"applyCompetingMassRejection", true, "apply competing mass rejection"}; + Configurable applyLifetimeSelection{"applyLifetimeSelection", true, "apply lifetime selection"}; + Configurable applyEtaDaughterSelection{"applyEtaDaughterSelection", true, "apply eta daughter selection"}; + Configurable doQAforSelectionVariables{"doQAforSelectionVariables", false, "enable QA plots"}; + } selectionFlags; + + struct : ConfigurableGroup { + Configurable yK0Selection{"yK0Selection", 0.5f, "rapidity selection for K0"}; + Configurable yLambdaSelection{"yLambdaSelection", 0.5f, "rapidity selection for Lambda"}; + Configurable dcaDaughterSelection{"dcaDaughterSelection", 1.0f, "DCA daughter selection"}; + Configurable cosPAngleSelection{"cosPAngleSelection", 0.97f, "cosine of pointing angle selection"}; + Configurable dcaDaughtersToPVSelection{"dcaDaughtersToPVSelection", 0.05f, "DCA daughters to primary vertex selection"}; + Configurable v0RadiusSelection{"v0RadiusSelection", 0.5f, "V0 radius selection"}; + Configurable armenterosSelection{"armenterosSelection", 0.2f, "Armenteros-Podolanski selection in [0]* alpha"}; + Configurable competingMassRejectionK0{"competingMassRejectionK0", 0.005f, "competing mass rejection for K0"}; + Configurable competingMassRejectionLambda{"competingMassRejectionLambda", 0.01f, "competing mass rejection for Lambda"}; + Configurable lifetimecutambda{"lifetimecutambda", 30, "lifetime cut for Lambda in cm"}; + Configurable lifetimecutak0{"lifetimecutak0", 20, "lifetime cut for K0 in cm"}; + Configurable etaDaughterSelection{"etaDaughterSelection", 0.8f, "eta daughter selection"}; + } selectionValues; + + uint16_t AppliedSelectionCheckMask; + double selectionCheck; + double selectionCheckPos; + static constexpr std::string_view kSelectionNames[] = {"DCAV0Daughters", "PointingAngle", "DCAtoPVNegDaughter", "DCAtoPVPosDaughter", "V0Radius", "ProperLifeTime"}; void init(InitContext&) { - histos.add("K0/hMassAllCandidates", "", kTH2D, {axisK0Mass, axisPt}); - histos.add("K0/hMassSelected", "", kTH2D, {axisK0Mass, axisPt}); + histos.add("K0/hMassAllCandidates", "", kTH2D, {histAxes.axisK0Mass, histAxes.axisPt}); + histos.add("K0/hMassSelected", "", kTH2D, {histAxes.axisK0Mass, histAxes.axisPt}); histos.add("K0/hSelections", "", kTH1D, {{10, 0, 10}}); histos.add("K0/hDCANegDaughter", "", kTH1D, {{200, -5, 5}}); histos.add("K0/hDCAPosDaughter", "", kTH1D, {{200, -5, 5}}); - histos.add("hPVz", "hPVz", kTH1F, {axisVertexZ}); + histos.add("hPVz", "hPVz", kTH1F, {histAxes.axisVertexZ}); + histos.add("hV0CandidateCounter", "hV0CandidateCounter", kTH1F, {{11, 0, 11}}); + histos.add("reconstructedCandidates/K0/hMass", "hMass", kTH2D, {histAxes.axisK0Mass, histAxes.axisPt}); + histos.add("reconstructedCandidates/Lambda/hMass", "hMass", kTH2D, {histAxes.axisLambdaMass, histAxes.axisPt}); + histos.add("reconstructedCandidates/hArmeterosBeforeAllSelections", "hArmeterosBeforeAllSelections", kTH2D, {{100, -1.0f, 1.0f}, {100, 0.0f, 0.25f}}); + histos.add("reconstructedCandidates/hArmeterosAfterAllSelections", "hArmeterosAfterAllSelections", kTH2D, {{100, -1.0f, 1.0f}, {100, 0.0f, 0.25f}}); + if (selectionFlags.doQAforSelectionVariables) { + if (!selectionFlags.applyDCAdaughtersToPVSelection) { + histos.add("reconstructedCandidates/K0/hDCAtoPVNegDaughter", "hDCAtoPVNegDaughter", kTH3D, {histAxes.axisK0Mass, histAxes.axisPt, histAxes.axisDCA}); + histos.add("reconstructedCandidates/K0/hDCAtoPVPosDaughter", "hDCAtoPVPosDaughter", kTH3D, {histAxes.axisK0Mass, histAxes.axisPt, histAxes.axisDCA}); + histos.add("reconstructedCandidates/Lambda/hDCAtoPVNegDaughter", "hDCAtoPVNegDaughter", kTH3D, {histAxes.axisLambdaMass, histAxes.axisPt, histAxes.axisDCA}); + histos.add("reconstructedCandidates/Lambda/hDCAtoPVPosDaughter", "hDCAtoPVPosDaughter", kTH3D, {histAxes.axisLambdaMass, histAxes.axisPt, histAxes.axisDCA}); + } + if (!selectionFlags.applyV0RadiusSelection) { + histos.add("reconstructedCandidates/K0/hV0Radius", "hV0Radius", kTH3D, {histAxes.axisK0Mass, histAxes.axisPt, histAxes.axisV0Radius}); + histos.add("reconstructedCandidates/Lambda/hV0Radius", "hV0Radius", kTH3D, {histAxes.axisLambdaMass, histAxes.axisPt, histAxes.axisV0Radius}); + } + if (!selectionFlags.applyDCAdaughterSelection) { + histos.add("reconstructedCandidates/K0/hDCAV0Daughters", "hDCAV0Daughters", kTH3D, {histAxes.axisK0Mass, histAxes.axisPt, histAxes.axisDCAV0Daughters}); + histos.add("reconstructedCandidates/Lambda/hDCAV0Daughters", "hDCAV0Daughters", kTH3D, {histAxes.axisLambdaMass, histAxes.axisPt, histAxes.axisDCAV0Daughters}); + } + if (!selectionFlags.applyCosOfPAngleSelection) { + histos.add("reconstructedCandidates/K0/hPointingAngle", "hPointingAngle", kTH3D, {histAxes.axisK0Mass, histAxes.axisPt, histAxes.axisPointingAngle}); + histos.add("reconstructedCandidates/Lambda/hPointingAngle", "hPointingAngle", kTH3D, {histAxes.axisLambdaMass, histAxes.axisPt, histAxes.axisPointingAngle}); + } + if (!selectionFlags.applyLifetimeSelection) { + histos.add("reconstructedCandidates/K0/hProperLifeTime", "hProperLifeTime", kTH3D, {histAxes.axisK0Mass, histAxes.axisPt, histAxes.axisProperLifeTime}); + histos.add("reconstructedCandidates/Lambda/hProperLifeTime", "hProperLifeTime", kTH3D, {histAxes.axisLambdaMass, histAxes.axisPt, histAxes.axisProperLifeTime}); + } + } + histos.addClone("reconstructedCandidates/Lambda/", "reconstructedCandidates/AntiLambda/"); + + AppliedSelectionCheckMask = 0; + if (!selectionFlags.applyDCAdaughterSelection) + SETBIT(AppliedSelectionCheckMask, 0); + if (!selectionFlags.applyCosOfPAngleSelection) + SETBIT(AppliedSelectionCheckMask, 1); + if (!selectionFlags.applyDCAdaughtersToPVSelection) { + SETBIT(AppliedSelectionCheckMask, 2); + SETBIT(AppliedSelectionCheckMask, 3); + } + if (!selectionFlags.applyV0RadiusSelection) + SETBIT(AppliedSelectionCheckMask, 4); + if (!selectionFlags.applyLifetimeSelection) + SETBIT(AppliedSelectionCheckMask, 5); } long int nEvents = 0; - void process(aod::Collisions const& collisions, aod::McCollisions const& mcCollisions, aod::UpgradeV0s const& v0Recos, alice3tracks const&) + void processAllFindableCandidates(aod::Collisions const& collisions, aod::McCollisions const& mcCollisions, aod::UpgradeV0s const& v0Recos, alice3tracks const&) { LOG(info) << "Event processed " << nEvents++ << " :" << collisions.size() << " " << mcCollisions.size(); for (const auto& collision : collisions) { @@ -108,7 +205,118 @@ struct alice3strangeness { } } } - PROCESS_SWITCH(alice3strangeness, process, "", true); + void processFoundV0Candidates(aod::Collision const& collision, fullV0Candidates const& v0Candidates, alice3tracks const&) + { + // if(collision.lutConfigId()!=idGeometry) + // return; + for (auto const& v0 : v0Candidates) { + bool isK0 = (v0.mK0Short() - o2::constants::physics::MassK0Short) < 0.3; + bool isLambda = (v0.mLambda() - o2::constants::physics::MassLambda0) < 0.2; + bool isAntiLambda = (v0.mAntiLambda() - o2::constants::physics::MassLambda0) < 0.2; + + histos.fill(HIST("reconstructedCandidates/hArmeterosBeforeAllSelections"), v0.alpha(), v0.qtarm()); + histos.fill(HIST("hV0CandidateCounter"), 1); + if (selectionFlags.applyRapiditySelection) { + if (isK0 && std::abs(v0.yK0Short()) < selectionValues.yK0Selection) + continue; + if ((isLambda || isAntiLambda) && std::abs(v0.yLambda()) < selectionValues.yLambdaSelection) + continue; + } + histos.fill(HIST("hV0CandidateCounter"), 2); + if (selectionFlags.applyDCAdaughterSelection) { + if (std::abs(v0.dcaV0Daughters()) > selectionValues.dcaDaughterSelection) + continue; + } else { + selectionCheck = v0.dcaV0Daughters(); + } + histos.fill(HIST("hV0CandidateCounter"), 3); + if (selectionFlags.applyCosOfPAngleSelection) { + if (v0.cosPA() < selectionValues.cosPAngleSelection) + continue; + } else { + selectionCheck = std::acos(v0.cosPA()); + } + histos.fill(HIST("hV0CandidateCounter"), 4); + if (selectionFlags.applyDCAdaughtersToPVSelection) { + if ((std::abs(v0.dcaNegToPV()) < selectionValues.dcaDaughtersToPVSelection) || + (std::abs(v0.dcaPosToPV()) < selectionValues.dcaDaughtersToPVSelection)) + continue; + } else { + selectionCheckPos = std::abs(v0.dcaPosToPV()); + selectionCheck = std::abs(v0.dcaNegToPV()); + } + histos.fill(HIST("hV0CandidateCounter"), 5); + if (selectionFlags.applyV0RadiusSelection) { + if (v0.v0radius() < selectionValues.v0RadiusSelection) + continue; + } else { + selectionCheck = v0.v0radius(); + } + histos.fill(HIST("hV0CandidateCounter"), 6); + if (isK0) { + if (selectionFlags.applyArmenterosSelection) { + if (v0.qtarm() < selectionValues.armenterosSelection * std::abs(v0.alpha())) + continue; + } + } + histos.fill(HIST("hV0CandidateCounter"), 7); + if (isK0 && selectionFlags.applyCompetingMassRejection) { + if (std::abs(v0.mLambda() - o2::constants::physics::MassLambda0) < selectionValues.competingMassRejectionK0) + continue; + if (std::abs(v0.mAntiLambda() - o2::constants::physics::MassLambda0) < selectionValues.competingMassRejectionK0) + continue; + } + if ((isLambda || isAntiLambda) && selectionFlags.applyCompetingMassRejection) { + if (std::abs(v0.mK0Short() - o2::constants::physics::MassK0Short) < selectionValues.competingMassRejectionLambda) + continue; + } + histos.fill(HIST("hV0CandidateCounter"), 8); + if (selectionFlags.applyLifetimeSelection) { + if (isK0 && v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassK0Short < selectionValues.lifetimecutak0) + continue; + if ((isLambda || isAntiLambda) && v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassLambda0 < selectionValues.lifetimecutambda) + continue; + } else { + if (isK0) + selectionCheck = v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassK0Short; + else + selectionCheck = v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassLambda0; + } + histos.fill(HIST("hV0CandidateCounter"), 9); + if (selectionFlags.applyEtaDaughterSelection) { + auto posTrack = v0.template posTrack_as(); + auto negTrack = v0.template negTrack_as(); + if (std::abs(posTrack.eta()) > selectionValues.etaDaughterSelection || std::abs(negTrack.eta()) > selectionValues.etaDaughterSelection) + continue; + } + histos.fill(HIST("hV0CandidateCounter"), 10); + + histos.fill(HIST("reconstructedCandidates/hArmeterosAfterAllSelections"), v0.alpha(), v0.qtarm()); + if (selectionFlags.doQAforSelectionVariables) { + static_for<0, 5>([&](auto i) { + constexpr int In = i.value; + if (TESTBIT(AppliedSelectionCheckMask, In)) { + if (In == 3) + selectionCheck = selectionCheckPos; + if (isK0) + histos.fill(HIST("reconstructedCandidates/K0/h") + HIST(kSelectionNames[In]), v0.mK0Short(), v0.pt(), selectionCheck); + if (isLambda) + histos.fill(HIST("reconstructedCandidates/Lambda/h") + HIST(kSelectionNames[In]), v0.mLambda(), v0.pt(), selectionCheck); + if (isAntiLambda) + histos.fill(HIST("reconstructedCandidates/AntiLambda/h") + HIST(kSelectionNames[In]), v0.mAntiLambda(), v0.pt(), selectionCheck); + } + }); + } + if (isK0) + histos.fill(HIST("reconstructedCandidates/K0/hMass"), v0.mK0Short(), v0.pt()); + if (isLambda) + histos.fill(HIST("reconstructedCandidates/Lambda/hMass"), v0.mLambda(), v0.pt()); + if (isAntiLambda) + histos.fill(HIST("reconstructedCandidates/AntiLambda/hMass"), v0.mAntiLambda(), v0.pt()); + } + } + PROCESS_SWITCH(alice3strangeness, processAllFindableCandidates, "", false); + PROCESS_SWITCH(alice3strangeness, processFoundV0Candidates, "", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)