Skip to content

[PWGUD] Added flags and occupancy in MC table #18760

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
ekryshen merged 2 commits into from
Nov 15, 2026
+111 −113
Merged

[PWGUD] Added flags and occupancy in MC table #18760

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
33f3a50
Added flags and occupancy in MC table
ekryshen Nov 15, 2026
7ff91b4
Megalinter fix: renamed flags to bcFlags
ekryshen Nov 15, 2026
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
224 changes: 111 additions & 113 deletions PWGUD/TableProducer/upcCandProducerBarrel.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -209,6 +209,8 @@ DECLARE_SOA_TABLE(UPCBarrelCands, "AOD", "UPCBARRELCANDS",
DECLARE_SOA_TABLE(UPCBarrelMcColls, "AOD", "UPCBMCCOLL",
upc_cand_prod_bar::TfId,
upc_cand_prod_bar::GlobalBC,
upc_cand_prod_bar::Flags,
upc_cand_prod_bar::OccupancyFT0C,
upc_cand_prod_bar::PosX,
upc_cand_prod_bar::PosY,
upc_cand_prod_bar::PosZ,
Expand Down Expand Up @@ -343,8 +345,6 @@ struct UpcCandProducerBarrel {
registry.add("hTF", "TF counts", HistType::kTH1D, {axisTFs});
registry.add("hTVX", "TVX counts", HistType::kTH2D, {axisBCs, axisTFs});
registry.add("hTVXRCT", "TVX counts after RCT selection", HistType::kTH2D, {axisBCs, axisTFs});
registry.add("hTVXupc", "TVX counts in upc ROFs", HistType::kTH2D, {axisBCs, axisTFs});
registry.add("hTVXRCTupc", "TVX counts after RCT selection in upc ROFs", HistType::kTH2D, {axisBCs, axisTFs});

registry.add("hVetoT00", "no RCT selection", HistType::kTH3D, {axisBCs, axisTFs, axisVeto});
registry.add("hVetoTV0", "no RCT selection", HistType::kTH3D, {axisBCs, axisTFs, axisVeto});
Expand Down Expand Up @@ -429,81 +429,37 @@ struct UpcCandProducerBarrel {
}
LOGP(info, "Starting new DF. Number of TFs = {}", tfIDs.size());

std::unordered_map<int64_t, SignalRef> signalParticleRefs;
if constexpr (isMC) {
std::unordered_set<int64_t> ambiguousTrackIds;
ambiguousTrackIds.reserve(ambiguousTracks->size());
for (const auto& ambiguousTrack : *ambiguousTracks) {
ambiguousTrackIds.insert(ambiguousTrack.trackId());
}
std::vector<int64_t> bcsWithFT0, bcsWithFV0, bcsWithFDD;
bcsWithFT0.reserve(ft0s.size());
bcsWithFV0.reserve(fv0s.size());
bcsWithFDD.reserve(fdds.size());

std::vector<int64_t> bcsWithMCCollision;
bcsWithMCCollision.reserve(mcCollisions->size());
// Use collisions from every source and retain same-BC pile-up
for (const auto& mcCollision : *mcCollisions) {
bcsWithMCCollision.push_back(static_cast<int64_t>(mcCollision.template bc_as<BCsWithSels>().globalBC()));
LOGP(debug, "Filling FT0 activity arrays");
for (const auto& ft0 : ft0s) {
const float sumAmpFT0A = ft0.sumAmpA();
const float sumAmpFT0C = ft0.sumAmpC();
if ((ft0.timeA() < MaxFITTime && sumAmpFT0A >= minAmpFT0A) || (ft0.timeC() < MaxFITTime && sumAmpFT0C >= minAmpFT0C)) {
bcsWithFT0.push_back(ft0.bc_as<BCsWithSels>().globalBC());
}
const auto minDistanceMC = buildMinimumDistanceMap(tfIDs, bcSOR, nBCsPerTF, bcsWithMCCollision, maxStoredDistance, true);

for (const auto& mcCollision : *mcCollisions) {
const auto bc = mcCollision.template bc_as<BCsWithSels>();
const auto globalBC = bc.globalBC();
const auto iTF = localTFIndex[bc.globalIndex()];
const auto bcInTF = (static_cast<int64_t>(globalBC) - bcSOR) % nBCsPerTF;
const auto particles = mcParticles->sliceBy(particlesPerMCCollision, mcCollision.globalIndex());
const bool isSignal = mcCollision.getSourceId() == signalSourceId.value;
}

std::vector<float> px, py, pz;
std::vector<int32_t> pdg, recoCollisionIds;
std::vector<uint8_t> recoFlags;
int32_t mult = 0;
// Output rows follow input MC collisions, making this the compact collision label
const auto newMCCollisionId = static_cast<int32_t>(selectedMCCollisions.lastIndex() + 1);
for (const auto& particle : particles) {
mult += particle.isPhysicalPrimary();
// Background rows keep only collision metadata and physical-primary multiplicity
if (!isSignal) {
continue;
}
const auto newMCParticleId = static_cast<int16_t>(pdg.size());
px.push_back(particle.px());
py.push_back(particle.py());
pz.push_back(particle.pz());
pdg.push_back(particle.pdgCode());
LOGP(debug, "Filling FV0 activity arrays");
for (const auto& fv0 : fv0s) {
const float sumAmpFV0A = std::accumulate(fv0.amplitude().begin(), fv0.amplitude().end(), 0.f);
if (fv0.time() < MaxFITTime && sumAmpFV0A >= minAmpFV0A) {
bcsWithFV0.push_back(fv0.bc_as<BCsWithSels>().globalBC());
}
}

uint8_t flags = particle.isPhysicalPrimary() ? kPhysicalPrimary : kNoFlags;
int32_t recoCollisionId = NoRecoCollision;
int nRecoTracks = 0;
// Inspect only reconstructed tracks matched to this signal particle
for (const auto& track : tracks.sliceBy(tracksPerMCParticle, particle.globalIndex())) {
if ((track.mcMask() & MCTrackFakeMask) != 0u) {
continue;
}
++nRecoTracks;
flags |= kHasRecoTrack;
flags |= track.isPVContributor() ? kHasPVContributor : kNoFlags;
flags |= ambiguousTrackIds.contains(track.globalIndex()) ? kHasAmbiguousTrack : kNoFlags;
if (track.has_collision()) {
recoCollisionId = track.collisionId();
flags |= collisions.iteratorAt(recoCollisionId).numContrib() > nTracks ? kHighMultiplicityCollision : kNoFlags;
}
}
// A single collision label is meaningful only for exactly one reconstructed track
if (nRecoTracks > 1) {
recoCollisionId = MultipleRecoTracks;
}
recoFlags.push_back(flags);
recoCollisionIds.push_back(recoCollisionId);
// Remap original particle IDs to output index for candidate labels
signalParticleRefs.emplace(particle.globalIndex(), SignalRef{newMCCollisionId, newMCParticleId});
}
selectedMCCollisions(tfIDs[iTF], globalBC, mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(),
px, py, pz, pdg, mult, minDistanceMC[iTF][bcInTF],
recoFlags, recoCollisionIds, isSignal);
LOGP(debug, "Filling FDD activity arrays");
for (const auto& fdd : fdds) {
const float sumAmpFDDA = std::accumulate(std::begin(fdd.chargeA()), std::end(fdd.chargeA()), 0.f);
const float sumAmpFDDC = std::accumulate(std::begin(fdd.chargeC()), std::end(fdd.chargeC()), 0.f);
if ((fdd.timeA() < MaxFITTime && sumAmpFDDA >= minAmpFDDA) || (fdd.timeC() < MaxFITTime && sumAmpFDDC >= minAmpFDDC)) {
bcsWithFDD.push_back(fdd.bc_as<BCsWithSels>().globalBC());
}
}

// collect number of all collisions and upc-mode collisions per ROF
LOGP(debug, "Collect number of all collisions and upc-mode collisions per ROF");
std::vector<std::map<int32_t, float>> vMultTracks(tfIDs.size(), std::map<int32_t, float>());
std::vector<std::vector<int>> vCollisionsPerROF(tfIDs.size(), std::vector<int>(nRofPerTF, 0));
Expand All @@ -524,16 +480,7 @@ struct UpcCandProducerBarrel {
auto hTVX = registry.get<TH2>(HIST("hTVX"));
auto hTVXRCT = registry.get<TH2>(HIST("hTVXRCT"));

auto hTVXupc = registry.get<TH2>(HIST("hTVXupc"));
auto hTVXRCTupc = registry.get<TH2>(HIST("hTVXRCTupc"));

std::vector<int64_t> bcsWithFT0, bcsWithFV0, bcsWithFDD;
bcsWithFT0.reserve(ft0s.size());
bcsWithFV0.reserve(fv0s.size());
bcsWithFDD.reserve(fdds.size());

LOGP(debug, "Starting lumi accounting");
// lumi accounting
std::vector<std::map<int32_t, float>> vMultFT0C(tfIDs.size(), std::map<int32_t, float>());
for (const auto& ft0 : ft0s) {
const auto bc = ft0.bc_as<BCsWithSels>();
Expand All @@ -542,8 +489,6 @@ struct UpcCandProducerBarrel {
const int32_t bcInTF = (globalBC - bcSOR) % nBCsPerTF;
const auto iTF = localTFIndex[bc.globalIndex()];
const auto tfID = tfIDs[iTF];
const int32_t rofInTF = getRofInTF(bcInTF);
bool upcMode = vCollisionsPerROF[iTF][rofInTF] == 0 || vUpcModePerROF[iTF][rofInTF] > 0;

// Multiplicity for FT0C-based occupancy estimates
if (ft0.timeC() < MaxFITTime && ft0.sumAmpC() >= minAmpFT0C) {
Expand All @@ -559,44 +504,12 @@ struct UpcCandProducerBarrel {
}

hTVX->Fill(bcInOrbit, tfID);
if (upcMode) {
hTVXupc->Fill(bcInOrbit, tfID);
}

if (!bc.selection_bit(aod::evsel::kNoTimeFrameBorder) || !rctChecker(bc)) {
continue;
}

hTVXRCT->Fill(bcInOrbit, tfID);
if (upcMode) {
hTVXRCTupc->Fill(bcInOrbit, tfID);
}
}

LOGP(debug, "Filling FT0 activity arrays");
for (const auto& ft0 : ft0s) {
const float sumAmpFT0A = ft0.sumAmpA();
const float sumAmpFT0C = ft0.sumAmpC();
if ((ft0.timeA() < MaxFITTime && sumAmpFT0A >= minAmpFT0A) || (ft0.timeC() < MaxFITTime && sumAmpFT0C >= minAmpFT0C)) {
bcsWithFT0.push_back(ft0.bc_as<BCsWithSels>().globalBC());
}
}

LOGP(debug, "Filling FV0 activity arrays");
for (const auto& fv0 : fv0s) {
const float sumAmpFV0A = std::accumulate(fv0.amplitude().begin(), fv0.amplitude().end(), 0.f);
if (fv0.time() < MaxFITTime && sumAmpFV0A >= minAmpFV0A) {
bcsWithFV0.push_back(fv0.bc_as<BCsWithSels>().globalBC());
}
}

LOGP(debug, "Filling FDD activity arrays");
for (const auto& fdd : fdds) {
const float sumAmpFDDA = std::accumulate(std::begin(fdd.chargeA()), std::end(fdd.chargeA()), 0.f);
const float sumAmpFDDC = std::accumulate(std::begin(fdd.chargeC()), std::end(fdd.chargeC()), 0.f);
if ((fdd.timeA() < MaxFITTime && sumAmpFDDA >= minAmpFDDA) || (fdd.timeC() < MaxFITTime && sumAmpFDDC >= minAmpFDDC)) {
bcsWithFDD.push_back(fdd.bc_as<BCsWithSels>().globalBC());
}
}

LOGP(debug, "Starting occupancy calculation");
Expand All @@ -620,6 +533,91 @@ struct UpcCandProducerBarrel {
}
}

std::unordered_map<int64_t, SignalRef> signalParticleRefs;
if constexpr (isMC) {
LOGP(debug, "Collecting MC info. Number of MC collisions = {}", mcCollisions->size());
std::unordered_set<int64_t> ambiguousTrackIds;
ambiguousTrackIds.reserve(ambiguousTracks->size());
for (const auto& ambiguousTrack : *ambiguousTracks) {
ambiguousTrackIds.insert(ambiguousTrack.trackId());
}

std::vector<int64_t> bcsWithMCCollision;
bcsWithMCCollision.reserve(mcCollisions->size());
// Use collisions from every source and retain same-BC pile-up
for (const auto& mcCollision : *mcCollisions) {
bcsWithMCCollision.push_back(static_cast<int64_t>(mcCollision.template bc_as<BCsWithSels>().globalBC()));
}
const auto minDistanceMC = buildMinimumDistanceMap(tfIDs, bcSOR, nBCsPerTF, bcsWithMCCollision, maxStoredDistance, true);

for (const auto& mcCollision : *mcCollisions) {
const auto bc = mcCollision.template bc_as<BCsWithSels>();
const int64_t globalBC = bc.globalBC();
const auto iTF = localTFIndex[bc.globalIndex()];
const int64_t bcInTF = (globalBC - bcSOR) % nBCsPerTF;
LOGP(debug, "setting flags for MC collision");
const int32_t rofInTF = getRofInTF(bcInTF);
const bool upcMode = vCollisionsPerROF[iTF][rofInTF] == 0 || vUpcModePerROF[iTF][rofInTF] > 0;
const bool noPileupInRof = vCollisionsPerROF[iTF][rofInTF] <= 1;
uint8_t bcFlags = 0;
bcFlags |= upcMode << 0;
bcFlags |= noPileupInRof << 1;
const auto occupancyFT0C = vOccupancyFT0C[iTF][bcInTF];

const auto particles = mcParticles->sliceBy(particlesPerMCCollision, mcCollision.globalIndex());
const bool isSignal = mcCollision.getSourceId() == signalSourceId.value;
LOGP(debug, "sourceId = {}", mcCollision.getSourceId());
std::vector<float> px, py, pz;
std::vector<int32_t> pdg, recoCollisionIds;
std::vector<uint8_t> recoFlags;
int32_t mult = 0;
// Output rows follow input MC collisions, making this the compact collision label
const auto newMCCollisionId = static_cast<int32_t>(selectedMCCollisions.lastIndex() + 1);
for (const auto& particle : particles) {
mult += particle.isPhysicalPrimary();
// Background rows keep only collision metadata and physical-primary multiplicity
if (!isSignal) {
continue;
}
const auto newMCParticleId = static_cast<int16_t>(pdg.size());
px.push_back(particle.px());
py.push_back(particle.py());
pz.push_back(particle.pz());
pdg.push_back(particle.pdgCode());

uint8_t flags = particle.isPhysicalPrimary() ? kPhysicalPrimary : kNoFlags;
int32_t recoCollisionId = NoRecoCollision;
int nRecoTracks = 0;
// Inspect only reconstructed tracks matched to this signal particle
for (const auto& track : tracks.sliceBy(tracksPerMCParticle, particle.globalIndex())) {
if ((track.mcMask() & MCTrackFakeMask) != 0u) {
continue;
}
++nRecoTracks;
flags |= kHasRecoTrack;
flags |= track.isPVContributor() ? kHasPVContributor : kNoFlags;
flags |= ambiguousTrackIds.contains(track.globalIndex()) ? kHasAmbiguousTrack : kNoFlags;
if (track.has_collision()) {
recoCollisionId = track.collisionId();
flags |= collisions.iteratorAt(recoCollisionId).numContrib() > nTracks ? kHighMultiplicityCollision : kNoFlags;
}
}
// A single collision label is meaningful only for exactly one reconstructed track
if (nRecoTracks > 1) {
recoCollisionId = MultipleRecoTracks;
}
recoFlags.push_back(flags);
recoCollisionIds.push_back(recoCollisionId);
// Remap original particle IDs to output index for candidate labels
signalParticleRefs.emplace(particle.globalIndex(), SignalRef{newMCCollisionId, newMCParticleId});
}
selectedMCCollisions(tfIDs[iTF], globalBC, bcFlags, occupancyFT0C,
mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(),
px, py, pz, pdg, mult, minDistanceMC[iTF][bcInTF],
recoFlags, recoCollisionIds, isSignal);
}
}

LOGP(debug, "Building minimum distance maps");
const uint32_t scanDistance = std::max(maxStoredDistance.value, vetoBCWindow.value);
const DistanceMap minDistanceFT0 = buildMinimumDistanceMap(tfIDs, bcSOR, nBCsPerTF, bcsWithFT0, scanDistance);
Expand Down
Loading