Skip to content

catch also tracks of particles of other MC collisions #1343

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
jgrosseo merged 1 commit into from
Oct 10, 2022
Merged

catch also tracks of particles of other MC collisions #1343

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
98 changes: 56 additions & 42 deletions PWGCF/TableProducer/filterCorrelations.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ struct FilterCF {
O2_DEFINE_CONFIGURABLE(cfgCutMCPt, float, 0.5f, "Minimal pT for particles")
O2_DEFINE_CONFIGURABLE(cfgCutMCEta, float, 0.8f, "Eta range for particles")
O2_DEFINE_CONFIGURABLE(cfgVerbosity, int, 1, "Verbosity level (0 = major, 1 = per collision)")
O2_DEFINE_CONFIGURABLE(cfgTrigger, int, 7, "Trigger choice: (0 = none, 7 = sel7, 8 = sel8)")

// Filters and input definitions
Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex;
Expand All @@ -74,8 +75,14 @@ struct FilterCF {
template <typename TCollision>
bool keepCollision(TCollision& collision)
{
// TODO needs Run 3 adjustement
return collision.alias()[kINT7] && collision.sel7();
if (cfgTrigger == 0) {
return true;
} else if (cfgTrigger == 7) {
return collision.alias()[kINT7] && collision.sel7();
} else if (cfgTrigger == 8) {
return collision.sel8();
}
return false;
}

void processData(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::CFMultiplicities>>::iterator const& collision, aod::BCsWithTimestamps const&, soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection>> const& tracks)
Expand Down Expand Up @@ -189,27 +196,25 @@ struct FilterCF {
}
PROCESS_SWITCH(FilterCF, processMC2, "Process MC: MC part", false);

void processMC(soa::Filtered<aod::McCollisions>::iterator const& mcCollision, aod::McParticles const& particles,
soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSels, aod::CFMultiplicities>> const& collisions,
// NOTE not filtering collisions here because in that case there can be tracks referring to MC particles which are not part of the selected MC collisions
Preslice<aod::McParticles> perMcCollision = aod::mcparticle::mcCollisionId;
void processMC(aod::McCollisions const& mcCollisions, aod::McParticles const& allParticles,
soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSels, aod::CFMultiplicities> const& allCollisions,
soa::Filtered<soa::Join<aod::Tracks, aod::McTrackLabels, aod::TrackSelection>> const& tracks,
aod::BCsWithTimestamps const&)
{
if (cfgVerbosity > 0) {
LOGF(info, "processMC: Particles for MC collision: %d | Vertex: %.1f", particles.size(), mcCollision.posZ());
}

bool* reconstructed = new bool[particles.size()];
int* mcParticleLabels = new int[particles.size()];
for (int i = 0; i < particles.size(); i++) {
bool* reconstructed = new bool[allParticles.size()];
int* mcParticleLabels = new int[allParticles.size()];
for (int i = 0; i < allParticles.size(); i++) {
reconstructed[i] = false;
mcParticleLabels[i] = -1;
}

// PASS 1 on collisions: check which particles are kept
for (auto& collision : collisions) {
for (auto& collision : allCollisions) {
auto groupedTracks = tracks.sliceBy(perCollision, collision.globalIndex());
if (cfgVerbosity > 0) {
LOGF(info, "processMC: Tracks for collision: %d | Vertex: %.1f (%d) | INT7: %d", groupedTracks.size(), collision.posZ(), collision.flags(), collision.sel7());
LOGF(info, "processMC: Tracks for collision %d: %d | Vertex: %.1f (%d) | INT7: %d", collision.globalIndex(), groupedTracks.size(), collision.posZ(), collision.flags(), collision.sel7());
}

if (!keepCollision(collision)) {
Expand All @@ -218,55 +223,64 @@ struct FilterCF {

for (auto& track : groupedTracks) {
if (track.has_mcParticle()) {
reconstructed[track.mcParticleId() - particles.begin().globalIndex()] = true;
reconstructed[track.mcParticleId()] = true;
}
}
}

// Store selected MC particles and MC collisions
int multiplicity = 0;
for (auto& particle : particles) {
int8_t sign = 0;
TParticlePDG* pdgparticle = pdg->GetParticle(particle.pdgCode());
if (pdgparticle != nullptr) {
sign = (pdgparticle->Charge() > 0) ? 1.0 : ((pdgparticle->Charge() < 0) ? -1.0 : 0.0);
}
bool primary = particle.isPhysicalPrimary() && sign != 0 && std::abs(particle.eta()) < cfgCutMCEta && particle.pt() > cfgCutMCPt;
if (primary) {
multiplicity++;
for (auto& mcCollision : mcCollisions) {
auto particles = allParticles.sliceBy(perMcCollision, mcCollision.globalIndex());

if (cfgVerbosity > 0) {
LOGF(info, "processMC: Particles for MC collision %d: %d | Vertex: %.1f", mcCollision.globalIndex(), particles.size(), mcCollision.posZ());
}
if (reconstructed[particle.index()] || primary) {
// keep particle

// use highest bit to flag if it is reconstructed
uint8_t flags = particle.flags() & ~aod::cfmcparticle::kReconstructed; // clear bit in case of clashes in the future
if (reconstructed[particle.index()]) {
flags |= aod::cfmcparticle::kReconstructed;
// Store selected MC particles and MC collisions
int multiplicity = 0;
for (auto& particle : particles) {
int8_t sign = 0;
TParticlePDG* pdgparticle = pdg->GetParticle(particle.pdgCode());
if (pdgparticle != nullptr) {
sign = (pdgparticle->Charge() > 0) ? 1.0 : ((pdgparticle->Charge() < 0) ? -1.0 : 0.0);
}
bool primary = particle.isPhysicalPrimary() && sign != 0 && std::abs(particle.eta()) < cfgCutMCEta && particle.pt() > cfgCutMCPt;
if (primary) {
multiplicity++;
}
if (reconstructed[particle.globalIndex()] || primary) {
// keep particle

// use highest bit to flag if it is reconstructed
uint8_t flags = particle.flags() & ~aod::cfmcparticle::kReconstructed; // clear bit in case of clashes in the future
if (reconstructed[particle.globalIndex()]) {
flags |= aod::cfmcparticle::kReconstructed;
}

// NOTE using "outputMcCollisions.lastIndex()+1" here to allow filling of outputMcCollisions *after* the loop
outputMcParticles(outputMcCollisions.lastIndex() + 1, truncateFloatFraction(particle.pt(), FLOAT_PRECISION), truncateFloatFraction(particle.eta(), FLOAT_PRECISION),
truncateFloatFraction(particle.phi(), FLOAT_PRECISION), sign, particle.pdgCode(), flags);
// NOTE using "outputMcCollisions.lastIndex()+1" here to allow filling of outputMcCollisions *after* the loop
outputMcParticles(outputMcCollisions.lastIndex() + 1, truncateFloatFraction(particle.pt(), FLOAT_PRECISION), truncateFloatFraction(particle.eta(), FLOAT_PRECISION),
truncateFloatFraction(particle.phi(), FLOAT_PRECISION), sign, particle.pdgCode(), flags);

// relabeling array
mcParticleLabels[particle.index()] = outputMcParticles.lastIndex();
// relabeling array
mcParticleLabels[particle.globalIndex()] = outputMcParticles.lastIndex();
}
}
outputMcCollisions(mcCollision.posZ(), multiplicity);
}
outputMcCollisions(mcCollision.posZ(), multiplicity);

// PASS 2 on collisions: store collisions and tracks
for (auto& collision : collisions) {
for (auto& collision : allCollisions) {
auto groupedTracks = tracks.sliceBy(perCollision, collision.globalIndex());
if (cfgVerbosity > 0) {
LOGF(info, "processMC: Tracks for collision: %d | Vertex: %.1f (%d) | INT7: %d", groupedTracks.size(), collision.posZ(), collision.flags(), collision.sel7());
LOGF(info, "processMC: Tracks for collision %d: %d | Vertex: %.1f (%d) | INT7: %d", collision.globalIndex(), groupedTracks.size(), collision.posZ(), collision.flags(), collision.sel7());
}

if (!keepCollision(collision)) {
continue;
}

auto bc = collision.bc_as<aod::BCsWithTimestamps>();
outputCollisions(outputMcCollisions.lastIndex(), bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp());
// NOTE works only when we store all MC collisions (as we do here)
outputCollisions(collision.mcCollisionId(), bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp());

for (auto& track : groupedTracks) {
uint8_t trackType = 0;
Expand All @@ -278,9 +292,9 @@ struct FilterCF {

int mcParticleId = track.mcParticleId();
if (mcParticleId >= 0) {
mcParticleId = mcParticleLabels[track.mcParticleId() - particles.begin().globalIndex()];
mcParticleId = mcParticleLabels[track.mcParticleId()];
if (mcParticleId < 0) {
LOGP(fatal, "processMC: Track {} is referring to a MC particle which we do not store {} {}", track.index(), track.mcParticleId(), mcParticleId);
LOGP(fatal, "processMC: Track {} is referring to a MC particle which we do not store {} {} (reco flag {})", track.index(), track.mcParticleId(), mcParticleId, reconstructed[track.mcParticleId()]);
}
}
outputTracks(outputCollisions.lastIndex(), mcParticleId,
Expand Down