Skip to content

Commit 04d537f

Browse files
a-m-andrushkoa-m-andrushko
committed
Add new MCTruth analysis
1 parent 5c1113a commit 04d537f

File tree

2 files changed

+246
-142
lines changed

2 files changed

+246
-142
lines changed

‎PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx‎

Lines changed: 154 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -74,24 +74,16 @@ using namespace o2::constants::physics;
7474
namespace o2::aod
7575
{
7676

77-
using FemtoFullCollision =
78-
soa::Join<aod::Collisions, aod::EvSels, aod::Mults>::iterator;
79-
using FemtoFullCollisionCentPP =
80-
soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Ms, aod::Mults>::iterator;
81-
using FemtoFullCollisionCentRun2 =
82-
soa::Join<aod::Collisions, aod::EvSels, aod::CentRun2V0Ms, aod::Mults>::iterator;
83-
using FemtoFullCollisionCentRun3 =
84-
soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::Mults>::iterator;
77+
using FemtoFullCollision = soa::Join<aod::Collisions, aod::EvSels, aod::Mults>::iterator;
78+
using FemtoFullCollisionCentPP = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Ms, aod::Mults>::iterator;
79+
using FemtoFullCollisionCentRun2 = soa::Join<aod::Collisions, aod::EvSels, aod::CentRun2V0Ms, aod::Mults>::iterator;
80+
using FemtoFullCollisionCentRun3 = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::Mults>::iterator;
8581
using FemtoFullCollisionMC = soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::McCollisionLabels>::iterator;
86-
using FemtoFullCollisionCentRun3MCs =
87-
soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::Mults, aod::McCollisionLabels>;
88-
using FemtoFullCollisionCentRun3MC =
89-
soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::Mults, aod::McCollisionLabels>::iterator;
90-
using FemtoFullTracks =
91-
soa::Join<aod::FullTracks, aod::TracksDCA, aod::TOFSignal, aod::pidTPCEl, aod::TrackSelection,
92-
aod::pidTPCMu, aod::pidTPCPi, aod::pidTPCKa, aod::pidTPCPr,
93-
aod::pidTPCDe, aod::pidTOFEl, aod::pidTOFMu, aod::pidTOFPi,
94-
aod::pidTOFKa, aod::pidTOFPr, aod::pidTOFDe>;
82+
using FemtoFullCollisionCentRun3MCs = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::Mults, aod::McCollisionLabels>;
83+
using FemtoFullCollisionCentRun3MC = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::Mults, aod::McCollisionLabels>::iterator;
84+
using FemtoFullTracks = soa::Join<aod::FullTracks, aod::TracksDCA, aod::TOFSignal, aod::TrackSelection,
85+
aod::pidTPCEl, aod::pidTPCMu, aod::pidTPCPi, aod::pidTPCKa, aod::pidTPCPr, aod::pidTPCDe,
86+
aod::pidTOFEl, aod::pidTOFMu, aod::pidTOFPi, aod::pidTOFKa, aod::pidTOFPr, aod::pidTOFDe>;
9587

9688
// using FilteredFullV0s = soa::Filtered<aod::V0Datas>; /// predefined Join
9789
// table for o2::aod::V0s = soa::Join<o2::aod::TransientV0s, o2::aod::StoredV0s>
@@ -232,6 +224,8 @@ struct FemtoUniverseProducerTask {
232224
Configurable<bool> confV0RejectKaons{"confV0RejectKaons", false, "Switch to reject kaons"};
233225
Configurable<float> confV0InvKaonMassLowLimit{"confV0InvKaonMassLowLimit", 0.48, "Lower limit of the V0 invariant mass for Kaon rejection"};
234226
Configurable<float> confV0InvKaonMassUpLimit{"confV0InvKaonMassUpLimit", 0.515, "Upper limit of the V0 invariant mass for Kaon rejection"};
227+
228+
Configurable<std::vector<int>> confV0PDGMCTruth{"confV0PDGMCTruth", std::vector<int>{2212, -211, 3122}, "PDG codes of V0 daughters and mother, the order must be as follows -- positive daughter, negative daughter, mother"};
235229
} ConfV0Selection;
236230

237231
struct : o2::framework::ConfigurableGroup {
@@ -875,6 +869,37 @@ struct FemtoUniverseProducerTask {
875869
}
876870
}
877871

872+
template <typename MCParticleType>
873+
void fillMCTruthParticle(MCParticleType const& mcparticle, o2::aod::femtouniverseparticle::ParticleType fdparttype)
874+
{
875+
auto pdgCode = mcparticle.pdgCode();
876+
int particleOrigin = 99;
877+
878+
// Determine particle origin
879+
if (std::abs(pdgCode) == std::abs(confPDGCodePartOne.value) || std::abs(pdgCode) == std::abs(confPDGCodePartTwo.value)) {
880+
if (mcparticle.isPhysicalPrimary()) {
881+
particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kPrimary;
882+
} else {
883+
auto mothers = mcparticle.template mothers_as<aod::McParticles>();
884+
if (!mothers.empty()) {
885+
auto mother = mothers.front();
886+
if (mother.producedByGenerator()) {
887+
particleOrigin = checkDaughterType(fdparttype, mother.pdgCode());
888+
} else {
889+
particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kMaterial;
890+
}
891+
}
892+
}
893+
} else {
894+
particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kFake;
895+
}
896+
897+
// Fill MC companion tables
898+
outputPartsMC(particleOrigin, pdgCode, mcparticle.pt(), mcparticle.eta(), mcparticle.phi());
899+
fillDebugParticleMC(mcparticle);
900+
outputPartsMCLabels(outputPartsMC.lastIndex());
901+
}
902+
878903
template <typename ParticleType>
879904
void fillMCParticlePhi(ParticleType const& kaon1, ParticleType const& kaon2)
880905
{
@@ -1244,6 +1269,77 @@ struct FemtoUniverseProducerTask {
12441269
}
12451270
}
12461271

1272+
template <typename MCParticlesType>
1273+
void fillV0MCTruth(MCParticlesType const& mcParticles)
1274+
{
1275+
for (const auto& mc : mcParticles) { // Loop over all MC Truth particles
1276+
if (mc.pdgCode() != ConfV0Selection.confV0PDGMCTruth.value[2])
1277+
continue; // Artificially single out V0s
1278+
1279+
auto daughters = mc.template daughters_as<aod::McParticles>(); // Access daughters (no differentiation of signs, it needs to be checked separately)
1280+
1281+
bool foundPos = false, foundNeg = false;
1282+
1283+
std::vector<int> childIDs = {0, 0};
1284+
int rowPos = 0;
1285+
int rowNeg = 0;
1286+
1287+
for (auto const& d : daughters) { // Loop over daughters
1288+
if (d.pdgCode() == ConfV0Selection.confV0PDGMCTruth.value[0]) { // Check for a positive child
1289+
foundPos = true;
1290+
1291+
outputParts(outputCollision.lastIndex(),
1292+
d.pt(),
1293+
d.eta(),
1294+
d.phi(),
1295+
aod::femtouniverseparticle::ParticleType::kMCTruthTrack,
1296+
0,
1297+
0,
1298+
d.pdgCode(),
1299+
childIDs, // {0, 0}
1300+
0,
1301+
0);
1302+
rowPos = outputParts.lastIndex();
1303+
fillMCTruthParticle(d, aod::femtouniverseparticle::ParticleType::kMCTruthTrack);
1304+
} else if (d.pdgCode() == ConfV0Selection.confV0PDGMCTruth.value[1]) { // Check for a negative child
1305+
foundNeg = true;
1306+
1307+
outputParts(outputCollision.lastIndex(),
1308+
d.pt(),
1309+
d.eta(),
1310+
d.phi(),
1311+
aod::femtouniverseparticle::ParticleType::kMCTruthTrack,
1312+
0,
1313+
0,
1314+
d.pdgCode(),
1315+
childIDs, // {0, 0}
1316+
0,
1317+
0);
1318+
rowNeg = outputParts.lastIndex();
1319+
fillMCTruthParticle(d, aod::femtouniverseparticle::ParticleType::kMCTruthTrack);
1320+
}
1321+
}
1322+
1323+
if (!foundPos || !foundNeg)
1324+
continue;
1325+
1326+
childIDs[0] = rowPos;
1327+
childIDs[1] = rowNeg;
1328+
outputParts(outputCollision.lastIndex(),
1329+
mc.pt(),
1330+
mc.eta(),
1331+
mc.phi(),
1332+
aod::femtouniverseparticle::ParticleType::kMCTruthTrack,
1333+
0,
1334+
0,
1335+
mc.pdgCode(),
1336+
childIDs,
1337+
0,
1338+
0);
1339+
fillMCTruthParticle(mc, aod::femtouniverseparticle::ParticleType::kMCTruthTrack);
1340+
}
1341+
}
1342+
12471343
template <bool isMC, typename CollisionType, typename CascadeType, typename TrackType>
12481344
void fillCascade(CollisionType const& col, CascadeType const& fullCascades, TrackType const&)
12491345
{
@@ -2429,47 +2525,6 @@ struct FemtoUniverseProducerTask {
24292525
}
24302526
PROCESS_SWITCH(FemtoUniverseProducerTask, processTruthAndFullMCV0, "Provide both MC truth and reco for tracks and V0s", false);
24312527

2432-
void processTruthAndFullMCCentRun3V0(
2433-
aod::McCollisions const& mccols,
2434-
aod::McParticles const& mcParticles,
2435-
aod::FemtoFullCollisionCentRun3MCs const& collisions,
2436-
soa::Filtered<soa::Join<aod::FemtoFullTracks, aod::McTrackLabels>> const& tracks,
2437-
soa::Join<o2::aod::V0Datas, aod::McV0Labels> const& fullV0s,
2438-
aod::BCsWithTimestamps const&)
2439-
{
2440-
2441-
// recos
2442-
std::set<int> recoMcIds;
2443-
for (const auto& col : collisions) {
2444-
auto groupedTracks = tracks.sliceBy(perCollisionTracks, col.globalIndex());
2445-
auto groupedV0Parts = fullV0s.sliceBy(perCollisionV0s, col.globalIndex());
2446-
getMagneticFieldTesla(col.bc_as<aod::BCsWithTimestamps>());
2447-
const auto colcheck = fillCollisionsCentRun3<false>(col);
2448-
if (colcheck) {
2449-
fillTracks<true>(groupedTracks);
2450-
fillV0<true>(col, groupedV0Parts, groupedTracks);
2451-
}
2452-
for (const auto& track : groupedTracks) {
2453-
if (trackCuts.isSelectedMinimal(track))
2454-
recoMcIds.insert(track.mcParticleId());
2455-
}
2456-
}
2457-
2458-
// truth
2459-
for (const auto& mccol : mccols) {
2460-
auto groupedCollisions = collisions.sliceBy(recoCollsPerMCCollCentPbPb, mccol.globalIndex());
2461-
for (const auto& col : groupedCollisions) {
2462-
const auto colcheck = fillMCTruthCollisionsCentRun3(col); // fills the reco collisions for mc collision
2463-
if (colcheck) {
2464-
auto groupedMCParticles = mcParticles.sliceBy(perMCCollision, mccol.globalIndex());
2465-
outputCollExtra(1.0, 1.0);
2466-
fillParticles<decltype(groupedMCParticles), true, true>(groupedMCParticles, recoMcIds); // fills mc particles
2467-
}
2468-
}
2469-
}
2470-
}
2471-
PROCESS_SWITCH(FemtoUniverseProducerTask, processTruthAndFullMCCentRun3V0, "Provide both MC truth and reco for tracks and V0s with centrality", false);
2472-
24732528
Preslice<soa::Join<o2::aod::CascDatas, aod::McCascLabels>> perCollisionCascs = aod::track::collisionId;
24742529
void processTruthAndFullMCCasc(
24752530
aod::McCollisions const& mccols,
@@ -2673,6 +2728,47 @@ struct FemtoUniverseProducerTask {
26732728
}
26742729
PROCESS_SWITCH(FemtoUniverseProducerTask, processV0CentRun3Data, "Provide experimental data for Run 3 with centrality for track track", false);
26752730

2731+
void processTruthAndFullMCCentRun3V0(
2732+
aod::McCollisions const& mccols,
2733+
aod::McParticles const& mcParticles,
2734+
aod::FemtoFullCollisionCentRun3MCs const& collisions,
2735+
soa::Filtered<soa::Join<aod::FemtoFullTracks, aod::McTrackLabels>> const& tracks,
2736+
soa::Join<o2::aod::V0Datas, aod::McV0Labels> const& fullV0s,
2737+
aod::BCsWithTimestamps const&)
2738+
{
2739+
2740+
// MCReco
2741+
std::set<int> recoMcIds;
2742+
for (const auto& col : collisions) { // loop over collisions
2743+
auto groupedTracks = tracks.sliceBy(perCollisionTracks, col.globalIndex()); // slicing for tracks
2744+
auto groupedV0Parts = fullV0s.sliceBy(perCollisionV0s, col.globalIndex()); // slicing for V0
2745+
getMagneticFieldTesla(col.bc_as<aod::BCsWithTimestamps>());
2746+
const auto colcheck = fillCollisionsCentRun3<false>(col);
2747+
if (colcheck) {
2748+
fillTracks<true>(groupedTracks);
2749+
fillV0<true>(col, groupedV0Parts, groupedTracks);
2750+
}
2751+
for (const auto& track : groupedTracks) {
2752+
if (trackCuts.isSelectedMinimal(track))
2753+
recoMcIds.insert(track.mcParticleId());
2754+
}
2755+
}
2756+
2757+
// MCTruth
2758+
for (const auto& mccol : mccols) {
2759+
auto groupedCollisions = collisions.sliceBy(recoCollsPerMCCollCentPbPb, mccol.globalIndex()); // slicing for MC collisions
2760+
auto groupedMCParticles = mcParticles.sliceBy(perMCCollision, mccol.globalIndex()); // slicing for MC particles
2761+
for (const auto& col : groupedCollisions) {
2762+
const auto colcheck = fillMCTruthCollisionsCentRun3(col);
2763+
if (colcheck) {
2764+
outputCollExtra(1.0, 1.0);
2765+
fillV0MCTruth(groupedMCParticles); // fills MC V0s and its daughters
2766+
}
2767+
}
2768+
}
2769+
}
2770+
PROCESS_SWITCH(FemtoUniverseProducerTask, processTruthAndFullMCCentRun3V0, "Provide both MC truth and reco for tracks and V0s with centrality", false);
2771+
26762772
void processCascadeCentRun3Data(aod::FemtoFullCollisionCentRun3 const& col,
26772773
aod::BCsWithTimestamps const&,
26782774
soa::Filtered<aod::FemtoFullTracks> const& tracks,

0 commit comments

Comments
 (0)