diff --git a/PWGLF/DataModel/Vtx3BodyTables.h b/PWGLF/DataModel/Vtx3BodyTables.h index 0cc4880dae3..8859d10ec61 100644 --- a/PWGLF/DataModel/Vtx3BodyTables.h +++ b/PWGLF/DataModel/Vtx3BodyTables.h @@ -25,6 +25,7 @@ DECLARE_SOA_INDEX_COLUMN_FULL(Track0, track0, int, Tracks, "_0"); //! DECLARE_SOA_INDEX_COLUMN_FULL(Track1, track1, int, Tracks, "_1"); //! DECLARE_SOA_INDEX_COLUMN_FULL(Track2, track2, int, Tracks, "_2"); //! DECLARE_SOA_INDEX_COLUMN(Collision, collision); //! +DECLARE_SOA_INDEX_COLUMN(Decay3Body, decay3body); //! // General 3 body Vtx properties: position, momentum DECLARE_SOA_COLUMN(PxTrack0, pxtrack0, float); //! track0 px at min @@ -109,7 +110,7 @@ DECLARE_SOA_EXPRESSION_COLUMN(Pz, pz, //! 3 body vtx pz } // namespace vtx3body DECLARE_SOA_TABLE_FULL(StoredVtx3BodyDatas, "Vtx3BodyDatas", "AOD", "Vtx3BodyDATA", //! - o2::soa::Index<>, vtx3body::Track0Id, vtx3body::Track1Id, vtx3body::Track2Id, vtx3body::CollisionId, + o2::soa::Index<>, vtx3body::Track0Id, vtx3body::Track1Id, vtx3body::Track2Id, vtx3body::CollisionId, vtx3body::Decay3BodyId, vtx3body::X, vtx3body::Y, vtx3body::Z, vtx3body::PxTrack0, vtx3body::PyTrack0, vtx3body::PzTrack0, vtx3body::PxTrack1, vtx3body::PyTrack1, vtx3body::PzTrack1, @@ -152,11 +153,31 @@ namespace vtx3body DECLARE_SOA_INDEX_COLUMN(Vtx3BodyData, vtx3BodyData); //! Index to Vtx3BodyData entry } -DECLARE_SOA_TABLE(Vtx3BodyDataLink, "AOD", "VtxDATALINK", //! Joinable table with Decay3bodys which links to Vtx3BodyData which is not produced for all entries +DECLARE_SOA_TABLE(Decay3BodyDataLink, "AOD", "DECAY3BODYLINK", //! Joinable table with Decay3bodys which links to Vtx3BodyData which is not produced for all entries vtx3body::Vtx3BodyDataId); -using Vtxs3BodyLinked = soa::Join; -using Vtx3BodyLinked = Vtxs3BodyLinked::iterator; +using Decay3BodysLinked = soa::Join; +using Decay3BodyLinked = Decay3BodysLinked::iterator; + +// Definition of labels for Vtx3BodyDatas +namespace mcvtx3bodylabel +{ +DECLARE_SOA_INDEX_COLUMN(McParticle, mcParticle); //! MC particle for Vtx3BodyDatas +} // namespace mcvtx3bodylabel + +DECLARE_SOA_TABLE(McVtx3BodyLabels, "AOD", "MCVTXLABEL", //! Table joinable with Vtx3BodyData containing the MC labels + mcvtx3bodylabel::McParticleId); +using McVtx3BodyLabel = McVtx3BodyLabels::iterator; + +// Definition of labels for Decay3Bodys // Full table, joinable with Decay3Bodys (CAUTION: NOT WITH Vtx3BodyDATA) +namespace mcfullvtx3bodylabel +{ +DECLARE_SOA_INDEX_COLUMN(McParticle, mcParticle); //! MC particle for Decay3Bodys +} // namespace mcfullvtx3bodylabel + +DECLARE_SOA_TABLE(McFullVtx3BodyLabels, "AOD", "MCFULLVTXLABEL", //! Table joinable with Decay3Bodys + mcfullvtx3bodylabel::McParticleId); +using McFullVtx3BodyLabel = McFullVtx3BodyLabels::iterator; } // namespace o2::aod #endif // PWGLF_DATAMODEL_VTX3BODYTABLES_H_ diff --git a/PWGLF/TableProducer/hypertriton3bodybuilder.cxx b/PWGLF/TableProducer/hypertriton3bodybuilder.cxx index 941b657e2f7..04c6370c477 100644 --- a/PWGLF/TableProducer/hypertriton3bodybuilder.cxx +++ b/PWGLF/TableProducer/hypertriton3bodybuilder.cxx @@ -49,7 +49,7 @@ using MyTracksIU = FullTracksExtIU; // in case requested using LabeledTracks = soa::Join; -struct hypertriton3bodybuilder { +struct hypertriton3bodyBuilder { Produces vtx3bodydata; Service ccdb; @@ -188,7 +188,7 @@ struct hypertriton3bodybuilder { initCCDB(bc); registry.fill(HIST("hEventCounter"), 0.5); - for (auto& vtx3body : decay3bodys) { // FIXME: turn into combination(...) + for (auto& vtx3body : decay3bodys) { registry.fill(HIST("hVtx3BodyCounter"), 0.5); @@ -264,7 +264,7 @@ struct hypertriton3bodybuilder { registry.fill(HIST("hVtx3BodyCounter"), 6.5); vtx3bodydata( - t0.globalIndex(), t1.globalIndex(), t2.globalIndex(), collision.globalIndex(), + t0.globalIndex(), t1.globalIndex(), t2.globalIndex(), collision.globalIndex(), vtx3body.globalIndex(), pos[0], pos[1], pos[2], p0[0], p0[1], p0[2], p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], fitter3body.getChi2AtPCACandidate(), @@ -273,9 +273,30 @@ struct hypertriton3bodybuilder { } }; +struct hypertriton3bodyDataLinkBuilder { + Produces vtxdataLink; + + void init(InitContext const&) {} + + void process(aod::Decay3Bodys const& decay3bodytable, aod::Vtx3BodyDatas const& vtxdatatable) + { + std::vector lIndices; + lIndices.reserve(decay3bodytable.size()); + for (int ii = 0; ii < decay3bodytable.size(); ii++) + lIndices[ii] = -1; + for (auto& vtxdata : vtxdatatable) { + lIndices[vtxdata.decay3bodyId()] = vtxdata.globalIndex(); + } + for (int ii = 0; ii < decay3bodytable.size(); ii++) { + vtxdataLink(lIndices[ii]); + } + } +}; + struct hypertriton3bodyLabelBuilder { - Produces vtxlabels; + Produces vtxlabels; + Produces vtxfulllabels; // for bookkeeping purposes: how many V0s come from same mother etc HistogramRegistry registry{ @@ -314,9 +335,15 @@ struct hypertriton3bodyLabelBuilder { } PROCESS_SWITCH(hypertriton3bodyLabelBuilder, processDoNotBuildLabels, "Do not produce MC label tables", true); - void processBuildLabels(aod::Collisions::iterator const& collision, aod::Decay3Bodys const& decay3bodys, LabeledTracks const&, aod::McParticles const& particlesMC) + void processBuildLabels(aod::Decay3BodysLinked const& decay3bodys, aod::Vtx3BodyDatas const& vtx3bodydatas, LabeledTracks const&, aod::McParticles const& particlesMC) { - for (auto& vtx : decay3bodys) { + std::vector lIndices; + lIndices.reserve(vtx3bodydatas.size()); + for (int ii = 0; ii < vtx3bodydatas.size(); ii++) { + lIndices[ii] = -1; + } + + for (auto& decay3body : decay3bodys) { int lLabel = -1; int lPDG = -1; @@ -324,22 +351,22 @@ struct hypertriton3bodyLabelBuilder { double MClifetime = -1; bool is3bodyDecay = false; // all considered V0s - auto lTrack0 = vtx.track0_as(); - auto lTrack1 = vtx.track1_as(); - auto lTrack2 = vtx.track2_as(); + auto lTrack0 = decay3body.track0_as(); + auto lTrack1 = decay3body.track1_as(); + auto lTrack2 = decay3body.track2_as(); registry.fill(HIST("hLabelCounter"), 0.5); // Association check // There might be smarter ways of doing this in the future if (!lTrack0.has_mcParticle() || !lTrack1.has_mcParticle() || !lTrack2.has_mcParticle()) { - vtxlabels(-1); + vtxfulllabels(-1); continue; } auto lMCTrack0 = lTrack0.mcParticle_as(); auto lMCTrack1 = lTrack1.mcParticle_as(); auto lMCTrack2 = lTrack2.mcParticle_as(); if (!lMCTrack0.has_mothers() || !lMCTrack1.has_mothers() || !lMCTrack2.has_mothers()) { - vtxlabels(-1); + vtxfulllabels(-1); continue; } @@ -350,14 +377,14 @@ struct hypertriton3bodyLabelBuilder { lLabel = lMother1.globalIndex(); lPt = lMother1.pt(); lPDG = lMother1.pdgCode(); - MClifetime = RecoDecay::sqrtSumOfSquares(lMCTrack2.vx() - lMother2.vx(), lMCTrack2.vy() - lMother2.vy(), lMCTrack2.vz() - lMother2.vz()) * 2.991 / lMother2.p(); + MClifetime = RecoDecay::sqrtSumOfSquares(lMCTrack2.vx() - lMother2.vx(), lMCTrack2.vy() - lMother2.vy(), lMCTrack2.vz() - lMother2.vz()) * o2::constants::physics::MassHyperTriton / lMother2.p(); is3bodyDecay = true; // vtxs with the same mother } } } } // end association check if (!is3bodyDecay) { - vtxlabels(-1); + vtxfulllabels(-1); continue; } registry.fill(HIST("hLabelCounter"), 1.5); @@ -365,7 +392,7 @@ struct hypertriton3bodyLabelBuilder { // Intended for cross-checks only // N.B. no rapidity cut! if (lPDG == 1010010030 && lMCTrack0.pdgCode() == 2212 && lMCTrack1.pdgCode() == -211 && lMCTrack2.pdgCode() == 1000010020) { - double hypertritonMCMass = RecoDecay::m(array{array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, array{RecoDecay::getMassPDG(kProton), RecoDecay::getMassPDG(kPiPlus), 1.87561}); + double hypertritonMCMass = RecoDecay::m(array{array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged, o2::constants::physics::MassDeuteron}); registry.fill(HIST("hHypertritonCounter"), 0.5); registry.fill(HIST("hHypertriton"), lPt); registry.fill(HIST("hHypertritonMass"), hypertritonMCMass); @@ -384,7 +411,7 @@ struct hypertriton3bodyLabelBuilder { } } if (lPDG == -1010010030 && lMCTrack0.pdgCode() == 211 && lMCTrack1.pdgCode() == -2212 && lMCTrack2.pdgCode() == -1000010020) { - double antiHypertritonMCMass = RecoDecay::m(array{array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, array{RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kProton), 1.87561}); + double antiHypertritonMCMass = RecoDecay::m(array{array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton, o2::constants::physics::MassDeuteron}); registry.fill(HIST("hHypertritonCounter"), 2.5); registry.fill(HIST("hAntiHypertriton"), lPt); registry.fill(HIST("hAntiHypertritonMass"), antiHypertritonMCMass); @@ -403,15 +430,20 @@ struct hypertriton3bodyLabelBuilder { } } - // Construct label table (note: this will be joinable with V0Datas) - vtxlabels( - lLabel); + // Construct label table + vtxfulllabels(lLabel); + if (decay3body.vtx3BodyDataId() != -1) { + lIndices[decay3body.vtx3BodyDataId()] = lLabel; + } + } + for (int ii = 0; ii < vtx3bodydatas.size(); ii++) { + vtxlabels(lIndices[ii]); } } PROCESS_SWITCH(hypertriton3bodyLabelBuilder, processBuildLabels, "Produce MC label tables", false); }; -struct hypertriton3bodyinitializer { +struct hypertriton3bodyInitializer { Spawns vtx3bodydatas; void init(InitContext const&) {} }; @@ -419,8 +451,9 @@ struct hypertriton3bodyinitializer { WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), }; } diff --git a/PWGLF/TableProducer/hypertriton3bodyfinder.cxx b/PWGLF/TableProducer/hypertriton3bodyfinder.cxx index c5b78d68056..622d7d50e0f 100644 --- a/PWGLF/TableProducer/hypertriton3bodyfinder.cxx +++ b/PWGLF/TableProducer/hypertriton3bodyfinder.cxx @@ -496,10 +496,10 @@ struct hypertriton3bodyFinder { float p2V0 = pt2V0 + pV0[2] * pV0[2], ptV0 = std::sqrt(pt2V0); // apply mass selections // float p2Pos = pP[0] * pP[0] + pP[1] * pP[1] + pP[2] * pP[2], p2Neg = pN[0] * pN[0] + pN[1] * pN[1] + pN[2] * pN[2]; - float massV0LambdaHyp = RecoDecay::m(array{array{pP[0], pP[1], pP[2]}, array{pN[0], pN[1], pN[2]}}, array{RecoDecay::getMassPDG(kProton), RecoDecay::getMassPDG(kPiPlus)}); - float massV0AntiLambdaHyp = RecoDecay::m(array{array{pP[0], pP[1], pP[2]}, array{pN[0], pN[1], pN[2]}}, array{RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kProton)}); + float massV0LambdaHyp = RecoDecay::m(array{array{pP[0], pP[1], pP[2]}, array{pN[0], pN[1], pN[2]}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}); + float massV0AntiLambdaHyp = RecoDecay::m(array{array{pP[0], pP[1], pP[2]}, array{pN[0], pN[1], pN[2]}}, array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton}); float massMargin = 20 * (0.001 * (1. + 0.5 * ptV0)) + 0.07; - if (massV0LambdaHyp - RecoDecay::getMassPDG(kLambda0) > massMargin && massV0AntiLambdaHyp - RecoDecay::getMassPDG(kLambda0) > massMargin) { + if (massV0LambdaHyp - o2::constants::physics::MassLambda > massMargin && massV0AntiLambdaHyp - o2::constants::physics::MassLambda > massMargin) { continue; } registry.fill(HIST("hV0Counter"), 7.5); @@ -642,7 +642,7 @@ struct hypertriton3bodyFinder { // Fix: Reconstruction H3L DCA Check vtx3bodydata( - t0.globalIndex(), t1.globalIndex(), t2.globalIndex(), collision.globalIndex(), + t0.globalIndex(), t1.globalIndex(), t2.globalIndex(), collision.globalIndex(), 0, vertexXYZ[0], vertexXYZ[1], vertexXYZ[2], p0[0], p0[1], p0[2], p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], fitter3body.getChi2AtPCACandidate(), diff --git a/PWGLF/Tasks/hypertriton3bodyMCQA.cxx b/PWGLF/Tasks/hypertriton3bodyMCQA.cxx index 22887116c97..05afc8b3cb5 100644 --- a/PWGLF/Tasks/hypertriton3bodyMCQA.cxx +++ b/PWGLF/Tasks/hypertriton3bodyMCQA.cxx @@ -71,10 +71,10 @@ struct hypertriton3bodyQaMc { registry.add("hMassHypertriton", "hMassHypertriton", {HistType::kTH1F, {massAxis}}); registry.add("hMassAntiHypertriton", "hMassAntiHypertriton", {HistType::kTH1F, {massAxis}}); } - void process(aod::Collision const& collision, aod::Vtx3BodyDatas const& vtx3BodyDatas, MyTracks const& tracks) + void process(aod::Collision const& collision, aod::Vtx3BodyDatas const& vtx3bodydatas, MyTracks const& tracks) { - for (auto& vtx : vtx3BodyDatas) { + for (auto& vtx : vtx3bodydatas) { registry.fill(HIST("hVtxRadius"), vtx.vtxradius()); registry.fill(HIST("hVtxCosPA"), vtx.vtxcosPA(collision.posX(), collision.posY(), collision.posZ())); registry.fill(HIST("hDCATrack0ToPV"), vtx.dcatrack0topv()); @@ -100,9 +100,9 @@ struct hypertriton3bodyAnalysisMc { {"hSelectedCandidatesCounter", "hSelectedCandidatesCounter", {HistType::kTH1F, {{9, 0.0f, 9.0f}}}}, {"hSelectedTrueHypertritonCounter", "hSelectedTrueHypertritonCounter", {HistType::kTH1F, {{9, 0.0f, 9.0f}}}}, {"hTestCounter", "hTestCounter", {HistType::kTH1F, {{9, 0.0f, 9.0f}}}}, - {"hMassHypertriton", "hMassHypertriton", {HistType::kTH1F, {{40, 2.95f, 3.05f}}}}, - {"hMassAntiHypertriton", "hMassAntiHypertriton", {HistType::kTH1F, {{40, 2.95f, 3.05f}}}}, - {"hMassHypertritonTotal", "hMassHypertritonTotal", {HistType::kTH1F, {{120, 2.9f, 3.2f}}}}, + {"hMassHypertriton", "hMassHypertriton", {HistType::kTH1F, {{100, 2.95f, 3.05f}}}}, + {"hMassAntiHypertriton", "hMassAntiHypertriton", {HistType::kTH1F, {{100, 2.95f, 3.05f}}}}, + {"hMassHypertritonTotal", "hMassHypertritonTotal", {HistType::kTH1F, {{300, 2.9f, 3.2f}}}}, {"hPtProton", "hPtProton", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, {"hPtAntiPion", "hPtAntiPion", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, {"hPtDeuteron", "hPtDeuteron", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, @@ -168,7 +168,7 @@ struct hypertriton3bodyAnalysisMc { Configurable minDeuteronPt{"minDeuteronPt", 0.6, "minDeuteronPt"}; Configurable maxDeuteronPt{"maxDeuteronPt", 10, "maxDeuteronPt"}; - void process(soa::Join::iterator const& collision, aod::Vtx3BodyDatas const& vtx3BodyDatas, MyTracks const& tracks, aod::McParticles const& particlesMC) + void process(soa::Join::iterator const& collision, aod::Vtx3BodyDatas const& vtx3bodydatas, MyTracks const& tracks, aod::McParticles const& particlesMC) { registry.fill(HIST("hSelectedEventCounter"), 0.5); if (eventSelection && !collision.sel8()) { @@ -176,7 +176,7 @@ struct hypertriton3bodyAnalysisMc { } registry.fill(HIST("hSelectedEventCounter"), 1.5); - for (auto& vtx : vtx3BodyDatas) { + for (auto& vtx : vtx3bodydatas) { // couldn't filter cosPA and radius variables (dynamic columns) int lLabel = -1; @@ -202,7 +202,7 @@ struct hypertriton3bodyAnalysisMc { if ((lPDG == 1010010030 && lMCTrack0.pdgCode() == 2212 && lMCTrack1.pdgCode() == -211 && lMCTrack2.pdgCode() == 1000010020) || (lPDG == -1010010030 && lMCTrack0.pdgCode() == 211 && lMCTrack1.pdgCode() == -2212 && lMCTrack2.pdgCode() == -1000010020)) { isFromHypertriton = true; - MClifetime = RecoDecay::sqrtSumOfSquares(lMCTrack2.vx() - lMother2.vx(), lMCTrack2.vy() - lMother2.vy(), lMCTrack2.vz() - lMother2.vz()) * 2.991 / lMother2.p(); + MClifetime = RecoDecay::sqrtSumOfSquares(lMCTrack2.vx() - lMother2.vx(), lMCTrack2.vy() - lMother2.vy(), lMCTrack2.vz() - lMother2.vz()) * o2::constants::physics::MassHyperTriton / lMother2.p(); } } } @@ -253,6 +253,11 @@ struct hypertriton3bodyAnalysisMc { registry.fill(HIST("hSelectedTrueHypertritonCounter"), 5.5); } + // 3sigma region for Dalitz plot + double mcsigma = 0.0018; // from local MC + double lowerlimit = o2::constants::physics::MassHyperTriton - 3 * mcsigma; + double upperlimit = o2::constants::physics::MassHyperTriton + 3 * mcsigma; + // Hypertriton if (TMath::Abs(vtx.track0_as().tpcNSigmaPr()) < TpcPidNsigmaCut && TMath::Abs(vtx.track1_as().tpcNSigmaPi()) < TpcPidNsigmaCut && TMath::Abs(vtx.track2_as().tpcNSigmaDe()) < TpcPidNsigmaCut) { registry.fill(HIST("hSelectedCandidatesCounter"), 6.5); @@ -290,6 +295,11 @@ struct hypertriton3bodyAnalysisMc { registry.fill(HIST("hMassHypertritonTotal"), vtx.mHypertriton()); registry.fill(HIST("h3dMassHypertriton"), 0., vtx.pt(), vtx.mHypertriton()); // collision.centV0M() instead of 0. once available registry.fill(HIST("h3dTotalHypertriton"), ct, vtx.pt(), vtx.mHypertriton()); + if (vtx.mHypertriton() > lowerlimit && vtx.mHypertriton() < upperlimit) { + registry.fill(HIST("hDalitz"), RecoDecay::m2(array{array{vtx.pxtrack0(), vtx.pytrack0(), vtx.pztrack0()}, array{vtx.pxtrack2(), vtx.pytrack2(), vtx.pztrack2()}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassDeuteron}), + RecoDecay::m2(array{array{vtx.pxtrack0(), vtx.pytrack0(), vtx.pztrack0()}, array{vtx.pxtrack1(), vtx.pytrack1(), vtx.pztrack1()}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged})); + } + if (isFromHypertriton) { registry.fill(HIST("h3dTotalTrueHypertriton"), MClifetime, lPt, vtx.mHypertriton()); } @@ -319,7 +329,7 @@ struct hypertriton3bodyAnalysisMc { } } - if (vtx.track0pt() > minPionPt && vtx.track0pt() < maxPionPt && vtx.track1pt() > minProtonPt && vtx.track1pt() < maxProtonPt && vtx.track2pt() > minDeuteronPt && vtx.track2pt() < maxDeuteronPt) { + if (vtx.track1pt() > minPionPt && vtx.track1pt() < maxPionPt && vtx.track0pt() > minProtonPt && vtx.track0pt() < maxProtonPt && vtx.track2pt() > minDeuteronPt && vtx.track2pt() < maxDeuteronPt) { registry.fill(HIST("hSelectedCandidatesCounter"), 7.5); if (isFromHypertriton) { registry.fill(HIST("hSelectedTrueHypertritonCounter"), 7.5); @@ -330,13 +340,18 @@ struct hypertriton3bodyAnalysisMc { registry.fill(HIST("hSelectedTrueHypertritonCounter"), 8.5); } - registry.fill(HIST("hPtAntiProton"), vtx.track0pt()); - registry.fill(HIST("hPtPion"), vtx.track1pt()); + registry.fill(HIST("hPtAntiProton"), vtx.track1pt()); + registry.fill(HIST("hPtPion"), vtx.track0pt()); registry.fill(HIST("hPtAntiDeuteron"), vtx.track2pt()); registry.fill(HIST("hMassAntiHypertriton"), vtx.mAntiHypertriton()); registry.fill(HIST("hMassHypertritonTotal"), vtx.mAntiHypertriton()); registry.fill(HIST("h3dMassAntiHypertriton"), 0., vtx.pt(), vtx.mAntiHypertriton()); // collision.centV0M() instead of 0. once available registry.fill(HIST("h3dTotalHypertriton"), ct, vtx.pt(), vtx.mAntiHypertriton()); + if (vtx.mAntiHypertriton() > lowerlimit && vtx.mAntiHypertriton() < upperlimit) { + registry.fill(HIST("hDalitz"), RecoDecay::m2(array{array{vtx.pxtrack1(), vtx.pytrack1(), vtx.pztrack1()}, array{vtx.pxtrack2(), vtx.pytrack2(), vtx.pztrack2()}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassDeuteron}), + RecoDecay::m2(array{array{vtx.pxtrack1(), vtx.pytrack1(), vtx.pztrack1()}, array{vtx.pxtrack0(), vtx.pytrack0(), vtx.pztrack0()}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged})); + } + if (isFromHypertriton) { registry.fill(HIST("h3dTotalTrueHypertriton"), MClifetime, lPt, vtx.mAntiHypertriton()); } @@ -351,13 +366,55 @@ struct hypertriton3bodyAnalysisMc { // PROCESS_SWITCH(hypertriton3bodyAnalysis, processRun3, "Process Run 3 data", true); }; +// check vtx3body with mclabels +struct hypertriton3bodyLabelCheck { + HistogramRegistry registry{ + "registry", + { + {"hLabeledVtxCounter", "hLabeledVtxCounter", {HistType::kTH1F, {{2, 0.0f, 2.0f}}}}, + {"hMassTruesH3L", "hMassTruesH3L", {HistType::kTH1F, {{100, 2.95f, 3.05f}}}}, + {"hMassTruesH3LMatter", "hMassTruesH3LMatter", {HistType::kTH1F, {{100, 2.95f, 3.05f}}}}, + {"hMassTruesH3LAntiMatter", "hMassTruesH3LAntiMatter", {HistType::kTH1F, {{100, 2.95f, 3.05f}}}}, + }, + }; + + void init(InitContext const&) + { + registry.get(HIST("hLabeledVtxCounter"))->GetXaxis()->SetBinLabel(1, "Readin"); + registry.get(HIST("hLabeledVtxCounter"))->GetXaxis()->SetBinLabel(2, "TrueMCH3L"); + } + + Configurable eventSelection{"eventSelection", true, "event selection"}; + + void process(soa::Join::iterator const& collision, soa::Join const& vtx3bodydatas, MyTracks const& tracks, aod::McParticles const& particlesMC) + { + if (eventSelection && !collision.sel8()) { + return; + } + + for (auto& vtx : vtx3bodydatas) { + registry.fill(HIST("hLabeledVtxCounter"), 0.5); + if (vtx.mcParticleId() != -1) { + registry.fill(HIST("hLabeledVtxCounter"), 1.5); + auto mcparticle = vtx.mcParticle_as(); + if (mcparticle.pdgCode() > 0) { + registry.fill(HIST("hMassTruesH3L"), vtx.mHypertriton()); + registry.fill(HIST("hMassTruesH3LMatter"), vtx.mHypertriton()); + } else { + registry.fill(HIST("hMassTruesH3L"), vtx.mAntiHypertriton()); + registry.fill(HIST("hMassTruesH3LAntiMatter"), vtx.mAntiHypertriton()); + } + } + } + } +}; + // check the properties of daughters candidates and true daughters struct hypertriton3bodyTrackMcinfo { // Basic checks HistogramRegistry registry{ "registry", { - {"hTotalCollCounter", "hTotalCollCounter", {HistType::kTH1F, {{1, 0.0f, 1.0f}}}}, {"hParticleCount", "hParticleCount", {HistType::kTH1F, {{7, 0.0f, 7.0f}}}}, {"hParticleCount2", "hParticleCount2", {HistType::kTH1F, {{7, 0.0f, 7.0f}}}}, // for tpcncls > 70 @@ -820,11 +877,12 @@ struct hypertriton3bodyMcParticleCount { {"hMcHypertritonDecayCount", "hMcHypertritonDecayCount", {HistType::kTH1F, {{3, 0.0f, 3.0f}}}}, {"hMcHypertritonPt", "hMcHypertritonPt", {HistType::kTH1F, {{300, 0.0f, 15.0f}}}}, {"hMcHypertritonLifeTime", "hMcHypertritonLifeTime", {HistType::kTH1F, {{500, 0.0f, 500.0f}}}}, + {"hMcHypertritonDaughterPionCount", "hMcHypertritonDaughterPionCount", {HistType::kTH1F, {{10, 0.0f, 10.0f}}}}, {"hMcProtonPt", "hMcProtonPt", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, {"hMcPionPt", "hMcPionPt", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, {"hMcDeuteronPt", "hMcDeuteronPt", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, - {"hMcRecoInvMass", "hMcRecoInvMass", {HistType::kTH1F, {{320, 2.8, 5.0f}}}}, + {"hMcRecoInvMass", "hMcRecoInvMass", {HistType::kTH1F, {{100, 2.95, 3.05f}}}}, }, }; @@ -953,18 +1011,17 @@ struct hypertriton3bodyMcParticleCount { registry.fill(HIST("hMcHypertritonDecayCount"), 0.5); registry.fill(HIST("hMcHypertritonDecayCount"), 2.5); } - double MClifetime = RecoDecay::sqrtSumOfSquares(dauDeuteronPos[0] - mcparticle.vx(), dauDeuteronPos[1] - mcparticle.vy(), dauDeuteronPos[2] - mcparticle.vz()) * 2.991 / mcparticle.p(); + double MClifetime = RecoDecay::sqrtSumOfSquares(dauDeuteronPos[0] - mcparticle.vx(), dauDeuteronPos[1] - mcparticle.vy(), dauDeuteronPos[2] - mcparticle.vz()) * o2::constants::physics::MassHyperTriton / mcparticle.p(); if ((haveProton && haveAntiPion && haveDeuteron && mcparticle.pdgCode() == 1010010030) || (haveAntiProton && havePion && haveAntiDeuteron && mcparticle.pdgCode() == -1010010030)) { - registry.fill(HIST("hMcRecoInvMass"), RecoDecay::m(array{array{dauProtonMom[0], dauProtonMom[1], dauProtonMom[2]}, array{dauPionMom[0], dauPionMom[1], dauPionMom[2]}, array{dauDeuteronMom[0], dauDeuteronMom[1], dauDeuteronMom[2]}}, array{RecoDecay::getMassPDG(kProton), RecoDecay::getMassPDG(kPiPlus), 1.87561})); + registry.fill(HIST("hMcRecoInvMass"), RecoDecay::m(array{array{dauProtonMom[0], dauProtonMom[1], dauProtonMom[2]}, array{dauPionMom[0], dauPionMom[1], dauPionMom[2]}, array{dauDeuteronMom[0], dauDeuteronMom[1], dauDeuteronMom[2]}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged, o2::constants::physics::MassDeuteron})); registry.fill(HIST("h3dMCDecayedHypertriton"), mcparticle.y(), mcparticle.pt(), MClifetime); } int daughterPionCount = 0; for (auto& mcparticleDaughter : mcparticle.daughters_as()) { - if (std::abs(mcparticleDaughter.pdgCode()) != 211) { - break; + if (std::abs(mcparticleDaughter.pdgCode()) == 211) { + daughterPionCount++; } - daughterPionCount++; } registry.fill(HIST("hMcHypertritonLifeTime"), MClifetime); @@ -1014,6 +1071,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), diff --git a/PWGLF/Tasks/hypertriton3bodyanalysis.cxx b/PWGLF/Tasks/hypertriton3bodyanalysis.cxx index 252f3dd4891..ca2b304b33d 100644 --- a/PWGLF/Tasks/hypertriton3bodyanalysis.cxx +++ b/PWGLF/Tasks/hypertriton3bodyanalysis.cxx @@ -70,9 +70,9 @@ struct hypertriton3bodyQa { registry.add("hMassHypertriton", "hMassHypertriton", {HistType::kTH1F, {massAxis}}); registry.add("hMassAntiHypertriton", "hMassAntiHypertriton", {HistType::kTH1F, {massAxis}}); } - void process(aod::Collision const& collision, aod::Vtx3BodyDatas const& vtx3BodyDatas, MyTracks const& tracks) + void process(aod::Collision const& collision, aod::Vtx3BodyDatas const& vtx3bodydatas, MyTracks const& tracks) { - for (auto& vtx : vtx3BodyDatas) { + for (auto& vtx : vtx3bodydatas) { registry.fill(HIST("hVtxRadius"), vtx.vtxradius()); registry.fill(HIST("hVtxCosPA"), vtx.vtxcosPA(collision.posX(), collision.posY(), collision.posZ())); registry.fill(HIST("hDCATrack0ToPV"), vtx.dcatrack0topv()); @@ -96,21 +96,28 @@ struct hypertriton3bodyAnalysis { { {"hSelectedEventCounter", "hSelectedEventCounter", {HistType::kTH1F, {{2, 0.0f, 2.0f}}}}, {"hSelectedCandidatesCounter", "hSelectedCandidatesCounter", {HistType::kTH1F, {{9, 0.0f, 9.0f}}}}, - {"hMassHypertriton", "hMassHypertriton", {HistType::kTH1F, {{40, 2.95f, 3.05f}}}}, - {"hMassAntiHypertriton", "hMassAntiHypertriton", {HistType::kTH1F, {{40, 2.95f, 3.05f}}}}, - {"hMassHypertritonTotal", "hMassHypertritonTotal", {HistType::kTH1F, {{120, 2.9f, 3.2f}}}}, + {"hMassHypertriton", "hMassHypertriton", {HistType::kTH1F, {{100, 2.95f, 3.05f}}}}, + {"hMassAntiHypertriton", "hMassAntiHypertriton", {HistType::kTH1F, {{100, 2.95f, 3.05f}}}}, + {"hMassHypertritonTotal", "hMassHypertritonTotal", {HistType::kTH1F, {{300, 2.9f, 3.2f}}}}, {"hPtProton", "hPtProton", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, {"hPtAntiPion", "hPtAntiPion", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, {"hPtDeuteron", "hPtDeuteron", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, {"hPtAntiProton", "hPtAntiProton", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, {"hPtPion", "hPtPion", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, {"hPtAntiDeuteron", "hPtAntiDeuteron", {HistType::kTH1F, {{200, 0.0f, 10.0f}}}}, - {"h3dMassHypertriton", "h3dMassHypertriton", {HistType::kTH3F, {{20, 0.0f, 100.0f, "Cent (%)"}, {200, 0.0f, 10.0f, "#it{p}_{T} (GeV/c)"}, {40, 2.95f, 3.05f, "Inv. Mass (GeV/c^{2})"}}}}, + {"h3dMassHypertriton", "h3dMassHypertriton", {HistType::kTH3F, {{20, 0.0f, 100.0f, "Cent (%)"}, {200, 0.0f, 10.0f, "#it{p}_{T} (GeV/c)"}, {100, 2.95f, 3.05f, "Inv. Mass (GeV/c^{2})"}}}}, {"hTPCPIDProton", "hTPCPIDProton", {HistType::kTH1F, {{240, -6.0f, 6.0f}}}}, {"hTPCPIDPion", "hTPCPIDPion", {HistType::kTH1F, {{240, -6.0f, 6.0f}}}}, {"hTPCPIDDeuteron", "hTPCPIDDeuteron", {HistType::kTH1F, {{240, -6.0f, 6.0f}}}}, - {"h3dMassAntiHypertriton", "h3dMassAntiHypertriton", {HistType::kTH3F, {{20, 0.0f, 100.0f, "Cent (%)"}, {200, 0.0f, 10.0f, "#it{p}_{T} (GeV/c)"}, {40, 2.95f, 3.05f, "Inv. Mass (GeV/c^{2})"}}}}, - {"h3dTotalHypertriton", "h3dTotalHypertriton", {HistType::kTH3F, {{50, 0, 50, "ct(cm)"}, {200, 0.0f, 10.0f, "#it{p}_{T} (GeV/c)"}, {40, 2.95f, 3.05f, "Inv. Mass (GeV/c^{2})"}}}}, + {"hProtonTPCBB", "hProtonTPCBB", {HistType::kTH2F, {{160, -8.0f, 8.0f, "p/z(GeV/c)"}, {200, 0.0f, 1000.0f, "TPCSignal"}}}}, + {"hPionTPCBB", "hPionTPCBB", {HistType::kTH2F, {{160, -8.0f, 8.0f, "p/z(GeV/c)"}, {200, 0.0f, 1000.0f, "TPCSignal"}}}}, + {"hDeuteronTPCBB", "hDeuteronTPCBB", {HistType::kTH2F, {{160, -8.0f, 8.0f, "p/z(GeV/c)"}, {200, 0.0f, 1000.0f, "TPCSignal"}}}}, + {"hProtonTPCVsPt", "hProtonTPCVsPt", {HistType::kTH2F, {{50, 0.0f, 5.0f, "#it{p}_{T} (GeV/c)"}, {240, -6.0f, 6.0f, "TPC n#sigma"}}}}, + {"hPionTPCVsPt", "hPionTPCVsPt", {HistType::kTH2F, {{20, 0.0f, 2.0f, "#it{p}_{T} (GeV/c)"}, {240, -6.0f, 6.0f, "TPC n#sigma"}}}}, + {"hDeuteronTPCVsPt", "hDeuteronTPCVsPt", {HistType::kTH2F, {{80, 0.0f, 8.0f, "#it{p}_{T} (GeV/c)"}, {240, -6.0f, 6.0f, "TPC n#sigma"}}}}, + {"h3dMassAntiHypertriton", "h3dMassAntiHypertriton", {HistType::kTH3F, {{20, 0.0f, 100.0f, "Cent (%)"}, {200, 0.0f, 10.0f, "#it{p}_{T} (GeV/c)"}, {100, 2.95f, 3.05f, "Inv. Mass (GeV/c^{2})"}}}}, + {"h3dTotalHypertriton", "h3dTotalHypertriton", {HistType::kTH3F, {{50, 0, 50, "ct(cm)"}, {200, 0.0f, 10.0f, "#it{p}_{T} (GeV/c)"}, {100, 2.95f, 3.05f, "Inv. Mass (GeV/c^{2})"}}}}, + {"hDalitz", "hDalitz", {HistType::kTH2F, {{120, 7.85, 8.45, "M^{2}(dp) (GeV^{2}/c^{4})"}, {60, 1.1, 1.4, "M^{2}(p#pi) (GeV^{2}/c^{4})"}}}}, }, }; @@ -159,7 +166,7 @@ struct hypertriton3bodyAnalysis { // Filter dcaFilterV0 = aod::vtx.ata::dcaV0daughters < dcavtx.au; - void process(soa::Join::iterator const& collision, aod::Vtx3BodyDatas const& vtx3BodyDatas, MyTracks const& tracks) + void process(soa::Join::iterator const& collision, aod::Vtx3BodyDatas const& vtx3bodydatas, MyTracks const& tracks) { registry.fill(HIST("hSelectedEventCounter"), 0.5); if (eventSelection && !collision.sel8()) { @@ -167,8 +174,7 @@ struct hypertriton3bodyAnalysis { } registry.fill(HIST("hSelectedEventCounter"), 1.5); - for (auto& vtx : vtx3BodyDatas) { - // couldn't filter cosPA and radius variables (dynamic columns) + for (auto& vtx : vtx3bodydatas) { registry.fill(HIST("hSelectedCandidatesCounter"), 0.5); if (vtx.vtxcosPA(collision.posX(), collision.posY(), collision.posZ()) < vtxcospa) { continue; @@ -192,6 +198,11 @@ struct hypertriton3bodyAnalysis { } registry.fill(HIST("hSelectedCandidatesCounter"), 5.5); + // 3sigma region for Dalitz plot + double mcsigma = 0.0018; // from local MC + double lowerlimit = o2::constants::physics::MassHyperTriton - 3 * mcsigma; + double upperlimit = o2::constants::physics::MassHyperTriton + 3 * mcsigma; + // Hypertriton if (TMath::Abs(vtx.track0_as().tpcNSigmaPr()) < TpcPidNsigmaCut && TMath::Abs(vtx.track1_as().tpcNSigmaPi()) < TpcPidNsigmaCut && TMath::Abs(vtx.track2_as().tpcNSigmaDe()) < TpcPidNsigmaCut) { @@ -209,10 +220,21 @@ struct hypertriton3bodyAnalysis { registry.fill(HIST("hTPCPIDProton"), vtx.track0_as().tpcNSigmaPr()); registry.fill(HIST("hTPCPIDPion"), vtx.track1_as().tpcNSigmaPi()); registry.fill(HIST("hTPCPIDDeuteron"), vtx.track2_as().tpcNSigmaDe()); + registry.fill(HIST("hProtonTPCBB"), vtx.track0_as().p(), vtx.track0_as().tpcSignal()); + registry.fill(HIST("hPionTPCBB"), -vtx.track1_as().p(), vtx.track1_as().tpcSignal()); + registry.fill(HIST("hDeuteronTPCBB"), vtx.track2_as().p(), vtx.track0_as().tpcSignal()); + registry.fill(HIST("hProtonTPCVsPt"), vtx.track0pt(), vtx.track0_as().tpcNSigmaPr()); + registry.fill(HIST("hPionTPCVsPt"), vtx.track1pt(), vtx.track1_as().tpcNSigmaPi()); + registry.fill(HIST("hDeuteronTPCVsPt"), vtx.track2pt(), vtx.track2_as().tpcNSigmaDe()); registry.fill(HIST("hMassHypertriton"), vtx.mHypertriton()); registry.fill(HIST("hMassHypertritonTotal"), vtx.mHypertriton()); registry.fill(HIST("h3dMassHypertriton"), 0., vtx.pt(), vtx.mHypertriton()); // collision.centV0M() instead of 0. once available registry.fill(HIST("h3dTotalHypertriton"), ct, vtx.pt(), vtx.mHypertriton()); + if (vtx.mHypertriton() > lowerlimit && vtx.mHypertriton() < upperlimit) { + registry.fill(HIST("hDalitz"), RecoDecay::m2(array{array{vtx.pxtrack0(), vtx.pytrack0(), vtx.pztrack0()}, array{vtx.pxtrack2(), vtx.pytrack2(), vtx.pztrack2()}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassDeuteron}), + RecoDecay::m2(array{array{vtx.pxtrack0(), vtx.pytrack0(), vtx.pztrack0()}, array{vtx.pxtrack1(), vtx.pytrack1(), vtx.pztrack1()}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged})); + } + if (saveDcaHist == 1) { registry.fill(HIST("h3dMassHypertritonDca"), vtx.dcaVtxdaughters(), vtx.pt(), vtx.mHypertriton()); } @@ -236,10 +258,21 @@ struct hypertriton3bodyAnalysis { registry.fill(HIST("hTPCPIDProton"), vtx.track1_as().tpcNSigmaPr()); registry.fill(HIST("hTPCPIDPion"), vtx.track0_as().tpcNSigmaPi()); registry.fill(HIST("hTPCPIDDeuteron"), vtx.track2_as().tpcNSigmaDe()); + registry.fill(HIST("hProtonTPCBB"), -vtx.track1_as().p(), vtx.track1_as().tpcSignal()); + registry.fill(HIST("hPionTPCBB"), vtx.track0_as().p(), vtx.track0_as().tpcSignal()); + registry.fill(HIST("hDeuteronTPCBB"), -vtx.track2_as().p(), vtx.track0_as().tpcSignal()); + registry.fill(HIST("hProtonTPCVsPt"), vtx.track1pt(), vtx.track1_as().tpcNSigmaPr()); + registry.fill(HIST("hPionTPCVsPt"), vtx.track0pt(), vtx.track0_as().tpcNSigmaPi()); + registry.fill(HIST("hDeuteronTPCVsPt"), vtx.track2pt(), vtx.track2_as().tpcNSigmaDe()); registry.fill(HIST("hMassAntiHypertriton"), vtx.mAntiHypertriton()); registry.fill(HIST("hMassHypertritonTotal"), vtx.mAntiHypertriton()); registry.fill(HIST("h3dMassAntiHypertriton"), 0., vtx.pt(), vtx.mAntiHypertriton()); // collision.centV0M() instead of 0. once available registry.fill(HIST("h3dTotalHypertriton"), ct, vtx.pt(), vtx.mAntiHypertriton()); + if (vtx.mAntiHypertriton() > lowerlimit && vtx.mAntiHypertriton() < upperlimit) { + registry.fill(HIST("hDalitz"), RecoDecay::m2(array{array{vtx.pxtrack1(), vtx.pytrack1(), vtx.pztrack1()}, array{vtx.pxtrack2(), vtx.pytrack2(), vtx.pztrack2()}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassDeuteron}), + RecoDecay::m2(array{array{vtx.pxtrack1(), vtx.pytrack1(), vtx.pztrack1()}, array{vtx.pxtrack0(), vtx.pytrack0(), vtx.pztrack0()}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged})); + } + if (saveDcaHist == 1) { registry.fill(HIST("h3dMassAntiHypertritonDca"), vtx.dcaVtxdaughters(), vtx.pt(), vtx.mAntiHypertriton()); } @@ -248,11 +281,8 @@ struct hypertriton3bodyAnalysis { } } } - - // PROCESS_SWITCH(hypertriton3bodyAnalysis, process, "Process Run 3 data", true); }; -// check the properties of daughters candidates and true daughters WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{