Skip to content

Commit 2d3e217

Browse files
ddobrigkalibuildromainschotter
authored
[PWGLF] Handle kink decays in association for V0/casc labeler (#8461)
Co-authored-by: ALICE Builder <alibuild@users.noreply.github.com> Co-authored-by: romainschotter <romain.schotter@gmail.com>
1 parent 8b2d8c1 commit 2d3e217

File tree

3 files changed

+173
-109
lines changed

3 files changed

+173
-109
lines changed

‎PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx‎

Lines changed: 76 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,8 @@ struct cascademcbuilder {
6161
Configurable<bool> addGeneratedOmegaMinus{"addGeneratedOmegaMinus", false, "add CascMCCore entry for generated, not-recoed OmegaMinus"};
6262
Configurable<bool> addGeneratedOmegaPlus{"addGeneratedOmegaPlus", false, "add CascMCCore entry for generated, not-recoed OmegaPlus"};
6363

64+
Configurable<bool> treatPiToMuDecays{"treatPiToMuDecays", true, "if true, will correctly capture pi -> mu and V0 label will still point to originating V0 decay in those cases. Nota bene: prong info will still be for the muon!"};
65+
6466
Configurable<float> rapidityWindow{"rapidityWindow", 0.5, "rapidity window to save non-recoed candidates"};
6567

6668
//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
@@ -92,6 +94,36 @@ struct cascademcbuilder {
9294
mcCascinfo thisInfo;
9395
//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
9496

97+
// kink handling
98+
template <typename mcpart>
99+
int getOriginatingParticle(mcpart const& part, int& indexForPositionOfDecay)
100+
{
101+
int returnValue = -1;
102+
if (part.has_mothers()) {
103+
auto const& motherList = part.template mothers_as<aod::McParticles>();
104+
if (motherList.size() == 1) {
105+
for (const auto& mother : motherList) {
106+
if (std::abs(part.pdgCode()) == 13 && treatPiToMuDecays) {
107+
// muon decay, de-ref mother twice
108+
if (mother.has_mothers()) {
109+
auto grandMotherList = mother.template mothers_as<aod::McParticles>();
110+
if (grandMotherList.size() == 1) {
111+
for (const auto& grandMother : grandMotherList) {
112+
returnValue = grandMother.globalIndex();
113+
indexForPositionOfDecay = mother.globalIndex(); // for V0 decay position: grab muon
114+
}
115+
}
116+
}
117+
} else {
118+
returnValue = mother.globalIndex();
119+
indexForPositionOfDecay = part.globalIndex();
120+
}
121+
}
122+
}
123+
}
124+
return returnValue;
125+
}
126+
95127
template <typename TCascadeTable, typename TMCParticleTable>
96128
void generateCascadeMCinfo(TCascadeTable cascTable, TMCParticleTable mcParticles)
97129
{
@@ -147,51 +179,51 @@ struct cascademcbuilder {
147179
thisInfo.processNegative = lMCNegTrack.getProcess();
148180
thisInfo.processBachelor = lMCBachTrack.getProcess();
149181

150-
// Step 1: check if the mother is the same, go up a level
151-
if (lMCNegTrack.has_mothers() && lMCPosTrack.has_mothers()) {
152-
for (auto& lNegMother : lMCNegTrack.template mothers_as<aod::McParticles>()) {
153-
for (auto& lPosMother : lMCPosTrack.template mothers_as<aod::McParticles>()) {
154-
if (lNegMother == lPosMother) {
155-
// acquire information
156-
thisInfo.lxyz[0] = lMCPosTrack.vx();
157-
thisInfo.lxyz[1] = lMCPosTrack.vy();
158-
thisInfo.lxyz[2] = lMCPosTrack.vz();
159-
thisInfo.pdgCodeV0 = lNegMother.pdgCode();
182+
// Step 0: treat pi -> mu + antineutrino
183+
// if present, de-reference original V0 correctly and provide label to original object
184+
// NOTA BENE: the prong info will still correspond to a muon, treat carefully!
185+
int negOriginating = -1, posOriginating = -1, bachOriginating = -1;
186+
int particleForLambdaDecayPositionIdx = -1, particleForCascadeDecayPositionIdx = -1;
187+
negOriginating = getOriginatingParticle(lMCNegTrack, particleForLambdaDecayPositionIdx);
188+
posOriginating = getOriginatingParticle(lMCPosTrack, particleForLambdaDecayPositionIdx);
189+
bachOriginating = getOriginatingParticle(lMCBachTrack, particleForCascadeDecayPositionIdx);
190+
191+
if (negOriginating > -1 && negOriginating == posOriginating) {
192+
auto originatingV0 = mcParticles.rawIteratorAt(negOriginating);
193+
auto particleForLambdaDecayPosition = mcParticles.rawIteratorAt(particleForLambdaDecayPositionIdx);
194+
195+
thisInfo.label = originatingV0.globalIndex();
196+
thisInfo.lxyz[0] = particleForLambdaDecayPosition.vx();
197+
thisInfo.lxyz[1] = particleForLambdaDecayPosition.vy();
198+
thisInfo.lxyz[2] = particleForLambdaDecayPosition.vz();
199+
200+
if (originatingV0.has_mothers()) {
201+
for (auto& lV0Mother : originatingV0.template mothers_as<aod::McParticles>()) {
202+
if (lV0Mother.globalIndex() == bachOriginating) { // found mother particle
203+
thisInfo.label = lV0Mother.globalIndex();
204+
205+
if (lV0Mother.has_mcCollision()) {
206+
thisInfo.mcCollision = lV0Mother.mcCollisionId(); // save this reference, please
207+
}
160208

161-
// if we got to this level, it means the mother particle exists and is the same
162-
// now we have to go one level up and compare to the bachelor mother too
163-
if (lNegMother.has_mothers() && lMCBachTrack.has_mothers()) {
164-
for (auto& lV0Mother : lNegMother.template mothers_as<aod::McParticles>()) {
165-
for (auto& lBachMother : lMCBachTrack.template mothers_as<aod::McParticles>()) {
166-
if (lV0Mother == lBachMother) {
167-
thisInfo.label = lV0Mother.globalIndex();
168-
169-
if (lV0Mother.has_mcCollision()) {
170-
thisInfo.mcCollision = lV0Mother.mcCollisionId(); // save this reference, please
171-
}
172-
173-
thisInfo.pdgCode = lV0Mother.pdgCode();
174-
thisInfo.isPhysicalPrimary = lV0Mother.isPhysicalPrimary();
175-
thisInfo.xyz[0] = lMCBachTrack.vx();
176-
thisInfo.xyz[1] = lMCBachTrack.vy();
177-
thisInfo.xyz[2] = lMCBachTrack.vz();
178-
thisInfo.momentum[0] = lV0Mother.px();
179-
thisInfo.momentum[1] = lV0Mother.py();
180-
thisInfo.momentum[2] = lV0Mother.pz();
181-
if (lV0Mother.has_mothers()) {
182-
for (auto& lV0GrandMother : lV0Mother.template mothers_as<aod::McParticles>()) {
183-
thisInfo.pdgCodeMother = lV0GrandMother.pdgCode();
184-
thisInfo.motherLabel = lV0GrandMother.globalIndex();
185-
}
186-
}
187-
}
188-
}
189-
} // end conditional V0-bach pair
190-
} // end has mothers
191-
} // end neg = pos mother conditional
192-
}
193-
} // end loop neg/pos mothers
194-
} // end conditional of mothers existing
209+
thisInfo.pdgCode = lV0Mother.pdgCode();
210+
thisInfo.isPhysicalPrimary = lV0Mother.isPhysicalPrimary();
211+
thisInfo.xyz[0] = originatingV0.vx();
212+
thisInfo.xyz[1] = originatingV0.vy();
213+
thisInfo.xyz[2] = originatingV0.vz();
214+
thisInfo.momentum[0] = lV0Mother.px();
215+
thisInfo.momentum[1] = lV0Mother.py();
216+
thisInfo.momentum[2] = lV0Mother.pz();
217+
if (lV0Mother.has_mothers()) {
218+
for (auto& lV0GrandMother : lV0Mother.template mothers_as<aod::McParticles>()) {
219+
thisInfo.pdgCodeMother = lV0GrandMother.pdgCode();
220+
thisInfo.motherLabel = lV0GrandMother.globalIndex();
221+
}
222+
}
223+
}
224+
} // end v0 mother loop
225+
} // end has_mothers check for V0
226+
} // end conditional of pos/neg originating being the same
195227
} // end association check
196228
// Construct label table (note: this will be joinable with CascDatas)
197229
casclabels(

‎PWGLF/TableProducer/Strangeness/lambdakzeromcbuilder.cxx‎

Lines changed: 89 additions & 61 deletions
Original file line numberDiff line numberDiff line change
@@ -55,6 +55,8 @@ struct lambdakzeromcbuilder {
5555
Configurable<bool> addGeneratedAntiLambda{"addGeneratedAntiLambda", false, "add V0MCCore entry for generated, not-recoed AntiLambda"};
5656
Configurable<bool> addGeneratedGamma{"addGeneratedGamma", false, "add V0MCCore entry for generated, not-recoed Gamma"};
5757

58+
Configurable<bool> treatPiToMuDecays{"treatPiToMuDecays", true, "if true, will correctly capture pi -> mu and V0 label will still point to originating V0 decay in those cases. Nota bene: prong info will still be for the muon!"};
59+
5860
Configurable<float> rapidityWindow{"rapidityWindow", 0.5, "rapidity window to save non-recoed candidates"};
5961

6062
HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
@@ -111,15 +113,38 @@ struct lambdakzeromcbuilder {
111113
std::array<float, 3> posP;
112114
std::array<float, 3> negP;
113115
std::array<float, 3> momentum;
114-
uint64_t packedMcParticleIndices;
115116
};
116117
mcV0info thisInfo;
117118
//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
118119

119-
// prong index combiner
120-
uint64_t combineProngIndices(uint32_t low, uint32_t high)
120+
// kink handling
121+
template <typename mcpart>
122+
int getOriginatingParticle(mcpart const& part, int& indexForPositionOfDecay)
121123
{
122-
return (((uint64_t)high) << 32) | ((uint64_t)low);
124+
int returnValue = -1;
125+
if (part.has_mothers()) {
126+
auto const& motherList = part.template mothers_as<aod::McParticles>();
127+
if (motherList.size() == 1) {
128+
for (const auto& mother : motherList) {
129+
if (std::abs(part.pdgCode()) == 13 && treatPiToMuDecays) {
130+
// muon decay, de-ref mother twice
131+
if (mother.has_mothers()) {
132+
auto grandMotherList = mother.template mothers_as<aod::McParticles>();
133+
if (grandMotherList.size() == 1) {
134+
for (const auto& grandMother : grandMotherList) {
135+
returnValue = grandMother.globalIndex();
136+
indexForPositionOfDecay = mother.globalIndex(); // for V0 decay position: grab muon
137+
}
138+
}
139+
}
140+
} else {
141+
returnValue = mother.globalIndex();
142+
indexForPositionOfDecay = part.globalIndex();
143+
}
144+
}
145+
}
146+
}
147+
return returnValue;
123148
}
124149

125150
//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
@@ -131,7 +156,6 @@ struct lambdakzeromcbuilder {
131156
std::vector<bool> mcParticleIsReco(mcParticles.size(), false); // mc Particle not recoed by V0s
132157

133158
for (auto& v0 : v0table) {
134-
thisInfo.packedMcParticleIndices = 0; // not de-referenced properly yet
135159
thisInfo.label = -1;
136160
thisInfo.motherLabel = -1;
137161
thisInfo.pdgCode = 0;
@@ -152,7 +176,6 @@ struct lambdakzeromcbuilder {
152176
auto lMCNegTrack = lNegTrack.mcParticle_as<aod::McParticles>();
153177
auto lMCPosTrack = lPosTrack.mcParticle_as<aod::McParticles>();
154178

155-
thisInfo.packedMcParticleIndices = combineProngIndices(lPosTrack.mcParticleId(), lNegTrack.mcParticleId());
156179
thisInfo.pdgCodePositive = lMCPosTrack.pdgCode();
157180
thisInfo.pdgCodeNegative = lMCNegTrack.pdgCode();
158181
thisInfo.processPositive = lMCPosTrack.getProcess();
@@ -163,65 +186,71 @@ struct lambdakzeromcbuilder {
163186
thisInfo.negP[0] = lMCNegTrack.px();
164187
thisInfo.negP[1] = lMCNegTrack.py();
165188
thisInfo.negP[2] = lMCNegTrack.pz();
166-
if (lMCNegTrack.has_mothers() && lMCPosTrack.has_mothers()) {
167-
for (auto& lNegMother : lMCNegTrack.mothers_as<aod::McParticles>()) {
168-
for (auto& lPosMother : lMCPosTrack.mothers_as<aod::McParticles>()) {
169-
if (lNegMother.globalIndex() == lPosMother.globalIndex()) {
170-
thisInfo.label = lNegMother.globalIndex();
171-
thisInfo.xyz[0] = lMCPosTrack.vx();
172-
thisInfo.xyz[1] = lMCPosTrack.vy();
173-
thisInfo.xyz[2] = lMCPosTrack.vz();
174-
175-
// MC pos. and neg. daughters are the same! Looking for replacement...
176-
if (lMCPosTrack.globalIndex() == lMCNegTrack.globalIndex()) {
177-
auto const& daughters = lNegMother.daughters_as<aod::McParticles>();
178-
for (auto& ldau : daughters) {
179-
// check if the candidate originate from a decay
180-
// if not, this is not a suitable candidate for one of the decay daughters
181-
if (ldau.getProcess() != 4) // see TMCProcess.h
182-
continue;
183-
184-
if (lMCPosTrack.pdgCode() < 0 && ldau.pdgCode() > 0) { // the positive track needs to be changed
185-
thisInfo.pdgCodePositive = ldau.pdgCode();
186-
thisInfo.processPositive = ldau.getProcess();
187-
thisInfo.posP[0] = ldau.px();
188-
thisInfo.posP[1] = ldau.py();
189-
thisInfo.posP[2] = ldau.pz();
190-
thisInfo.xyz[0] = ldau.vx();
191-
thisInfo.xyz[1] = ldau.vy();
192-
thisInfo.xyz[2] = ldau.vz();
193-
}
194-
if (lMCNegTrack.pdgCode() > 0 && ldau.pdgCode() < 0) { // the negative track needs to be changed
195-
thisInfo.pdgCodeNegative = ldau.pdgCode();
196-
thisInfo.processNegative = ldau.getProcess();
197-
thisInfo.negP[0] = ldau.px();
198-
thisInfo.negP[1] = ldau.py();
199-
thisInfo.negP[2] = ldau.pz();
200-
}
201-
}
202-
}
203189

204-
if (lNegMother.has_mcCollision()) {
205-
thisInfo.mcCollision = lNegMother.mcCollisionId(); // save this reference, please
206-
}
190+
// check for pi -> mu + antineutrino decay
191+
// if present, de-reference original V0 correctly and provide label to original object
192+
// NOTA BENE: the prong info will still correspond to a muon, treat carefully!
193+
int negOriginating = -1, posOriginating = -1, particleForDecayPositionIdx = -1;
194+
negOriginating = getOriginatingParticle(lMCNegTrack, particleForDecayPositionIdx);
195+
posOriginating = getOriginatingParticle(lMCPosTrack, particleForDecayPositionIdx);
196+
197+
if (negOriginating > -1 && negOriginating == posOriginating) {
198+
auto originatingV0 = mcParticles.rawIteratorAt(negOriginating);
199+
auto particleForDecayPosition = mcParticles.rawIteratorAt(particleForDecayPositionIdx);
200+
201+
thisInfo.label = originatingV0.globalIndex();
202+
thisInfo.xyz[0] = particleForDecayPosition.vx();
203+
thisInfo.xyz[1] = particleForDecayPosition.vy();
204+
thisInfo.xyz[2] = particleForDecayPosition.vz();
205+
206+
// MC pos. and neg. daughters are the same! Looking for replacement...
207+
// if (lMCPosTrack.globalIndex() == lMCNegTrack.globalIndex()) {
208+
// auto const& daughters = lNegMother.daughters_as<aod::McParticles>();
209+
// for (auto& ldau : daughters) {
210+
// // check if the candidate originates from a decay
211+
// // if not, this is not a suitable candidate for one of the decay daughters
212+
// if (ldau.getProcess() != 4) // see TMCProcess.h
213+
// continue;
214+
215+
// if (lMCPosTrack.pdgCode() < 0 && ldau.pdgCode() > 0) { // the positive track needs to be changed
216+
// thisInfo.pdgCodePositive = ldau.pdgCode();
217+
// thisInfo.processPositive = ldau.getProcess();
218+
// thisInfo.posP[0] = ldau.px();
219+
// thisInfo.posP[1] = ldau.py();
220+
// thisInfo.posP[2] = ldau.pz();
221+
// thisInfo.xyz[0] = ldau.vx();
222+
// thisInfo.xyz[1] = ldau.vy();
223+
// thisInfo.xyz[2] = ldau.vz();
224+
// }
225+
// if (lMCNegTrack.pdgCode() > 0 && ldau.pdgCode() < 0) { // the negative track needs to be changed
226+
// thisInfo.pdgCodeNegative = ldau.pdgCode();
227+
// thisInfo.processNegative = ldau.getProcess();
228+
// thisInfo.negP[0] = ldau.px();
229+
// thisInfo.negP[1] = ldau.py();
230+
// thisInfo.negP[2] = ldau.pz();
231+
// }
232+
// }
233+
// }
234+
235+
if (originatingV0.has_mcCollision()) {
236+
thisInfo.mcCollision = originatingV0.mcCollisionId(); // save this reference, please
237+
}
207238

208-
// acquire information
209-
thisInfo.pdgCode = lNegMother.pdgCode();
210-
thisInfo.isPhysicalPrimary = lNegMother.isPhysicalPrimary();
211-
thisInfo.momentum[0] = lNegMother.px();
212-
thisInfo.momentum[1] = lNegMother.py();
213-
thisInfo.momentum[2] = lNegMother.pz();
214-
215-
if (lNegMother.has_mothers()) {
216-
for (auto& lNegGrandMother : lNegMother.mothers_as<aod::McParticles>()) {
217-
thisInfo.pdgCodeMother = lNegGrandMother.pdgCode();
218-
thisInfo.motherLabel = lNegGrandMother.globalIndex();
219-
}
220-
}
221-
}
239+
// acquire information
240+
thisInfo.pdgCode = originatingV0.pdgCode();
241+
thisInfo.isPhysicalPrimary = originatingV0.isPhysicalPrimary();
242+
thisInfo.momentum[0] = originatingV0.px();
243+
thisInfo.momentum[1] = originatingV0.py();
244+
thisInfo.momentum[2] = originatingV0.pz();
245+
246+
if (originatingV0.has_mothers()) {
247+
for (auto& lV0Mother : originatingV0.mothers_as<aod::McParticles>()) {
248+
thisInfo.pdgCodeMother = lV0Mother.pdgCode();
249+
thisInfo.motherLabel = lV0Mother.globalIndex();
222250
}
223251
}
224252
}
253+
225254
} // end association check
226255
// Construct label table (note: this will be joinable with V0Datas!)
227256
v0labels(
@@ -308,7 +337,6 @@ struct lambdakzeromcbuilder {
308337
if (populateV0MCCoresAsymmetric) {
309338
// first step: add any un-recoed v0mmcores that were requested
310339
for (auto& mcParticle : mcParticles) {
311-
thisInfo.packedMcParticleIndices = 0;
312340
thisInfo.label = -1;
313341
thisInfo.motherLabel = -1;
314342
thisInfo.pdgCode = 0;

‎PWGLF/Tasks/Strangeness/derivedlambdakzeroanalysis.cxx‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -148,6 +148,7 @@ struct derivedlambdakzeroanalysis {
148148

149149
// for MC
150150
Configurable<bool> doMCAssociation{"doMCAssociation", true, "if MC, do MC association"};
151+
Configurable<bool> doTreatPiToMuon{"doTreatPiToMuon", false, "Take pi decay into muon into account in MC"};
151152
Configurable<bool> doCollisionAssociationQA{"doCollisionAssociationQA", true, "check collision association"};
152153

153154
// Machine learning evaluation for pre-selection and corresponding information generation
@@ -816,19 +817,22 @@ struct derivedlambdakzeroanalysis {
816817
// precalculate this information so that a check is one mask operation, not many
817818
{
818819
uint64_t bitMap = 0;
819-
// check for specific particle species
820+
bool isPositiveProton = v0.pdgCodePositive() == 2212;
821+
bool isPositivePion = v0.pdgCodePositive() == 211 || (doTreatPiToMuon && v0.pdgCodePositive() == -13);
822+
bool isNegativeProton = v0.pdgCodeNegative() == -2212;
823+
bool isNegativePion = v0.pdgCodeNegative() == -211 || (doTreatPiToMuon && v0.pdgCodeNegative() == 13);
820824

821-
if (v0.pdgCode() == 310 && v0.pdgCodePositive() == 211 && v0.pdgCodeNegative() == -211) {
825+
if (v0.pdgCode() == 310 && isPositivePion && isNegativePion) {
822826
bitset(bitMap, selConsiderK0Short);
823827
if (v0.isPhysicalPrimary())
824828
bitset(bitMap, selPhysPrimK0Short);
825829
}
826-
if (v0.pdgCode() == 3122 && v0.pdgCodePositive() == 2212 && v0.pdgCodeNegative() == -211) {
830+
if (v0.pdgCode() == 3122 && isPositiveProton && isNegativePion) {
827831
bitset(bitMap, selConsiderLambda);
828832
if (v0.isPhysicalPrimary())
829833
bitset(bitMap, selPhysPrimLambda);
830834
}
831-
if (v0.pdgCode() == -3122 && v0.pdgCodePositive() == 211 && v0.pdgCodeNegative() == -2212) {
835+
if (v0.pdgCode() == -3122 && isPositivePion && isNegativeProton) {
832836
bitset(bitMap, selConsiderAntiLambda);
833837
if (v0.isPhysicalPrimary())
834838
bitset(bitMap, selPhysPrimAntiLambda);

0 commit comments

Comments
 (0)