From a3737b7f9fcd878ccab576ac143b660cc955b0cc Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Sun, 5 May 2024 12:14:12 +0200 Subject: [PATCH] PWGLF: add found tag in findable exercise (#5991) * PWGLF: add found tag in findable exercise * Fix default * Many bugfixes * Please consider the following formatting changes (#288) * Changes in QA counting * Please consider the following formatting changes (#289) * Megalinter fix * Please consider the following formatting changes (#290) * Asymmetric V0MCCore build mode for duplicate handling * Extra bookkeeping * Please consider the following formatting changes (#291) * Fix bugs * Please consider the following formatting changes (#292) * V0Cores vs V0MCCores: always interlink * Please consider the following formatting changes (#293) --------- Co-authored-by: ALICE Builder --- PWGLF/DataModel/LFStrangenessTables.h | 39 +++ .../Strangeness/cascadebuilder.cxx | 224 ++++++++++++------ .../Strangeness/cascademcbuilder.cxx | 96 ++++---- .../Strangeness/cascademcfinder.cxx | 19 +- .../Strangeness/lambdakzerobuilder.cxx | 66 +++++- .../Strangeness/lambdakzeromcbuilder.cxx | 156 +++++++++--- .../Strangeness/lambdakzeromcfinder.cxx | 7 +- .../Strangeness/strangederivedbuilder.cxx | 85 +++++++ 8 files changed, 531 insertions(+), 161 deletions(-) diff --git a/PWGLF/DataModel/LFStrangenessTables.h b/PWGLF/DataModel/LFStrangenessTables.h index afd7faa3182..d7e0628d274 100644 --- a/PWGLF/DataModel/LFStrangenessTables.h +++ b/PWGLF/DataModel/LFStrangenessTables.h @@ -693,16 +693,40 @@ namespace v0data DECLARE_SOA_INDEX_COLUMN(V0Data, v0Data); //! Index to V0Data entry DECLARE_SOA_INDEX_COLUMN(V0fCData, v0fCData); //! Index to V0Data entry DECLARE_SOA_INDEX_COLUMN_FULL(V0MC, v0MC, int, V0MCCores, "_MC"); //! +DECLARE_SOA_INDEX_COLUMN(V0MCCore, v0MCCore); } // namespace v0data DECLARE_SOA_TABLE(V0DataLink, "AOD", "V0DATALINK", //! Joinable table with V0s which links to V0Data which is not produced for all entries o2::soa::Index<>, v0data::V0DataId, v0data::V0fCDataId); DECLARE_SOA_TABLE(V0MCRefs, "AOD", "V0MCREF", //! index table when using AO2Ds o2::soa::Index<>, v0data::V0MCId); +DECLARE_SOA_TABLE(V0CoreMCLabels, "AOD", "V0COREMCLABEL", //! optional table to refer to V0MCCores if not joinable + o2::soa::Index<>, v0data::V0MCCoreId); using V0sLinked = soa::Join; using V0Linked = V0sLinked::iterator; +namespace v0data +{ +DECLARE_SOA_COLUMN(IsFound, isFound, bool); //! is this FindableV0 actually in the V0s table? +} + +// Major bypass for simultaneous found vs findable study +DECLARE_SOA_TABLE(FindableV0s, "AOD", "FindableV0", //! Will store findable + o2::soa::Index<>, v0::CollisionId, + v0::PosTrackId, v0::NegTrackId, + v0::V0Type, + v0::IsStandardV0, + v0::IsPhotonV0, + v0::IsCollinearV0, + o2::soa::Marker<1>); + +DECLARE_SOA_TABLE(V0FoundTags, "AOD", "V0FoundTag", //! found or not? + v0data::IsFound); + +using FindableV0sLinked = soa::Join; +using FindableV0Linked = FindableV0sLinked::iterator; + // helper for building namespace v0tag { @@ -1239,6 +1263,21 @@ using CascadeLinked = CascadesLinked::iterator; using KFCascadesLinked = soa::Join; using KFCascadeLinked = KFCascadesLinked::iterator; +namespace cascdata +{ +DECLARE_SOA_INDEX_COLUMN(FindableV0, findableV0); //! V0 index +DECLARE_SOA_COLUMN(IsFound, isFound, bool); //! is this FindableCascade actually in the Cascades table? +} // namespace cascdata + +DECLARE_SOA_TABLE(FindableCascades, "AOD", "FINDABLECASCS", //! Run 3 cascade table + o2::soa::Index<>, cascade::CollisionId, cascdata::FindableV0Id, cascade::BachelorId, o2::soa::Marker<1>); + +DECLARE_SOA_TABLE(CascFoundTags, "AOD", "CascFoundTag", //! found or not? + cascdata::IsFound); + +using FindableCascadesLinked = soa::Join; +using FindableCascadeLinked = FindableCascadesLinked::iterator; + namespace casctag { DECLARE_SOA_COLUMN(IsInteresting, isInteresting, bool); //! will this be built or not? diff --git a/PWGLF/TableProducer/Strangeness/cascadebuilder.cxx b/PWGLF/TableProducer/Strangeness/cascadebuilder.cxx index f85d3e2d34f..b981070d9bb 100644 --- a/PWGLF/TableProducer/Strangeness/cascadebuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/cascadebuilder.cxx @@ -99,6 +99,7 @@ using TracksExtraWithPIDandLabels = soa::Join; using V0fCfull = soa::Join; using TaggedCascades = soa::Join; +using TaggedFindableCascades = soa::Join; // For MC association in pre-selection using LabeledTracksExtra = soa::Join; @@ -251,6 +252,7 @@ struct cascadeBuilder { // For manual sliceBy Preslice perCollision = o2::aod::cascade::collisionId; + Preslice perCollisionFindable = o2::aod::cascade::collisionId; Preslice perCascade = o2::aod::strangenesstracking::cascadeId; // Define o2 fitter, 2-prong, active memory (no need to redefine per event) @@ -533,8 +535,8 @@ struct cascadeBuilder { lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(ccdbConfigurations.lutPath)); } - if (doprocessRun2 == false && doprocessRun3 == false && doprocessRun3withStrangenessTracking == false && doprocessRun3withKFParticle == false) { - LOGF(fatal, "Neither processRun2 nor processRun3 nor processRun3withstrangenesstracking enabled. Please choose one!"); + if (doprocessRun2 == false && doprocessRun3 == false && doprocessRun3withStrangenessTracking == false && doprocessRun3withKFParticle == false && doprocessFindableRun3 == false) { + LOGF(fatal, "Neither processRun2 nor processRun3 nor processRun3withstrangenesstracking nor processFindableRun3 enabled. Please choose one!"); } if (doprocessRun2 == true && doprocessRun3 == true) { LOGF(fatal, "Cannot enable processRun2 and processRun3 at the same time. Please choose one."); @@ -1522,79 +1524,101 @@ struct cascadeBuilder { return true; } + template + void processCascadeCandidate(TV0Index const& v0index, TCascade const& cascade) + { + bool validCascadeCandidate = false; + if (v0index.has_v0Data()) { + // this V0 passed both standard V0 and cascade V0 selections + auto v0row = v0index.template v0Data_as(); + validCascadeCandidate = buildCascadeCandidate(cascade, v0row); + } else if (v0index.has_v0fCData()) { + // this V0 passes only V0-for-cascade selections, use that instead + auto v0row = v0index.template v0fCData_as(); + validCascadeCandidate = buildCascadeCandidate(cascade, v0row); + } else { + return; // this was inadequately linked, should not happen + } + if (!validCascadeCandidate) + return; // doesn't pass cascade selections + + // round the DCA variables to a certain precision if asked + if (roundDCAVariables) + roundCascadeCandidateVariables(); + + cascidx(/*cascadecandidate.v0Id, */ cascade.globalIndex(), + cascadecandidate.positiveId, cascadecandidate.negativeId, + cascadecandidate.bachelorId, cascade.collisionId()); + cascdata(cascadecandidate.charge, cascadecandidate.mXi, cascadecandidate.mOmega, + cascadecandidate.pos[0], cascadecandidate.pos[1], cascadecandidate.pos[2], + cascadecandidate.v0pos[0], cascadecandidate.v0pos[1], cascadecandidate.v0pos[2], + cascadecandidate.v0mompos[0], cascadecandidate.v0mompos[1], cascadecandidate.v0mompos[2], + cascadecandidate.v0momneg[0], cascadecandidate.v0momneg[1], cascadecandidate.v0momneg[2], + cascadecandidate.bachP[0], cascadecandidate.bachP[1], cascadecandidate.bachP[2], + cascadecandidate.bachP[0] + cascadecandidate.v0mompos[0] + cascadecandidate.v0momneg[0], // <--- redundant but ok + cascadecandidate.bachP[1] + cascadecandidate.v0mompos[1] + cascadecandidate.v0momneg[1], // <--- redundant but ok + cascadecandidate.bachP[2] + cascadecandidate.v0mompos[2] + cascadecandidate.v0momneg[2], // <--- redundant but ok + cascadecandidate.v0dcadau, cascadecandidate.dcacascdau, + cascadecandidate.v0dcapostopv, cascadecandidate.v0dcanegtopv, + cascadecandidate.bachDCAxy, cascadecandidate.cascDCAxy, cascadecandidate.cascDCAz); // <--- no corresponding stratrack information available + if (createCascTrackXs) { + cascTrackXs(cascadecandidate.positiveX, cascadecandidate.negativeX, cascadecandidate.bachelorX); + } + cascbb(cascadecandidate.bachBaryonCosPA, cascadecandidate.bachBaryonDCAxyToPV); + + // populate cascade covariance matrices if required by any other task + if (createCascCovMats) { + // Calculate position covariance matrix + auto covVtxV = fitter.calcPCACovMatrix(0); + // std::array positionCovariance; + float positionCovariance[6]; + positionCovariance[0] = covVtxV(0, 0); + positionCovariance[1] = covVtxV(1, 0); + positionCovariance[2] = covVtxV(1, 1); + positionCovariance[3] = covVtxV(2, 0); + positionCovariance[4] = covVtxV(2, 1); + positionCovariance[5] = covVtxV(2, 2); + // store momentum covariance matrix + std::array covTv0 = {0.}; + std::array covTbachelor = {0.}; + // std::array momentumCovariance; + float momentumCovariance[6]; + lV0Track.getCovXYZPxPyPzGlo(covTv0); + lBachelorTrack.getCovXYZPxPyPzGlo(covTbachelor); + constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int i = 0; i < 6; i++) { + momentumCovariance[i] = covTv0[MomInd[i]] + covTbachelor[MomInd[i]]; + } + casccovs(positionCovariance, momentumCovariance); + } + } + template void buildStrangenessTables(TCascTable const& cascades) { statisticsRegistry.eventCounter++; - for (auto& cascade : cascades) { // de-reference from V0 pool, either specific for cascades or general // use templatizing to avoid code duplication - bool validCascadeCandidate = false; - auto v0index = cascade.template v0_as(); - if (v0index.has_v0Data()) { - // this V0 passed both standard V0 and cascade V0 selections - auto v0row = v0index.template v0Data_as(); - validCascadeCandidate = buildCascadeCandidate(cascade, v0row); - } else if (v0index.has_v0fCData()) { - // this V0 passes only V0-for-cascade selections, use that instead - auto v0row = v0index.template v0fCData_as(); - validCascadeCandidate = buildCascadeCandidate(cascade, v0row); - } else { - continue; // this was inadequately linked, should not happen - } - if (!validCascadeCandidate) - continue; // doesn't pass cascade selections - // round the DCA variables to a certain precision if asked - if (roundDCAVariables) - roundCascadeCandidateVariables(); + auto v0index = cascade.template v0_as(); + processCascadeCandidate(v0index, cascade); + } + // En masse filling at end of process call + fillHistos(); + resetHistos(); + } - cascidx(/*cascadecandidate.v0Id, */ cascade.globalIndex(), - cascadecandidate.positiveId, cascadecandidate.negativeId, - cascadecandidate.bachelorId, cascade.collisionId()); - cascdata(cascadecandidate.charge, cascadecandidate.mXi, cascadecandidate.mOmega, - cascadecandidate.pos[0], cascadecandidate.pos[1], cascadecandidate.pos[2], - cascadecandidate.v0pos[0], cascadecandidate.v0pos[1], cascadecandidate.v0pos[2], - cascadecandidate.v0mompos[0], cascadecandidate.v0mompos[1], cascadecandidate.v0mompos[2], - cascadecandidate.v0momneg[0], cascadecandidate.v0momneg[1], cascadecandidate.v0momneg[2], - cascadecandidate.bachP[0], cascadecandidate.bachP[1], cascadecandidate.bachP[2], - cascadecandidate.bachP[0] + cascadecandidate.v0mompos[0] + cascadecandidate.v0momneg[0], // <--- redundant but ok - cascadecandidate.bachP[1] + cascadecandidate.v0mompos[1] + cascadecandidate.v0momneg[1], // <--- redundant but ok - cascadecandidate.bachP[2] + cascadecandidate.v0mompos[2] + cascadecandidate.v0momneg[2], // <--- redundant but ok - cascadecandidate.v0dcadau, cascadecandidate.dcacascdau, - cascadecandidate.v0dcapostopv, cascadecandidate.v0dcanegtopv, - cascadecandidate.bachDCAxy, cascadecandidate.cascDCAxy, cascadecandidate.cascDCAz); // <--- no corresponding stratrack information available - if (createCascTrackXs) { - cascTrackXs(cascadecandidate.positiveX, cascadecandidate.negativeX, cascadecandidate.bachelorX); - } - cascbb(cascadecandidate.bachBaryonCosPA, cascadecandidate.bachBaryonDCAxyToPV); + template + void buildFindableStrangenessTables(TCascTable const& cascades) + { + statisticsRegistry.eventCounter++; + for (auto& cascade : cascades) { + // de-reference from V0 pool, either specific for cascades or general + // use templatizing to avoid code duplication - // populate cascade covariance matrices if required by any other task - if (createCascCovMats) { - // Calculate position covariance matrix - auto covVtxV = fitter.calcPCACovMatrix(0); - // std::array positionCovariance; - float positionCovariance[6]; - positionCovariance[0] = covVtxV(0, 0); - positionCovariance[1] = covVtxV(1, 0); - positionCovariance[2] = covVtxV(1, 1); - positionCovariance[3] = covVtxV(2, 0); - positionCovariance[4] = covVtxV(2, 1); - positionCovariance[5] = covVtxV(2, 2); - // store momentum covariance matrix - std::array covTv0 = {0.}; - std::array covTbachelor = {0.}; - // std::array momentumCovariance; - float momentumCovariance[6]; - lV0Track.getCovXYZPxPyPzGlo(covTv0); - lBachelorTrack.getCovXYZPxPyPzGlo(covTbachelor); - constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component - for (int i = 0; i < 6; i++) { - momentumCovariance[i] = covTv0[MomInd[i]] + covTbachelor[MomInd[i]]; - } - casccovs(positionCovariance, momentumCovariance); - } + auto v0index = cascade.template findableV0_as(); + processCascadeCandidate(v0index, cascade); } // En masse filling at end of process call fillHistos(); @@ -1890,6 +1914,20 @@ struct cascadeBuilder { } PROCESS_SWITCH(cascadeBuilder, processRun3, "Produce Run 3 cascade tables", true); + void processFindableRun3(aod::Collisions const& collisions, aod::FindableV0sLinked const&, V0full const&, soa::Filtered const& cascades, FullTracksExtIU const&, aod::BCsWithTimestamps const&) + { + for (const auto& collision : collisions) { + // Fire up CCDB + auto bc = collision.bc_as(); + initCCDB(bc); + // Do analysis with collision-grouped V0s, retain full collision information + const uint64_t collIdx = collision.globalIndex(); + auto CascadeTable_thisCollision = cascades.sliceBy(perCollisionFindable, collIdx); + buildFindableStrangenessTables(CascadeTable_thisCollision); + } + } + PROCESS_SWITCH(cascadeBuilder, processFindableRun3, "Produce Run 3 findable cascade tables", false); + void processRun3withKFParticle(aod::Collisions const& collisions, soa::Filtered const& cascades, FullTracksExtIU const&, aod::BCsWithTimestamps const&, aod::V0s const&) { for (const auto& collision : collisions) { @@ -1970,6 +2008,11 @@ struct cascadePreselector { void init(InitContext const&) { + // check settings and stop if not viable + if (doprocessBuildAll == false && doprocessBuildMCAssociated == false && doprocessBuildValiddEdx == false && doprocessBuildValiddEdxMCAssociated == false && doprocessBuildFindable == false) { + LOGF(fatal, "No processBuild function enabled. Please choose one."); + } + auto h = histos.add("hPreselectorStatistics", "hPreselectorStatistics", kTH1D, {{5, -0.5f, 4.5f}}); h->GetXaxis()->SetBinLabel(1, "All"); h->GetXaxis()->SetBinLabel(2, "Tracks OK"); @@ -2015,11 +2058,9 @@ struct cascadePreselector { //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* /// function to check track quality - template - void checkTrackQuality(TCascadeObject const& lCascadeCandidate, uint16_t& maskElement, bool passdEdx = false) + template + void checkTrackQuality(TCascadeObject const& lCascadeCandidate, TV0Type const& v0, uint16_t& maskElement, bool passdEdx = false) { - auto v0 = lCascadeCandidate.template v0_as(); - // Acquire all three daughter tracks, please auto lBachTrack = lCascadeCandidate.template bachelor_as(); auto lNegTrack = v0.template negTrack_as(); @@ -2198,7 +2239,8 @@ struct cascadePreselector { { initializeMasks(cascades.size()); for (auto& casc : cascades) { - checkTrackQuality(casc, selectionMask[casc.globalIndex()], true); + auto v0 = casc.v0(); + checkTrackQuality(casc, v0, selectionMask[casc.globalIndex()], true); } if (!doprocessSkipCascadesNotUsedInTrackedCascades) checkAndFinalize(); @@ -2209,7 +2251,8 @@ struct cascadePreselector { initializeMasks(cascades.size()); for (auto& casc : cascades) { checkPDG(casc, selectionMask[casc.globalIndex()]); - checkTrackQuality(casc, selectionMask[casc.globalIndex()], true); + auto v0 = casc.v0(); + checkTrackQuality(casc, v0, selectionMask[casc.globalIndex()], true); } // end cascades loop if (!doprocessSkipCascadesNotUsedInTrackedCascades) checkAndFinalize(); @@ -2220,7 +2263,8 @@ struct cascadePreselector { initializeMasks(cascades.size()); for (auto& casc : cascades) { checkdEdx(casc, selectionMask[casc.globalIndex()]); - checkTrackQuality(casc, selectionMask[casc.globalIndex()]); + auto v0 = casc.v0(); + checkTrackQuality(casc, v0, selectionMask[casc.globalIndex()]); } if (!doprocessSkipCascadesNotUsedInTrackedCascades) checkAndFinalize(); @@ -2232,12 +2276,25 @@ struct cascadePreselector { for (auto& casc : cascades) { checkPDG(casc, selectionMask[casc.globalIndex()]); checkdEdx(casc, selectionMask[casc.globalIndex()]); - checkTrackQuality(casc, selectionMask[casc.globalIndex()]); + auto v0 = casc.v0(); + checkTrackQuality(casc, v0, selectionMask[casc.globalIndex()]); } if (!doprocessSkipCascadesNotUsedInTrackedCascades) checkAndFinalize(); } //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + /// This process function ensures that all findable cascades are built. + /// Not to be used with processSkip. + void processBuildFindable(aod::FindableCascades const& cascades, aod::FindableV0s const&, aod::V0Datas const&, aod::TracksExtra const&) + { + initializeMasks(cascades.size()); + for (auto& casc : cascades) { + auto v0 = casc.findableV0(); + checkTrackQuality(casc, v0, selectionMask[casc.globalIndex()], true); + } + checkAndFinalize(); + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* /// This process function checks for the use of Cascades in strangeness tracked cascades /// They are then marked appropriately; the user could then operate /// the cascadebuilder to construct only those Cascades. @@ -2254,6 +2311,7 @@ struct cascadePreselector { PROCESS_SWITCH(cascadePreselector, processBuildMCAssociated, "Switch to build MC-associated cascades", false); PROCESS_SWITCH(cascadePreselector, processBuildValiddEdx, "Switch to build cascades with dE/dx preselection", false); PROCESS_SWITCH(cascadePreselector, processBuildValiddEdxMCAssociated, "Switch to build MC-associated cascades with dE/dx preselection", false); + PROCESS_SWITCH(cascadePreselector, processBuildFindable, "Switch to build all findable cascades", false); /// skipper option (choose in addition to a processBuild if you like) PROCESS_SWITCH(cascadePreselector, processSkipCascadesNotUsedInTrackedCascades, "Switch to skip cascades not used in cascade tracking", false); //*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<* @@ -2273,7 +2331,7 @@ struct cascadeLinkBuilder { void init(InitContext const&) {} // build Cascade -> CascData link table - void process(aod::Cascades const& casctable, aod::CascDatas const& cascdatatable) + void processFound(aod::Cascades const& casctable, aod::CascDatas const& cascdatatable) { std::vector lIndices; lIndices.reserve(casctable.size()); @@ -2286,6 +2344,24 @@ struct cascadeLinkBuilder { cascdataLink(lIndices[ii]); } } + + // build Cascade -> CascData link table + void processFindable(aod::FindableCascades const& casctable, aod::CascDatas const& cascdatatable) + { + std::vector lIndices; + lIndices.reserve(casctable.size()); + for (int ii = 0; ii < casctable.size(); ii++) + lIndices[ii] = -1; + for (auto& cascdata : cascdatatable) { + lIndices[cascdata.cascadeId()] = cascdata.globalIndex(); + } + for (int ii = 0; ii < casctable.size(); ii++) { + cascdataLink(lIndices[ii]); + } + } + + PROCESS_SWITCH(cascadeLinkBuilder, processFound, "process found Cascades (default)", true); + PROCESS_SWITCH(cascadeLinkBuilder, processFindable, "process findable Cascades", false); }; struct kfcascadeLinkBuilder { diff --git a/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx b/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx index c19436a0193..a138744f43c 100644 --- a/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx @@ -55,11 +55,10 @@ struct cascademcbuilder { void init(InitContext const&) {} - //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - // build cascade labels - void processCascades(aod::CascDatas const& casctable, aod::V0sLinked const&, aod::V0Datas const& /*v0table*/, aod::McTrackLabels const&, aod::McParticles const&) + template + void generateCascadeMCinfo(TCascadeTable cascTable) { - for (auto& casc : casctable) { + for (auto& casc : cascTable) { int pdgCode = -1, pdgCodeMother = -1; int pdgCodePositive = -1, pdgCodeNegative = -1, pdgCodeBachelor = -1, pdgCodeV0 = -1; bool isPhysicalPrimary = false; @@ -72,16 +71,16 @@ struct cascademcbuilder { int lLabel = -1, lMotherLabel = -1; // Acquire all three daughter tracks, please - auto lBachTrack = casc.bachelor_as(); - auto lNegTrack = casc.negTrack_as(); - auto lPosTrack = casc.posTrack_as(); + auto lBachTrack = casc.template bachelor_as(); + auto lNegTrack = casc.template negTrack_as(); + auto lPosTrack = casc.template posTrack_as(); // Association check // There might be smarter ways of doing this in the future if (lNegTrack.has_mcParticle() && lPosTrack.has_mcParticle() && lBachTrack.has_mcParticle()) { - auto lMCBachTrack = lBachTrack.mcParticle_as(); - auto lMCNegTrack = lNegTrack.mcParticle_as(); - auto lMCPosTrack = lPosTrack.mcParticle_as(); + auto lMCBachTrack = lBachTrack.template mcParticle_as(); + auto lMCNegTrack = lNegTrack.template mcParticle_as(); + auto lMCPosTrack = lPosTrack.template mcParticle_as(); pdgCodePositive = lMCPosTrack.pdgCode(); pdgCodeNegative = lMCNegTrack.pdgCode(); @@ -98,8 +97,8 @@ struct cascademcbuilder { // Step 1: check if the mother is the same, go up a level if (lMCNegTrack.has_mothers() && lMCPosTrack.has_mothers()) { - for (auto& lNegMother : lMCNegTrack.mothers_as()) { - for (auto& lPosMother : lMCPosTrack.mothers_as()) { + for (auto& lNegMother : lMCNegTrack.template mothers_as()) { + for (auto& lPosMother : lMCPosTrack.template mothers_as()) { if (lNegMother == lPosMother) { // acquire information xlmc = lMCPosTrack.vx(); @@ -109,27 +108,29 @@ struct cascademcbuilder { // if we got to this level, it means the mother particle exists and is the same // now we have to go one level up and compare to the bachelor mother too - for (auto& lV0Mother : lNegMother.mothers_as()) { - for (auto& lBachMother : lMCBachTrack.mothers_as()) { - if (lV0Mother == lBachMother) { - lLabel = lV0Mother.globalIndex(); - pdgCode = lV0Mother.pdgCode(); - isPhysicalPrimary = lV0Mother.isPhysicalPrimary(); - xmc = lMCBachTrack.vx(); - ymc = lMCBachTrack.vy(); - zmc = lMCBachTrack.vz(); - px = lV0Mother.px(); - py = lV0Mother.py(); - pz = lV0Mother.pz(); - if (lV0Mother.has_mothers()) { - for (auto& lV0GrandMother : lV0Mother.mothers_as()) { - pdgCodeMother = lV0GrandMother.pdgCode(); - lMotherLabel = lV0GrandMother.globalIndex(); + if (lNegMother.has_mothers() && lMCBachTrack.has_mothers()) { + for (auto& lV0Mother : lNegMother.template mothers_as()) { + for (auto& lBachMother : lMCBachTrack.template mothers_as()) { + if (lV0Mother == lBachMother) { + lLabel = lV0Mother.globalIndex(); + pdgCode = lV0Mother.pdgCode(); + isPhysicalPrimary = lV0Mother.isPhysicalPrimary(); + xmc = lMCBachTrack.vx(); + ymc = lMCBachTrack.vy(); + zmc = lMCBachTrack.vz(); + px = lV0Mother.px(); + py = lV0Mother.py(); + pz = lV0Mother.pz(); + if (lV0Mother.has_mothers()) { + for (auto& lV0GrandMother : lV0Mother.template mothers_as()) { + pdgCodeMother = lV0GrandMother.pdgCode(); + lMotherLabel = lV0GrandMother.globalIndex(); + } } } } - } - } // end conditional V0-bach pair + } // end conditional V0-bach pair + } // end has mothers } // end neg = pos mother conditional } } // end loop neg/pos mothers @@ -151,6 +152,20 @@ struct cascademcbuilder { } // end casctable loop } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + // build cascade labels + void processCascades(aod::CascDatas const& casctable, aod::V0sLinked const&, aod::V0Datas const& /*v0table*/, aod::McTrackLabels const&, aod::McParticles const&) + { + generateCascadeMCinfo(casctable); + } + + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + // build findable cascade labels + void processFindableCascades(aod::CascDatas const& casctable, aod::FindableV0sLinked const&, aod::V0Datas const& /*v0table*/, aod::McTrackLabels const&, aod::McParticles const&) + { + generateCascadeMCinfo(casctable); + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* // build kf cascade labels void processKFCascades(aod::KFCascDatas const& casctable, aod::V0s const&, aod::McTrackLabels const&, aod::McParticles const&) @@ -159,25 +174,25 @@ struct cascademcbuilder { int lLabel = -1; // Acquire all three daughter tracks, please - auto lBachTrack = casc.bachelor_as(); - auto lNegTrack = casc.negTrack_as(); - auto lPosTrack = casc.posTrack_as(); + auto lBachTrack = casc.template bachelor_as(); + auto lNegTrack = casc.template negTrack_as(); + auto lPosTrack = casc.template posTrack_as(); // Association check // There might be smarter ways of doing this in the future if (lNegTrack.has_mcParticle() && lPosTrack.has_mcParticle() && lBachTrack.has_mcParticle()) { - auto lMCBachTrack = lBachTrack.mcParticle_as(); - auto lMCNegTrack = lNegTrack.mcParticle_as(); - auto lMCPosTrack = lPosTrack.mcParticle_as(); + auto lMCBachTrack = lBachTrack.template mcParticle_as(); + auto lMCNegTrack = lNegTrack.template mcParticle_as(); + auto lMCPosTrack = lPosTrack.template mcParticle_as(); // Step 1: check if the mother is the same, go up a level if (lMCNegTrack.has_mothers() && lMCPosTrack.has_mothers()) { - for (auto& lNegMother : lMCNegTrack.mothers_as()) { - for (auto& lPosMother : lMCPosTrack.mothers_as()) { + for (auto& lNegMother : lMCNegTrack.template mothers_as()) { + for (auto& lPosMother : lMCPosTrack.template mothers_as()) { if (lNegMother == lPosMother) { // if we got to this level, it means the mother particle exists and is the same // now we have to go one level up and compare to the bachelor mother too - for (auto& lV0Mother : lNegMother.mothers_as()) { - for (auto& lBachMother : lMCBachTrack.mothers_as()) { + for (auto& lV0Mother : lNegMother.template mothers_as()) { + for (auto& lBachMother : lMCBachTrack.template mothers_as()) { if (lV0Mother == lBachMother) { lLabel = lV0Mother.globalIndex(); } @@ -289,6 +304,7 @@ struct cascademcbuilder { } PROCESS_SWITCH(cascademcbuilder, processCascades, "Produce regular cascade label tables", true); + PROCESS_SWITCH(cascademcbuilder, processFindableCascades, "Produce findable cascade label tables", false); PROCESS_SWITCH(cascademcbuilder, processKFCascades, "Produce KF cascade label tables", false); PROCESS_SWITCH(cascademcbuilder, processTrackedCascades, "Produce tracked cascade label tables", false); PROCESS_SWITCH(cascademcbuilder, processBBTags, "Produce cascade bach-baryon correlation tags", true); diff --git a/PWGLF/TableProducer/Strangeness/cascademcfinder.cxx b/PWGLF/TableProducer/Strangeness/cascademcfinder.cxx index 46c4933cdab..8d26f0ea33c 100644 --- a/PWGLF/TableProducer/Strangeness/cascademcfinder.cxx +++ b/PWGLF/TableProducer/Strangeness/cascademcfinder.cxx @@ -24,6 +24,11 @@ // david.dobrigkeit.chinellato@cern.ch // +#include +#include +#include +#include + #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" #include "Framework/AnalysisDataModel.h" @@ -40,7 +45,6 @@ #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Centrality.h" #include "DataFormatsParameters/GRPObject.h" -#include "DataFormatsParameters/GRPObject.h" #include "DataFormatsParameters/GRPMagField.h" #include "CCDB/BasicCCDBManager.h" @@ -49,12 +53,8 @@ #include #include #include -#include #include #include -#include -#include -#include using namespace o2; using namespace o2::framework; @@ -62,11 +62,12 @@ using namespace o2::framework::expressions; using std::array; using namespace ROOT::Math; +// WARNING: the cascade findable uses findable V0s as well using LabeledTracks = soa::Join; -using LabeledFullV0s = soa::Join; +using LabeledFullV0s = soa::Join; struct cascademcfinder { - Produces cascades; + Produces cascades; HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -89,8 +90,8 @@ struct cascademcfinder { void init(InitContext&) { // initialize histograms - const AxisSpec axisNTimesCollRecoed{(int)10, -0.5f, +9.5f, ""}; - const AxisSpec axisPt{(int)100, +0.0f, +10.0f, "p_{T} (GeV/c)"}; + const AxisSpec axisNTimesCollRecoed{static_cast(10), -0.5f, +9.5f, ""}; + const AxisSpec axisPt{static_cast(100), +0.0f, +10.0f, "p_{T} (GeV/c)"}; histos.add("hNTimesCollRecoed", "hNTimesCollRecoed", kTH1F, {axisNTimesCollRecoed}); diff --git a/PWGLF/TableProducer/Strangeness/lambdakzerobuilder.cxx b/PWGLF/TableProducer/Strangeness/lambdakzerobuilder.cxx index 67304d45ceb..0c6a0b843a2 100644 --- a/PWGLF/TableProducer/Strangeness/lambdakzerobuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/lambdakzerobuilder.cxx @@ -90,6 +90,7 @@ using TracksExtraWithPIDandLabels = soa::Join; +using TaggedFindableV0s = soa::Join; // For MC association in pre-selection using LabeledTracksExtra = soa::Join; @@ -528,12 +529,18 @@ struct lambdakzeroBuilder { LOGF(info, "LUT load done!"); } - if (doprocessRun2 == false && doprocessRun3 == false) { - LOGF(fatal, "Neither processRun2 nor processRun3 enabled. Please choose one."); + if (doprocessRun2 == false && doprocessRun3 == false && doprocessFindableRun3 == false) { + LOGF(fatal, "Neither processRun2 nor processRun3 nor processFindableRun3 enabled. Please choose one."); } if (doprocessRun2 == true && doprocessRun3 == true) { LOGF(fatal, "Cannot enable processRun2 and processRun3 at the same time. Please choose one."); } + if (doprocessRun2 == true && doprocessFindableRun3 == true) { + LOGF(fatal, "Cannot enable processRun2 and processFindableRun3 at the same time. Please choose one."); + } + if (doprocessRun3 == true && doprocessFindableRun3 == true) { + LOGF(fatal, "Cannot enable processRun3 and processFindableRun3 at the same time. Please choose one."); + } if (d_UseAutodetectMode) { double loosest_v0cospa = 100; @@ -1283,6 +1290,20 @@ struct lambdakzeroBuilder { buildStrangenessTables(V0s); } PROCESS_SWITCH(lambdakzeroBuilder, processRun3, "Produce Run 3 V0 tables", true); + + void processFindableRun3(aod::Collisions const& collisions, soa::Filtered const& V0s, FullTracksExtIU const&, aod::BCsWithTimestamps const& bcs) + { + statisticsRegistry.eventCounter += collisions.size(); + // Fire up CCDB + auto bc = collisions.size() ? collisions.begin().bc_as() : bcs.begin(); + if (!bcs.size()) { + LOGF(warn, "No BC found, skipping this DF."); + return; + } + initCCDB(bc); + buildStrangenessTables(V0s); + } + PROCESS_SWITCH(lambdakzeroBuilder, processFindableRun3, "Produce Run 3 V0 tables with all findable candidates", false); }; //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* @@ -1342,6 +1363,16 @@ struct lambdakzeroPreselector { void init(InitContext const&) { + // check settings and stop if not viable + if (doprocessBuildAll == false && doprocessBuildMCAssociated == false && doprocessBuildValiddEdx == false && doprocessBuildValiddEdxMCAssociated == false && doprocessBuildFindable == false) { + LOGF(fatal, "No processBuild function enabled. Please choose one."); + } + LOGF(info, "Process function processBuildAll status: %i", static_cast(doprocessBuildAll)); + LOGF(info, "Process function doprocessBuildMCAssociated status: %i", static_cast(doprocessBuildMCAssociated)); + LOGF(info, "Process function doprocessBuildValiddEdx status: %i", static_cast(doprocessBuildValiddEdx)); + LOGF(info, "Process function doprocessBuildValiddEdxMCAssociated status: %i", static_cast(doprocessBuildValiddEdxMCAssociated)); + LOGF(info, "Process function doprocessBuildFindable status: %i", static_cast(doprocessBuildFindable)); + auto h = histos.add("hPreselectorStatistics", "hPreselectorStatistics", kTH1D, {{6, -0.5f, 5.5f}}); h->GetXaxis()->SetBinLabel(1, "All"); h->GetXaxis()->SetBinLabel(2, "Tracks OK"); @@ -1581,6 +1612,16 @@ struct lambdakzeroPreselector { checkAndFinalize(); } //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + /// This process function ensures that all findable V0s are built. + void processBuildFindable(aod::FindableV0s const& v0table, aod::TracksExtra const&) + { + initializeMasks(v0table.size()); + for (auto const& v0 : v0table) { + checkTrackQuality(v0, selectionMask[v0.globalIndex()], true); + } + checkAndFinalize(); + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* /// This process function checks for the use of V0s in cascades /// They are then marked appropriately; the user could then operate /// the lambdakzerobuilder to construct only those V0s. @@ -1609,6 +1650,7 @@ struct lambdakzeroPreselector { PROCESS_SWITCH(lambdakzeroPreselector, processBuildMCAssociated, "Switch to build MC-associated V0s", false); PROCESS_SWITCH(lambdakzeroPreselector, processBuildValiddEdx, "Switch to build V0s with dE/dx preselection", false); PROCESS_SWITCH(lambdakzeroPreselector, processBuildValiddEdxMCAssociated, "Switch to build MC-associated V0s with dE/dx preselection", false); + PROCESS_SWITCH(lambdakzeroPreselector, processBuildFindable, "Switch to build findable V0s. Requires lambdakzeromcfinder", false); /// skippers options (choose one in addition to a processBuild if you like) PROCESS_SWITCH(lambdakzeroPreselector, processSkipV0sNotUsedInCascades, "skip all V0s not used in cascades", false); PROCESS_SWITCH(lambdakzeroPreselector, processSkipV0sNotUsedInTrackedCascades, "skip all V0s not used in tracked cascades", false); @@ -1623,7 +1665,7 @@ struct lambdakzeroV0DataLinkBuilder { //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* // build V0 -> V0Data link table - void process(aod::V0s const& v0table, aod::V0Datas const& v0datatable, aod::V0fCDatas const& v0fcdatatable) + void processFound(aod::V0s const& v0table, aod::V0Datas const& v0datatable, aod::V0fCDatas const& v0fcdatatable) { std::vector lIndices, lfCIndices; lIndices.reserve(v0table.size()); @@ -1643,6 +1685,24 @@ struct lambdakzeroV0DataLinkBuilder { } } //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + // build V0Findable -> V0Data link table + void processFindable(aod::FindableV0s const& v0table, aod::V0Datas const& v0datatable) + { + std::vector lIndices; + lIndices.reserve(v0table.size()); + for (int ii = 0; ii < v0table.size(); ii++) { + lIndices[ii] = -1; + } + for (auto& v0data : v0datatable) { + lIndices[v0data.v0Id()] = v0data.globalIndex(); + } + for (int ii = 0; ii < v0table.size(); ii++) { + v0dataLink(lIndices[ii], -1); + } + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + PROCESS_SWITCH(lambdakzeroV0DataLinkBuilder, processFound, "process found V0s (default)", true); + PROCESS_SWITCH(lambdakzeroV0DataLinkBuilder, processFindable, "process findable V0s", false); }; // Extends the v0data table with expression columns diff --git a/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx b/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx index 2df904a6274..82918346104 100644 --- a/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx @@ -44,54 +44,98 @@ using std::array; struct lambdakzeromcbuilder { Produces v0labels; // MC labels for V0s Produces v0mccores; // optionally aggregate information from MC side for posterior analysis (derived data) + Produces v0CoreMCLabels; // interlink V0Cores -> V0MCCores in asymmetric mode - Configurable populateV0MCCores{"populateV0MCCores", false, "populate V0MCCores table for derived data analysis"}; + Configurable populateV0MCCoresSymmetric{"populateV0MCCoresSymmetric", false, "populate V0MCCores table for derived data analysis, keep V0MCCores joinable with V0Cores"}; + Configurable populateV0MCCoresAsymmetric{"populateV0MCCoresAsymmetric", false, "populate V0MCCores table for derived data analysis, create V0Cores -> V0MCCores interlink. Saves only labeled V0s."}; - void init(InitContext const&) {} + HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + void init(InitContext const&) + { + if (populateV0MCCoresAsymmetric) { + LOGF(info, "Asymmetric V0MCCores filling enabled!"); + } + if (populateV0MCCoresSymmetric) { + LOGF(info, "Symmetric V0MCCores filling enabled!"); + } + if (populateV0MCCoresAsymmetric && populateV0MCCoresSymmetric) { + LOGF(fatal, "Error in configuration: please select only one out of populateV0MCCoresAsymmetric and populateV0MCCoresSymmetric! Crashing!"); + } + + // for storing basic statistics + auto h = histos.add("hBuildingStatistics", "hBuildingStatistics", kTH1F, {{4, -0.5, 3.5f}}); + h->GetXaxis()->SetBinLabel(1, "V0Cores population"); + h->GetXaxis()->SetBinLabel(2, "V0MCCores population"); + h->GetXaxis()->SetBinLabel(3, "x check: duplicates"); + h->GetXaxis()->SetBinLabel(4, "x check: unique"); + } + + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + // Helper struct to contain V0MCCore information prior to filling + struct mcV0info { + int label; + int motherLabel; + int pdgCode; + int pdgCodeMother; + int pdgCodePositive; + int pdgCodeNegative; + bool isPhysicalPrimary; + std::array xyz; + std::array posP; + std::array negP; + uint64_t packedMcParticleIndices; + }; + mcV0info thisInfo; + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + + // prong index combiner + uint64_t combineProngIndices(uint32_t low, uint32_t high) + { + return (((uint64_t)high) << 32) | ((uint64_t)low); + } //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* // build V0 labels void process(aod::V0Datas const& v0table, aod::McTrackLabels const&, aod::McParticles const& /*particlesMC*/) { - for (auto& v0 : v0table) { - int lLabel = -1, lMotherLabel = -1; - int pdgCode = -1, pdgCodeMother = -1, pdgCodePositive = -1, pdgCodeNegative = -1; - bool isPhysicalPrimary = false; - float xmc = -999.0f, ymc = -999.0f, zmc = -999.0f; - float pxposmc = -999.0f, pyposmc = -999.0f, pzposmc = -999.0f; - float pxnegmc = -999.0f, pynegmc = -999.0f, pznegmc = -999.0f; + // to be used if using the populateV0MCCoresAsymmetric mode, kept empty otherwise + std::vector mcV0infos; // V0MCCore information + for (auto& v0 : v0table) { + thisInfo.packedMcParticleIndices = 0; // not de-referenced properly yet auto lNegTrack = v0.negTrack_as(); auto lPosTrack = v0.posTrack_as(); // Association check // There might be smarter ways of doing this in the future if (lNegTrack.has_mcParticle() && lPosTrack.has_mcParticle()) { + thisInfo.packedMcParticleIndices = combineProngIndices(lPosTrack.mcParticleId(), lNegTrack.mcParticleId()); auto lMCNegTrack = lNegTrack.mcParticle_as(); auto lMCPosTrack = lPosTrack.mcParticle_as(); - pdgCodePositive = lMCPosTrack.pdgCode(); - pdgCodeNegative = lMCNegTrack.pdgCode(); - pxposmc = lMCPosTrack.px(); - pyposmc = lMCPosTrack.py(); - pzposmc = lMCPosTrack.pz(); - pxnegmc = lMCNegTrack.px(); - pynegmc = lMCNegTrack.py(); - pznegmc = lMCNegTrack.pz(); + thisInfo.pdgCodePositive = lMCPosTrack.pdgCode(); + thisInfo.pdgCodeNegative = lMCNegTrack.pdgCode(); + thisInfo.posP[0] = lMCPosTrack.px(); + thisInfo.posP[1] = lMCPosTrack.py(); + thisInfo.posP[2] = lMCPosTrack.pz(); + thisInfo.negP[0] = lMCNegTrack.px(); + thisInfo.negP[1] = lMCNegTrack.py(); + thisInfo.negP[2] = lMCNegTrack.pz(); if (lMCNegTrack.has_mothers() && lMCPosTrack.has_mothers()) { for (auto& lNegMother : lMCNegTrack.mothers_as()) { for (auto& lPosMother : lMCPosTrack.mothers_as()) { if (lNegMother.globalIndex() == lPosMother.globalIndex()) { - lLabel = lNegMother.globalIndex(); + thisInfo.label = lNegMother.globalIndex(); // acquire information - xmc = lMCPosTrack.vx(); - ymc = lMCPosTrack.vy(); - zmc = lMCPosTrack.vz(); - pdgCode = lNegMother.pdgCode(); - isPhysicalPrimary = lNegMother.isPhysicalPrimary(); + thisInfo.xyz[0] = lMCPosTrack.vx(); + thisInfo.xyz[1] = lMCPosTrack.vy(); + thisInfo.xyz[2] = lMCPosTrack.vz(); + thisInfo.pdgCode = lNegMother.pdgCode(); + thisInfo.isPhysicalPrimary = lNegMother.isPhysicalPrimary(); if (lNegMother.has_mothers()) { for (auto& lNegGrandMother : lNegMother.mothers_as()) { - pdgCodeMother = lNegGrandMother.pdgCode(); - lMotherLabel = lNegGrandMother.globalIndex(); + thisInfo.pdgCodeMother = lNegGrandMother.pdgCode(); + thisInfo.motherLabel = lNegGrandMother.globalIndex(); } } } @@ -101,15 +145,65 @@ struct lambdakzeromcbuilder { } // end association check // Construct label table (note: this will be joinable with V0Datas!) v0labels( - lLabel, lMotherLabel); - if (populateV0MCCores) { + thisInfo.label, thisInfo.motherLabel); + + // ---] Symmetric populate [--- + // in this approach, V0Cores will be joinable with V0MCCores. + // this is the most pedagogical approach, but it is also more limited + // and it might use more disk space unnecessarily. + if (populateV0MCCoresSymmetric) { v0mccores( - lLabel, pdgCode, pdgCodeMother, pdgCodePositive, pdgCodeNegative, - isPhysicalPrimary, xmc, ymc, zmc, - pxposmc, pyposmc, pzposmc, - pxnegmc, pynegmc, pznegmc); + thisInfo.label, thisInfo.pdgCode, + thisInfo.pdgCodeMother, thisInfo.pdgCodePositive, thisInfo.pdgCodeNegative, + thisInfo.isPhysicalPrimary, thisInfo.xyz[0], thisInfo.xyz[1], thisInfo.xyz[2], + thisInfo.posP[0], thisInfo.posP[1], thisInfo.posP[2], + thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2]); + + // n.b. placing the interlink index here allows for the writing of + // code that is agnostic with respect to the joinability of + // V0Cores and V0MCCores (always dereference -> safe) + v0CoreMCLabels(v0.globalIndex()); // interlink index + } + // ---] Asymmetric populate [--- + // in this approach, V0Cores will NOT be joinable with V0MCCores. + // an additional reference to V0MCCore that IS joinable with V0Cores + // will be provided to the user. + if (populateV0MCCoresAsymmetric) { + int thisV0MCCoreIndex = -1; + // step 1: check if this element is already provided in the table + // using the packedIndices variable calculated above + for (uint32_t ii = 0; ii < mcV0infos.size(); ii++) { + if (thisInfo.packedMcParticleIndices == mcV0infos[ii].packedMcParticleIndices && mcV0infos[ii].packedMcParticleIndices > 0) { + thisV0MCCoreIndex = ii; + histos.fill(HIST("hBuildingStatistics"), 2.0f); // found + break; // this exists already in list + } + } + if (thisV0MCCoreIndex < 0) { + // this V0MCCore does not exist yet. Create it and reference it + histos.fill(HIST("hBuildingStatistics"), 3.0f); // new + thisV0MCCoreIndex = mcV0infos.size(); + mcV0infos.push_back(thisInfo); + } + v0CoreMCLabels(thisV0MCCoreIndex); // interlink index + } + } + + // now populate V0MCCores if in asymmetric mode + if (populateV0MCCoresAsymmetric) { + for (auto info : mcV0infos) { + v0mccores( + info.label, info.pdgCode, + info.pdgCodeMother, info.pdgCodePositive, info.pdgCodeNegative, + info.isPhysicalPrimary, info.xyz[0], info.xyz[1], info.xyz[2], + info.posP[0], info.posP[1], info.posP[2], + info.negP[0], info.negP[1], info.negP[2]); } } + + // collect operating parameters + histos.fill(HIST("hBuildingStatistics"), 0.0f, v0table.size()); + histos.fill(HIST("hBuildingStatistics"), 1.0f, mcV0infos.size()); } }; diff --git a/PWGLF/TableProducer/Strangeness/lambdakzeromcfinder.cxx b/PWGLF/TableProducer/Strangeness/lambdakzeromcfinder.cxx index a7d898bc0f2..5e0c7494282 100644 --- a/PWGLF/TableProducer/Strangeness/lambdakzeromcfinder.cxx +++ b/PWGLF/TableProducer/Strangeness/lambdakzeromcfinder.cxx @@ -68,7 +68,7 @@ using LabeledTracks = soa::Join; struct lambdakzeromcfinder { - Produces v0; + Produces v0; Produces fullv0labels; HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -178,9 +178,8 @@ struct lambdakzeromcfinder { { int nPosReco = 0; int nNegReco = 0; - const int maxReco = 20; - int trackIndexPositive[maxReco]; - int trackIndexNegative[maxReco]; + int trackIndexPositive[20]; + int trackIndexNegative[20]; int positivePdg = 211; int negativePdg = -211; diff --git a/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx b/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx index f88e72a8f2b..fe0880e7031 100644 --- a/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx @@ -125,6 +125,11 @@ struct strangederivedbuilder { Produces geOmegaMinus; Produces geOmegaPlus; + //__________________________________________________ + // Found tags for findable exercise + Produces v0FoundTags; + Produces cascFoundTags; + // histogram registry for bookkeeping HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -213,6 +218,16 @@ struct strangederivedbuilder { histos.add("h2dNVerticesVsCentrality", "h2dNVerticesVsCentrality", kTH2D, {axisCentrality, axisNVertices}); + if (doprocessV0FoundTags || doprocessCascFoundTags) { + auto h = histos.add("hFoundTagsCounters", "hFoundTagsCounters", kTH1D, {{6, -0.5f, 5.5f}}); + h->GetXaxis()->SetBinLabel(1, "Found V0s"); + h->GetXaxis()->SetBinLabel(2, "Findable V0s"); + h->GetXaxis()->SetBinLabel(3, "Findable & found V0s"); + h->GetXaxis()->SetBinLabel(4, "Found Cascades"); + h->GetXaxis()->SetBinLabel(5, "Findable Cascades"); + h->GetXaxis()->SetBinLabel(6, "Findable & found Cascades"); + } + // for QA and test purposes auto hRawCentrality = histos.add("hRawCentrality", "hRawCentrality", kTH1F, {axisRawCentrality}); @@ -783,6 +798,70 @@ struct strangederivedbuilder { StraFV0AQVs(collision.qvecFV0ARe(), collision.qvecFV0AIm(), collision.sumAmplFV0A()); } + uint64_t combineProngIndices(uint32_t low, uint32_t high) + { + return (((uint64_t)high) << 32) | ((uint64_t)low); + } + + void processV0FoundTags(aod::V0s const& foundV0s, aod::V0Datas const& findableV0s, aod::FindableV0s const& /* added to avoid troubles */) + { + histos.fill(HIST("hFoundTagsCounters"), 0.0f, foundV0s.size()); + histos.fill(HIST("hFoundTagsCounters"), 1.0f, findableV0s.size()); + + // pack the found V0s in a long long + std::vector foundV0sPacked; + foundV0sPacked.reserve(foundV0s.size()); + for (auto const& foundV0 : foundV0s) { + foundV0sPacked[foundV0.globalIndex()] = combineProngIndices(foundV0.posTrackId(), foundV0.negTrackId()); + } + + bool hasBeenFound = false; + for (auto const& findableV0 : findableV0s) { + uint64_t indexPack = combineProngIndices(findableV0.posTrackId(), findableV0.negTrackId()); + for (uint32_t ic = 0; ic < foundV0s.size(); ic++) { + if (indexPack == foundV0sPacked[ic]) { + hasBeenFound = true; + histos.fill(HIST("hFoundTagsCounters"), 2.0f); + break; + } + } + v0FoundTags(hasBeenFound); + } + } + + using uint128_t = __uint128_t; + uint128_t combineProngIndices128(uint32_t pos, uint32_t neg, uint32_t bach) + { + return (((uint128_t)pos) << 64) | (((uint128_t)neg) << 32) | ((uint128_t)bach); + } + + void processCascFoundTags(aod::Cascades const& foundCascades, aod::CascDatas const& findableCascades, aod::V0s const&, aod::FindableCascades const& /* added to avoid troubles */) + { + histos.fill(HIST("hFoundTagsCounters"), 3.0f, foundCascades.size()); + histos.fill(HIST("hFoundTagsCounters"), 4.0f, findableCascades.size()); + + // pack the found V0s in a long long + std::vector foundCascadesPacked; + foundCascadesPacked.reserve(foundCascades.size()); + for (auto const& foundCascade : foundCascades) { + auto v0 = foundCascade.v0(); + foundCascadesPacked[foundCascade.globalIndex()] = combineProngIndices128(v0.posTrackId(), v0.negTrackId(), foundCascade.bachelorId()); + } + + bool hasBeenFound = false; + for (auto const& findableCascade : findableCascades) { + uint128_t indexPack = combineProngIndices128(findableCascade.posTrackId(), findableCascade.negTrackId(), findableCascade.bachelorId()); + for (uint32_t ic = 0; ic < foundCascades.size(); ic++) { + if (indexPack == foundCascadesPacked[ic]) { + hasBeenFound = true; + histos.fill(HIST("hFoundTagsCounters"), 5.0f); + break; + } + } + cascFoundTags(hasBeenFound); + } + } + PROCESS_SWITCH(strangederivedbuilder, processCollisionsV0sOnly, "Produce collisions (V0s only)", true); PROCESS_SWITCH(strangederivedbuilder, processCollisions, "Produce collisions (V0s + casc)", true); PROCESS_SWITCH(strangederivedbuilder, processCollisionsMC, "Produce collisions (V0s + casc)", false); @@ -795,11 +874,17 @@ struct strangederivedbuilder { PROCESS_SWITCH(strangederivedbuilder, processPureSimulation, "Produce pure simulated information", true); PROCESS_SWITCH(strangederivedbuilder, processReconstructedSimulation, "Produce reco-ed simulated information", true); PROCESS_SWITCH(strangederivedbuilder, processBinnedGenerated, "Produce binned generated information", false); + + // event plane information PROCESS_SWITCH(strangederivedbuilder, processFT0AQVectors, "Produce FT0A Q-vectors table", false); PROCESS_SWITCH(strangederivedbuilder, processFT0CQVectors, "Produce FT0C Q-vectors table", false); PROCESS_SWITCH(strangederivedbuilder, processFT0CQVectorsLF, "Produce FT0C Q-vectors table using LF temporary calibration", false); PROCESS_SWITCH(strangederivedbuilder, processFT0MQVectors, "Produce FT0M Q-vectors table", false); PROCESS_SWITCH(strangederivedbuilder, processFV0AQVectors, "Produce FV0A Q-vectors table", false); + + // dedicated findable functionality + PROCESS_SWITCH(strangederivedbuilder, processV0FoundTags, "Produce FoundV0Tags for findable exercise", false); + PROCESS_SWITCH(strangederivedbuilder, processCascFoundTags, "Produce FoundCascTags for findable exercise", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)