Skip to content

PWGLF: mc core for all gen #6565

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
njacazio merged 2 commits into from
Jun 19, 2024
Merged

PWGLF: mc core for all gen #6565

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b07803e
PWGLF: mc core for all gen
ddobrigk Jun 17, 2024
8db1480
Please consider the following formatting changes (#308)
alibuild Jun 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions PWGLF/DataModel/LFStrangenessTables.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>);

Expand Down
247 changes: 194 additions & 53 deletions PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,25 +50,70 @@ struct cascademcbuilder {
Produces<aod::McTraCascLabels> tracasclabels; // MC labels for tracked cascades
Produces<aod::McCascBBTags> bbtags; // bb tags (inv structure tagging)
Produces<aod::CascMCCores> cascmccores; // optionally aggregate information from MC side for posterior analysis (derived data)
Produces<aod::CascCoreMCLabels> cascCoreMClabels; // optionally aggregate information from MC side for posterior analysis (derived data)
Produces<aod::CascMCCollRefs> cascmccollrefs; // references MC collisions from MC cascades

Configurable<bool> populateCascMCCores{"populateCascMCCores", false, "populate CascMCCores table for derived data analysis"};
Configurable<bool> populateCascMCCoresSymmetric{"populateCascMCCoresSymmetric", false, "populate CascMCCores table for derived data analysis, keep CascMCCores joinable with CascCores"};
Configurable<bool> populateCascMCCoresAsymmetric{"populateCascMCCoresAsymmetric", false, "populate CascMCCores table for derived data analysis, create CascCores -> CascMCCores interlink. Saves only labeled Cascades."};

Configurable<bool> addGeneratedXiMinus{"addGeneratedXiMinus", false, "add CascMCCore entry for generated, not-recoed XiMinus"};
Configurable<bool> addGeneratedXiPlus{"addGeneratedXiPlus", false, "add CascMCCore entry for generated, not-recoed XiPlus"};
Configurable<bool> addGeneratedOmegaMinus{"addGeneratedOmegaMinus", false, "add CascMCCore entry for generated, not-recoed OmegaMinus"};
Configurable<bool> addGeneratedOmegaPlus{"addGeneratedOmegaPlus", false, "add CascMCCore entry for generated, not-recoed OmegaPlus"};

Configurable<float> 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<float, 3> xyz;
std::array<float, 3> lxyz;
std::array<float, 3> posP;
std::array<float, 3> negP;
std::array<float, 3> bachP;
std::array<float, 3> momentum;
int mcParticlePositive;
int mcParticleNegative;
int mcParticleBachelor;
};
mcCascinfo thisInfo;
//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*

void init(InitContext const&) {}

template <typename TCascadeTable>
void generateCascadeMCinfo(TCascadeTable cascTable)
template <typename TCascadeTable, typename TMCParticleTable>
void generateCascadeMCinfo(TCascadeTable cascTable, TMCParticleTable mcParticles)
{

// to be used if using the asymmetric mode, kept empty otherwise
std::vector<mcCascinfo> mcCascinfos; // V0MCCore information
std::vector<bool> 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<aod::McTrackLabels>();
Expand All @@ -82,49 +127,57 @@ struct cascademcbuilder {
auto lMCNegTrack = lNegTrack.template mcParticle_as<aod::McParticles>();
auto lMCPosTrack = lPosTrack.template mcParticle_as<aod::McParticles>();

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()) {
for (auto& lNegMother : lMCNegTrack.template mothers_as<aod::McParticles>()) {
for (auto& lPosMother : lMCPosTrack.template mothers_as<aod::McParticles>()) {
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
if (lNegMother.has_mothers() && lMCBachTrack.has_mothers()) {
for (auto& lV0Mother : lNegMother.template mothers_as<aod::McParticles>()) {
for (auto& lBachMother : lMCBachTrack.template mothers_as<aod::McParticles>()) {
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<aod::McParticles>()) {
pdgCodeMother = lV0GrandMother.pdgCode();
lMotherLabel = lV0GrandMother.globalIndex();
thisInfo.pdgCodeMother = lV0GrandMother.pdgCode();
thisInfo.motherLabel = lV0GrandMother.globalIndex();
}
}
}
Expand All @@ -138,32 +191,120 @@ 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);
}

//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
Expand Down
Loading