Skip to content

Commit 76134d5

Browse files
morgmattmorgmatt
committed
Removed nested loop from processFindable method, substituted with index map structure. Added decayradius detailed analysis of the filters
1 parent d8f1e57 commit 76134d5

File tree

1 file changed

+184
-112
lines changed

1 file changed

+184
-112
lines changed

‎PWGLF/TableProducer/Strangeness/sigmaminustask.cxx‎

Lines changed: 184 additions & 112 deletions
Original file line numberDiff line numberDiff line change
@@ -18,11 +18,9 @@
1818
#include "Common/DataModel/EventSelection.h"
1919
#include "Common/DataModel/PIDResponse.h"
2020

21-
#include "Common/Core/trackUtilities.h"
2221
#include "Framework/AnalysisTask.h"
2322
#include "Framework/runDataProcessing.h"
2423
#include "ReconstructionDataFormats/PID.h"
25-
#include "ReconstructionDataFormats/Track.h"
2624

2725
using namespace o2;
2826
using namespace o2::framework;
@@ -66,13 +64,13 @@ struct sigmaminustask {
6664
const AxisSpec xiMassAxis{100, 1.2, 1.6, "m_{#Xi} (GeV/#it{c}^{2})"};
6765
const AxisSpec pdgAxis{10001, -5000, 5000, "PDG code"};
6866
const AxisSpec vertexZAxis{100, -15., 15., "vrtx_{Z} [cm]"};
69-
const AxisSpec dcaMothAxis{100, 0, 1, "DCA [cm]"};
67+
const AxisSpec dcaMothAxis{100, 0, 0.2, "DCA [cm]"};
7068
const AxisSpec dcaDaugAxis{200, 0, 20, "DCA [cm]"};
7169
const AxisSpec radiusAxis{100, -1, 40, "Decay radius [cm]"};
7270

73-
const AxisSpec ptResolutionAxis{100, -2, 2, "(#it{p}_{T}^{rec} - #it{p}_{T}^{gen}) / #it{p}_{T}^{gen}"};
74-
const AxisSpec massResolutionAxis{100, -2, 2, "(m_{rec} - m_{gen}) / m_{gen}"};
75-
const AxisSpec radiusResolutionAxis{100, -2, 2, "(r_{rec} - r_{gen}) / r_{gen}"};
71+
const AxisSpec ptResolutionAxis{100, -1.0, 1.0, "(#it{p}_{T}^{rec} - #it{p}_{T}^{gen}) / #it{p}_{T}^{gen}"};
72+
const AxisSpec massResolutionAxis{100, -0.5, 0.5, "(m_{rec} - m_{gen}) / m_{gen}"};
73+
const AxisSpec radiusResolutionAxis{100, -0.5, 0.5, "(r_{rec} - r_{gen}) / r_{gen}"};
7674

7775
const AxisSpec boolAxis{2, -0.5, 1.5, "Boolean value (0=false, 1=true)"};
7876
const AxisSpec filtersAxis{10, -0.5, 9.5, "Filter index"};
@@ -93,7 +91,8 @@ struct sigmaminustask {
9391

9492
rSigmaMinus.add("h2MassResolution", "h2MassResolution", {HistType::kTH2F, {ptAxis, massResolutionAxis}});
9593
rSigmaMinus.add("h2PtResolution", "h2PtResolution", {HistType::kTH2F, {ptAxis, ptResolutionAxis}});
96-
94+
rSigmaMinus.add("h2RadiusResolution", "h2RadiusResolution", {HistType::kTH2F, {ptAxis, radiusResolutionAxis}});
95+
9796
rSigmaMinus.add("h2NSigmaTOFPiPt", "h2NSigmaTOFPiPt", {HistType::kTH2F, {ptAxis, nSigmaPiAxis}});
9897
rSigmaMinus.add("h2NSigmaTOFPrPt", "h2NSigmaTOFPrPt", {HistType::kTH2F, {ptAxis, nSigmaPrAxis}});
9998

@@ -108,9 +107,8 @@ struct sigmaminustask {
108107
// Add findable Sigma histograms
109108
rFindable.add("h2MassPtFindableAll", "h2MassPtFindableAll", {HistType::kTH2F, {ptAxis, sigmaMassAxis}});
110109
rFindable.add("hFilterIndex", "hFilterIndex", {HistType::kTH1F, {filtersAxis}});
111-
rFindable.add("h2MCRadiusFilterIndex", "h2RadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}});
110+
rFindable.add("h2MCRadiusFilterIndex", "h2MCRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}});
112111
rFindable.add("h2RecRadiusFilterIndex", "h2RecRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}});
113-
114112
}
115113
}
116114

@@ -293,114 +291,188 @@ struct sigmaminustask {
293291

294292
PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false);
295293

296-
void processFindable(TracksFull const& tracks, aod::McTrackLabels const& trackLabelsMC, aod::KinkCands const& kinkCands, aod::McParticles const&, CollisionsFullMC const&)
297-
{
298-
for (const auto& motherTrack : tracks) {
299-
// Check if mother is Sigma in MC
300-
auto mcLabelMother = trackLabelsMC.rawIteratorAt(motherTrack.globalIndex());
301-
if (!mcLabelMother.has_mcParticle()) {
294+
void processFindable(aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC,
295+
TracksFull const& tracks, CollisionsFullMC const& mcCollisions){
296+
// A - generated findable track pairs map: mcMother.globalIndex() -> (motherTrack.globalIndex(), daughterTrack.globalIndex())
297+
std::unordered_map<int64_t, std::pair<int64_t, int64_t>> allCandsIndices;
298+
299+
for (const auto& track : tracks) {
300+
auto mcLabel = trackLabelsMC.rawIteratorAt(track.globalIndex());
301+
if (!mcLabel.has_mcParticle()) {
302302
continue;
303303
}
304-
auto mcMother = mcLabelMother.mcParticle_as<aod::McParticles>();
305-
if (std::abs(mcMother.pdgCode()) != 3112 && std::abs(mcMother.pdgCode()) != 3222) {
304+
auto mcParticle = mcLabel.mcParticle_as<aod::McParticles>();
305+
306+
if (mcParticle.has_daughters() && (std::abs(mcParticle.pdgCode()) == 3112 || std::abs(mcParticle.pdgCode()) == 3222)) {
307+
allCandsIndices[mcParticle.globalIndex()] = {track.globalIndex(), -1};
308+
}
309+
}
310+
311+
for (const auto& track : tracks) {
312+
auto mcLabel = trackLabelsMC.rawIteratorAt(track.globalIndex());
313+
if (!mcLabel.has_mcParticle()) {
306314
continue;
307315
}
308-
309-
for (const auto& daughterTrack : tracks) {
310-
// Check if daughter is pi/proton
311-
auto mcLabelDaughter = trackLabelsMC.rawIteratorAt(daughterTrack.globalIndex());
312-
if (!mcLabelDaughter.has_mcParticle()) {
313-
continue;
314-
}
315-
auto mcDaughter = mcLabelDaughter.mcParticle_as<aod::McParticles>();
316-
if (std::abs(mcDaughter.pdgCode()) != 211 && std::abs(mcDaughter.pdgCode()) != 2212) {
317-
continue;
318-
}
319-
320-
// Verify the MC mother-daughter relationship
321-
bool isValidPair = false;
322-
if (mcDaughter.has_mothers()) {
323-
for (const auto& mother : mcDaughter.mothers_as<aod::McParticles>()) {
324-
if (mother.globalIndex() == mcMother.globalIndex()) {
325-
isValidPair = true;
326-
break;
327-
}
316+
auto mcParticle = mcLabel.mcParticle_as<aod::McParticles>();
317+
318+
if (mcParticle.has_mothers() && (std::abs(mcParticle.pdgCode()) == 211 || std::abs(mcParticle.pdgCode()) == 2212)) {
319+
for (const auto& mother : mcParticle.mothers_as<aod::McParticles>()) {
320+
auto it = allCandsIndices.find(mother.globalIndex());
321+
if (it != allCandsIndices.end()) {
322+
it->second.second = track.globalIndex();
323+
break;
328324
}
329325
}
330-
if (!isValidPair) {
331-
continue;
332-
}
326+
}
327+
}
333328

334-
float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p());
335-
int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1;
336-
rSigmaMinus.fill(HIST("h2MassPtFindable"), sigmaSign * motherTrack.pt(), mcMass);
337-
338-
// Define filter index and progressively apply kinkbuilder cuts to track pairs
339-
int filterIndex = 0;
340-
rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex);
341-
342-
// 1 - tracks with right ITS, TPC, TOF signals
343-
if (motherTrack.has_collision() && motherTrack.hasITS() && !motherTrack.hasTPC() && !motherTrack.hasTOF() &&
344-
daughterTrack.hasITS() && daughterTrack.hasTPC()) {
345-
filterIndex += 1;
346-
rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex);
347-
} else {
348-
continue;
349-
}
329+
// B - reconstructed kinkcands map: mcMother.globalIndex() -> kinkCand.globalIndex()
330+
std::unordered_map<int64_t, int64_t> findableToKinkCand;
331+
for (const auto& kinkCand : kinkCands) {
332+
auto motherTrack = kinkCand.trackMoth_as<TracksFull>();
333+
auto daughterTrack = kinkCand.trackDaug_as<TracksFull>();
334+
335+
auto mcLabMoth = trackLabelsMC.rawIteratorAt(motherTrack.globalIndex());
336+
auto mcLabDaug = trackLabelsMC.rawIteratorAt(daughterTrack.globalIndex());
337+
if (!mcLabMoth.has_mcParticle() || !mcLabDaug.has_mcParticle()) {
338+
continue;
339+
}
340+
auto mcMother = mcLabMoth.mcParticle_as<aod::McParticles>();
341+
auto mcDaughter = mcLabDaug.mcParticle_as<aod::McParticles>();
342+
343+
if (std::abs(mcMother.pdgCode()) != 3112 && std::abs(mcMother.pdgCode()) != 3222) {
344+
continue;
345+
}
346+
if (std::abs(mcDaughter.pdgCode()) != 211 && std::abs(mcDaughter.pdgCode()) != 2212) {
347+
continue;
348+
}
349+
350+
auto findableIt = allCandsIndices.find(mcMother.globalIndex());
351+
if (findableIt != allCandsIndices.end() &&
352+
findableIt->second.first == motherTrack.globalIndex() &&
353+
findableIt->second.second == daughterTrack.globalIndex()) {
354+
355+
findableToKinkCand[mcMother.globalIndex()] = kinkCand.globalIndex();
356+
}
357+
}
358+
359+
360+
// C - loop on valid pairs for findable analysis
361+
for (const auto& [mcMotherIndex, trackIndices] : allCandsIndices) {
362+
if (trackIndices.second == -1 || trackIndices.first == -1) {
363+
continue;
364+
}
365+
// Retrieve mother and daughter tracks and mcParticles
366+
auto motherTrack = tracks.rawIteratorAt(trackIndices.first);
367+
auto daughterTrack = tracks.rawIteratorAt(trackIndices.second);
368+
auto mcLabMoth = trackLabelsMC.rawIteratorAt(motherTrack.globalIndex());
369+
auto mcLabDaug = trackLabelsMC.rawIteratorAt(daughterTrack.globalIndex());
370+
if (!mcLabMoth.has_mcParticle() || !mcLabDaug.has_mcParticle()) {
371+
continue;
372+
}
373+
auto mcMother = mcLabMoth.mcParticle_as<aod::McParticles>();
374+
auto mcDaughter = mcLabDaug.mcParticle_as<aod::McParticles>();
375+
376+
// Compute mass and radii
377+
float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p());
378+
int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1;
379+
float mcRadius = std::sqrt((mcMother.vx() - mcDaughter.vx()) * (mcMother.vx() - mcDaughter.vx()) +
380+
(mcMother.vy() - mcDaughter.vy()) * (mcMother.vy() - mcDaughter.vy()));
381+
float recRadius = -1.0;
382+
if (findableToKinkCand.find(mcMother.globalIndex()) != findableToKinkCand.end()) {
383+
auto kinkCand = kinkCands.rawIteratorAt(findableToKinkCand[mcMother.globalIndex()]);
384+
recRadius = std::sqrt(kinkCand.xDecVtx() * kinkCand.xDecVtx() + kinkCand.yDecVtx() * kinkCand.yDecVtx());
385+
}
386+
387+
rFindable.fill(HIST("h2MassPtFindableAll"), sigmaSign * mcMother.pt(), mcMass);
388+
389+
// Define filter index and progressively apply kinkbuilder cuts to track pairs
390+
int filterIndex = 0;
391+
rFindable.fill(HIST("hFilterIndex"), filterIndex);
392+
rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius);
393+
rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius);
394+
395+
// 1 - tracks with right ITS, TPC, TOF signals
396+
if (motherTrack.has_collision() && motherTrack.hasITS() && !motherTrack.hasTPC() && !motherTrack.hasTOF() &&
397+
daughterTrack.hasITS() && daughterTrack.hasTPC()) {
398+
filterIndex += 1;
399+
rFindable.fill(HIST("hFilterIndex"), filterIndex);
400+
rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius);
401+
rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius);
402+
} else {
403+
continue;
404+
}
350405

351-
// 2, 3 - mother track ITS properties
352-
if (motherTrack.itsNCls() < 6 &&
406+
// 2, 3 - mother track ITS properties
407+
if (motherTrack.itsNCls() < 6 &&
353408
motherTrack.itsNClsInnerBarrel() == 3 && motherTrack.itsChi2NCl() < 36) {
354-
filterIndex += 1;
355-
rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex);
356-
}
357-
if (filterIndex > 1 && motherTrack.pt() > 0.5) {
358-
filterIndex += 1;
359-
rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex);
360-
} else {
361-
continue;
362-
}
409+
filterIndex += 1;
410+
rFindable.fill(HIST("hFilterIndex"), filterIndex);
411+
rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius);
412+
rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius);
413+
} else {
414+
continue;
415+
}
363416

364-
// 4 - daughter track ITS+TPC properties
365-
if (daughterTrack.itsNClsInnerBarrel() == 0 && daughterTrack.itsNCls() < 4 &&
366-
daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > 80) {
367-
filterIndex += 1;
368-
rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex);
369-
} else {
370-
continue;
371-
}
372-
373-
// 5 - geometric cuts: eta
374-
if (std::abs(motherTrack.eta()) < 1.0 && std::abs(daughterTrack.eta()) < 1.0) {
375-
filterIndex += 1;
376-
rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex);
377-
} else {
378-
continue;
379-
}
380-
// 6 - geometric cuts: phi difference
381-
if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < 50.0) {
382-
filterIndex += 1;
383-
rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex);
384-
} else {
385-
continue;
386-
}
387-
// 7 - geometric cuts: z difference
388-
o2::track::TrackParCov trackParCovMoth = getTrackParCov(motherTrack);
389-
o2::track::TrackParCov trackParCovDaug = getTrackParCov(daughterTrack);
390-
if (std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()) < 20.0) {
391-
filterIndex += 1;
392-
rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex);
393-
} else {
394-
continue;
395-
}
396-
// 8 - collision selection
397-
auto collision = motherTrack.template collision_as<CollisionsFullMC>();
398-
if (!(std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) {
399-
filterIndex += 1;
400-
rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex);
401-
}
402-
417+
if (motherTrack.pt() > 0.5) {
418+
filterIndex += 1;
419+
rFindable.fill(HIST("hFilterIndex"), filterIndex);
420+
rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius);
421+
rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius);
422+
} else {
423+
continue;
424+
}
425+
426+
// 4 - daughter track ITS+TPC properties
427+
if (daughterTrack.itsNClsInnerBarrel() == 0 && daughterTrack.itsNCls() < 4 &&
428+
daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > 80) {
429+
filterIndex += 1;
430+
rFindable.fill(HIST("hFilterIndex"), filterIndex);
431+
rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius);
432+
rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius);
433+
} else {
434+
continue;
435+
}
436+
437+
// 5 - geometric cuts: eta
438+
if (std::abs(motherTrack.eta()) < 1.0 && std::abs(daughterTrack.eta()) < 1.0) {
439+
filterIndex += 1;
440+
rFindable.fill(HIST("hFilterIndex"), filterIndex);
441+
rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius);
442+
rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius);
443+
} else {
444+
continue;
403445
}
446+
447+
// 6 - geometric cuts: phi difference
448+
if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < 50.0) {
449+
filterIndex += 1;
450+
rFindable.fill(HIST("hFilterIndex"), filterIndex);
451+
rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius);
452+
rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius);
453+
} else {
454+
continue;
455+
}
456+
457+
// 7 - radius cut
458+
if (recRadius > 19.6213f) {
459+
filterIndex += 1;
460+
rFindable.fill(HIST("hFilterIndex"), filterIndex);
461+
rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius);
462+
rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius);
463+
} else {
464+
continue;
465+
}
466+
467+
// 8 - collision selection
468+
auto collision = motherTrack.template collision_as<CollisionsFullMC>();
469+
if (!(std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) {
470+
filterIndex += 1;
471+
rFindable.fill(HIST("hFilterIndex"), filterIndex);
472+
rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius);
473+
rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius);
474+
}
475+
404476
}
405477
}
406478

@@ -414,10 +486,10 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
414486
}
415487

416488
// Next steps:
417-
// 0. Resolution histograms should have relative values, not absolute OK
418-
// 1. New h2 with genRadius (recRadius) vs FilterIndex
419-
// 2. Get recRadius through a map on kinkCands, put a negative value if the candidate is not reconstructed
420-
// 2.1 Consider adding step in filters with the cuts on radius
421-
// 2.2 Add h2 of radius resolution vs pt
422-
// 3. Rewrite the findable method using maps to avoid the nested loop
489+
// 0. Resolution histograms should have relative values, not absolute: OK
490+
// 1. New h2 with genRadius (recRadius) vs FilterIndex: OK
491+
// 2. Get recRadius through a map on kinkCands, put a negative value if the candidate is not reconstructed: OK
492+
// 2.1 Consider adding step in filters with the cuts on radius: OK
493+
// 2.2 Add h2 of radius resolution vs pt: OK
494+
// 3. Rewrite the findable method using maps to avoid the nested loop: OK
423495
// 4. For generated h2, avoid mass axis, use a bool axis to easily distinguish sigma minus and sigma plus

0 commit comments

Comments
 (0)