From b07803e28fda80882b8cfabe4440519bbe22d8ff Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Mon, 17 Jun 2024 18:26:09 +0200 Subject: [PATCH 1/2] PWGLF: mc core for all gen --- PWGLF/DataModel/LFStrangenessTables.h | 7 + .../Strangeness/cascademcbuilder.cxx | 249 ++++++++++++++---- .../Strangeness/lambdakzeromcbuilder.cxx | 72 ++++- .../Strangeness/strangederivedbuilder.cxx | 28 +- 4 files changed, 275 insertions(+), 81 deletions(-) diff --git a/PWGLF/DataModel/LFStrangenessTables.h b/PWGLF/DataModel/LFStrangenessTables.h index 03a916d1b06..9baca729fb4 100644 --- a/PWGLF/DataModel/LFStrangenessTables.h +++ b/PWGLF/DataModel/LFStrangenessTables.h @@ -1190,6 +1190,13 @@ DECLARE_SOA_TABLE(CascMCCores, "AOD", "CASCMCCORE", //! bachelor-baryon correlat cascdata::PxBachMC, cascdata::PyBachMC, cascdata::PzBachMC, cascdata::PxMC, cascdata::PyMC, cascdata::PzMC); +namespace cascdata +{ +DECLARE_SOA_INDEX_COLUMN(CascMCCore, cascMCCore); //! Index to CascMCCore entry +} + +DECLARE_SOA_TABLE(CascCoreMCLabels, "AOD", "CASCCOREMCLABEL", //! optional table to refer to CascMCCores if not joinable + o2::soa::Index<>, cascdata::CascMCCoreId); DECLARE_SOA_TABLE(CascMCCollRefs, "AOD", "CASCMCCOLLREF", //! refers MC candidate back to proper MC Collision o2::soa::Index<>, cascdata::StraMCCollisionId, o2::soa::Marker<3>); diff --git a/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx b/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx index a138744f43c..e57a9959a82 100644 --- a/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx @@ -50,25 +50,70 @@ struct cascademcbuilder { Produces tracasclabels; // MC labels for tracked cascades Produces bbtags; // bb tags (inv structure tagging) Produces cascmccores; // optionally aggregate information from MC side for posterior analysis (derived data) + Produces cascCoreMClabels; // optionally aggregate information from MC side for posterior analysis (derived data) + Produces cascmccollrefs; // references MC collisions from MC cascades - Configurable populateCascMCCores{"populateCascMCCores", false, "populate CascMCCores table for derived data analysis"}; + Configurable populateCascMCCoresSymmetric{"populateCascMCCoresSymmetric", false, "populate CascMCCores table for derived data analysis, keep CascMCCores joinable with CascCores"}; + Configurable populateCascMCCoresAsymmetric{"populateCascMCCoresAsymmetric", false, "populate CascMCCores table for derived data analysis, create CascCores -> CascMCCores interlink. Saves only labeled Cascades."}; + + Configurable addGeneratedXiMinus{"addGeneratedXiMinus", false, "add CascMCCore entry for generated, not-recoed XiMinus"}; + Configurable addGeneratedXiPlus{"addGeneratedXiPlus", false, "add CascMCCore entry for generated, not-recoed XiPlus"}; + Configurable addGeneratedOmegaMinus{"addGeneratedOmegaMinus", false, "add CascMCCore entry for generated, not-recoed OmegaMinus"}; + Configurable addGeneratedOmegaPlus{"addGeneratedOmegaPlus", false, "add CascMCCore entry for generated, not-recoed OmegaPlus"}; + + Configurable rapidityWindow{"rapidityWindow", 0.5, "rapidity window to save non-recoed candidates"}; + + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + // Helper struct to contain CascMCCore information prior to filling + struct mcCascinfo { + int label; + int motherLabel; + int mcCollision; + int pdgCode; + int pdgCodeMother; + int pdgCodeV0; + int pdgCodePositive; + int pdgCodeNegative; + int pdgCodeBachelor; + bool isPhysicalPrimary; + std::array xyz; + std::array lxyz; + std::array posP; + std::array negP; + std::array bachP; + std::array momentum; + int mcParticlePositive; + int mcParticleNegative; + int mcParticleBachelor; + }; + mcCascinfo thisInfo; + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* void init(InitContext const&) {} - template - void generateCascadeMCinfo(TCascadeTable cascTable) + template + void generateCascadeMCinfo(TCascadeTable cascTable, TMCParticleTable mcParticles) { + + // to be used if using the asymmetric mode, kept empty otherwise + std::vector mcCascinfos; // V0MCCore information + std::vector mcParticleIsReco(mcParticles.size(), false); // mc Particle not recoed by V0s + for (auto& casc : cascTable) { - int pdgCode = -1, pdgCodeMother = -1; - int pdgCodePositive = -1, pdgCodeNegative = -1, pdgCodeBachelor = -1, pdgCodeV0 = -1; - bool isPhysicalPrimary = false; - float xmc = -999.0f, ymc = -999.0f, zmc = -999.0f; - float xlmc = -999.0f, ylmc = -999.0f, zlmc = -999.0f; - float pxposmc = -999.0f, pyposmc = -999.0f, pzposmc = -999.0f; - float pxnegmc = -999.0f, pynegmc = -999.0f, pznegmc = -999.0f; - float pxbachmc = -999.0f, pybachmc = -999.0f, pzbachmc = -999.0f; - float px = -999.0f, py = -999.0f, pz = -999.0f; - int lLabel = -1, lMotherLabel = -1; + thisInfo.pdgCode = -1, thisInfo.pdgCodeMother = -1; + thisInfo.pdgCodePositive = -1, thisInfo.pdgCodeNegative = -1; + thisInfo.pdgCodeBachelor = -1, thisInfo.pdgCodeV0 = -1; + thisInfo.isPhysicalPrimary = false; + thisInfo.xyz[0] = -999.0f, thisInfo.xyz[1] = -999.0f, thisInfo.xyz[2] = -999.0f; + thisInfo.lxyz[0] = -999.0f, thisInfo.lxyz[1] = -999.0f, thisInfo.lxyz[2] = -999.0f; + thisInfo.posP[0] = -999.0f, thisInfo.posP[1] = -999.0f, thisInfo.posP[2] = -999.0f; + thisInfo.negP[0] = -999.0f, thisInfo.negP[1] = -999.0f, thisInfo.negP[2] = -999.0f; + thisInfo.bachP[0] = -999.0f, thisInfo.bachP[1] = -999.0f, thisInfo.bachP[2] = -999.0f; + thisInfo.momentum[0] = -999.0f, thisInfo.momentum[1] = -999.0f, thisInfo.momentum[2] = -999.0f; + thisInfo.label = -1, thisInfo.motherLabel = -1; + thisInfo.mcParticlePositive=-1; + thisInfo.mcParticleNegative=-1; + thisInfo.mcParticleBachelor=-1; // Acquire all three daughter tracks, please auto lBachTrack = casc.template bachelor_as(); @@ -82,18 +127,21 @@ struct cascademcbuilder { auto lMCNegTrack = lNegTrack.template mcParticle_as(); auto lMCPosTrack = lPosTrack.template mcParticle_as(); - pdgCodePositive = lMCPosTrack.pdgCode(); - pdgCodeNegative = lMCNegTrack.pdgCode(); - pdgCodeBachelor = lMCBachTrack.pdgCode(); - pxposmc = lMCPosTrack.px(); - pyposmc = lMCPosTrack.py(); - pzposmc = lMCPosTrack.pz(); - pxnegmc = lMCNegTrack.px(); - pynegmc = lMCNegTrack.py(); - pznegmc = lMCNegTrack.pz(); - pxbachmc = lMCBachTrack.px(); - pybachmc = lMCBachTrack.py(); - pzbachmc = lMCBachTrack.pz(); + thisInfo.mcParticlePositive = lMCPosTrack.globalIndex(); + thisInfo.mcParticleNegative = lMCNegTrack.globalIndex(); + thisInfo.mcParticleBachelor = lMCBachTrack.globalIndex(); + thisInfo.pdgCodePositive = lMCPosTrack.pdgCode(); + thisInfo.pdgCodeNegative = lMCNegTrack.pdgCode(); + thisInfo.pdgCodeBachelor = lMCBachTrack.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(); + thisInfo.bachP[0] = lMCBachTrack.px(); + thisInfo.bachP[1] = lMCBachTrack.py(); + thisInfo.bachP[2] = lMCBachTrack.pz(); // Step 1: check if the mother is the same, go up a level if (lMCNegTrack.has_mothers() && lMCPosTrack.has_mothers()) { @@ -101,10 +149,10 @@ struct cascademcbuilder { for (auto& lPosMother : lMCPosTrack.template mothers_as()) { if (lNegMother == lPosMother) { // acquire information - xlmc = lMCPosTrack.vx(); - ylmc = lMCPosTrack.vy(); - zlmc = lMCPosTrack.vz(); - pdgCodeV0 = lNegMother.pdgCode(); + thisInfo.lxyz[0] = lMCPosTrack.vx(); + thisInfo.lxyz[1] = lMCPosTrack.vy(); + thisInfo.lxyz[2] = lMCPosTrack.vz(); + thisInfo.pdgCodeV0 = lNegMother.pdgCode(); // 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 @@ -112,19 +160,24 @@ struct cascademcbuilder { 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(); + thisInfo.label = lV0Mother.globalIndex(); + + if(lV0Mother.has_mcCollision()){ + thisInfo.mcCollision = lV0Mother.mcCollisionId(); // save this reference, please + } + + thisInfo.pdgCode = lV0Mother.pdgCode(); + thisInfo.isPhysicalPrimary = lV0Mother.isPhysicalPrimary(); + thisInfo.xyz[0] = lMCBachTrack.vx(); + thisInfo.xyz[1] = lMCBachTrack.vy(); + thisInfo.xyz[2] = lMCBachTrack.vz(); + thisInfo.momentum[0] = lV0Mother.px(); + thisInfo.momentum[1] = lV0Mother.py(); + thisInfo.momentum[2] = lV0Mother.pz(); if (lV0Mother.has_mothers()) { for (auto& lV0GrandMother : lV0Mother.template mothers_as()) { - pdgCodeMother = lV0GrandMother.pdgCode(); - lMotherLabel = lV0GrandMother.globalIndex(); + thisInfo.pdgCodeMother = lV0GrandMother.pdgCode(); + thisInfo.motherLabel = lV0GrandMother.globalIndex(); } } } @@ -138,32 +191,122 @@ struct cascademcbuilder { } // end association check // Construct label table (note: this will be joinable with CascDatas) casclabels( - lLabel, lMotherLabel); - if (populateCascMCCores) { + thisInfo.label, thisInfo.motherLabel); + + // Mark mcParticle as recoed (no searching necessary afterwards) + if(thisInfo.motherLabel>-1){ + mcParticleIsReco[thisInfo.motherLabel] = true; + } + + if (populateCascMCCoresSymmetric) { cascmccores( - pdgCode, pdgCodeMother, pdgCodeV0, isPhysicalPrimary, - pdgCodePositive, pdgCodeNegative, pdgCodeBachelor, - xmc, ymc, zmc, xlmc, ylmc, zlmc, - pxposmc, pyposmc, pzposmc, - pxnegmc, pynegmc, pznegmc, - pxbachmc, pybachmc, pzbachmc, - px, py, pz); + thisInfo.pdgCode, thisInfo.pdgCodeMother, thisInfo.pdgCodeV0, thisInfo.isPhysicalPrimary, + thisInfo.pdgCodePositive, thisInfo.pdgCodeNegative, thisInfo.pdgCodeBachelor, + thisInfo.xyz[0], thisInfo.xyz[1], thisInfo.xyz[2], + thisInfo.lxyz[0], thisInfo.lxyz[1], thisInfo.lxyz[2], + thisInfo.posP[0], thisInfo.posP[1], thisInfo.posP[2], + thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2], + thisInfo.bachP[0], thisInfo.bachP[1], thisInfo.bachP[2], + thisInfo.momentum[0], thisInfo.momentum[1], thisInfo.momentum[2]); + cascmccollrefs(thisInfo.mcCollision); + } + + if (populateCascMCCoresAsymmetric) { + int thisCascMCCoreIndex = -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 < mcCascinfos.size(); ii++) { + if ( + thisInfo.mcParticlePositive == mcCascinfos[ii].mcParticlePositive && mcCascinfos[ii].mcParticlePositive > 0 && + thisInfo.mcParticleNegative == mcCascinfos[ii].mcParticleNegative && mcCascinfos[ii].mcParticleNegative > 0 && + thisInfo.mcParticleBachelor == mcCascinfos[ii].mcParticleBachelor && mcCascinfos[ii].mcParticleBachelor > 0 + ) { + thisCascMCCoreIndex = ii; + break; // this exists already in list + } + } + if (thisCascMCCoreIndex < 0) { + // this CascMCCore does not exist yet. Create it and reference it + thisCascMCCoreIndex = mcCascinfos.size(); + mcCascinfos.push_back(thisInfo); + } + cascCoreMClabels(thisCascMCCoreIndex); // interlink: reconstructed -> MC index } } // end casctable loop + + // now populate V0MCCores if in asymmetric mode + if (populateCascMCCoresAsymmetric) { + // first step: add any un-recoed v0mmcores that were requested + for (auto & mcParticle : mcParticles) { + thisInfo.pdgCode = -1, thisInfo.pdgCodeMother = -1; + thisInfo.pdgCodePositive = -1, thisInfo.pdgCodeNegative = -1; + thisInfo.pdgCodeBachelor = -1, thisInfo.pdgCodeV0 = -1; + thisInfo.isPhysicalPrimary = false; + thisInfo.xyz[0] = 0.0f, thisInfo.xyz[1] = 0.0f, thisInfo.xyz[2] = 0.0f; + thisInfo.lxyz[0] = 0.0f, thisInfo.lxyz[1] = 0.0f, thisInfo.lxyz[2] = 0.0f; + thisInfo.posP[0] = 0.0f, thisInfo.posP[1] = 0.0f, thisInfo.posP[2] = 0.0f; + thisInfo.negP[0] = 0.0f, thisInfo.negP[1] = 0.0f, thisInfo.negP[2] = 0.0f; + thisInfo.bachP[0] = 0.0f, thisInfo.bachP[1] = 0.0f, thisInfo.bachP[2] = 0.0f; + thisInfo.momentum[0] = 0.0f, thisInfo.momentum[1] = 0.0f, thisInfo.momentum[2] = 0.0f; + thisInfo.label = -1, thisInfo.motherLabel = -1; + thisInfo.mcParticlePositive=-1; + thisInfo.mcParticleNegative=-1; + thisInfo.mcParticleBachelor=-1; + + if(mcParticleIsReco[mcParticle.globalIndex()] == true) + continue; // skip if already created in list + + if(TMath::Abs(mcParticle.y()) > rapidityWindow) + continue; // skip outside midrapidity + + if( + (addGeneratedXiMinus && mcParticle.pdgCode() == 3312) || + (addGeneratedXiPlus && mcParticle.pdgCode() == -3312) || + (addGeneratedOmegaMinus && mcParticle.pdgCode() == 3334) || + (addGeneratedOmegaPlus && mcParticle.pdgCode() == -3334) + ){ + thisInfo.pdgCode = mcParticle.pdgCode(); + thisInfo.isPhysicalPrimary = mcParticle.isPhysicalPrimary(); + + if(mcParticle.has_mcCollision()){ + thisInfo.mcCollision = mcParticle.mcCollisionId(); // save this reference, please + } + thisInfo.momentum[0] = mcParticle.px(); + thisInfo.momentum[1] = mcParticle.py(); + thisInfo.momentum[2] = mcParticle.pz(); + + // if I got here, it means this MC particle was not recoed and is of interest. Add it please + mcCascinfos.push_back(thisInfo); + } + } + + for (auto thisInfo : mcCascinfos) { + cascmccores( // a lot of the info below will be compressed in case of not-recoed MC (good!) + thisInfo.pdgCode, thisInfo.pdgCodeMother, thisInfo.pdgCodeV0, thisInfo.isPhysicalPrimary, + thisInfo.pdgCodePositive, thisInfo.pdgCodeNegative, thisInfo.pdgCodeBachelor, + thisInfo.xyz[0], thisInfo.xyz[1], thisInfo.xyz[2], + thisInfo.lxyz[0], thisInfo.lxyz[1], thisInfo.lxyz[2], + thisInfo.posP[0], thisInfo.posP[1], thisInfo.posP[2], + thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2], + thisInfo.bachP[0], thisInfo.bachP[1], thisInfo.bachP[2], + thisInfo.momentum[0], thisInfo.momentum[1], thisInfo.momentum[2]); + cascmccollrefs(thisInfo.mcCollision); + } + } } //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* // build cascade labels - void processCascades(aod::CascDatas const& casctable, aod::V0sLinked const&, aod::V0Datas const& /*v0table*/, aod::McTrackLabels const&, aod::McParticles const&) + void processCascades(aod::CascDatas const& casctable, aod::V0sLinked const&, aod::V0Datas const& /*v0table*/, aod::McTrackLabels const&, aod::McParticles const& mcParticles) { - generateCascadeMCinfo(casctable); + generateCascadeMCinfo(casctable, mcParticles); } //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* // build findable cascade labels - void processFindableCascades(aod::CascDatas const& casctable, aod::FindableV0sLinked const&, aod::V0Datas const& /*v0table*/, aod::McTrackLabels const&, aod::McParticles const&) + void processFindableCascades(aod::CascDatas const& casctable, aod::FindableV0sLinked const&, aod::V0Datas const& /*v0table*/, aod::McTrackLabels const&, aod::McParticles const& mcParticles) { - generateCascadeMCinfo(casctable); + generateCascadeMCinfo(casctable, mcParticles); } //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* diff --git a/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx b/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx index 82bda1aea0b..728e5a958c5 100644 --- a/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx @@ -45,10 +45,18 @@ 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 + Produces v0mccollref; // references collisions from V0MCCores 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."}; + Configurable addGeneratedK0Short{"addGeneratedK0Short", false, "add V0MCCore entry for generated, not-recoed K0Short"}; + Configurable addGeneratedLambda{"addGeneratedLambda", false, "add V0MCCore entry for generated, not-recoed Lambda"}; + Configurable addGeneratedAntiLambda{"addGeneratedAntiLambda", false, "add V0MCCore entry for generated, not-recoed AntiLambda"}; + Configurable addGeneratedGamma{"addGeneratedGamma", false, "add V0MCCore entry for generated, not-recoed Gamma"}; + + Configurable rapidityWindow{"rapidityWindow", 0.5, "rapidity window to save non-recoed candidates"}; + HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; void init(InitContext const&) @@ -80,6 +88,7 @@ struct lambdakzeromcbuilder { int pdgCodeMother; int pdgCodePositive; int pdgCodeNegative; + int mcCollision; bool isPhysicalPrimary; std::array xyz; std::array posP; @@ -97,10 +106,11 @@ struct lambdakzeromcbuilder { //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* // build V0 labels - void process(aod::V0Datas const& v0table, aod::McTrackLabels const&, aod::McParticles const& /*particlesMC*/) + void process(aod::V0Datas const& v0table, aod::McTrackLabels const&, aod::McParticles const& mcParticles) { // to be used if using the populateV0MCCoresAsymmetric mode, kept empty otherwise std::vector mcV0infos; // V0MCCore information + std::vector mcParticleIsReco(mcParticles.size(), false); // mc Particle not recoed by V0s for (auto& v0 : v0table) { thisInfo.packedMcParticleIndices = 0; // not de-referenced properly yet @@ -110,6 +120,10 @@ struct lambdakzeromcbuilder { thisInfo.pdgCodeMother = 0; thisInfo.pdgCodePositive = 0; thisInfo.pdgCodeNegative = 0; + thisInfo.mcCollision = -1; + thisInfo.xyz[0] = thisInfo.xyz[1] = thisInfo.xyz[2] = 0.0f; + thisInfo.posP[0] = thisInfo.posP[1] = thisInfo.posP[2] = 0.0f; + thisInfo.negP[0] = thisInfo.negP[1] = thisInfo.negP[2] = 0.0f; auto lNegTrack = v0.negTrack_as(); auto lPosTrack = v0.posTrack_as(); @@ -132,6 +146,11 @@ struct lambdakzeromcbuilder { for (auto& lPosMother : lMCPosTrack.mothers_as()) { if (lNegMother.globalIndex() == lPosMother.globalIndex()) { thisInfo.label = lNegMother.globalIndex(); + + if(lNegMother.has_mcCollision()){ + thisInfo.mcCollision = lNegMother.mcCollisionId(); // save this reference, please + } + // acquire information thisInfo.xyz[0] = lMCPosTrack.vx(); thisInfo.xyz[1] = lMCPosTrack.vy(); @@ -153,6 +172,11 @@ struct lambdakzeromcbuilder { v0labels( thisInfo.label, thisInfo.motherLabel); + // Mark mcParticle as recoed (no searching necessary afterwards) + if(thisInfo.motherLabel>-1){ + mcParticleIsReco[thisInfo.motherLabel] = true; + } + // ---] Symmetric populate [--- // in this approach, V0Cores will be joinable with V0MCCores. // this is the most pedagogical approach, but it is also more limited @@ -164,6 +188,7 @@ struct lambdakzeromcbuilder { 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]); + v0mccollref(thisInfo.mcCollision); // n.b. placing the interlink index here allows for the writing of // code that is agnostic with respect to the joinability of @@ -197,6 +222,50 @@ struct lambdakzeromcbuilder { // now populate V0MCCores if in asymmetric mode if (populateV0MCCoresAsymmetric) { + // first step: add any un-recoed v0mmcores that were requested + for (auto & mcParticle : mcParticles) { + thisInfo.packedMcParticleIndices = 0; + thisInfo.label = -1; + thisInfo.motherLabel = -1; + thisInfo.pdgCode = 0; + thisInfo.pdgCodeMother = 0; + thisInfo.pdgCodePositive = 0; + thisInfo.pdgCodeNegative = 0; + thisInfo.mcCollision = -1; + thisInfo.xyz[0] = thisInfo.xyz[1] = thisInfo.xyz[2] = 0.0f; + thisInfo.posP[0] = thisInfo.posP[1] = thisInfo.posP[2] = 0.0f; + thisInfo.negP[0] = thisInfo.negP[1] = thisInfo.negP[2] = 0.0f; + + if(mcParticleIsReco[mcParticle.globalIndex()] == true) + continue; // skip if already created in list + + if(TMath::Abs(mcParticle.y()) > rapidityWindow) + continue; // skip outside midrapidity + + if( + (addGeneratedK0Short && mcParticle.pdgCode() == 310) || + (addGeneratedLambda && mcParticle.pdgCode() == 3122) || + (addGeneratedAntiLambda && mcParticle.pdgCode() == -3122) || + (addGeneratedGamma && mcParticle.pdgCode() == 22) + ){ + thisInfo.pdgCode = mcParticle.pdgCode(); + thisInfo.isPhysicalPrimary = mcParticle.isPhysicalPrimary(); + + if(mcParticle.has_mcCollision()){ + thisInfo.mcCollision = mcParticle.mcCollisionId(); // save this reference, please + } + + // guarantee compressibility: keep momentum entirely in the positive prong + // WARNING: THIS IS AS ARBITRARY AS IT GETS but should be ok + thisInfo.posP[0] = mcParticle.px(); + thisInfo.posP[1] = mcParticle.py(); + thisInfo.posP[2] = mcParticle.pz(); + + // if I got here, it means this MC particle was not recoed and is of interest. Add it please + mcV0infos.push_back(thisInfo); + } + } + for (auto info : mcV0infos) { v0mccores( info.label, info.pdgCode, @@ -204,6 +273,7 @@ struct lambdakzeromcbuilder { 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]); + v0mccollref(thisInfo.mcCollision); } } diff --git a/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx b/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx index 4db5528e1b8..351eac60497 100644 --- a/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx @@ -77,9 +77,7 @@ struct strangederivedbuilder { Produces strangeEvSels; // characterises collisions / sel8 selection Produces strangeStamps; // provides timestamps, run numbers Produces v0collref; // references collisions from V0s - Produces v0mccollref; // references collisions from V0s Produces casccollref; // references collisions from cascades - Produces cascmccollref; // references collisions from V0s Produces kfcasccollref; // references collisions from KF cascades Produces tracasccollref; // references collisions from tracked cascades @@ -366,7 +364,7 @@ struct strangederivedbuilder { tracasccollref(TraCascadeCollIndices[casc.globalIndex()]); } - void processCollisionsMC(soa::Join const& collisions, soa::Join const& V0s, soa::Join const& Cascades, aod::KFCascDatas const& KFCascades, aod::TraCascDatas const& TraCascades, aod::BCsWithTimestamps const&, soa::Join const& mcCollisions, aod::McParticles const&) + void processCollisionsMC(soa::Join const& collisions, soa::Join const& V0s, soa::Join const& V0MCCores, soa::Join const& Cascades, aod::KFCascDatas const& KFCascades, aod::TraCascDatas const& TraCascades, aod::BCsWithTimestamps const&, soa::Join const& mcCollisions, aod::McParticles const&) { // create collision indices beforehand std::vector V0CollIndices(V0s.size(), -1); // index -1: no collision @@ -441,38 +439,14 @@ struct strangederivedbuilder { KFCascadeCollIndices[casc.globalIndex()] = strangeColl.lastIndex(); for (const auto& casc : TraCascTable_thisColl) TraCascadeCollIndices[casc.globalIndex()] = strangeColl.lastIndex(); - - // populate MC collision references - for (const auto& v0 : V0Table_thisColl) { - uint32_t indMCColl = -1; - if (v0.has_mcParticle()) { - auto mcParticle = v0.mcParticle(); - if (mcParticle.has_mcCollision()) { - indMCColl = mcParticle.mcCollisionId(); - } - } - V0MCCollIndices[v0.globalIndex()] = indMCColl; - } - for (const auto& casc : CascTable_thisColl) { - uint32_t indMCColl = -1; - if (casc.has_mcParticle()) { - auto mcParticle = casc.mcParticle(); - if (mcParticle.has_mcCollision()) { - indMCColl = mcParticle.mcCollisionId(); - } - } - CascadeMCCollIndices[casc.globalIndex()] = indMCColl; - } } // populate references, including those that might not be assigned for (const auto& v0 : V0s) { v0collref(V0CollIndices[v0.globalIndex()]); - v0mccollref(V0MCCollIndices[v0.globalIndex()]); } for (const auto& casc : Cascades) { casccollref(CascadeCollIndices[casc.globalIndex()]); - cascmccollref(CascadeMCCollIndices[casc.globalIndex()]); } for (const auto& casc : KFCascades) kfcasccollref(KFCascadeCollIndices[casc.globalIndex()]); From 8db1480e76c2fe1229e338c232cb3ce37338a3d5 Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Tue, 18 Jun 2024 15:06:56 +0200 Subject: [PATCH 2/2] Please consider the following formatting changes (#308) --- .../Strangeness/cascademcbuilder.cxx | 62 +++++++++---------- .../Strangeness/lambdakzeromcbuilder.cxx | 31 +++++----- 2 files changed, 45 insertions(+), 48 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx b/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx index e57a9959a82..4084afdcae4 100644 --- a/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx @@ -96,7 +96,7 @@ struct cascademcbuilder { { // to be used if using the asymmetric mode, kept empty otherwise - std::vector mcCascinfos; // V0MCCore information + std::vector mcCascinfos; // V0MCCore information std::vector mcParticleIsReco(mcParticles.size(), false); // mc Particle not recoed by V0s for (auto& casc : cascTable) { @@ -111,9 +111,9 @@ struct cascademcbuilder { thisInfo.bachP[0] = -999.0f, thisInfo.bachP[1] = -999.0f, thisInfo.bachP[2] = -999.0f; thisInfo.momentum[0] = -999.0f, thisInfo.momentum[1] = -999.0f, thisInfo.momentum[2] = -999.0f; thisInfo.label = -1, thisInfo.motherLabel = -1; - thisInfo.mcParticlePositive=-1; - thisInfo.mcParticleNegative=-1; - thisInfo.mcParticleBachelor=-1; + thisInfo.mcParticlePositive = -1; + thisInfo.mcParticleNegative = -1; + thisInfo.mcParticleBachelor = -1; // Acquire all three daughter tracks, please auto lBachTrack = casc.template bachelor_as(); @@ -162,7 +162,7 @@ struct cascademcbuilder { if (lV0Mother == lBachMother) { thisInfo.label = lV0Mother.globalIndex(); - if(lV0Mother.has_mcCollision()){ + if (lV0Mother.has_mcCollision()) { thisInfo.mcCollision = lV0Mother.mcCollisionId(); // save this reference, please } @@ -194,19 +194,19 @@ struct cascademcbuilder { thisInfo.label, thisInfo.motherLabel); // Mark mcParticle as recoed (no searching necessary afterwards) - if(thisInfo.motherLabel>-1){ - mcParticleIsReco[thisInfo.motherLabel] = true; + if (thisInfo.motherLabel > -1) { + mcParticleIsReco[thisInfo.motherLabel] = true; } if (populateCascMCCoresSymmetric) { cascmccores( thisInfo.pdgCode, thisInfo.pdgCodeMother, thisInfo.pdgCodeV0, thisInfo.isPhysicalPrimary, thisInfo.pdgCodePositive, thisInfo.pdgCodeNegative, thisInfo.pdgCodeBachelor, - thisInfo.xyz[0], thisInfo.xyz[1], thisInfo.xyz[2], + thisInfo.xyz[0], thisInfo.xyz[1], thisInfo.xyz[2], thisInfo.lxyz[0], thisInfo.lxyz[1], thisInfo.lxyz[2], - thisInfo.posP[0], thisInfo.posP[1], thisInfo.posP[2], - thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2], - thisInfo.bachP[0], thisInfo.bachP[1], thisInfo.bachP[2], + thisInfo.posP[0], thisInfo.posP[1], thisInfo.posP[2], + thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2], + thisInfo.bachP[0], thisInfo.bachP[1], thisInfo.bachP[2], thisInfo.momentum[0], thisInfo.momentum[1], thisInfo.momentum[2]); cascmccollrefs(thisInfo.mcCollision); } @@ -219,10 +219,9 @@ struct cascademcbuilder { if ( thisInfo.mcParticlePositive == mcCascinfos[ii].mcParticlePositive && mcCascinfos[ii].mcParticlePositive > 0 && thisInfo.mcParticleNegative == mcCascinfos[ii].mcParticleNegative && mcCascinfos[ii].mcParticleNegative > 0 && - thisInfo.mcParticleBachelor == mcCascinfos[ii].mcParticleBachelor && mcCascinfos[ii].mcParticleBachelor > 0 - ) { + thisInfo.mcParticleBachelor == mcCascinfos[ii].mcParticleBachelor && mcCascinfos[ii].mcParticleBachelor > 0) { thisCascMCCoreIndex = ii; - break; // this exists already in list + break; // this exists already in list } } if (thisCascMCCoreIndex < 0) { @@ -236,8 +235,8 @@ struct cascademcbuilder { // now populate V0MCCores if in asymmetric mode if (populateCascMCCoresAsymmetric) { - // first step: add any un-recoed v0mmcores that were requested - for (auto & mcParticle : mcParticles) { + // first step: add any un-recoed v0mmcores that were requested + for (auto& mcParticle : mcParticles) { thisInfo.pdgCode = -1, thisInfo.pdgCodeMother = -1; thisInfo.pdgCodePositive = -1, thisInfo.pdgCodeNegative = -1; thisInfo.pdgCodeBachelor = -1, thisInfo.pdgCodeV0 = -1; @@ -249,32 +248,31 @@ struct cascademcbuilder { thisInfo.bachP[0] = 0.0f, thisInfo.bachP[1] = 0.0f, thisInfo.bachP[2] = 0.0f; thisInfo.momentum[0] = 0.0f, thisInfo.momentum[1] = 0.0f, thisInfo.momentum[2] = 0.0f; thisInfo.label = -1, thisInfo.motherLabel = -1; - thisInfo.mcParticlePositive=-1; - thisInfo.mcParticleNegative=-1; - thisInfo.mcParticleBachelor=-1; + thisInfo.mcParticlePositive = -1; + thisInfo.mcParticleNegative = -1; + thisInfo.mcParticleBachelor = -1; - if(mcParticleIsReco[mcParticle.globalIndex()] == true) + if (mcParticleIsReco[mcParticle.globalIndex()] == true) continue; // skip if already created in list - - if(TMath::Abs(mcParticle.y()) > rapidityWindow) + + if (TMath::Abs(mcParticle.y()) > rapidityWindow) continue; // skip outside midrapidity - if( - (addGeneratedXiMinus && mcParticle.pdgCode() == 3312) || + if ( + (addGeneratedXiMinus && mcParticle.pdgCode() == 3312) || (addGeneratedXiPlus && mcParticle.pdgCode() == -3312) || (addGeneratedOmegaMinus && mcParticle.pdgCode() == 3334) || - (addGeneratedOmegaPlus && mcParticle.pdgCode() == -3334) - ){ + (addGeneratedOmegaPlus && mcParticle.pdgCode() == -3334)) { thisInfo.pdgCode = mcParticle.pdgCode(); thisInfo.isPhysicalPrimary = mcParticle.isPhysicalPrimary(); - if(mcParticle.has_mcCollision()){ + if (mcParticle.has_mcCollision()) { thisInfo.mcCollision = mcParticle.mcCollisionId(); // save this reference, please } thisInfo.momentum[0] = mcParticle.px(); thisInfo.momentum[1] = mcParticle.py(); thisInfo.momentum[2] = mcParticle.pz(); - + // if I got here, it means this MC particle was not recoed and is of interest. Add it please mcCascinfos.push_back(thisInfo); } @@ -284,11 +282,11 @@ struct cascademcbuilder { cascmccores( // a lot of the info below will be compressed in case of not-recoed MC (good!) thisInfo.pdgCode, thisInfo.pdgCodeMother, thisInfo.pdgCodeV0, thisInfo.isPhysicalPrimary, thisInfo.pdgCodePositive, thisInfo.pdgCodeNegative, thisInfo.pdgCodeBachelor, - thisInfo.xyz[0], thisInfo.xyz[1], thisInfo.xyz[2], + thisInfo.xyz[0], thisInfo.xyz[1], thisInfo.xyz[2], thisInfo.lxyz[0], thisInfo.lxyz[1], thisInfo.lxyz[2], - thisInfo.posP[0], thisInfo.posP[1], thisInfo.posP[2], - thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2], - thisInfo.bachP[0], thisInfo.bachP[1], thisInfo.bachP[2], + thisInfo.posP[0], thisInfo.posP[1], thisInfo.posP[2], + thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2], + thisInfo.bachP[0], thisInfo.bachP[1], thisInfo.bachP[2], thisInfo.momentum[0], thisInfo.momentum[1], thisInfo.momentum[2]); cascmccollrefs(thisInfo.mcCollision); } diff --git a/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx b/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx index 728e5a958c5..8fe57bc144f 100644 --- a/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx @@ -45,7 +45,7 @@ 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 - Produces v0mccollref; // references collisions from V0MCCores + Produces v0mccollref; // references collisions from V0MCCores 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."}; @@ -147,7 +147,7 @@ struct lambdakzeromcbuilder { if (lNegMother.globalIndex() == lPosMother.globalIndex()) { thisInfo.label = lNegMother.globalIndex(); - if(lNegMother.has_mcCollision()){ + if (lNegMother.has_mcCollision()) { thisInfo.mcCollision = lNegMother.mcCollisionId(); // save this reference, please } @@ -173,8 +173,8 @@ struct lambdakzeromcbuilder { thisInfo.label, thisInfo.motherLabel); // Mark mcParticle as recoed (no searching necessary afterwards) - if(thisInfo.motherLabel>-1){ - mcParticleIsReco[thisInfo.motherLabel] = true; + if (thisInfo.motherLabel > -1) { + mcParticleIsReco[thisInfo.motherLabel] = true; } // ---] Symmetric populate [--- @@ -222,8 +222,8 @@ struct lambdakzeromcbuilder { // now populate V0MCCores if in asymmetric mode if (populateV0MCCoresAsymmetric) { - // first step: add any un-recoed v0mmcores that were requested - for (auto & mcParticle : mcParticles) { + // first step: add any un-recoed v0mmcores that were requested + for (auto& mcParticle : mcParticles) { thisInfo.packedMcParticleIndices = 0; thisInfo.label = -1; thisInfo.motherLabel = -1; @@ -236,31 +236,30 @@ struct lambdakzeromcbuilder { thisInfo.posP[0] = thisInfo.posP[1] = thisInfo.posP[2] = 0.0f; thisInfo.negP[0] = thisInfo.negP[1] = thisInfo.negP[2] = 0.0f; - if(mcParticleIsReco[mcParticle.globalIndex()] == true) + if (mcParticleIsReco[mcParticle.globalIndex()] == true) continue; // skip if already created in list - - if(TMath::Abs(mcParticle.y()) > rapidityWindow) + + if (TMath::Abs(mcParticle.y()) > rapidityWindow) continue; // skip outside midrapidity - if( - (addGeneratedK0Short && mcParticle.pdgCode() == 310) || + if ( + (addGeneratedK0Short && mcParticle.pdgCode() == 310) || (addGeneratedLambda && mcParticle.pdgCode() == 3122) || (addGeneratedAntiLambda && mcParticle.pdgCode() == -3122) || - (addGeneratedGamma && mcParticle.pdgCode() == 22) - ){ + (addGeneratedGamma && mcParticle.pdgCode() == 22)) { thisInfo.pdgCode = mcParticle.pdgCode(); thisInfo.isPhysicalPrimary = mcParticle.isPhysicalPrimary(); - if(mcParticle.has_mcCollision()){ + if (mcParticle.has_mcCollision()) { thisInfo.mcCollision = mcParticle.mcCollisionId(); // save this reference, please } - // guarantee compressibility: keep momentum entirely in the positive prong + // guarantee compressibility: keep momentum entirely in the positive prong // WARNING: THIS IS AS ARBITRARY AS IT GETS but should be ok thisInfo.posP[0] = mcParticle.px(); thisInfo.posP[1] = mcParticle.py(); thisInfo.posP[2] = mcParticle.pz(); - + // if I got here, it means this MC particle was not recoed and is of interest. Add it please mcV0infos.push_back(thisInfo); }