Skip to content

Commit 1674243

Browse files
ddobrigkalibuild
ddobrigk
and
alibuild
authored
PWGLF: mc core for all gen (#6565)
* PWGLF: mc core for all gen * Please consider the following formatting changes (#308) --------- Co-authored-by: ALICE Builder <alibuild@users.noreply.github.com>
1 parent bcefd65 commit 1674243

4 files changed

Lines changed: 272 additions & 81 deletions

File tree

‎PWGLF/DataModel/LFStrangenessTables.h‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1190,6 +1190,13 @@ DECLARE_SOA_TABLE(CascMCCores, "AOD", "CASCMCCORE", //! bachelor-baryon correlat
11901190
cascdata::PxBachMC, cascdata::PyBachMC, cascdata::PzBachMC,
11911191
cascdata::PxMC, cascdata::PyMC, cascdata::PzMC);
11921192

1193+
namespace cascdata
1194+
{
1195+
DECLARE_SOA_INDEX_COLUMN(CascMCCore, cascMCCore); //! Index to CascMCCore entry
1196+
}
1197+
1198+
DECLARE_SOA_TABLE(CascCoreMCLabels, "AOD", "CASCCOREMCLABEL", //! optional table to refer to CascMCCores if not joinable
1199+
o2::soa::Index<>, cascdata::CascMCCoreId);
11931200
DECLARE_SOA_TABLE(CascMCCollRefs, "AOD", "CASCMCCOLLREF", //! refers MC candidate back to proper MC Collision
11941201
o2::soa::Index<>, cascdata::StraMCCollisionId, o2::soa::Marker<3>);
11951202

‎PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx‎

Lines changed: 194 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -50,25 +50,70 @@ struct cascademcbuilder {
5050
Produces<aod::McTraCascLabels> tracasclabels; // MC labels for tracked cascades
5151
Produces<aod::McCascBBTags> bbtags; // bb tags (inv structure tagging)
5252
Produces<aod::CascMCCores> cascmccores; // optionally aggregate information from MC side for posterior analysis (derived data)
53+
Produces<aod::CascCoreMCLabels> cascCoreMClabels; // optionally aggregate information from MC side for posterior analysis (derived data)
54+
Produces<aod::CascMCCollRefs> cascmccollrefs; // references MC collisions from MC cascades
5355

54-
Configurable<bool> populateCascMCCores{"populateCascMCCores", false, "populate CascMCCores table for derived data analysis"};
56+
Configurable<bool> populateCascMCCoresSymmetric{"populateCascMCCoresSymmetric", false, "populate CascMCCores table for derived data analysis, keep CascMCCores joinable with CascCores"};
57+
Configurable<bool> populateCascMCCoresAsymmetric{"populateCascMCCoresAsymmetric", false, "populate CascMCCores table for derived data analysis, create CascCores -> CascMCCores interlink. Saves only labeled Cascades."};
58+
59+
Configurable<bool> addGeneratedXiMinus{"addGeneratedXiMinus", false, "add CascMCCore entry for generated, not-recoed XiMinus"};
60+
Configurable<bool> addGeneratedXiPlus{"addGeneratedXiPlus", false, "add CascMCCore entry for generated, not-recoed XiPlus"};
61+
Configurable<bool> addGeneratedOmegaMinus{"addGeneratedOmegaMinus", false, "add CascMCCore entry for generated, not-recoed OmegaMinus"};
62+
Configurable<bool> addGeneratedOmegaPlus{"addGeneratedOmegaPlus", false, "add CascMCCore entry for generated, not-recoed OmegaPlus"};
63+
64+
Configurable<float> rapidityWindow{"rapidityWindow", 0.5, "rapidity window to save non-recoed candidates"};
65+
66+
//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
67+
// Helper struct to contain CascMCCore information prior to filling
68+
struct mcCascinfo {
69+
int label;
70+
int motherLabel;
71+
int mcCollision;
72+
int pdgCode;
73+
int pdgCodeMother;
74+
int pdgCodeV0;
75+
int pdgCodePositive;
76+
int pdgCodeNegative;
77+
int pdgCodeBachelor;
78+
bool isPhysicalPrimary;
79+
std::array<float, 3> xyz;
80+
std::array<float, 3> lxyz;
81+
std::array<float, 3> posP;
82+
std::array<float, 3> negP;
83+
std::array<float, 3> bachP;
84+
std::array<float, 3> momentum;
85+
int mcParticlePositive;
86+
int mcParticleNegative;
87+
int mcParticleBachelor;
88+
};
89+
mcCascinfo thisInfo;
90+
//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
5591

5692
void init(InitContext const&) {}
5793

58-
template <typename TCascadeTable>
59-
void generateCascadeMCinfo(TCascadeTable cascTable)
94+
template <typename TCascadeTable, typename TMCParticleTable>
95+
void generateCascadeMCinfo(TCascadeTable cascTable, TMCParticleTable mcParticles)
6096
{
97+
98+
// to be used if using the asymmetric mode, kept empty otherwise
99+
std::vector<mcCascinfo> mcCascinfos; // V0MCCore information
100+
std::vector<bool> mcParticleIsReco(mcParticles.size(), false); // mc Particle not recoed by V0s
101+
61102
for (auto& casc : cascTable) {
62-
int pdgCode = -1, pdgCodeMother = -1;
63-
int pdgCodePositive = -1, pdgCodeNegative = -1, pdgCodeBachelor = -1, pdgCodeV0 = -1;
64-
bool isPhysicalPrimary = false;
65-
float xmc = -999.0f, ymc = -999.0f, zmc = -999.0f;
66-
float xlmc = -999.0f, ylmc = -999.0f, zlmc = -999.0f;
67-
float pxposmc = -999.0f, pyposmc = -999.0f, pzposmc = -999.0f;
68-
float pxnegmc = -999.0f, pynegmc = -999.0f, pznegmc = -999.0f;
69-
float pxbachmc = -999.0f, pybachmc = -999.0f, pzbachmc = -999.0f;
70-
float px = -999.0f, py = -999.0f, pz = -999.0f;
71-
int lLabel = -1, lMotherLabel = -1;
103+
thisInfo.pdgCode = -1, thisInfo.pdgCodeMother = -1;
104+
thisInfo.pdgCodePositive = -1, thisInfo.pdgCodeNegative = -1;
105+
thisInfo.pdgCodeBachelor = -1, thisInfo.pdgCodeV0 = -1;
106+
thisInfo.isPhysicalPrimary = false;
107+
thisInfo.xyz[0] = -999.0f, thisInfo.xyz[1] = -999.0f, thisInfo.xyz[2] = -999.0f;
108+
thisInfo.lxyz[0] = -999.0f, thisInfo.lxyz[1] = -999.0f, thisInfo.lxyz[2] = -999.0f;
109+
thisInfo.posP[0] = -999.0f, thisInfo.posP[1] = -999.0f, thisInfo.posP[2] = -999.0f;
110+
thisInfo.negP[0] = -999.0f, thisInfo.negP[1] = -999.0f, thisInfo.negP[2] = -999.0f;
111+
thisInfo.bachP[0] = -999.0f, thisInfo.bachP[1] = -999.0f, thisInfo.bachP[2] = -999.0f;
112+
thisInfo.momentum[0] = -999.0f, thisInfo.momentum[1] = -999.0f, thisInfo.momentum[2] = -999.0f;
113+
thisInfo.label = -1, thisInfo.motherLabel = -1;
114+
thisInfo.mcParticlePositive = -1;
115+
thisInfo.mcParticleNegative = -1;
116+
thisInfo.mcParticleBachelor = -1;
72117

73118
// Acquire all three daughter tracks, please
74119
auto lBachTrack = casc.template bachelor_as<aod::McTrackLabels>();
@@ -82,49 +127,57 @@ struct cascademcbuilder {
82127
auto lMCNegTrack = lNegTrack.template mcParticle_as<aod::McParticles>();
83128
auto lMCPosTrack = lPosTrack.template mcParticle_as<aod::McParticles>();
84129

85-
pdgCodePositive = lMCPosTrack.pdgCode();
86-
pdgCodeNegative = lMCNegTrack.pdgCode();
87-
pdgCodeBachelor = lMCBachTrack.pdgCode();
88-
pxposmc = lMCPosTrack.px();
89-
pyposmc = lMCPosTrack.py();
90-
pzposmc = lMCPosTrack.pz();
91-
pxnegmc = lMCNegTrack.px();
92-
pynegmc = lMCNegTrack.py();
93-
pznegmc = lMCNegTrack.pz();
94-
pxbachmc = lMCBachTrack.px();
95-
pybachmc = lMCBachTrack.py();
96-
pzbachmc = lMCBachTrack.pz();
130+
thisInfo.mcParticlePositive = lMCPosTrack.globalIndex();
131+
thisInfo.mcParticleNegative = lMCNegTrack.globalIndex();
132+
thisInfo.mcParticleBachelor = lMCBachTrack.globalIndex();
133+
thisInfo.pdgCodePositive = lMCPosTrack.pdgCode();
134+
thisInfo.pdgCodeNegative = lMCNegTrack.pdgCode();
135+
thisInfo.pdgCodeBachelor = lMCBachTrack.pdgCode();
136+
thisInfo.posP[0] = lMCPosTrack.px();
137+
thisInfo.posP[1] = lMCPosTrack.py();
138+
thisInfo.posP[2] = lMCPosTrack.pz();
139+
thisInfo.negP[0] = lMCNegTrack.px();
140+
thisInfo.negP[1] = lMCNegTrack.py();
141+
thisInfo.negP[2] = lMCNegTrack.pz();
142+
thisInfo.bachP[0] = lMCBachTrack.px();
143+
thisInfo.bachP[1] = lMCBachTrack.py();
144+
thisInfo.bachP[2] = lMCBachTrack.pz();
97145

98146
// Step 1: check if the mother is the same, go up a level
99147
if (lMCNegTrack.has_mothers() && lMCPosTrack.has_mothers()) {
100148
for (auto& lNegMother : lMCNegTrack.template mothers_as<aod::McParticles>()) {
101149
for (auto& lPosMother : lMCPosTrack.template mothers_as<aod::McParticles>()) {
102150
if (lNegMother == lPosMother) {
103151
// acquire information
104-
xlmc = lMCPosTrack.vx();
105-
ylmc = lMCPosTrack.vy();
106-
zlmc = lMCPosTrack.vz();
107-
pdgCodeV0 = lNegMother.pdgCode();
152+
thisInfo.lxyz[0] = lMCPosTrack.vx();
153+
thisInfo.lxyz[1] = lMCPosTrack.vy();
154+
thisInfo.lxyz[2] = lMCPosTrack.vz();
155+
thisInfo.pdgCodeV0 = lNegMother.pdgCode();
108156

109157
// if we got to this level, it means the mother particle exists and is the same
110158
// now we have to go one level up and compare to the bachelor mother too
111159
if (lNegMother.has_mothers() && lMCBachTrack.has_mothers()) {
112160
for (auto& lV0Mother : lNegMother.template mothers_as<aod::McParticles>()) {
113161
for (auto& lBachMother : lMCBachTrack.template mothers_as<aod::McParticles>()) {
114162
if (lV0Mother == lBachMother) {
115-
lLabel = lV0Mother.globalIndex();
116-
pdgCode = lV0Mother.pdgCode();
117-
isPhysicalPrimary = lV0Mother.isPhysicalPrimary();
118-
xmc = lMCBachTrack.vx();
119-
ymc = lMCBachTrack.vy();
120-
zmc = lMCBachTrack.vz();
121-
px = lV0Mother.px();
122-
py = lV0Mother.py();
123-
pz = lV0Mother.pz();
163+
thisInfo.label = lV0Mother.globalIndex();
164+
165+
if (lV0Mother.has_mcCollision()) {
166+
thisInfo.mcCollision = lV0Mother.mcCollisionId(); // save this reference, please
167+
}
168+
169+
thisInfo.pdgCode = lV0Mother.pdgCode();
170+
thisInfo.isPhysicalPrimary = lV0Mother.isPhysicalPrimary();
171+
thisInfo.xyz[0] = lMCBachTrack.vx();
172+
thisInfo.xyz[1] = lMCBachTrack.vy();
173+
thisInfo.xyz[2] = lMCBachTrack.vz();
174+
thisInfo.momentum[0] = lV0Mother.px();
175+
thisInfo.momentum[1] = lV0Mother.py();
176+
thisInfo.momentum[2] = lV0Mother.pz();
124177
if (lV0Mother.has_mothers()) {
125178
for (auto& lV0GrandMother : lV0Mother.template mothers_as<aod::McParticles>()) {
126-
pdgCodeMother = lV0GrandMother.pdgCode();
127-
lMotherLabel = lV0GrandMother.globalIndex();
179+
thisInfo.pdgCodeMother = lV0GrandMother.pdgCode();
180+
thisInfo.motherLabel = lV0GrandMother.globalIndex();
128181
}
129182
}
130183
}
@@ -138,32 +191,120 @@ struct cascademcbuilder {
138191
} // end association check
139192
// Construct label table (note: this will be joinable with CascDatas)
140193
casclabels(
141-
lLabel, lMotherLabel);
142-
if (populateCascMCCores) {
194+
thisInfo.label, thisInfo.motherLabel);
195+
196+
// Mark mcParticle as recoed (no searching necessary afterwards)
197+
if (thisInfo.motherLabel > -1) {
198+
mcParticleIsReco[thisInfo.motherLabel] = true;
199+
}
200+
201+
if (populateCascMCCoresSymmetric) {
143202
cascmccores(
144-
pdgCode, pdgCodeMother, pdgCodeV0, isPhysicalPrimary,
145-
pdgCodePositive, pdgCodeNegative, pdgCodeBachelor,
146-
xmc, ymc, zmc, xlmc, ylmc, zlmc,
147-
pxposmc, pyposmc, pzposmc,
148-
pxnegmc, pynegmc, pznegmc,
149-
pxbachmc, pybachmc, pzbachmc,
150-
px, py, pz);
203+
thisInfo.pdgCode, thisInfo.pdgCodeMother, thisInfo.pdgCodeV0, thisInfo.isPhysicalPrimary,
204+
thisInfo.pdgCodePositive, thisInfo.pdgCodeNegative, thisInfo.pdgCodeBachelor,
205+
thisInfo.xyz[0], thisInfo.xyz[1], thisInfo.xyz[2],
206+
thisInfo.lxyz[0], thisInfo.lxyz[1], thisInfo.lxyz[2],
207+
thisInfo.posP[0], thisInfo.posP[1], thisInfo.posP[2],
208+
thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2],
209+
thisInfo.bachP[0], thisInfo.bachP[1], thisInfo.bachP[2],
210+
thisInfo.momentum[0], thisInfo.momentum[1], thisInfo.momentum[2]);
211+
cascmccollrefs(thisInfo.mcCollision);
212+
}
213+
214+
if (populateCascMCCoresAsymmetric) {
215+
int thisCascMCCoreIndex = -1;
216+
// step 1: check if this element is already provided in the table
217+
// using the packedIndices variable calculated above
218+
for (uint32_t ii = 0; ii < mcCascinfos.size(); ii++) {
219+
if (
220+
thisInfo.mcParticlePositive == mcCascinfos[ii].mcParticlePositive && mcCascinfos[ii].mcParticlePositive > 0 &&
221+
thisInfo.mcParticleNegative == mcCascinfos[ii].mcParticleNegative && mcCascinfos[ii].mcParticleNegative > 0 &&
222+
thisInfo.mcParticleBachelor == mcCascinfos[ii].mcParticleBachelor && mcCascinfos[ii].mcParticleBachelor > 0) {
223+
thisCascMCCoreIndex = ii;
224+
break; // this exists already in list
225+
}
226+
}
227+
if (thisCascMCCoreIndex < 0) {
228+
// this CascMCCore does not exist yet. Create it and reference it
229+
thisCascMCCoreIndex = mcCascinfos.size();
230+
mcCascinfos.push_back(thisInfo);
231+
}
232+
cascCoreMClabels(thisCascMCCoreIndex); // interlink: reconstructed -> MC index
151233
}
152234
} // end casctable loop
235+
236+
// now populate V0MCCores if in asymmetric mode
237+
if (populateCascMCCoresAsymmetric) {
238+
// first step: add any un-recoed v0mmcores that were requested
239+
for (auto& mcParticle : mcParticles) {
240+
thisInfo.pdgCode = -1, thisInfo.pdgCodeMother = -1;
241+
thisInfo.pdgCodePositive = -1, thisInfo.pdgCodeNegative = -1;
242+
thisInfo.pdgCodeBachelor = -1, thisInfo.pdgCodeV0 = -1;
243+
thisInfo.isPhysicalPrimary = false;
244+
thisInfo.xyz[0] = 0.0f, thisInfo.xyz[1] = 0.0f, thisInfo.xyz[2] = 0.0f;
245+
thisInfo.lxyz[0] = 0.0f, thisInfo.lxyz[1] = 0.0f, thisInfo.lxyz[2] = 0.0f;
246+
thisInfo.posP[0] = 0.0f, thisInfo.posP[1] = 0.0f, thisInfo.posP[2] = 0.0f;
247+
thisInfo.negP[0] = 0.0f, thisInfo.negP[1] = 0.0f, thisInfo.negP[2] = 0.0f;
248+
thisInfo.bachP[0] = 0.0f, thisInfo.bachP[1] = 0.0f, thisInfo.bachP[2] = 0.0f;
249+
thisInfo.momentum[0] = 0.0f, thisInfo.momentum[1] = 0.0f, thisInfo.momentum[2] = 0.0f;
250+
thisInfo.label = -1, thisInfo.motherLabel = -1;
251+
thisInfo.mcParticlePositive = -1;
252+
thisInfo.mcParticleNegative = -1;
253+
thisInfo.mcParticleBachelor = -1;
254+
255+
if (mcParticleIsReco[mcParticle.globalIndex()] == true)
256+
continue; // skip if already created in list
257+
258+
if (TMath::Abs(mcParticle.y()) > rapidityWindow)
259+
continue; // skip outside midrapidity
260+
261+
if (
262+
(addGeneratedXiMinus && mcParticle.pdgCode() == 3312) ||
263+
(addGeneratedXiPlus && mcParticle.pdgCode() == -3312) ||
264+
(addGeneratedOmegaMinus && mcParticle.pdgCode() == 3334) ||
265+
(addGeneratedOmegaPlus && mcParticle.pdgCode() == -3334)) {
266+
thisInfo.pdgCode = mcParticle.pdgCode();
267+
thisInfo.isPhysicalPrimary = mcParticle.isPhysicalPrimary();
268+
269+
if (mcParticle.has_mcCollision()) {
270+
thisInfo.mcCollision = mcParticle.mcCollisionId(); // save this reference, please
271+
}
272+
thisInfo.momentum[0] = mcParticle.px();
273+
thisInfo.momentum[1] = mcParticle.py();
274+
thisInfo.momentum[2] = mcParticle.pz();
275+
276+
// if I got here, it means this MC particle was not recoed and is of interest. Add it please
277+
mcCascinfos.push_back(thisInfo);
278+
}
279+
}
280+
281+
for (auto thisInfo : mcCascinfos) {
282+
cascmccores( // a lot of the info below will be compressed in case of not-recoed MC (good!)
283+
thisInfo.pdgCode, thisInfo.pdgCodeMother, thisInfo.pdgCodeV0, thisInfo.isPhysicalPrimary,
284+
thisInfo.pdgCodePositive, thisInfo.pdgCodeNegative, thisInfo.pdgCodeBachelor,
285+
thisInfo.xyz[0], thisInfo.xyz[1], thisInfo.xyz[2],
286+
thisInfo.lxyz[0], thisInfo.lxyz[1], thisInfo.lxyz[2],
287+
thisInfo.posP[0], thisInfo.posP[1], thisInfo.posP[2],
288+
thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2],
289+
thisInfo.bachP[0], thisInfo.bachP[1], thisInfo.bachP[2],
290+
thisInfo.momentum[0], thisInfo.momentum[1], thisInfo.momentum[2]);
291+
cascmccollrefs(thisInfo.mcCollision);
292+
}
293+
}
153294
}
154295

155296
//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
156297
// build cascade labels
157-
void processCascades(aod::CascDatas const& casctable, aod::V0sLinked const&, aod::V0Datas const& /*v0table*/, aod::McTrackLabels const&, aod::McParticles const&)
298+
void processCascades(aod::CascDatas const& casctable, aod::V0sLinked const&, aod::V0Datas const& /*v0table*/, aod::McTrackLabels const&, aod::McParticles const& mcParticles)
158299
{
159-
generateCascadeMCinfo(casctable);
300+
generateCascadeMCinfo(casctable, mcParticles);
160301
}
161302

162303
//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
163304
// build findable cascade labels
164-
void processFindableCascades(aod::CascDatas const& casctable, aod::FindableV0sLinked const&, aod::V0Datas const& /*v0table*/, aod::McTrackLabels const&, aod::McParticles const&)
305+
void processFindableCascades(aod::CascDatas const& casctable, aod::FindableV0sLinked const&, aod::V0Datas const& /*v0table*/, aod::McTrackLabels const&, aod::McParticles const& mcParticles)
165306
{
166-
generateCascadeMCinfo(casctable);
307+
generateCascadeMCinfo(casctable, mcParticles);
167308
}
168309

169310
//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*

0 commit comments

Comments
 (0)