diff --git a/PWGCF/DataModel/CorrelationsDerived.h b/PWGCF/DataModel/CorrelationsDerived.h index ec0f2aa11c2..1cf6ec4292f 100644 --- a/PWGCF/DataModel/CorrelationsDerived.h +++ b/PWGCF/DataModel/CorrelationsDerived.h @@ -40,7 +40,7 @@ using CFMcCollisionWithExtra = CFMcCollisionsWithExtra::iterator; namespace cfmcparticle { -DECLARE_SOA_INDEX_COLUMN(CFMcCollision, cfMcCollision); //! Index to reduced MC collision +DECLARE_SOA_INDEX_COLUMN(CFMcCollision, cfMcCollision); //! Index to reduced MC collision; o2-linter: disable=name/o2-column (preserve the established derived-table API) DECLARE_SOA_COLUMN(Pt, pt, float); //! pT (GeV/c) DECLARE_SOA_COLUMN(Eta, eta, float); //! Pseudorapidity DECLARE_SOA_COLUMN(Phi, phi, float); //! Phi angle @@ -58,14 +58,26 @@ using CFMcParticle = CFMcParticles::iterator; namespace cfmultiplicity { DECLARE_SOA_COLUMN(Multiplicity, multiplicity, float); -} -DECLARE_SOA_TABLE(CFMultiplicities, "AOD", "CFMULTIPLICITY", cfmultiplicity::Multiplicity); +DECLARE_SOA_COLUMN(MultiplicityEstimator, multiplicityEstimator, uint8_t); //! Source used for the multiplicity value +enum EstimatorType : uint8_t { + Tracks, + FT0M, + FT0C, + FT0CVariant1, + FT0CVariant2, + FT0A, + CentNGlobal, + Run2V0M, + MCParticles, +}; +} // namespace cfmultiplicity +DECLARE_SOA_TABLE(CFMultiplicities, "AOD", "CFMULTIPLICITY", cfmultiplicity::Multiplicity, cfmultiplicity::MultiplicityEstimator); using CFMultiplicity = CFMultiplicities::iterator; namespace cfcollision { -DECLARE_SOA_INDEX_COLUMN(CFMcCollision, cfMcCollision); //! Index to reduced MC collision +DECLARE_SOA_INDEX_COLUMN(CFMcCollision, cfMcCollision); //! Index to reduced MC collision; o2-linter: disable=name/o2-column (preserve the established derived-table API) DECLARE_SOA_COLUMN(Multiplicity, multiplicity, float); //! Centrality/multiplicity value } // namespace cfcollision DECLARE_SOA_TABLE(CFCollisions, "AOD", "CFCOLLISION", //! Reduced collision table @@ -79,10 +91,20 @@ using CFCollLabel = CFCollLabels::iterator; using CFCollisionsWithLabel = soa::Join; using CFCollisionWithLabel = CFCollisionsWithLabel::iterator; +namespace cfcollisionextra +{ +DECLARE_SOA_COLUMN(MultiplicityCorrected, multiplicityCorrected, float); //! Efficiency-corrected track count +} // namespace cfcollisionextra +DECLARE_SOA_TABLE(CFCollisionsExtra, "AOD", "CFCOLLSEXTRA", //! Row-aligned extension of CFCollisions; filled only when multiplicity efficiency is configured + cfcollisionextra::MultiplicityCorrected); +using CFCollisionExtra = CFCollisionsExtra::iterator; +using CFCollisionsWithExtra = soa::Join; +using CFCollisionWithExtra = CFCollisionsWithExtra::iterator; + namespace cftrack { -DECLARE_SOA_INDEX_COLUMN(CFCollision, cfCollision); //! Index to collision -DECLARE_SOA_INDEX_COLUMN(CFMcParticle, cfMCParticle); //! Index to MC particle +DECLARE_SOA_INDEX_COLUMN(CFCollision, cfCollision); //! Index to collision; o2-linter: disable=name/o2-column (preserve the established derived-table API) +DECLARE_SOA_INDEX_COLUMN(CFMcParticle, cfMCParticle); //! Index to MC particle; o2-linter: disable=name/o2-column (preserve the established derived-table API) DECLARE_SOA_COLUMN(Pt, pt, float); //! pT (GeV/c) DECLARE_SOA_COLUMN(Eta, eta, float); //! Pseudorapidity DECLARE_SOA_COLUMN(Phi, phi, float); //! Phi angle @@ -147,8 +169,8 @@ using CFMcParticleRef = CFMcParticleRefs::iterator; namespace cf2prongtrack { -DECLARE_SOA_INDEX_COLUMN_FULL(CFTrackProng0, cfTrackProng0, int, CFTracks, "_0"); //! Index to prong 1 CFTrack -DECLARE_SOA_INDEX_COLUMN_FULL(CFTrackProng1, cfTrackProng1, int, CFTracks, "_1"); //! Index to prong 2 CFTrack +DECLARE_SOA_INDEX_COLUMN_FULL(CFTrackProng0, cfTrackProng0, int, CFTracks, "_0"); //! Index to prong 1 CFTrack; o2-linter: disable=name/o2-column (preserve the established derived-table API) +DECLARE_SOA_INDEX_COLUMN_FULL(CFTrackProng1, cfTrackProng1, int, CFTracks, "_1"); //! Index to prong 2 CFTrack; o2-linter: disable=name/o2-column (preserve the established derived-table API) DECLARE_SOA_COLUMN(Pt, pt, float); //! pT (GeV/c) DECLARE_SOA_COLUMN(Eta, eta, float); //! Pseudorapidity DECLARE_SOA_COLUMN(Phi, phi, float); //! Phi angle @@ -201,8 +223,8 @@ using CF2ProngTrackml = CF2ProngTrackmls::iterator; namespace cf2prongmcpart { -DECLARE_SOA_INDEX_COLUMN_FULL(CFParticleDaugh0, cfParticleDaugh0, int, CFMcParticles, "_0"); //! Index to prong 1 CFMcParticle -DECLARE_SOA_INDEX_COLUMN_FULL(CFParticleDaugh1, cfParticleDaugh1, int, CFMcParticles, "_1"); //! Index to prong 2 CFMcParticle +DECLARE_SOA_INDEX_COLUMN_FULL(CFParticleDaugh0, cfParticleDaugh0, int, CFMcParticles, "_0"); //! Index to prong 1 CFMcParticle; o2-linter: disable=name/o2-column (preserve the established derived-table API) +DECLARE_SOA_INDEX_COLUMN_FULL(CFParticleDaugh1, cfParticleDaugh1, int, CFMcParticles, "_1"); //! Index to prong 2 CFMcParticle; o2-linter: disable=name/o2-column (preserve the established derived-table API) DECLARE_SOA_COLUMN(Decay, decay, uint8_t); //! Particle decay and flags DECLARE_SOA_DYNAMIC_COLUMN(McDecay, mcDecay, [](uint8_t decay) -> uint8_t { return decay & 0x3f; }); //! MC particle decay enum ParticleDecayFlags { diff --git a/PWGCF/TableProducer/filterCorrelations.cxx b/PWGCF/TableProducer/filterCorrelations.cxx index 3063ff72605..2f79e838d54 100644 --- a/PWGCF/TableProducer/filterCorrelations.cxx +++ b/PWGCF/TableProducer/filterCorrelations.cxx @@ -9,6 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +// o2-linter: disable=name/workflow-file (historic file contains several table-producer tasks) #include "PWGCF/DataModel/CorrelationsDerived.h" #include "Common/CCDB/EventSelectionParams.h" @@ -21,6 +22,7 @@ #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" +#include #include #include #include @@ -35,14 +37,21 @@ #include #include +#include #include +#include #include #include #include +#include +#include #include #include // required for is_detected +#include +#include +#include #include #include @@ -58,6 +67,7 @@ using namespace o2::math_utils::detail; struct FilterCF { Service pdg; + Service ccdb; enum TrackSelectionCuts1 : uint8_t { kTrackSelected = BIT(0), @@ -103,6 +113,10 @@ struct FilterCF { O2_DEFINE_CONFIGURABLE(chi2peritscluster, float, 36, "maximum Chi2 / cluster for the ITS track segment") O2_DEFINE_CONFIGURABLE(cfgEstimatorBitMask, uint16_t, 0, "BitMask for multiplicity estimators to be included in the CFMultSet tables."); + O2_DEFINE_CONFIGURABLE(cfgEfficiencyMultiplicity, std::string, "", "Multiplicity efficiency (RecoAll / MC): CCDB path or local ROOT file with a 4D ccdb_object (eta, pT, multiplicity, z-vtx); empty disables CFCollisionsExtra output") + O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, int, 0, "0 = CCDB efficiency, 1 = local ROOT efficiency") + O2_DEFINE_CONFIGURABLE(cfgMultiplicityTrackBitMask, uint16_t, 0, "Required track-type bits for corrected multiplicity; match cfgTrackBitMask used to produce the efficiency (0 = all stored tracks)") + // Filters and input definitions Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex; Filter collisionVertexTypeFilter = (cfgCollisionFlags == 0) || ((aod::collision::flags & cfgCollisionFlags) == cfgCollisionFlags); @@ -119,6 +133,7 @@ struct FilterCF { HistogramRegistry registrytrackQA{"TrackQA", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; Produces outputCollisions; + Produces outputCollisionsExtra; Produces outputTracks; Produces outputMcCollisionLabels; @@ -135,12 +150,42 @@ struct FilterCF { Produces outputMultSets; std::vector multiplicities{}; + // Own local histograms independently of their input file. CCDB owns its objects. + std::unique_ptr localMultiplicityEfficiency; + static constexpr int MultiplicityEfficiencyDimensions = 4; + // persistent caches std::vector mcReconstructedCache; std::vector mcParticleLabelsCache; void init(InitContext&) { + if (!cfgEfficiencyMultiplicity.value.empty()) { + if (cfgLocalEfficiency != 0 && cfgLocalEfficiency != 1) { + LOGF(fatal, "cfgLocalEfficiency must be 0 (CCDB) or 1 (local ROOT file)"); + } + if (cfgMultiplicityTrackBitMask > std::numeric_limits::max()) { + LOGF(fatal, "cfgMultiplicityTrackBitMask must fit the 8-bit track type"); + } + if (cfgLocalEfficiency == 1) { + std::unique_ptr file(TFile::Open(cfgEfficiencyMultiplicity.value.c_str(), "READ")); + if (!file) { + LOGF(fatal, "Could not open multiplicity efficiency file %s", cfgEfficiencyMultiplicity.value.c_str()); + return; + } + if (file->IsZombie()) { + LOGF(fatal, "Multiplicity efficiency file %s is invalid", cfgEfficiencyMultiplicity.value.c_str()); + return; + } + auto* efficiency = dynamic_cast(file->Get("ccdb_object")); + validateMultiplicityEfficiency(efficiency); + localMultiplicityEfficiency.reset(static_cast(efficiency->Clone())); + } else { + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + } + } if (doprocessTrackQA) { registrytrackQA.add("zvtx", "Z Vertex position; posz (cm); Events", HistType::kTH1F, {{100, -12, 12}}); registrytrackQA.add("eta", "eta distribution; eta; arb. units", HistType::kTH1F, {{100, -2, 2}}); @@ -156,7 +201,7 @@ struct FilterCF { } template - bool keepCollision(TCollision& collision) + bool keepCollision(const TCollision& collision) { bool isMultSelected = false; if (collision.multiplicity() >= cfgMinMultiplicity) @@ -164,23 +209,23 @@ struct FilterCF { if (cfgTrigger == 0) { return true; - } else if (cfgTrigger == 7) { + } else if (cfgTrigger == 7) { // o2-linter: disable=magic-number (documented legacy trigger-selection code) return isMultSelected && collision.alias_bit(kINT7) && collision.sel7(); - } else if (cfgTrigger == 8) { + } else if (cfgTrigger == 8) { // o2-linter: disable=magic-number (documented legacy trigger-selection code) return isMultSelected && collision.sel8(); - } else if (cfgTrigger == 9) { // relevant only for Pb-Pb + } else if (cfgTrigger == 9) { // relevant only for Pb-Pb; o2-linter: disable=magic-number (documented legacy trigger-selection code) return isMultSelected && collision.sel8() && collision.selection_bit(aod::evsel::kNoSameBunchPileup) && collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) && collision.selection_bit(aod::evsel::kIsGoodITSLayersAll); - } else if (cfgTrigger == 10) { // TVX trigger only (sel8 selection before April, 2024) + } else if (cfgTrigger == 10) { // TVX trigger only (sel8 selection before April, 2024); o2-linter: disable=magic-number (documented legacy trigger-selection code) return isMultSelected && collision.selection_bit(aod::evsel::kIsTriggerTVX); - } else if (cfgTrigger == 11) { // sel8 selection for MC + } else if (cfgTrigger == 11) { // sel8 selection for MC; o2-linter: disable=magic-number (documented legacy trigger-selection code) return isMultSelected && collision.selection_bit(aod::evsel::kIsTriggerTVX) && collision.selection_bit(aod::evsel::kNoTimeFrameBorder); - } else if (cfgTrigger == 12) { // relevant only for Pb-Pb with occupancy cuts and rejection of the collisions which have other events nearby + } else if (cfgTrigger == 12) { // relevant only for Pb-Pb with occupancy cuts and rejection of nearby collisions; o2-linter: disable=magic-number (documented legacy trigger-selection code) int occupancy = collision.trackOccupancyInTimeRange(); if (occupancy >= cfgMinOcc && occupancy < cfgMaxOcc) return isMultSelected && collision.sel8() && collision.selection_bit(aod::evsel::kNoSameBunchPileup) && collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) && collision.selection_bit(aod::evsel::kNoCollInTimeRangeStandard) && collision.selection_bit(aod::evsel::kIsGoodITSLayersAll); else return false; - } else if (cfgTrigger == 13) { // relevant for pO/OO/NeNe --recommended by Physics Board on 27.01.2026 + } else if (cfgTrigger == 13) { // relevant for pO/OO/NeNe, recommended by Physics Board on 27.01.2026; o2-linter: disable=magic-number (documented legacy trigger-selection code) return isMultSelected && collision.sel8() && collision.selection_bit(aod::evsel::kNoSameBunchPileup) && collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV); } return false; @@ -193,18 +238,18 @@ struct FilterCF { { o2::aod::ITSResponse itsResponse; - if (ITSProtonselection && candidate.pt() <= 0.6 && !(itsResponse.nSigmaITS(candidate) > nsigmaCutITSProton)) { + if (ITSProtonselection && candidate.pt() <= 0.6 && !(itsResponse.nSigmaITS(candidate) > nsigmaCutITSProton)) { // o2-linter: disable=magic-number (established proton PID momentum boundary) return false; } - if (ITSProtonselection && candidate.pt() > 0.6 && candidate.pt() <= 0.8 && !(itsResponse.nSigmaITS(candidate) > nsigmaCutITSProton)) { + if (ITSProtonselection && candidate.pt() > 0.6 && candidate.pt() <= 0.8 && !(itsResponse.nSigmaITS(candidate) > nsigmaCutITSProton)) { // o2-linter: disable=magic-number (established proton PID momentum boundaries) return false; } if (candidate.hasTOF()) { - if (candidate.pt() < 0.7 && std::abs(candidate.tpcNSigmaPr()) < nsigmaCutTPCProton) { + if (candidate.pt() < 0.7 && std::abs(candidate.tpcNSigmaPr()) < nsigmaCutTPCProton) { // o2-linter: disable=magic-number (established proton PID momentum boundary) return true; } - if (candidate.p() >= 0.7 && std::abs(candidate.tpcNSigmaPr()) < nsigmaCutTPCProton && std::abs(candidate.tofNSigmaPr()) < nsigmaCutTOFProton) { + if (candidate.p() >= 0.7 && std::abs(candidate.tpcNSigmaPr()) < nsigmaCutTPCProton && std::abs(candidate.tofNSigmaPr()) < nsigmaCutTOFProton) { // o2-linter: disable=magic-number (established proton PID momentum boundary) return true; } } else { @@ -258,7 +303,7 @@ struct FilterCF { } } return trackType; - } else if (cfgTrackSelection == 2) { + } else if (cfgTrackSelection == 2) { // o2-linter: disable=magic-number (documented track-selection mode) uint8_t trackType = 0; if constexpr (HasProtonPID::value) { if (track.isGlobalTrack() && (track.itsNCls() >= itsnclusters) && (track.tpcNClsCrossedRows() >= tpcncrossedrows) && selectionPIDProton(track)) { @@ -280,6 +325,67 @@ struct FilterCF { return dcaXyConst + dcaXySlope / pt; // a + b/pT } + void validateMultiplicityEfficiency(const THn* efficiency) const + { + if (!efficiency || efficiency->GetNdimensions() != MultiplicityEfficiencyDimensions) { + LOGF(fatal, "Multiplicity efficiency from %s must be a 4D THn with axes (eta, pT, multiplicity, z-vtx)", cfgEfficiencyMultiplicity.value.c_str()); + } + } + + THn* loadMultiplicityEfficiency(uint64_t timestamp) + { + if (cfgLocalEfficiency == 1) { + return localMultiplicityEfficiency.get(); + } + // Query each collision so the manager can refresh its cache at validity boundaries. + auto* efficiency = ccdb->getForTimeStamp>(cfgEfficiencyMultiplicity.value, timestamp); + validateMultiplicityEfficiency(efficiency); + return efficiency; + } + + template + float getCorrectedMultiplicity(const TCollision& collision, const TTracks& tracks, uint64_t timestamp) + { + if (collision.multiplicityEstimator() != aod::cfmultiplicity::Tracks) { + LOGF(fatal, "Efficiency-corrected multiplicity requires MultiplicitySelector::processTracks, but estimator type %u was configured", static_cast(collision.multiplicityEstimator())); + } + auto* efficiency = loadMultiplicityEfficiency(timestamp); + double correctedMultiplicity = 0.; + for (const auto& track : tracks) { + // Match the tracks written by the corresponding data/MC producer path. + if constexpr (applyDCA) { + if (std::abs(track.dcaXY()) > getMaxDCAxy(track.pt()) || std::abs(track.dcaZ()) > dcazmax) { + continue; + } + } + const auto mask = static_cast(cfgMultiplicityTrackBitMask.value); + if (mask != 0 && (getTrackType(track) & mask) != mask) { + continue; + } + + // The map contains RecoAll / MC, not inverse-efficiency weights. + // Keep the original estimator as the map coordinate, including for centrality. + const std::array values{track.eta(), track.pt(), collision.multiplicity(), collision.posZ()}; + std::array bins{}; + for (int axis = 0; axis < MultiplicityEfficiencyDimensions; ++axis) { + auto* efficiencyAxis = efficiency->GetAxis(axis); + bins[axis] = efficiencyAxis->FindFixBin(values[axis]); + if (!std::isfinite(values[axis]) || bins[axis] < 1 || bins[axis] > efficiencyAxis->GetNbins()) { + LOGF(fatal, "Multiplicity efficiency from %s does not cover axis %d value %g", cfgEfficiencyMultiplicity.value.c_str(), axis, values[axis]); + } + } + const double eff = efficiency->GetBinContent(bins.data()); + if (!std::isfinite(eff) || eff <= 0.) { + LOGF(fatal, "Invalid multiplicity efficiency %g from %s at bins (%d, %d, %d, %d)", eff, cfgEfficiencyMultiplicity.value.c_str(), bins[0], bins[1], bins[2], bins[3]); + } + correctedMultiplicity += 1. / eff; + } + if (!std::isfinite(correctedMultiplicity) || correctedMultiplicity > std::numeric_limits::max()) { + LOGF(fatal, "Corrected multiplicity cannot be represented as a float: %g", correctedMultiplicity); + } + return static_cast(correctedMultiplicity); + } + template using HasMultTables = decltype(std::declval().multNTracksPV()); @@ -299,6 +405,9 @@ struct FilterCF { auto bc = collision.template bc_as(); outputCollisions(bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp()); + if (!cfgEfficiencyMultiplicity.value.empty()) { + outputCollisionsExtra(getCorrectedMultiplicity(collision, tracks, bc.timestamp())); + } if constexpr (std::experimental::is_detected::value) { multiplicities.clear(); @@ -317,7 +426,7 @@ struct FilterCF { if (cfgTransientTables) outputCollRefs(collision.globalIndex()); - for (auto& track : tracks) { + for (const auto& track : tracks) { float maxDCAxy = getMaxDCAxy(track.pt()); if ((std::abs(track.dcaXY()) > maxDCAxy) || (std::abs(track.dcaZ()) > dcazmax)) { continue; @@ -402,7 +511,7 @@ struct FilterCF { } // PASS 1 on collisions: check which particles are kept - for (auto& collision : allCollisions) { + for (const auto& collision : allCollisions) { auto groupedTracks = tracks.sliceBy(perCollision, collision.globalIndex()); if (cfgVerbosity > 0) { LOGF(info, "processMC: Tracks for collision %d: %d | Vertex: %.1f (%d) | INT7: %d", collision.globalIndex(), groupedTracks.size(), collision.posZ(), collision.flags(), collision.sel7()); @@ -412,14 +521,14 @@ struct FilterCF { continue; } - for (auto& track : groupedTracks) { + for (const auto& track : groupedTracks) { if (track.has_mcParticle()) { mcReconstructedCache[track.mcParticleId()] = true; } } } - for (auto& mcCollision : mcCollisions) { + for (const auto& mcCollision : mcCollisions) { auto particles = allParticles.sliceBy(perMcCollision, mcCollision.globalIndex()); if (cfgVerbosity > 0) { @@ -428,7 +537,7 @@ struct FilterCF { // Store selected MC particles and MC collisions int multiplicity = 0; - for (auto& particle : particles) { + for (const auto& particle : particles) { int8_t sign = 0; TParticlePDG* pdgparticle = pdg->GetParticle(particle.pdgCode()); if (pdgparticle != nullptr) { @@ -464,7 +573,7 @@ struct FilterCF { } // PASS 2 on collisions: store collisions and tracks - for (auto& collision : allCollisions) { + for (const auto& collision : allCollisions) { auto groupedTracks = tracks.sliceBy(perCollision, collision.globalIndex()); if (cfgVerbosity > 0) { LOGF(info, "processMC: Tracks for collision %d: %d | Vertex: %.1f (%d) | INT7: %d", collision.globalIndex(), groupedTracks.size(), collision.posZ(), collision.flags(), collision.sel7()); @@ -477,6 +586,9 @@ struct FilterCF { auto bc = collision.template bc_as(); // NOTE works only when we store all MC collisions (as we do here) outputCollisions(bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp()); + if (!cfgEfficiencyMultiplicity.value.empty()) { + outputCollisionsExtra(getCorrectedMultiplicity(collision, groupedTracks, bc.timestamp())); + } outputMcCollisionLabels(collision.mcCollisionId()); if constexpr (std::experimental::is_detected::value) { @@ -497,7 +609,7 @@ struct FilterCF { if (cfgTransientTables) outputCollRefs(collision.globalIndex()); - for (auto& track : groupedTracks) { + for (const auto& track : groupedTracks) { int mcParticleId = track.mcParticleId(); if (mcParticleId >= 0) { mcParticleId = mcParticleLabelsCache[track.mcParticleId()]; @@ -552,7 +664,7 @@ struct FilterCF { void processMCGen(McCollisionsWithHepMC::iterator const& mcCollision, aod::McParticles const& particles) { float multiplicity = 0.0f; - for (auto& particle : particles) { + for (const auto& particle : particles) { if (!particle.isPhysicalPrimary() || std::abs(particle.eta()) > cfgCutMCEta || particle.pt() < cfgCutMCPt) continue; int8_t sign = 0; @@ -617,69 +729,69 @@ struct MultiplicitySelector { void processTracks(aod::Collision const&, soa::Filtered> const& tracks) { - output(tracks.size()); + output(tracks.size(), aod::cfmultiplicity::Tracks); } PROCESS_SWITCH(MultiplicitySelector, processTracks, "Select track count as multiplicity", false); void processFT0M(aod::CentFT0Ms const& centralities) { - for (auto& c : centralities) { - output(c.centFT0M()); + for (const auto& c : centralities) { + output(c.centFT0M(), aod::cfmultiplicity::FT0M); } } PROCESS_SWITCH(MultiplicitySelector, processFT0M, "Select FT0M centrality as multiplicity", false); void processFT0C(aod::CentFT0Cs const& centralities) { - for (auto& c : centralities) { - output(c.centFT0C()); + for (const auto& c : centralities) { + output(c.centFT0C(), aod::cfmultiplicity::FT0C); } } PROCESS_SWITCH(MultiplicitySelector, processFT0C, "Select FT0C centrality as multiplicity", false); void processFT0CVariant1(aod::CentFT0CVariant1s const& centralities) { - for (auto& c : centralities) { - output(c.centFT0CVariant1()); + for (const auto& c : centralities) { + output(c.centFT0CVariant1(), aod::cfmultiplicity::FT0CVariant1); } } PROCESS_SWITCH(MultiplicitySelector, processFT0CVariant1, "Select FT0CVariant1 centrality as multiplicity", false); void processFT0CVariant2(aod::CentFT0CVariant2s const& centralities) { - for (auto& c : centralities) { - output(c.centFT0CVariant2()); + for (const auto& c : centralities) { + output(c.centFT0CVariant2(), aod::cfmultiplicity::FT0CVariant2); } } PROCESS_SWITCH(MultiplicitySelector, processFT0CVariant2, "Select FT0CVariant2 centrality as multiplicity", false); void processFT0A(aod::CentFT0As const& centralities) { - for (auto& c : centralities) { - output(c.centFT0A()); + for (const auto& c : centralities) { + output(c.centFT0A(), aod::cfmultiplicity::FT0A); } } PROCESS_SWITCH(MultiplicitySelector, processFT0A, "Select FT0A centrality as multiplicity", false); void processCentNGlobal(aod::CentNGlobals const& centralities) { - for (auto& c : centralities) { - output(c.centNGlobal()); + for (const auto& c : centralities) { + output(c.centNGlobal(), aod::cfmultiplicity::CentNGlobal); } } PROCESS_SWITCH(MultiplicitySelector, processCentNGlobal, "Select CentNGlobal centrality as multiplicity", false); void processRun2V0M(aod::CentRun2V0Ms const& centralities) { - for (auto& c : centralities) { - output(c.centRun2V0M()); + for (const auto& c : centralities) { + output(c.centRun2V0M(), aod::cfmultiplicity::Run2V0M); } } PROCESS_SWITCH(MultiplicitySelector, processRun2V0M, "Select V0M centrality as multiplicity", true); void processMCGen(aod::McCollision const&, aod::McParticles const& particles) { - output(particles.size()); + output(particles.size(), aod::cfmultiplicity::MCParticles); } PROCESS_SWITCH(MultiplicitySelector, processMCGen, "Select MC particle count as multiplicity", false); }; diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index 908ce848d7d..441f25dc85e 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -70,6 +70,7 @@ #include #include #include +#include #include #include @@ -99,7 +100,7 @@ struct TwoParticleCorrelationsMpi { ; Configurable cfgLocalEfficiency{"cfgLocalEfficiency", 0, "0 = OFF and 1 = ON for local efficiency"}; Configurable cfgDropStepRECO{"cfgDropStepRECO", false, "choice to drop step RECO if efficiency correction is used"}; - Configurable cfgCentBinsForMC{"cfgCentBinsForMC", 0, "0 = OFF and 1 = ON for data like multiplicity/centrality bins for MC steps"}; + Configurable cfgCentBinsForMC{"cfgCentBinsForMC", 0, "0 = generated multiplicity; 1 = reconstructed multiplicity and all associated collisions; 2 = reconstructed multiplicity and first associated collision only in processMCEfficiency"}; Configurable cfgTrackBitMask{"cfgTrackBitMask", 0, "BitMask for track selection systematics; refer to the enum TrackSelectionCuts in filtering task"}; Configurable cfgMultCorrelationsMask{"cfgMultCorrelationsMask", 0, "Selection bitmask for the multiplicity correlations. This should match the filter selection cfgEstimatorBitMask."}; Configurable cfgMultCutFormula{"cfgMultCutFormula", "", "Multiplicity correlations cut formula. A result greater than zero results in accepted event. Parameters: [cFT0C] FT0C centrality, [mFV0A] V0A multiplicity, [mGlob] global track multiplicity, [mPV] PV track multiplicity, [cFT0M] FT0M centrality"}; @@ -284,14 +285,27 @@ struct TwoParticleCorrelationsMpi { using AodTracks = soa::Filtered>; using DerivedCollisions = soa::Filtered; + using DerivedCollisionsCorrected = soa::Filtered; using DerivedTracks = soa::Filtered; + static constexpr int McEfficiencySingleRecoCollisionMode = 2; void init(o2::framework::InitContext&) { if (cfgUserAxis < NoUserAxis || cfgUserAxis > EventSeedAxis) { LOGF(fatal, "Unsupported cfgUserAxis=%d; use 0 (off), 1 (invariant mass), or 2 (event seed)", cfgUserAxis.value); } - if (doprocessMCSameDerived && (doprocessSameDerived || doprocessSameDerivedMultSet)) { + if (cfgCentBinsForMC < 0 || cfgCentBinsForMC > McEfficiencySingleRecoCollisionMode) { + LOGF(fatal, "Unsupported cfgCentBinsForMC=%d; use 0 (generated multiplicity), 1 (all reconstructed collisions), or 2 (first reconstructed collision only for efficiency)", cfgCentBinsForMC.value); + } + const int enabledDerivedSameProcesses = doprocessSameDerived + doprocessSameDerivedCorrected + doprocessSameDerivedMultSet + doprocessSameDerivedMultSetCorrected; + if (enabledDerivedSameProcesses > 1) { + LOGF(fatal, "Only one reconstructed derived same-event process can be enabled"); + } + const int enabledDerivedMixedProcesses = doprocessMixedDerived + doprocessMixedDerivedCorrected + doprocessMixedDerivedMultSet + doprocessMixedDerivedMultSetCorrected; + if (enabledDerivedMixedProcesses > 1) { + LOGF(fatal, "Only one reconstructed derived mixed-event process can be enabled"); + } + if (doprocessMCSameDerived && enabledDerivedSameProcesses > 0) { LOGF(fatal, "processMCSameDerived is mutually exclusive with the reconstructed derived same-event processes because it also fills those outputs"); } if (doprocessSameGenMC && doprocessMCSameDerived) { @@ -347,7 +361,7 @@ struct TwoParticleCorrelationsMpi { registry.add("yields", "multiplicity/centrality vs pT vs eta", {HistType::kTH3F, {{100, 0, 100, "/multiplicity/centrality"}, {40, 0, 20, "p_{T}"}, {100, -2, 2, "#eta"}}}); registry.add("etaphi", "multiplicity/centrality vs eta vs phi", {HistType::kTH3F, {{100, 0, 100, "multiplicity/centrality"}, {100, -2, 2, "#eta"}, {200, 0, o2::constants::math::TwoPI, "#varphi"}}}); - if (doprocessSameDerivedMultSet) { + if (doprocessSameDerivedMultSet || doprocessSameDerivedMultSetCorrected) { if (cfgMultCorrelationsMask == 0) { LOGF(fatal, "cfgMultCorrelationsMask can not be 0 when MultSet process functions are in use."); } @@ -370,6 +384,16 @@ struct TwoParticleCorrelationsMpi { registry.add("multCorrelations", "Multiplicity correlations", {HistType::kTHnSparseF, multAxes}); } registry.add("multiplicity", "event multiplicity", {HistType::kTH1F, {{1000, 0, 100, "/multiplicity/centrality"}}}); + if (doprocessMCEfficiency) { + registry.add("mcEfficiencyDiagnostics/reconstructedTracks", "selected reconstructed tracks per collision;N_{tracks};collisions", {HistType::kTH1F, {{1001, -0.5, 1000.5}}}); + registry.add("mcEfficiencyDiagnostics/matchedTracks", "MC-matched reconstructed tracks per collision;N_{matched};collisions", {HistType::kTH1F, {{1001, -0.5, 1000.5}}}); + registry.add("mcEfficiencyDiagnostics/uniqueMatchedParticles", "unique matched MC particles per collision;N_{unique labels};collisions", {HistType::kTH1F, {{1001, -0.5, 1000.5}}}); + registry.add("mcEfficiencyDiagnostics/duplicateMatchedTracks", "matched tracks beyond one per MC label;N_{matched}-N_{unique labels};collisions", {HistType::kTH1F, {{501, -0.5, 500.5}}}); + registry.add("mcEfficiencyDiagnostics/duplicateFraction", "fraction of matched tracks sharing an MC label;(N_{matched}-N_{unique labels})/N_{matched};collisions", {HistType::kTH1F, {{101, -0.005, 1.005}}}); + registry.add("mcEfficiencyDiagnostics/tracksPerMcParticle", "reconstructed tracks per matched MC-particle label;tracks per MC label;MC labels", {HistType::kTH1F, {{21, -0.5, 20.5}}}); + registry.add("mcEfficiencyDiagnostics/generatedVsUniqueMatchedPrimaries", "generated versus uniquely matched physical primaries;N_{generated primary};N_{unique matched primary}", {HistType::kTH2F, {{501, -0.5, 500.5}, {501, -0.5, 500.5}}}); + registry.add("mcEfficiencyDiagnostics/primaryTracksPerMcParticle", "reconstructed tracks per matched physical-primary label;tracks per primary MC label;MC labels", {HistType::kTH1F, {{21, -0.5, 20.5}}}); + } if (eventSeedEstimatorEnabled) { registry.add("eventSeedEstimator", "event-level template estimator", {HistType::kTHnSparseF, {{100, 0, 100, "multiplicity"}, {100, -0.5, 99.5, "N_{trig}"}, {200, 0, 20, "Y_{near}"}, {200, 0, 20, "Y_{away}"}, {200, 0, 100, "N_{uncorrelated seeds}"}}}); registry.add("eventSeedPairProbabilities", "summed pair probabilities", {HistType::kTH3F, {{200, 0, 200, "#Sigma P_{baseline}"}, {200, 0, 200, "#Sigma P_{near}"}, {200, 0, 200, "#Sigma P_{away}"}}}); @@ -555,6 +579,18 @@ struct TwoParticleCorrelationsMpi { template using HasMultSet = decltype(std::declval().multiplicities()); + template + using HasCorrectedMultiplicity = decltype(std::declval().multiplicityCorrected()); + + template + static float getAnalysisMultiplicity(const TCollision& collision) + { + if constexpr (std::experimental::is_detected::value) { + return collision.multiplicityCorrected(); + } + return collision.multiplicity(); + } + template void fillQA(const TCollision& collision, float multiplicity, const TTracks& tracks) { @@ -1670,22 +1706,22 @@ struct TwoParticleCorrelationsMpi { template void processSameDerivedT(CollType const& collision, TTracks1 const& tracks1, TTracks2 const& tracks2, const int* trueNMPI = nullptr) { - using BinningTypeDerived = ColumnBinningPolicy; - BinningTypeDerived configurableBinningDerived{{axisVertex, axisMultiplicity}, true}; // true is for 'ignore overflows' (true by default). Underflows and overflows will have bin -1. + auto getMultiplicity = [](const auto& col) { return getAnalysisMultiplicity(col); }; + using BinningTypeDerived = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getMultiplicity)>; + BinningTypeDerived configurableBinningDerived{{getMultiplicity}, {axisVertex, axisMultiplicity}, true}; // true is for 'ignore overflows' (true by default). Underflows and overflows will have bin -1. + const auto multiplicity = getAnalysisMultiplicity(collision); if (cfgVerbosity > 0) { - LOGF(info, "processSameDerivedT: Tracks for collision: %d/%d | Vertex: %.1f | Multiplicity/Centrality: %.1f", tracks1.size(), tracks2.size(), collision.posZ(), collision.multiplicity()); + LOGF(info, "processSameDerivedT: Tracks for collision: %d/%d | Vertex: %.1f | Multiplicity/Centrality: %.1f", tracks1.size(), tracks2.size(), collision.posZ(), multiplicity); } loadEfficiency(collision.timestamp()); loadCcdbYieldTemplates(collision.timestamp()); - const auto multiplicity = collision.multiplicity(); - int field = 0; if (cfgTwoTrackCut > 0) { field = getMagneticField(collision.timestamp()); } - int bin = configurableBinningDerived.getBin({collision.posZ(), collision.multiplicity()}); + int bin = configurableBinningDerived.getBin(std::tuple(collision.posZ(), multiplicity)); registry.fill(HIST("eventcount_same"), bin); registry.fill(HIST("trackcount_same"), bin, tracks1.size()); if constexpr (std::experimental::is_detected::value) { @@ -1748,6 +1784,12 @@ struct TwoParticleCorrelationsMpi { } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processSameDerived, "Process same event on derived data", false); + void processSameDerivedCorrected(DerivedCollisionsCorrected::iterator const& collision, soa::Filtered const& tracks) + { + processSameDerivedT(collision, tracks, tracks); + } + PROCESS_SWITCH(TwoParticleCorrelationsMpi, processSameDerivedCorrected, "Process same event on derived data with corrected multiplicity", false); + void processSameDerivedMultSet(soa::Filtered>::iterator const& collision, soa::Filtered const& tracks) { if (!passOutlier(collision)) { @@ -1757,6 +1799,15 @@ struct TwoParticleCorrelationsMpi { } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processSameDerivedMultSet, "Process same event on derived data with multiplicity sets", false); + void processSameDerivedMultSetCorrected(soa::Filtered>::iterator const& collision, soa::Filtered const& tracks) + { + if (!passOutlier(collision)) { + return; + } + processSameDerivedT(collision, tracks, tracks); + } + PROCESS_SWITCH(TwoParticleCorrelationsMpi, processSameDerivedMultSetCorrected, "Process same event on derived data with corrected multiplicity and multiplicity sets", false); + using BinningTypeAOD = ColumnBinningPolicy; void processMixedAOD(AodCollisions const& collisions, AodTracks const& tracks, aod::BCsWithTimestamps const&) { @@ -1810,7 +1861,7 @@ struct TwoParticleCorrelationsMpi { void processMixedDerivedT(CollType const& collisions, TrackTypes&&... tracks) { auto getMultiplicity = - [this](auto& col) { + [this](const auto& col) { if constexpr (std::experimental::is_detected::value) { if (!passOutlier(col)) { return -1.0f; @@ -1818,7 +1869,7 @@ struct TwoParticleCorrelationsMpi { } else { (void)this; // fix compile error on unused 'this' capture } - return col.multiplicity(); + return getAnalysisMultiplicity(col); }; using BinningTypeDerived = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getMultiplicity)>; @@ -1841,7 +1892,7 @@ struct TwoParticleCorrelationsMpi { } if (cfgVerbosity > 0) { - LOGF(info, "processMixedDerived: Mixed collisions bin: %d pair: [%d, %d] %d (%.3f, %.3f), %d (%.3f, %.3f)", bin, it.isNewWindow(), it.currentWindowNeighbours(), collision1.globalIndex(), collision1.posZ(), collision1.multiplicity(), collision2.globalIndex(), collision2.posZ(), collision2.multiplicity()); + LOGF(info, "processMixedDerived: Mixed collisions bin: %d pair: [%d, %d] %d (%.3f, %.3f), %d (%.3f, %.3f)", bin, it.isNewWindow(), it.currentWindowNeighbours(), collision1.globalIndex(), collision1.posZ(), multiplicity, collision2.globalIndex(), collision2.posZ(), getAnalysisMultiplicity(collision2)); } bool hasEfficiencyMixed = (cfg.mEfficiencyAssociated != nullptr || cfg.mEfficiencyTrigger != nullptr); @@ -1855,14 +1906,14 @@ struct TwoParticleCorrelationsMpi { if (cfgUserAxis == EventSeedAxis) { loadCcdbYieldTemplates(collision1.timestamp()); if constexpr (std::is_same_v, std::remove_cvref_t>) { - triggerEventSeed = getEventSeedUserAxisValue(collision1.multiplicity(), estimateEventSeedWithoutFilling(mixed, tracks1, collision1.multiplicity(), collision1.posZ(), field)); + triggerEventSeed = getEventSeedUserAxisValue(multiplicity, estimateEventSeedWithoutFilling(mixed, tracks1, multiplicity, collision1.posZ(), field)); } else { LOGF(fatal, "Event-seed user axis for mixed events requires the same trigger and associated track table so the trigger event can be estimated independently"); } } if (fillRecoMixed) { - fillContainerEvent(mixed, collision1.multiplicity(), CorrelationContainer::kCFStepReconstructed); + fillContainerEvent(mixed, multiplicity, CorrelationContainer::kCFStepReconstructed); } } @@ -1871,14 +1922,14 @@ struct TwoParticleCorrelationsMpi { registry.fill(HIST("eventcount_mixed"), bin); registry.fill(HIST("trackcount_mixed"), bin, tracks1.size(), tracks2.size()); if (fillRecoMixed) { - fillCorrelations(mixed, tracks1, tracks2, collision1.multiplicity(), collision1.posZ(), field, eventWeight, nullptr, nullptr, nullptr, triggerEventSeed); + fillCorrelations(mixed, tracks1, tracks2, multiplicity, collision1.posZ(), field, eventWeight, nullptr, nullptr, nullptr, triggerEventSeed); } if (hasEfficiencyMixed) { if (it.isNewWindow()) { - fillContainerEvent(mixed, collision1.multiplicity(), CorrelationContainer::kCFStepCorrected); + fillContainerEvent(mixed, multiplicity, CorrelationContainer::kCFStepCorrected); } - fillCorrelations(mixed, tracks1, tracks2, collision1.multiplicity(), collision1.posZ(), field, eventWeight, nullptr, nullptr, nullptr, triggerEventSeed); + fillCorrelations(mixed, tracks1, tracks2, multiplicity, collision1.posZ(), field, eventWeight, nullptr, nullptr, nullptr, triggerEventSeed); } } } @@ -1889,12 +1940,24 @@ struct TwoParticleCorrelationsMpi { } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMixedDerived, "Process mixed events on derived data", false); + void processMixedDerivedCorrected(DerivedCollisionsCorrected const& collisions, DerivedTracks const& tracks) + { + processMixedDerivedT(collisions, tracks); + } + PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMixedDerivedCorrected, "Process mixed events on derived data with corrected multiplicity", false); + void processMixedDerivedMultSet(soa::Filtered> const& collisions, DerivedTracks const& tracks) { processMixedDerivedT(collisions, tracks); } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMixedDerivedMultSet, "Process mixed events on derived data with multiplicity sets", false); + void processMixedDerivedMultSetCorrected(soa::Filtered> const& collisions, DerivedTracks const& tracks) + { + processMixedDerivedT(collisions, tracks); + } + PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMixedDerivedMultSetCorrected, "Process mixed events on derived data with corrected multiplicity and multiplicity sets", false); + int getSpecies(int pdgCode) { switch (pdgCode) { @@ -1926,22 +1989,36 @@ struct TwoParticleCorrelationsMpi { } auto multiplicity = mcCollision.multiplicity(); + const bool useSingleRecoCollision = cfgCentBinsForMC == McEfficiencySingleRecoCollisionMode; if (cfgCentBinsForMC > 0) { if (collisions.size() == 0) { return; } - for (const auto& collision : collisions) { - multiplicity = collision.multiplicity(); + if (useSingleRecoCollision) { + multiplicity = collisions.begin().multiplicity(); + } else { + for (const auto& collision : collisions) { + multiplicity = collision.multiplicity(); + } } } // Primaries + int generatedPrimaries = 0; for (const auto& mcParticle : mcParticles) { if (mcParticle.isPhysicalPrimary() && mcParticle.sign() != 0 && !(std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), mcParticle.pdgCode()) != cfgMcTriggerPDGs->end())) { + ++generatedPrimaries; same->getTrackHistEfficiency()->Fill(CorrelationContainer::MC, mcParticle.eta(), mcParticle.pt(), getSpecies(mcParticle.pdgCode()), multiplicity, mcCollision.posZ()); } } for (const auto& collision : collisions) { + if (useSingleRecoCollision && collision.globalIndex() != collisions.begin().globalIndex()) { + continue; + } auto groupedTracks = tracks.sliceBy(perCollision, collision.globalIndex()); + int reconstructedTracks = 0; + int matchedTracks = 0; + std::unordered_map tracksPerMcParticle; + std::unordered_map primaryTracksPerMcParticle; if (cfgVerbosity > 0) { LOGF(info, " Reconstructed collision at vtx-z = %f", collision.posZ()); LOGF(info, " which has %d tracks", groupedTracks.size()); @@ -1951,9 +2028,13 @@ struct TwoParticleCorrelationsMpi { if (cfgTrackBitMask > 0 && (track.trackType() & (uint8_t)cfgTrackBitMask) != (uint8_t)cfgTrackBitMask) { continue; } + ++reconstructedTracks; if (track.has_cfMCParticle()) { + ++matchedTracks; + ++tracksPerMcParticle[track.cfMCParticleId()]; const auto& mcParticle = track.cfMCParticle(); if (mcParticle.isPhysicalPrimary()) { + ++primaryTracksPerMcParticle[track.cfMCParticleId()]; same->getTrackHistEfficiency()->Fill(CorrelationContainer::RecoPrimaries, mcParticle.eta(), mcParticle.pt(), getSpecies(mcParticle.pdgCode()), multiplicity, mcCollision.posZ()); } same->getTrackHistEfficiency()->Fill(CorrelationContainer::RecoAll, mcParticle.eta(), mcParticle.pt(), getSpecies(mcParticle.pdgCode()), multiplicity, mcCollision.posZ()); @@ -1963,6 +2044,19 @@ struct TwoParticleCorrelationsMpi { same->getTrackHistEfficiency()->Fill(CorrelationContainer::Fake, track.eta(), track.pt(), 0, multiplicity, mcCollision.posZ()); } } + const int duplicateMatchedTracks = matchedTracks - static_cast(tracksPerMcParticle.size()); + registry.fill(HIST("mcEfficiencyDiagnostics/reconstructedTracks"), reconstructedTracks); + registry.fill(HIST("mcEfficiencyDiagnostics/matchedTracks"), matchedTracks); + registry.fill(HIST("mcEfficiencyDiagnostics/uniqueMatchedParticles"), tracksPerMcParticle.size()); + registry.fill(HIST("mcEfficiencyDiagnostics/duplicateMatchedTracks"), duplicateMatchedTracks); + registry.fill(HIST("mcEfficiencyDiagnostics/duplicateFraction"), matchedTracks > 0 ? static_cast(duplicateMatchedTracks) / matchedTracks : 0.f); + for (const auto& entry : tracksPerMcParticle) { + registry.fill(HIST("mcEfficiencyDiagnostics/tracksPerMcParticle"), entry.second); + } + for (const auto& entry : primaryTracksPerMcParticle) { + registry.fill(HIST("mcEfficiencyDiagnostics/primaryTracksPerMcParticle"), entry.second); + } + registry.fill(HIST("mcEfficiencyDiagnostics/generatedVsUniqueMatchedPrimaries"), generatedPrimaries, primaryTracksPerMcParticle.size()); } } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMCEfficiency, "MC: Extract efficiencies", false); @@ -1984,7 +2078,7 @@ struct TwoParticleCorrelationsMpi { } } - if (!(doprocessMCSameDerived || doprocessSameDerived || doprocessSameDerivedMultSet)) { + if (!(doprocessMCSameDerived || doprocessSameDerived || doprocessSameDerivedCorrected || doprocessSameDerivedMultSet || doprocessSameDerivedMultSetCorrected)) { if constexpr (std::experimental::is_detected::value) { fillQA(mcCollision, multiplicity, mcCollision.posZ(), mcParticles1, mcParticles2); } else {