diff --git a/PWGCF/TableProducer/filterCorrelations.cxx b/PWGCF/TableProducer/filterCorrelations.cxx index 92bb3c619c4..b5c1c37f950 100644 --- a/PWGCF/TableProducer/filterCorrelations.cxx +++ b/PWGCF/TableProducer/filterCorrelations.cxx @@ -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; @@ -74,8 +75,14 @@ struct FilterCF { template 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>::iterator const& collision, aod::BCsWithTimestamps const&, soa::Filtered> const& tracks) @@ -189,27 +196,25 @@ struct FilterCF { } PROCESS_SWITCH(FilterCF, processMC2, "Process MC: MC part", false); - void processMC(soa::Filtered::iterator const& mcCollision, aod::McParticles const& particles, - soa::SmallGroups> 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 perMcCollision = aod::mcparticle::mcCollisionId; + void processMC(aod::McCollisions const& mcCollisions, aod::McParticles const& allParticles, + soa::Join const& allCollisions, soa::Filtered> 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)) { @@ -218,47 +223,55 @@ 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)) { @@ -266,7 +279,8 @@ struct FilterCF { } auto bc = collision.bc_as(); - 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; @@ -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,