diff --git a/PWGCF/DataModel/CorrelationsDerived.h b/PWGCF/DataModel/CorrelationsDerived.h index ec0f2aa11c2..4401f9cb14d 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,31 +58,42 @@ using CFMcParticle = CFMcParticles::iterator; namespace cfmultiplicity { DECLARE_SOA_COLUMN(Multiplicity, multiplicity, float); -} +} // namespace cfmultiplicity DECLARE_SOA_TABLE(CFMultiplicities, "AOD", "CFMULTIPLICITY", cfmultiplicity::Multiplicity); using CFMultiplicity = CFMultiplicities::iterator; namespace cfcollision { -DECLARE_SOA_INDEX_COLUMN(CFMcCollision, cfMcCollision); //! Index to reduced MC collision -DECLARE_SOA_COLUMN(Multiplicity, multiplicity, float); //! Centrality/multiplicity value +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 +DECLARE_SOA_COLUMN(BestRecoCollision, bestRecoCollision, bool); //! Whether this is the best reconstructed collision for the associated MC collision (largest number of contributors) } // namespace cfcollision DECLARE_SOA_TABLE(CFCollisions, "AOD", "CFCOLLISION", //! Reduced collision table o2::soa::Index<>, bc::RunNumber, collision::PosZ, cfcollision::Multiplicity, timestamp::Timestamp); DECLARE_SOA_TABLE(CFCollLabels, "AOD", "CFCOLLLABEL", //! Labels for reduced collision table - cfcollision::CFMcCollisionId); + cfcollision::CFMcCollisionId, cfcollision::BestRecoCollision); using CFCollision = CFCollisions::iterator; 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 @@ -114,7 +125,6 @@ enum MultiplicityEstimators : uint8_t { MultNTracksGlobal = 0x8, CentFT0M = 0x10, }; - inline constexpr uint32_t NMultiplicityEstimators = __builtin_ctz(CentFT0M) + 1; } // namespace cfmultset @@ -147,8 +157,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 +211,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..8f246c14653 100644 --- a/PWGCF/TableProducer/filterCorrelations.cxx +++ b/PWGCF/TableProducer/filterCorrelations.cxx @@ -9,10 +9,12 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +// o2-linter: disable=name/workflow-file (file contains several table-producer tasks) #include "PWGCF/DataModel/CorrelationsDerived.h" #include "Common/CCDB/EventSelectionParams.h" #include "Common/CCDB/TriggerAliases.h" +#include "Common/Core/TableHelper.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" @@ -21,6 +23,8 @@ #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" +#include +#include #include #include #include @@ -35,29 +39,38 @@ #include #include +#include #include +#include #include #include #include +#include +#include +#include #include #include // required for is_detected +#include +#include +#include #include #include -#include - using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::math_utils::detail; -#define FLOAT_PRECISION 0xFFFFFFF0 -#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; +constexpr std::uint32_t kFloatPrecision = 0xFFFFFFF0u; +// clang-format off +#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME { #NAME, (DEFAULT), (HELP) } // NOLINT(bugprone-macro-parentheses) +// clang-format on struct FilterCF { - Service pdg; + Service pdg{}; + Service ccdb{}; enum TrackSelectionCuts1 : uint8_t { kTrackSelected = BIT(0), @@ -73,36 +86,40 @@ struct FilterCF { }; // Configuration - O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 7.0f, "Accepted z-vertex range") - O2_DEFINE_CONFIGURABLE(cfgCutPt, float, 0.5f, "Minimal pT for tracks") - O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks") - 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, 9 = sel8 + kNoSameBunchPileup + kIsGoodZvtxFT0vsPV, 10 = sel8 before April, 2024, 11 = sel8 for MC, 12 = sel8 with low occupancy cut, 13 = sel8 + kNoSameBunchPileup + kIsGoodITSLayersAll -- for OO/NeNe) ") - O2_DEFINE_CONFIGURABLE(cfgMinOcc, int, 0, "minimum occupancy selection") - O2_DEFINE_CONFIGURABLE(cfgMaxOcc, int, 3000, "maximum occupancy selection") - O2_DEFINE_CONFIGURABLE(cfgCollisionFlags, uint16_t, aod::collision::CollisionFlagsRun2::Run2VertexerTracks, "Request collision flags if non-zero (0 = off, 1 = Run2VertexerTracks)") - O2_DEFINE_CONFIGURABLE(cfgTransientTables, bool, false, "Output transient tables for collision and track IDs to enable successive filtering tasks") - O2_DEFINE_CONFIGURABLE(cfgTrackSelection, int, 0, "Type of track selection (0 = Run 2/3 without systematics | 1 = Run 3 with systematics | 2 = Run 3 with proton pid selection)") - O2_DEFINE_CONFIGURABLE(cfgMinMultiplicity, float, -1, "Minimum multiplicity considered for filtering (if value positive)") - O2_DEFINE_CONFIGURABLE(cfgMcSpecialPDGs, std::vector, {}, "Special MC PDG codes to include in the MC primary particle output (additional to charged particles). Empty = charged particles only.") // needed for some neutral particles - O2_DEFINE_CONFIGURABLE(nsigmaCutTPCProton, float, 3, "proton nsigma TPC") - O2_DEFINE_CONFIGURABLE(nsigmaCutTOFProton, float, 3, "proton nsigma TOF") - O2_DEFINE_CONFIGURABLE(ITSProtonselection, bool, false, "flag for ITS proton nsigma selection") - O2_DEFINE_CONFIGURABLE(nsigmaCutITSProton, float, 3, "proton nsigma ITS") - O2_DEFINE_CONFIGURABLE(dcaxymax, float, 999.f, "maximum dcaxy of tracks") - O2_DEFINE_CONFIGURABLE(dcazmax, float, 999.f, "maximum dcaz of tracks") - O2_DEFINE_CONFIGURABLE(enablePtDepDCAxy, bool, false, "Enable pT-dependent DCAxy cut: |DCAxy| < a + b/pT") - O2_DEFINE_CONFIGURABLE(dcaXyConst, float, 0.004f, "Constant term 'a' for pT-dependent DCAxy cut: |DCAxy| < a + b/pT (cm)") - O2_DEFINE_CONFIGURABLE(dcaXySlope, float, 0.013f, "Slope term 'b' for pT-dependent DCAxy cut: |DCAxy| < a + b/pT (cm x GeV/c)") - O2_DEFINE_CONFIGURABLE(itsnclusters, int, 5, "minimum number of ITS clusters for tracks") - O2_DEFINE_CONFIGURABLE(tpcncrossedrows, int, 80, "minimum number of TPC crossed rows for tracks") - O2_DEFINE_CONFIGURABLE(tpcnclusters, int, 50, "minimum number of TPC clusters found") - O2_DEFINE_CONFIGURABLE(chi2pertpccluster, float, 2.5, "maximum Chi2 / cluster for the TPC track segment") - O2_DEFINE_CONFIGURABLE(chi2peritscluster, float, 36, "maximum Chi2 / cluster for the ITS track segment") + O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 7.0f, "Accepted z-vertex range"); + O2_DEFINE_CONFIGURABLE(cfgCutPt, float, 0.5f, "Minimal pT for tracks"); + O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks"); + 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, 9 = sel8 + kNoSameBunchPileup + kIsGoodZvtxFT0vsPV, 10 = sel8 before April, 2024, 11 = sel8 for MC, 12 = sel8 with low occupancy cut, 13 = sel8 + kNoSameBunchPileup + kIsGoodITSLayersAll -- for OO/NeNe) "); + O2_DEFINE_CONFIGURABLE(cfgMinOcc, int, 0, "minimum occupancy selection"); + O2_DEFINE_CONFIGURABLE(cfgMaxOcc, int, 3000, "maximum occupancy selection"); + O2_DEFINE_CONFIGURABLE(cfgCollisionFlags, uint16_t, aod::collision::CollisionFlagsRun2::Run2VertexerTracks, "Request collision flags if non-zero (0 = off, 1 = Run2VertexerTracks)"); + O2_DEFINE_CONFIGURABLE(cfgTransientTables, bool, false, "Output transient tables for collision and track IDs to enable successive filtering tasks"); + O2_DEFINE_CONFIGURABLE(cfgTrackSelection, int, 0, "Type of track selection (0 = Run 2/3 without systematics | 1 = Run 3 with systematics | 2 = Run 3 with proton pid selection)"); + O2_DEFINE_CONFIGURABLE(cfgMinMultiplicity, float, -1, "Minimum multiplicity considered for filtering (if value positive)"); + O2_DEFINE_CONFIGURABLE(cfgMcSpecialPDGs, std::vector, std::vector{}, "Special MC PDG codes to include in the MC primary particle output (additional to charged particles). Empty = charged particles only."); // needed for some neutral particles + O2_DEFINE_CONFIGURABLE(nsigmaCutTPCProton, float, 3, "proton nsigma TPC"); + O2_DEFINE_CONFIGURABLE(nsigmaCutTOFProton, float, 3, "proton nsigma TOF"); + O2_DEFINE_CONFIGURABLE(ITSProtonselection, bool, false, "flag for ITS proton nsigma selection"); + O2_DEFINE_CONFIGURABLE(nsigmaCutITSProton, float, 3, "proton nsigma ITS"); + O2_DEFINE_CONFIGURABLE(dcaxymax, float, 999.f, "maximum dcaxy of tracks"); + O2_DEFINE_CONFIGURABLE(dcazmax, float, 999.f, "maximum dcaz of tracks"); + O2_DEFINE_CONFIGURABLE(enablePtDepDCAxy, bool, false, "Enable pT-dependent DCAxy cut: |DCAxy| < a + b/pT"); + O2_DEFINE_CONFIGURABLE(dcaXyConst, float, 0.004f, "Constant term 'a' for pT-dependent DCAxy cut: |DCAxy| < a + b/pT (cm)"); + O2_DEFINE_CONFIGURABLE(dcaXySlope, float, 0.013f, "Slope term 'b' for pT-dependent DCAxy cut: |DCAxy| < a + b/pT (cm x GeV/c)"); + O2_DEFINE_CONFIGURABLE(itsnclusters, int, 5, "minimum number of ITS clusters for tracks"); + O2_DEFINE_CONFIGURABLE(tpcncrossedrows, int, 80, "minimum number of TPC crossed rows for tracks"); + O2_DEFINE_CONFIGURABLE(tpcnclusters, int, 50, "minimum number of TPC clusters found"); + O2_DEFINE_CONFIGURABLE(chi2pertpccluster, float, 2.5, "maximum Chi2 / cluster for the TPC track segment"); + 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); @@ -114,11 +131,12 @@ struct FilterCF { Filter mcCollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVertex; OutputObj yields{TH3F("yields", "centrality vs pT vs eta", 100, 0, 100, 40, 0, 20, 100, -2, 2)}; - OutputObj etaphi{TH3F("etaphi", "centrality vs eta vs phi", 100, 0, 100, 100, -2, 2, 200, 0, 2 * M_PI)}; + OutputObj etaphi{TH3F("etaphi", "centrality vs eta vs phi", 100, 0, 100, 100, -2, 2, 200, 0, o2::constants::math::TwoPI)}; HistogramRegistry registrytrackQA{"TrackQA", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; Produces outputCollisions; + Produces outputCollisionsExtra; Produces outputTracks; Produces outputMcCollisionLabels; @@ -133,14 +151,55 @@ struct FilterCF { Produces outputMcParticleRefs; Produces outputMultSets; - std::vector multiplicities{}; + std::vector multiplicities; + + // Own local histograms independently of their input file. CCDB owns its objects. + std::unique_ptr localMultiplicityEfficiency; + THn* mEfficiency = nullptr; + static constexpr int MultiplicityEfficiencyDimensions = 4; // persistent caches std::vector mcReconstructedCache; std::vector mcParticleLabelsCache; - void init(InitContext&) + void init(InitContext& initContext) { + if (!cfgEfficiencyMultiplicity.value.empty()) { + bool processTracksEnabled = false; + if (!o2::common::core::getTaskOptionValue(initContext, "multiplicity-selector", "processTracks", processTracksEnabled, false)) { + LOGF(fatal, "Could not determine whether MultiplicitySelector::processTracks is enabled"); + } + if (!processTracksEnabled) { + LOGF(fatal, "Efficiency-corrected multiplicity requires MultiplicitySelector::processTracks"); + } + 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")); + 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()); + return; + } + localMultiplicityEfficiency.reset(dynamic_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,31 +215,38 @@ struct FilterCF { } template - bool keepCollision(TCollision& collision) + bool keepCollision(const TCollision& collision) { bool isMultSelected = false; - if (collision.multiplicity() >= cfgMinMultiplicity) + if (collision.multiplicity() >= cfgMinMultiplicity) { isMultSelected = true; - + } if (cfgTrigger == 0) { return true; - } else if (cfgTrigger == 7) { + } + 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) { + } + 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 + } + 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) + } + 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 + } + 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 + } + 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) + 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 + } + return false; + } + 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 +259,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 { @@ -233,11 +299,13 @@ struct FilterCF { if (cfgTrackSelection == 0) { if (track.isGlobalTrack()) { return 1; - } else if (track.isGlobalTrackSDD()) { + } + if (track.isGlobalTrackSDD()) { return 2; } return 0; - } else if (cfgTrackSelection == 1) { + } + if (cfgTrackSelection == 1) { uint8_t trackType = 0; if (track.isGlobalTrack()) { trackType |= kTrackSelected; @@ -258,7 +326,8 @@ struct FilterCF { } } return trackType; - } else if (cfgTrackSelection == 2) { + } + 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 +349,66 @@ struct FilterCF { return dcaXyConst + dcaXySlope / pt; // a + b/pT } + THn* loadMultiplicityEfficiency(uint64_t timestamp) + { + if (cfgLocalEfficiency == 1) { + return localMultiplicityEfficiency.get(); + } + if (!mEfficiency || !ccdb->isCachedObjectValid(cfgEfficiencyMultiplicity.value, timestamp)) { + mEfficiency = ccdb->getForTimeStamp>(cfgEfficiencyMultiplicity.value, timestamp); + } + if (!mEfficiency || mEfficiency->GetNdimensions() != MultiplicityEfficiencyDimensions) { + LOGF(fatal, "Multiplicity efficiency from %s must be a 4D THn with axes (eta, pT, multiplicity, z-vtx)", cfgEfficiencyMultiplicity.value.c_str()); + } + return mEfficiency; + } + + template + float getCorrectedMultiplicity(const TCollision& collision, const TTracks& tracks, uint64_t timestamp) + { + auto* efficiency = loadMultiplicityEfficiency(timestamp); + double correctedMultiplicity = 0.; + size_t skippedTracks = 0; + for (const auto& track : tracks) { + // Match the tracks written by the corresponding data/MC producer path. + if (!isTrackSelected(track, true)) { + 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()}; + const double eff = efficiency->GetBinContent(efficiency->GetBin(values.data())); + if (!std::isfinite(eff) || eff <= 0.) { + ++skippedTracks; + continue; + } + correctedMultiplicity += 1. / eff; + } + if (cfgVerbosity > 0 && skippedTracks > 0) { + LOGF(warning, "Skipped %zu tracks with invalid efficiency while correcting collision %lld", skippedTracks, static_cast(collision.globalIndex())); + } + 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 + bool isTrackSelected(const TTrack& track, bool checkTrackBitMask = false) + { + const float maxDCAxy = getMaxDCAxy(track.pt()); + if (std::abs(track.dcaXY()) > maxDCAxy || std::abs(track.dcaZ()) > dcazmax) { + return false; + } + if (checkTrackBitMask) { + const auto mask = static_cast(cfgMultiplicityTrackBitMask.value); + if (mask != 0 && (getTrackType(track) & mask) != mask) { + return false; + } + } + return true; + } + template using HasMultTables = decltype(std::declval().multNTracksPV()); @@ -299,33 +428,42 @@ 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(); - if (cfgEstimatorBitMask & aod::cfmultset::CentFT0C) + if (cfgEstimatorBitMask & aod::cfmultset::CentFT0C) { multiplicities.push_back(collision.centFT0C()); - if (cfgEstimatorBitMask & aod::cfmultset::MultFV0A) + } + if (cfgEstimatorBitMask & aod::cfmultset::MultFV0A) { multiplicities.push_back(collision.multFV0A()); - if (cfgEstimatorBitMask & aod::cfmultset::MultNTracksPV) + } + if (cfgEstimatorBitMask & aod::cfmultset::MultNTracksPV) { multiplicities.push_back(collision.multNTracksPV()); - if (cfgEstimatorBitMask & aod::cfmultset::MultNTracksGlobal) + } + if (cfgEstimatorBitMask & aod::cfmultset::MultNTracksGlobal) { multiplicities.push_back(collision.multNTracksGlobal()); - if (cfgEstimatorBitMask & aod::cfmultset::CentFT0M) + } + if (cfgEstimatorBitMask & aod::cfmultset::CentFT0M) { multiplicities.push_back(collision.centFT0M()); + } outputMultSets(multiplicities); } - if (cfgTransientTables) + if (cfgTransientTables) { outputCollRefs(collision.globalIndex()); - for (auto& track : tracks) { - float maxDCAxy = getMaxDCAxy(track.pt()); - if ((std::abs(track.dcaXY()) > maxDCAxy) || (std::abs(track.dcaZ()) > dcazmax)) { + } + for (const auto& track : tracks) { + if (!isTrackSelected(track)) { continue; } outputTracks(outputCollisions.lastIndex(), track.pt(), track.eta(), track.phi(), track.sign(), getTrackType(track)); - if (cfgTransientTables) + if (cfgTransientTables) { outputTrackRefs(collision.globalIndex(), track.globalIndex()); + } yields->Fill(collision.multiplicity(), track.pt(), track.eta()); etaphi->Fill(collision.multiplicity(), track.eta(), track.phi()); @@ -360,8 +498,7 @@ struct FilterCF { if (!track.isGlobalTrack()) { continue; // trackQA for global tracks only } - float maxDCAxy = getMaxDCAxy(track.pt()); - if ((std::abs(track.dcaXY()) > maxDCAxy) || (std::abs(track.dcaZ()) > dcazmax)) { + if (!isTrackSelected(track)) { continue; } registrytrackQA.fill(HIST("eta"), track.eta()); @@ -371,10 +508,12 @@ struct FilterCF { registrytrackQA.fill(HIST("tpcxrows"), track.tpcNClsCrossedRows()); registrytrackQA.fill(HIST("tpcnclst"), track.tpcNClsFound()); registrytrackQA.fill(HIST("itsnclst"), track.itsNCls()); - if (track.tpcNClsFound() > 0) + if (track.tpcNClsFound() > 0) { registrytrackQA.fill(HIST("chi2tpc"), track.tpcChi2NCl()); - if (track.itsNCls() > 0) + } + if (track.itsNCls() > 0) { registrytrackQA.fill(HIST("chi2its"), track.itsChi2NCl()); + } } } PROCESS_SWITCH(FilterCF, processTrackQA, "Process track QA", false); @@ -401,8 +540,11 @@ struct FilterCF { mcParticleLabelsCache.push_back(-1); } + std::vector bestRecoCollisionIndices(mcCollisions.size(), -1); + std::vector bestRecoCollisionNContrib(mcCollisions.size(), -1); + // 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 +554,20 @@ struct FilterCF { continue; } - for (auto& track : groupedTracks) { + const auto mcCollisionId = collision.mcCollisionId(); + if (mcCollisionId >= 0 && mcCollisionId < static_cast(bestRecoCollisionIndices.size()) && collision.numContrib() > bestRecoCollisionNContrib[mcCollisionId]) { + bestRecoCollisionNContrib[mcCollisionId] = collision.numContrib(); + bestRecoCollisionIndices[mcCollisionId] = collision.globalIndex(); + } + + 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,11 +576,11 @@ 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) { - sign = (pdgparticle->Charge() > 0) ? 1.0 : ((pdgparticle->Charge() < 0) ? -1.0 : 0.0); + sign = (pdgparticle->Charge() > 0) ? 1 : ((pdgparticle->Charge() < 0) ? -1 : 0); } bool special = !cfgMcSpecialPDGs->empty() && std::find(cfgMcSpecialPDGs->begin(), cfgMcSpecialPDGs->end(), particle.pdgCode()) != cfgMcSpecialPDGs->end(); @@ -450,10 +598,11 @@ struct FilterCF { } // 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); - if (cfgTransientTables) + outputMcParticles(outputMcCollisions.lastIndex() + 1, truncateFloatFraction(particle.pt(), kFloatPrecision), truncateFloatFraction(particle.eta(), kFloatPrecision), + truncateFloatFraction(particle.phi(), kFloatPrecision), sign, particle.pdgCode(), flags); + if (cfgTransientTables) { outputMcParticleRefs(outputMcCollisions.lastIndex() + 1, particle.globalIndex()); + } // relabeling array mcParticleLabelsCache[particle.globalIndex()] = outputMcParticles.lastIndex(); @@ -464,7 +613,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,27 +626,39 @@ 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()); - outputMcCollisionLabels(collision.mcCollisionId()); + if (!cfgEfficiencyMultiplicity.value.empty()) { + outputCollisionsExtra(getCorrectedMultiplicity(collision, groupedTracks, bc.timestamp())); + } + + const auto mcCollisionId = collision.mcCollisionId(); + const bool bestRecoCollision = mcCollisionId >= 0 && mcCollisionId < static_cast(bestRecoCollisionIndices.size()) && bestRecoCollisionIndices[mcCollisionId] == collision.globalIndex(); + outputMcCollisionLabels(mcCollisionId, bestRecoCollision); if constexpr (std::experimental::is_detected::value) { multiplicities.clear(); - if (cfgEstimatorBitMask & aod::cfmultset::CentFT0C) + if (cfgEstimatorBitMask & aod::cfmultset::CentFT0C) { multiplicities.push_back(collision.centFT0C()); - if (cfgEstimatorBitMask & aod::cfmultset::MultFV0A) + } + if (cfgEstimatorBitMask & aod::cfmultset::MultFV0A) { multiplicities.push_back(collision.multFV0A()); - if (cfgEstimatorBitMask & aod::cfmultset::MultNTracksPV) + } + if (cfgEstimatorBitMask & aod::cfmultset::MultNTracksPV) { multiplicities.push_back(collision.multNTracksPV()); - if (cfgEstimatorBitMask & aod::cfmultset::MultNTracksGlobal) + } + if (cfgEstimatorBitMask & aod::cfmultset::MultNTracksGlobal) { multiplicities.push_back(collision.multNTracksGlobal()); - if (cfgEstimatorBitMask & aod::cfmultset::CentFT0M) + } + if (cfgEstimatorBitMask & aod::cfmultset::CentFT0M) { multiplicities.push_back(collision.centFT0M()); + } outputMultSets(multiplicities); } - if (cfgTransientTables) + 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()]; @@ -507,8 +668,9 @@ struct FilterCF { } outputTracks(outputCollisions.lastIndex(), truncateFloatFraction(track.pt()), truncateFloatFraction(track.eta()), truncateFloatFraction(track.phi()), track.sign(), getTrackType(track)); outputTrackLabels(mcParticleId); - if (cfgTransientTables) + if (cfgTransientTables) { outputTrackRefs(collision.globalIndex(), track.globalIndex()); + } yields->Fill(collision.multiplicity(), track.pt(), track.eta()); etaphi->Fill(collision.multiplicity(), track.eta(), track.phi()); @@ -522,7 +684,7 @@ struct FilterCF { using McCollisionsWithHepMC = soa::Join; void processMC(McCollisionsWithHepMC const& mcCollisions, aod::McParticles const& allParticles, soa::Join const& allCollisions, - soa::Filtered> const& tracks, + soa::Filtered> const& tracks, aod::BCsWithTimestamps const& bcs) { processMCT(mcCollisions, allParticles, allCollisions, tracks, bcs); @@ -541,7 +703,7 @@ struct FilterCF { void processMCMults(McCollisionsWithHepMC const& mcCollisions, aod::McParticles const& allParticles, soa::Join const& allCollisions, - soa::Filtered> const& tracks, + soa::Filtered> const& tracks, aod::BCsWithTimestamps const& bcs) { processMCT(mcCollisions, allParticles, allCollisions, tracks, bcs); @@ -552,16 +714,20 @@ struct FilterCF { void processMCGen(McCollisionsWithHepMC::iterator const& mcCollision, aod::McParticles const& particles) { float multiplicity = 0.0f; - for (auto& particle : particles) { - if (!particle.isPhysicalPrimary() || std::abs(particle.eta()) > cfgCutMCEta || particle.pt() < cfgCutMCPt) + for (const auto& particle : particles) { + if (!particle.isPhysicalPrimary() || std::abs(particle.eta()) > cfgCutMCEta || particle.pt() < cfgCutMCPt) { continue; + } int8_t sign = 0; - if (TParticlePDG* pdgparticle = pdg->GetParticle(particle.pdgCode())) - if ((sign = pdgparticle->Charge()) != 0) + if (TParticlePDG* pdgparticle = pdg->GetParticle(particle.pdgCode())) { + sign = static_cast(pdgparticle->Charge()); + if (sign != 0) { multiplicity += 1.0f; - outputMcParticles(outputMcCollisions.lastIndex() + 1, truncateFloatFraction(particle.pt(), FLOAT_PRECISION), - truncateFloatFraction(particle.eta(), FLOAT_PRECISION), - truncateFloatFraction(particle.phi(), FLOAT_PRECISION), + } + } + outputMcParticles(outputMcCollisions.lastIndex() + 1, truncateFloatFraction(particle.pt(), kFloatPrecision), + truncateFloatFraction(particle.eta(), kFloatPrecision), + truncateFloatFraction(particle.phi(), kFloatPrecision), sign, particle.pdgCode(), particle.flags()); } outputMcCollisions(mcCollision.posZ(), multiplicity); @@ -573,8 +739,8 @@ struct FilterCF { struct MultiplicitySelector { Produces output; - O2_DEFINE_CONFIGURABLE(cfgCutPt, float, 0.5f, "Minimal pT for tracks") - O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks") + O2_DEFINE_CONFIGURABLE(cfgCutPt, float, 0.5f, "Minimal pT for tracks"); + O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks"); Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPt); Filter trackSelection = (requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true); @@ -623,7 +789,7 @@ struct MultiplicitySelector { void processFT0M(aod::CentFT0Ms const& centralities) { - for (auto& c : centralities) { + for (const auto& c : centralities) { output(c.centFT0M()); } } @@ -631,7 +797,7 @@ struct MultiplicitySelector { void processFT0C(aod::CentFT0Cs const& centralities) { - for (auto& c : centralities) { + for (const auto& c : centralities) { output(c.centFT0C()); } } @@ -639,7 +805,7 @@ struct MultiplicitySelector { void processFT0CVariant1(aod::CentFT0CVariant1s const& centralities) { - for (auto& c : centralities) { + for (const auto& c : centralities) { output(c.centFT0CVariant1()); } } @@ -647,7 +813,7 @@ struct MultiplicitySelector { void processFT0CVariant2(aod::CentFT0CVariant2s const& centralities) { - for (auto& c : centralities) { + for (const auto& c : centralities) { output(c.centFT0CVariant2()); } } @@ -655,7 +821,7 @@ struct MultiplicitySelector { void processFT0A(aod::CentFT0As const& centralities) { - for (auto& c : centralities) { + for (const auto& c : centralities) { output(c.centFT0A()); } } @@ -663,7 +829,7 @@ struct MultiplicitySelector { void processCentNGlobal(aod::CentNGlobals const& centralities) { - for (auto& c : centralities) { + for (const auto& c : centralities) { output(c.centNGlobal()); } } @@ -671,7 +837,7 @@ struct MultiplicitySelector { void processRun2V0M(aod::CentRun2V0Ms const& centralities) { - for (auto& c : centralities) { + for (const auto& c : centralities) { output(c.centRun2V0M()); } } diff --git a/PWGCF/Tasks/correlations.cxx b/PWGCF/Tasks/correlations.cxx index 76698f1d6c7..46d372b88ae 100644 --- a/PWGCF/Tasks/correlations.cxx +++ b/PWGCF/Tasks/correlations.cxx @@ -13,6 +13,8 @@ /// \brief task for the correlation calculations with CF-filtered tracks for O2 analysis /// \author Jan Fiete Grosse-Oetringhaus , Jasper Parkkila +// o2-linter: disable=name/workflow-file (preserve the established workflow filename and executable name) + #include "PWGCF/Core/CorrelationContainer.h" #include "PWGCF/Core/PairCuts.h" #include "PWGCF/DataModel/CorrelationsDerived.h" @@ -46,7 +48,6 @@ #include #include #include -#include #include #include @@ -72,7 +73,7 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace constants::math; -#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; +#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; // NOLINT(bugprone-macro-parentheses) // NOTE This is a nice idea but will again make it impossible to use subwagon configurations... // namespace o2::aod @@ -85,7 +86,7 @@ using namespace constants::math; // cfcorreff::Correction); // } // namespace o2::aod -static constexpr float kCfgPairCutDefaults[1][5] = {{-1, -1, -1, -1, -1}}; +static constexpr std::array, 1> kCfgPairCutDefaults = {{{-1, -1, -1, -1, -1}}}; // o2-linter: disable=name/constexpr-constant (preserve the established identifier) struct CorrelationTask { SliceCache cache; @@ -106,12 +107,13 @@ struct CorrelationTask { O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, int, 0, "0 = OFF and 1 = ON for local efficiency"); O2_DEFINE_CONFIGURABLE(cfgDropStepRECO, bool, false, "choice to drop step RECO if efficiency correction is used") O2_DEFINE_CONFIGURABLE(cfgCentBinsForMC, int, 0, "0 = OFF and 1 = ON for data like multiplicity/centrality bins for MC steps"); + O2_DEFINE_CONFIGURABLE(cfgRequireRecoCollision, int, 1, "0 = all generated collisions; 1 = only generated collisions with exactly 1 reconstructed collision; 2 = select reconstructed collision with largest number of contributors per generated collision") O2_DEFINE_CONFIGURABLE(cfgTrackBitMask, uint16_t, 0, "BitMask for track selection systematics; refer to the enum TrackSelectionCuts in filtering task"); O2_DEFINE_CONFIGURABLE(cfgMultCorrelationsMask, uint16_t, 0, "Selection bitmask for the multiplicity correlations. This should match the filter selection cfgEstimatorBitMask.") O2_DEFINE_CONFIGURABLE(cfgMultCutFormula, std::string, "", "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") // Suggested values: Photon: 0.004; K0 and Lambda: 0.005 - Configurable> cfgPairCut{"cfgPairCut", {kCfgPairCutDefaults[0], 5, {"Photon", "K0", "Lambda", "Phi", "Rho"}}, "Pair cuts on various particles"}; + Configurable> cfgPairCut{"cfgPairCut", {kCfgPairCutDefaults.front().data(), 5, {"Photon", "K0", "Lambda", "Phi", "Rho"}}, "Pair cuts on various particles"}; O2_DEFINE_CONFIGURABLE(cfgEfficiencyTrigger, std::string, "", "CCDB path to efficiency object for trigger particles") O2_DEFINE_CONFIGURABLE(cfgEfficiencyAssociated, std::string, "", "CCDB path to efficiency object for associated particles") @@ -175,7 +177,7 @@ struct CorrelationTask { std::vector p2indexCache; std::unique_ptr multCutFormula; - std::array multCutFormulaParamIndex; + std::array multCutFormulaParamIndex{}; struct Config { bool mPairCuts = false; @@ -187,7 +189,7 @@ struct CorrelationTask { HistogramRegistry registry{"registry"}; PairCuts mPairCuts; - Service ccdb; + Service ccdb{}; int mCachedRunNumber{-1}; // cached run number for magnetic field -- to avoid re-fetching the magnetic field for the same run, assuming that the magnetic field remains the same for the same run int mCachedMagField{0}; // cached magnetic field --reduces number of calls to the CCDB @@ -197,15 +199,27 @@ struct CorrelationTask { using DerivedCollisions = soa::Filtered; using DerivedTracks = soa::Filtered; + enum RecoCollisionSelection { + AllGeneratedCollisions = 0, + RequireOneRecoCollision, + RequireBestRecoCollision + }; + void init(o2::framework::InitContext&) { + if (cfgRequireRecoCollision < 0 || cfgRequireRecoCollision > RecoCollisionSelection::RequireBestRecoCollision) { + LOGF(fatal, "Unsupported cfgRequireRecoCollision=%d; use 0 (no reco. collision required), 1 (exactly one reco. collision required), 2 (at least one. reco, and best reco. collision selected)", cfgRequireRecoCollision.value); + } if (doprocessSame2ProngDerivedML || doprocessSame2Prong2ProngML || doprocessMixed2ProngDerivedML || doprocessMixed2Prong2ProngML || doprocessMCEfficiency2ProngML || doprocessMCReflection2ProngML) { - if (cfgPtDepMLbkg->empty() || cfgPtCentDepMLbkgSel->empty()) + if (cfgPtDepMLbkg->empty() || cfgPtCentDepMLbkgSel->empty()) { LOGF(fatal, "cfgPtDepMLbkg or cfgPtCentDepMLbkgSel can not be empty when ML 2-prong selections are used."); - if (cfgPtDepMLbkg->size() != cfgPtCentDepMLbkgSel->size()) + } + if (cfgPtDepMLbkg->size() != cfgPtCentDepMLbkgSel->size()) { LOGF(fatal, "cfgPtDepMLbkg and cfgPtCentDepMLbkgSel must be same size."); - if (!cfgPtCentDepMLpromptSel->empty() && cfgPtCentDepMLpromptSel->size() != cfgPtDepMLbkg->size()) + } + if (!cfgPtCentDepMLpromptSel->empty() && cfgPtCentDepMLpromptSel->size() != cfgPtDepMLbkg->size()) { LOGF(fatal, "cfgPtDepMLbkg and cfgPtCentDepMLpromptSel must be same size."); + } } 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"}}}); @@ -228,19 +242,25 @@ struct CorrelationTask { registry.add("invMassReflected", "2-prong invariant mass (GeV/c^2)", {HistType::kTH3F, {axisSpecMass, axisPtTrigger, axisMultiplicity}}); } if (doprocessSameDerivedMultSet) { - if (cfgMultCorrelationsMask == 0) + if (cfgMultCorrelationsMask == 0) { LOGF(fatal, "cfgMultCorrelationsMask can not be 0 when MultSet process functions are in use."); + } std::vector multAxes; - if (cfgMultCorrelationsMask & aod::cfmultset::CentFT0C) + if (cfgMultCorrelationsMask & aod::cfmultset::CentFT0C) { multAxes.emplace_back(axisMultCorrCent, "FT0C centrality"); - if (cfgMultCorrelationsMask & aod::cfmultset::MultFV0A) + } + if (cfgMultCorrelationsMask & aod::cfmultset::MultFV0A) { multAxes.emplace_back(axisMultCorrV0, "V0A multiplicity"); - if (cfgMultCorrelationsMask & aod::cfmultset::MultNTracksPV) + } + if (cfgMultCorrelationsMask & aod::cfmultset::MultNTracksPV) { multAxes.emplace_back(axisMultCorrMult, "Nch PV"); - if (cfgMultCorrelationsMask & aod::cfmultset::MultNTracksGlobal) + } + if (cfgMultCorrelationsMask & aod::cfmultset::MultNTracksGlobal) { multAxes.emplace_back(axisMultCorrMult, "Nch Global"); - if (cfgMultCorrelationsMask & aod::cfmultset::CentFT0M) + } + if (cfgMultCorrelationsMask & aod::cfmultset::CentFT0M) { multAxes.emplace_back(axisMultCorrCent, "FT0M centrality"); + } registry.add("multCorrelations", "Multiplicity correlations", {HistType::kTHnSparseF, multAxes}); } registry.add("multiplicity", "event multiplicity", {HistType::kTH1F, {{1000, 0, 100, "/multiplicity/centrality"}}}); @@ -306,10 +326,12 @@ struct CorrelationTask { userAxis.emplace_back(axisInvMass, "m (GeV/c^2)"); userMixingAxis.emplace_back(axisInvMass, "m (GeV/c^2)"); } - if (doprocessSame2Prong2Prong || doprocessSame2Prong2ProngML) + if (doprocessSame2Prong2Prong || doprocessSame2Prong2ProngML) { userAxis.emplace_back(axisInvMass, "m (GeV/c^2)"); - if (doprocessMixed2Prong2Prong || doprocessMixed2Prong2ProngML) + } + if (doprocessMixed2Prong2Prong || doprocessMixed2Prong2ProngML) { userMixingAxis.emplace_back(axisInvMass, "m (GeV/c^2)"); + } same.setObject(new CorrelationContainer("sameEvent", "sameEvent", corrAxis, effAxis, userAxis)); mixed.setObject(new CorrelationContainer("mixedEvent", "mixedEvent", corrAxis, effAxis, userMixingAxis)); @@ -317,12 +339,14 @@ struct CorrelationTask { same->setTrackEtaCut(cfgCutEta); mixed->setTrackEtaCut(cfgCutEta); - if (!cfgEfficiencyAssociated.value.empty()) + if (!cfgEfficiencyAssociated.value.empty()) { efficiencyAssociatedCache.reserve(512); + } if (doprocessMCEfficiency2Prong || doprocessMCEfficiency2ProngML || doprocessMCReflection2ProngML) { p2indexCache.reserve(16); - if (cfgMcTriggerPDGs->empty()) + if (cfgMcTriggerPDGs->empty()) { LOGF(fatal, "At least one PDG code in {} is to be selected to process 2-prong efficiency.", cfgMcTriggerPDGs.name); + } } // o2-ccdb-upload -p Users/jgrosseo/correlations/LHC15o -f /tmp/correction_2011_global.root -k correction @@ -387,8 +411,9 @@ struct CorrelationTask { { registry.fill(HIST("multiplicity"), multiplicity); if constexpr (std::experimental::is_detected::value) { - if (std::popcount(cfgMultCorrelationsMask.value) != static_cast(collision.multiplicities().size())) + if (std::popcount(cfgMultCorrelationsMask.value) != static_cast(collision.multiplicities().size())) { LOGF(fatal, "Multiplicity selections (cfgMultCorrelationsMask = 0x%x) do not match the size of the table column (%ld). The histogram filling relies on the preservation of order.", cfgMultCorrelationsMask.value, collision.multiplicities().size()); + } // need to convert to vec of doubles since THnSparse has no way to fill vec of floats directly std::vector v(collision.multiplicities().begin(), collision.multiplicities().end()); registry.get(HIST("multCorrelations")).get()->Fill(v.data()); @@ -409,21 +434,25 @@ struct CorrelationTask { { for (const auto& track1 : tracks1) { if constexpr (std::experimental::is_detected::value && std::experimental::is_detected::value) { - if (cfgDecayParticleMask != 0 && (cfgDecayParticleMask & (1u << static_cast(track1.decay()))) == 0u) + if (cfgDecayParticleMask != 0 && (cfgDecayParticleMask & (1u << static_cast(track1.decay()))) == 0u) { continue; + } if constexpr (std::experimental::is_detected::value) { - if (!passMLScore(track1)) + if (!passMLScore(track1)) { continue; + } } registry.fill(HIST("invMass"), track1.invMass(), track1.pt(), multiplicity, posZ); for (const auto& track2 : tracks2) { if constexpr (std::experimental::is_detected::value && std::experimental::is_detected::value) { if (doprocessSame2Prong2Prong || doprocessMixed2Prong2Prong || doprocessSame2Prong2ProngML || doprocessMixed2Prong2ProngML) { - if (cfgDecayParticleMask != 0 && (cfgDecayParticleMask & (1u << static_cast(track1.decay()))) == 0u) + if (cfgDecayParticleMask != 0 && (cfgDecayParticleMask & (1u << static_cast(track1.decay()))) == 0u) { continue; + } if constexpr (std::experimental::is_detected::value) { - if (!passMLScore(track2)) + if (!passMLScore(track2)) { continue; + } } if constexpr (std::experimental::is_detected::value) { @@ -452,12 +481,14 @@ struct CorrelationTask { } } // no shared prong for two mothers - if (cfgCorrelationMethod == 1 && track1.decay() != track2.decay()) + if (cfgCorrelationMethod == 1 && track1.decay() != track2.decay()) { continue; - if (cfgCorrelationMethod == 2 && track1.decay() == track2.decay()) + } + if (cfgCorrelationMethod == 2 && track1.decay() == track2.decay()) { // o2-linter: disable=magic-number (value is the established ddbar correlation-method mode) continue; + } registry.fill(HIST("invMassTwoPart"), track1.invMass(), track2.invMass(), track1.pt(), track2.pt(), multiplicity); - registry.fill(HIST("invMassTwoPartDPhi"), track1.invMass(), track2.invMass(), track1.pt(), track2.pt(), TVector2::Phi_0_2pi(track1.phi() - track2.phi() + TMath::Pi() / 2.0) - TMath::Pi() / 2.0); + registry.fill(HIST("invMassTwoPartDPhi"), track1.invMass(), track2.invMass(), track1.pt(), track2.pt(), TVector2::Phi_0_2pi(track1.phi() - track2.phi() + PIHalf) - PIHalf); if (std::abs(track1.phi() - track2.phi()) < constants::math::PI * 0.5) { registry.fill(HIST("invMassTwoPartDEta"), track1.invMass(), track2.invMass(), track1.pt(), track2.pt(), track1.eta() - track2.eta()); } @@ -466,8 +497,9 @@ struct CorrelationTask { } } if constexpr (std::experimental::is_detected::value) { - if (!cfgMcTriggerPDGs->empty() && std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), track1.pdgCode()) == cfgMcTriggerPDGs->end()) + if (!cfgMcTriggerPDGs->empty() && std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), track1.pdgCode()) == cfgMcTriggerPDGs->end()) { continue; + } } registry.fill(HIST("yieldsTrigger"), multiplicity, track1.pt(), track1.eta()); registry.fill(HIST("etaphiTrigger"), multiplicity, track1.eta(), track1.phi()); @@ -476,7 +508,7 @@ struct CorrelationTask { } template - bool fillCollisionAOD(TTarget target, TCollision collision, float multiplicity) + bool fillCollisionAOD(TTarget target, const TCollision& collision, float multiplicity) { target->fillEvent(multiplicity, CorrelationContainer::kCFStepAll); @@ -536,11 +568,13 @@ struct CorrelationTask { template bool passOutlier(CollType const& collision) { - if (cfgMultCutFormula.value.empty()) + if (cfgMultCutFormula.value.empty()) { return true; + } for (uint i = 0; i < aod::cfmultset::NMultiplicityEstimators; ++i) { - if ((cfgMultCorrelationsMask.value & (1u << i)) == 0 || multCutFormulaParamIndex[i] == ~0u) + if ((cfgMultCorrelationsMask.value & (1u << i)) == 0 || multCutFormulaParamIndex[i] == ~0u) { continue; + } auto estIndex = std::popcount(cfgMultCorrelationsMask.value & ((1u << i) - 1)); multCutFormula->SetParameter(multCutFormulaParamIndex[i], collision.multiplicities()[estIndex]); } @@ -550,8 +584,9 @@ struct CorrelationTask { template std::tuple getV0Rapidity(const T& track) { - if constexpr (!std::experimental::is_detected::value) + if constexpr (!std::experimental::is_detected::value) { return {false, 0.0f}; // no decay type, return dummy rapidity + } const auto decayType = track.decay(); float mass = 0.f; @@ -575,7 +610,7 @@ struct CorrelationTask { const float p2 = px * px + py * py + pz * pz; - const float E = std::sqrt(p2 + mass * mass); + const float E = std::sqrt(p2 + mass * mass); // o2-linter: disable=name/function-variable (E is the conventional symbol for particle energy) return {true, 0.5f * std::log((E + pz) / (E - pz))}; } @@ -606,45 +641,52 @@ struct CorrelationTask { if constexpr (std::experimental::is_detected::value) { // If the MC trigger particle is on the trigger PDG code list, we will accept them regardless of their charge. if (!cfgMcTriggerPDGs->empty()) { - if (std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), track1.pdgCode()) == cfgMcTriggerPDGs->end()) + if (std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), track1.pdgCode()) == cfgMcTriggerPDGs->end()) { continue; + } } else { // otherwise check the sign against the configuration if (cfgTriggerCharge != 0) { - if (cfgTriggerCharge * track1.sign() < 0) + if (cfgTriggerCharge * track1.sign() < 0) { continue; + } } else if (track1.sign() == 0) { continue; // reject neutral MC particles } } } else if constexpr (std::experimental::is_detected::value) { // Check reco objects that have the sign attribute. There are no neutrals to deal with. - if (cfgTriggerCharge != 0 && cfgTriggerCharge * track1.sign() < 0) + if (cfgTriggerCharge != 0 && cfgTriggerCharge * track1.sign() < 0) { continue; + } } if constexpr (std::experimental::is_detected::value) { - if (((track1.mcDecay() != aod::cf2prongtrack::D0ToPiK) && (track1.mcDecay() != aod::cf2prongtrack::D0barToKPiExclusive)) || (!cfgPtCentDepMLpromptSel->empty() && (track1.decay() & aod::cf2prongmcpart::Prompt) == 0)) + if (((track1.mcDecay() != aod::cf2prongtrack::D0ToPiK) && (track1.mcDecay() != aod::cf2prongtrack::D0barToKPiExclusive)) || (!cfgPtCentDepMLpromptSel->empty() && (track1.decay() & aod::cf2prongmcpart::Prompt) == 0)) { continue; + } } else if constexpr (std::experimental::is_detected::value) { if (cfgDecayParticleMask != 0 && (cfgDecayParticleMask & (1u << static_cast(track1.decay()))) == 0u) { continue; // skip particles that do not match the decay mask } if (cfgV0RapidityMax > 0) { auto [t, y] = getV0Rapidity(track1); - if (t && std::abs(y) > cfgV0RapidityMax) + if (t && std::abs(y) > cfgV0RapidityMax) { continue; // V0s are not allowed to be outside the rapidity range + } registry.fill(HIST("yvspt"), y, track1.pt()); } } if constexpr (std::experimental::is_detected::value) { - if (track1.cfParticleDaugh0Id() < 0 && track1.cfParticleDaugh1Id() < 0) + if (track1.cfParticleDaugh0Id() < 0 && track1.cfParticleDaugh1Id() < 0) { continue; // these we could not match + } } if constexpr (std::experimental::is_detected::value) { - if (!passMLScore(track1)) + if (!passMLScore(track1)) { continue; + } } // ML selection float triggerWeight = eventWeight; @@ -655,9 +697,9 @@ struct CorrelationTask { } if (cfgMassAxis) { - if constexpr (std::experimental::is_detected::value) + if constexpr (std::experimental::is_detected::value) { target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, track1.invMass(), triggerWeight); - else if constexpr (std::experimental::is_detected::value) { + } else if constexpr (std::experimental::is_detected::value) { // TParticlePDG *p = pdg->GetParticle(track1.pdgCode()); // target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, p->Mass(), triggerWeight); target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, 1.8, triggerWeight); @@ -699,26 +741,31 @@ struct CorrelationTask { } if constexpr (std::experimental::is_detected::value) { // skip those that are specifically chosen to be triggers - if (!cfgMcTriggerPDGs->empty() && std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), track2.pdgCode()) != cfgMcTriggerPDGs->end()) + if (!cfgMcTriggerPDGs->empty() && std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), track2.pdgCode()) != cfgMcTriggerPDGs->end()) { continue; // TODO: fix cases like MC D0-D0 + } } // Daughter track and particle checks if constexpr (std::experimental::is_detected::value) { - if (track2.globalIndex() == track1.cfTrackProng0Id()) // do not correlate daughter tracks of the same event + if (track2.globalIndex() == track1.cfTrackProng0Id()) { // do not correlate daughter tracks of the same event continue; + } } if constexpr (std::experimental::is_detected::value) { - if (track2.globalIndex() == track1.cfTrackProng1Id()) // do not correlate daughter tracks of the same event + if (track2.globalIndex() == track1.cfTrackProng1Id()) { // do not correlate daughter tracks of the same event continue; + } } if constexpr (std::experimental::is_detected::value) { - if (track2.globalIndex() == track1.cfParticleDaugh0Id()) // do not correlate daughter particles of the same event + if (track2.globalIndex() == track1.cfParticleDaugh0Id()) { // do not correlate daughter particles of the same event continue; + } } if constexpr (std::experimental::is_detected::value) { - if (track2.globalIndex() == track1.cfParticleDaugh1Id()) // do not correlate daughter particles of the same event + if (track2.globalIndex() == track1.cfParticleDaugh1Id()) { // do not correlate daughter particles of the same event continue; + } } if constexpr (step <= CorrelationContainer::kCFStepTracked && !std::experimental::is_detected::value) { @@ -729,8 +776,9 @@ struct CorrelationTask { // If decay attributes are found for the second track/particle, we assume 2p-2p correlation if constexpr (std::experimental::is_detected::value) { - if ((((track2.mcDecay()) != aod::cf2prongtrack::D0ToPiK) && ((track2.mcDecay()) != aod::cf2prongtrack::D0barToKPiExclusive)) || (!cfgPtCentDepMLpromptSel->empty() && (track2.decay() & aod::cf2prongmcpart::Prompt) == 0)) + if ((((track2.mcDecay()) != aod::cf2prongtrack::D0ToPiK) && ((track2.mcDecay()) != aod::cf2prongtrack::D0barToKPiExclusive)) || (!cfgPtCentDepMLpromptSel->empty() && (track2.decay() & aod::cf2prongmcpart::Prompt) == 0)) { continue; + } } else if constexpr (std::experimental::is_detected::value) { if (cfgDecayParticleMask != 0 && (cfgDecayParticleMask & (1u << static_cast(track2.decay()))) == 0u) { continue; // skip particles that do not match the decay mask @@ -746,10 +794,12 @@ struct CorrelationTask { } if constexpr (std::experimental::is_detected::value && std::experimental::is_detected::value) { - if (cfgCorrelationMethod == 1 && track1.decay() != track2.decay()) + if (cfgCorrelationMethod == 1 && track1.decay() != track2.decay()) { continue; - if (cfgCorrelationMethod == 2 && track1.decay() == track2.decay()) + } + if (cfgCorrelationMethod == 2 && track1.decay() == track2.decay()) { // o2-linter: disable=magic-number (value is the established ddbar correlation-method mode) continue; + } } if constexpr (std::experimental::is_detected::value) { @@ -786,8 +836,9 @@ struct CorrelationTask { if constexpr (std::experimental::is_detected::value) { // TODO: support for MC D0-D0 case if (cfgAssociatedCharge != 0) { - if (cfgAssociatedCharge * track2.sign() < 0) + if (cfgAssociatedCharge * track2.sign() < 0) { continue; + } } else if (track2.sign() == 0) { // mc particles come in neutrals, need to check explicitly continue; } @@ -822,16 +873,18 @@ struct CorrelationTask { float deltaPhi = RecoDecay::constrainAngle(track1.phi() - track2.phi(), -o2::constants::math::PIHalf); if constexpr (std::experimental::is_detected::value) { - if (!passMLScore(track2)) + if (!passMLScore(track2)) { continue; + } } // ML selection // last param is the weight if (cfgMassAxis && (doprocessSame2Prong2Prong || doprocessMixed2Prong2Prong || doprocessSame2Prong2ProngML || doprocessMixed2Prong2ProngML) && !(doprocessSame2ProngDerived || doprocessSame2ProngDerivedML || doprocessMixed2ProngDerived || doprocessMixed2ProngDerivedML || doprocessMixed2ProngDerivedMixedPhi)) { - if constexpr (std::experimental::is_detected::value && std::experimental::is_detected::value) + if constexpr (std::experimental::is_detected::value && std::experimental::is_detected::value) { target->getPairHist()->Fill(step, track1.eta() - track2.eta(), track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, track2.invMass(), track1.invMass(), associatedWeight); - else + } else { LOGF(fatal, "Can not fill mass axis without invMass column. \n no mass for two particles"); + } } else if (cfgMassAxis) { if constexpr (std::experimental::is_detected::value) { target->getPairHist()->Fill(step, track1.eta() - track2.eta(), track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, track1.invMass(), associatedWeight); @@ -862,7 +915,7 @@ struct CorrelationTask { if (cfg.mEfficiencyTrigger == nullptr) { LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgEfficiencyTrigger.value.c_str()); } - LOGF(info, "Loaded efficiency histogram for trigger particles from %s (%p)", cfgEfficiencyTrigger.value.c_str(), (void*)cfg.mEfficiencyTrigger); + LOGF(info, "Loaded efficiency histogram for trigger particles from %s (%p)", cfgEfficiencyTrigger.value.c_str(), static_cast(cfg.mEfficiencyTrigger)); } if (cfgEfficiencyAssociated.value.empty() == false) { if (cfgLocalEfficiency > 0) { @@ -874,7 +927,7 @@ struct CorrelationTask { if (cfg.mEfficiencyAssociated == nullptr) { LOGF(fatal, "Could not load efficiency histogram for associated particles from %s", cfgEfficiencyAssociated.value.c_str()); } - LOGF(info, "Loaded efficiency histogram for associated particles from %s (%p)", cfgEfficiencyAssociated.value.c_str(), (void*)cfg.mEfficiencyAssociated); + LOGF(info, "Loaded efficiency histogram for associated particles from %s (%p)", cfgEfficiencyAssociated.value.c_str(), static_cast(cfg.mEfficiencyAssociated)); } cfg.efficiencyLoaded = true; } @@ -932,10 +985,11 @@ struct CorrelationTask { int bin = configurableBinningDerived.getBin({collision.posZ(), collision.multiplicity()}); registry.fill(HIST("eventcount_same"), bin); registry.fill(HIST("trackcount_same"), bin, tracks1.size()); - if constexpr (std::experimental::is_detected::value) + if constexpr (std::experimental::is_detected::value) { fillQA(collision, multiplicity, collision.posZ(), tracks1, tracks2); - else + } else { fillQA(collision, multiplicity, tracks1); + } const bool hasEfficiency = (cfg.mEfficiencyAssociated != nullptr || cfg.mEfficiencyTrigger != nullptr); const bool fillReco = !(cfgDropStepRECO && hasEfficiency); @@ -958,8 +1012,9 @@ struct CorrelationTask { void processSameDerivedMultSet(soa::Filtered>::iterator const& collision, soa::Filtered const& tracks) { - if (!passOutlier(collision)) + if (!passOutlier(collision)) { return; + } processSameDerivedT(collision, tracks, tracks); } PROCESS_SWITCH(CorrelationTask, processSameDerivedMultSet, "Process same event on derived data with multiplicity sets", false); @@ -1037,8 +1092,9 @@ struct CorrelationTask { auto getMultiplicity = [this](auto& col) { if constexpr (std::experimental::is_detected::value) { - if (!passOutlier(col)) + if (!passOutlier(col)) { return -1.0f; + } } else { (void)this; // fix compile error on unused 'this' capture } @@ -1212,10 +1268,11 @@ struct CorrelationTask { case -2212: return 2; } - if (std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), pdgCode) != cfgMcTriggerPDGs->end()) + if (std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), pdgCode) != cfgMcTriggerPDGs->end()) { return 4; // NOTE - if changed, the number in processMCEfficiency2Prong needs to be changed too since we skip the getSpecies call - else // The efficiency histogram is hardcoded to contain 5 species. Anything special will have the 4th slot. + } else { // The efficiency histogram is hardcoded to contain 5 species. Anything special will have the 4th slot. return 3; + } } // NOTE SmallGroups includes soa::Filtered always @@ -1226,13 +1283,30 @@ struct CorrelationTask { LOGF(info, "MC collision at vtx-z = %f with %d mc particles and %d reconstructed collisions", mcCollision.posZ(), mcParticles.size(), collisions.size()); } + // Select reconstructed collisions as specified by cfgRequireRecoCollision auto multiplicity = mcCollision.multiplicity(); - if (cfgCentBinsForMC > 0) { + if (cfgRequireRecoCollision > 0) { if (collisions.size() == 0) { return; } - for (const auto& collision : collisions) { - multiplicity = collision.multiplicity(); + if (cfgRequireRecoCollision == RecoCollisionSelection::RequireOneRecoCollision) { // cfgRequireRecoCollision == 1 + if (collisions.size() != 1) { + return; + } + multiplicity = collisions.begin().multiplicity(); + } + if (cfgRequireRecoCollision == RecoCollisionSelection::RequireBestRecoCollision) { // cfgRequireRecoCollision == 2 + bool foundBestCollision = false; + for (const auto& collision : collisions) { + if (collision.bestRecoCollision()) { + multiplicity = collision.multiplicity(); + foundBestCollision = true; + break; + } + } + if (!foundBestCollision) { + return; + } } } // Primaries @@ -1241,7 +1315,11 @@ struct CorrelationTask { same->getTrackHistEfficiency()->Fill(CorrelationContainer::MC, mcParticle.eta(), mcParticle.pt(), getSpecies(mcParticle.pdgCode()), multiplicity, mcCollision.posZ()); } } + const bool useBestCollision = cfgRequireRecoCollision == RecoCollisionSelection::AllGeneratedCollisions || cfgRequireRecoCollision == RecoCollisionSelection::RequireBestRecoCollision; // For cfgRequireRecoCollision == 0 still need to reject split vertices for (const auto& collision : collisions) { + if (useBestCollision && !collision.bestRecoCollision()) { + continue; + } auto groupedTracks = tracks.sliceBy(perCollision, collision.globalIndex()); if (cfgVerbosity > 0) { LOGF(info, " Reconstructed collision at vtx-z = %f", collision.posZ()); @@ -1249,12 +1327,14 @@ struct CorrelationTask { } for (const auto& track : groupedTracks) { - if (cfgTrackBitMask > 0 && (track.trackType() & (uint8_t)cfgTrackBitMask) != (uint8_t)cfgTrackBitMask) + if (cfgTrackBitMask > 0 && (track.trackType() & (uint8_t)cfgTrackBitMask) != (uint8_t)cfgTrackBitMask) { continue; + } if (track.has_cfMCParticle()) { const auto& mcParticle = track.cfMCParticle(); - if ((doprocessMCEfficiency2Prong || doprocessMCEfficiency2ProngML) && std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), mcParticle.pdgCode()) != cfgMcTriggerPDGs->end()) + if ((doprocessMCEfficiency2Prong || doprocessMCEfficiency2ProngML) && std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), mcParticle.pdgCode()) != cfgMcTriggerPDGs->end()) { continue; // properly booked by the 2Prong efficiency function, ignore here + } if (mcParticle.isPhysicalPrimary()) { same->getTrackHistEfficiency()->Fill(CorrelationContainer::RecoPrimaries, mcParticle.eta(), mcParticle.pt(), getSpecies(mcParticle.pdgCode()), multiplicity, mcCollision.posZ()); } @@ -1270,7 +1350,7 @@ struct CorrelationTask { PROCESS_SWITCH(CorrelationTask, processMCEfficiency, "MC: Extract efficiencies", false); template - void processMCEfficiency2ProngT(soa::Filtered::iterator const& mcCollision, soa::Join const& mcParticles, soa::SmallGroups const& collisions, aod::CFTracksWithLabel const&, p2type const& p2tracks, Preslice& perCollision2Prong) + void processMCEfficiency2ProngT(soa::Filtered::iterator const& mcCollision, soa::Join const& mcParticles, soa::SmallGroups const& collisions, aod::CFTracksWithLabel const&, p2type const& p2tracks, Preslice& collision2ProngPreslice) { auto multiplicity = mcCollision.multiplicity(); if (cfgCentBinsForMC > 0) { @@ -1285,68 +1365,81 @@ struct CorrelationTask { p2indexCache.clear(); for (const auto& mcParticle : mcParticles) { if (std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), mcParticle.pdgCode()) != cfgMcTriggerPDGs->end()) { - if ((mcParticle.mcDecay() != aod::cf2prongtrack::D0ToPiK) && (mcParticle.mcDecay() != aod::cf2prongtrack::D0barToKPiExclusive)) + if ((mcParticle.mcDecay() != aod::cf2prongtrack::D0ToPiK) && (mcParticle.mcDecay() != aod::cf2prongtrack::D0barToKPiExclusive)) { continue; // wrong decay channel - if (mcParticle.cfParticleDaugh0Id() < 0 && mcParticle.cfParticleDaugh1Id() < 0) + } + if (mcParticle.cfParticleDaugh0Id() < 0 && mcParticle.cfParticleDaugh1Id() < 0) { continue; // daughters not found + } if constexpr (!reflectionSpec) { - if (cfgPtCentDepMLpromptSel->empty() || (mcParticle.decay() & aod::cf2prongmcpart::Prompt) != 0) + if (cfgPtCentDepMLpromptSel->empty() || (mcParticle.decay() & aod::cf2prongmcpart::Prompt) != 0) { same->getTrackHistEfficiency()->Fill(CorrelationContainer::MC, mcParticle.eta(), mcParticle.pt(), 4, multiplicity, mcCollision.posZ()); + } } p2indexCache.push_back(mcParticle.globalIndex()); } } for (const auto& collision : collisions) { - auto grouped2ProngTracks = p2tracks.sliceBy(perCollision2Prong, collision.globalIndex()); + auto grouped2ProngTracks = p2tracks.sliceBy(collision2ProngPreslice, collision.globalIndex()); for (const auto& p2track : grouped2ProngTracks) { if constexpr (!reflectionSpec) { - if (cfgDecayParticleMask != 0 && (cfgDecayParticleMask & (1u << static_cast(p2track.decay()))) == 0u) + if (cfgDecayParticleMask != 0 && (cfgDecayParticleMask & (1u << static_cast(p2track.decay()))) == 0u) { continue; + } } // Check if the mc particles of the prongs are found. if constexpr (std::experimental::is_detected::value) { - if (!passMLScore(p2track)) + if (!passMLScore(p2track)) { continue; + } } - if constexpr (!reflectionSpec) + if constexpr (!reflectionSpec) { same->getTrackHistEfficiency()->Fill(CorrelationContainer::RecoAll, p2track.eta(), p2track.pt(), 4, multiplicity, mcCollision.posZ()); + } auto fillMC2p = [&](const aod::CFTracksWithLabel::iterator& p) -> bool { - if (!p.has_cfMCParticle()) + if (!p.has_cfMCParticle()) { return false; + } auto m = std::find_if(p2indexCache.begin(), p2indexCache.end(), [&](const auto& t) -> bool { const auto& mcParticle = mcParticles.iteratorAt(t - mcParticles.begin().globalIndex()); return (p.cfMCParticleId() == mcParticle.cfParticleDaugh0Id() || p.cfMCParticleId() == mcParticle.cfParticleDaugh1Id()); }); - if (m == p2indexCache.end()) + if (m == p2indexCache.end()) { return false; + } const auto& mcParticle = mcParticles.iteratorAt(*m - mcParticles.begin().globalIndex()); if constexpr (!reflectionSpec) { - if (!cfgPtCentDepMLpromptSel->empty() && (mcParticle.decay() & aod::cf2prongmcpart::Prompt) == 0) + if (!cfgPtCentDepMLpromptSel->empty() && (mcParticle.decay() & aod::cf2prongmcpart::Prompt) == 0) { return true; // a valid candidate but not a prompt + } same->getTrackHistEfficiency()->Fill(CorrelationContainer::RecoPrimaries, mcParticle.eta(), mcParticle.pt(), 4, multiplicity, mcCollision.posZ()); } else { if ((mcParticle.mcDecay() == aod::cf2prongtrack::D0barToKPiExclusive && (p2track.decay() == aod::cf2prongtrack::D0barToKPiExclusive || p2track.decay() == aod::cf2prongtrack::D0barToKPi)) || - (mcParticle.mcDecay() == aod::cf2prongtrack::D0ToPiK && p2track.decay() == aod::cf2prongtrack::D0ToPiK)) + (mcParticle.mcDecay() == aod::cf2prongtrack::D0ToPiK && p2track.decay() == aod::cf2prongtrack::D0ToPiK)) { registry.fill(HIST("invMassSignal"), p2track.invMass(), p2track.pt(), multiplicity); - else // one particle may be filled into both histograms through duplicates + } else { // one particle may be filled into both histograms through duplicates registry.fill(HIST("invMassReflected"), p2track.invMass(), p2track.pt(), multiplicity); + } } return true; }; if (p2track.has_cfTrackProng0()) { // - if (const auto& p0 = p2track.template cfTrackProng0_as(); fillMC2p(p0)) + if (const auto& p0 = p2track.template cfTrackProng0_as(); fillMC2p(p0)) { continue; + } } if (p2track.has_cfTrackProng1()) { - if (const auto& p1 = p2track.template cfTrackProng1_as(); fillMC2p(p1)) + if (const auto& p1 = p2track.template cfTrackProng1_as(); fillMC2p(p1)) { continue; + } } // fake track - if constexpr (!reflectionSpec) + if constexpr (!reflectionSpec) { same->getTrackHistEfficiency()->Fill(CorrelationContainer::Fake, p2track.eta(), p2track.pt(), 4, multiplicity, mcCollision.posZ()); + } } } } @@ -1389,10 +1482,11 @@ struct CorrelationTask { } if (!(doprocessSameDerived || doprocessSameDerivedMultSet || doprocessSame2ProngDerived || doprocessSame2ProngDerivedML || doprocessSame2Prong2Prong || doprocessSame2Prong2ProngML)) { - if constexpr (std::experimental::is_detected::value) + if constexpr (std::experimental::is_detected::value) { fillQA(mcCollision, multiplicity, mcCollision.posZ(), mcParticles1, mcParticles2); - else + } else { fillQA(mcCollision, multiplicity, mcParticles1); + } } same->fillEvent(multiplicity, CorrelationContainer::kCFStepAll); @@ -1439,11 +1533,13 @@ struct CorrelationTask { bool useMCMultiplicity = (cfgCentBinsForMC == 0); auto getMultiplicity = [&collisions, &useMCMultiplicity, this](auto& col) { - if (useMCMultiplicity) + if (useMCMultiplicity) { return col.multiplicity(); + } auto groupedCollisions = collisions.sliceBy(collisionPerMCCollision, col.globalIndex()); - if (groupedCollisions.size() == 0) + if (groupedCollisions.size() == 0) { return -1.0f; + } return groupedCollisions.begin().multiplicity(); }; diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index 908ce848d7d..c17eb1de1b7 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -39,7 +39,6 @@ #include #include #include -#include #include #include @@ -99,7 +98,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"}; 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,6 +283,9 @@ struct TwoParticleCorrelationsMpi { using AodTracks = soa::Filtered>; using DerivedCollisions = soa::Filtered; + using DerivedCollisionsCorrected = soa::Filtered; + using DerivedCollisionsMultSet = soa::Filtered>; + using DerivedCollisionsMultSetCorrected = soa::Filtered>; using DerivedTracks = soa::Filtered; void init(o2::framework::InitContext&) @@ -291,7 +293,18 @@ struct TwoParticleCorrelationsMpi { 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 > 1) { + LOGF(fatal, "Unsupported cfgCentBinsForMC=%d; use 0 (generated multiplicity), 1 (reconstructed multiplicity and all associated collisions)", cfgCentBinsForMC.value); + } + const int enabledDerivedSameProcesses = static_cast(doprocessSameDerived) + static_cast(doprocessSameDerivedCorrected) + static_cast(doprocessSameDerivedMultSet) + static_cast(doprocessSameDerivedMultSetCorrected); + if (enabledDerivedSameProcesses > 1) { + LOGF(fatal, "Only one reconstructed derived same-event process can be enabled"); + } + const int enabledDerivedMixedProcesses = static_cast(doprocessMixedDerived) + static_cast(doprocessMixedDerivedCorrected) + static_cast(doprocessMixedDerivedMultSet) + static_cast(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 +360,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."); } @@ -369,7 +382,7 @@ struct TwoParticleCorrelationsMpi { } registry.add("multCorrelations", "Multiplicity correlations", {HistType::kTHnSparseF, multAxes}); } - registry.add("multiplicity", "event multiplicity", {HistType::kTH1F, {{1000, 0, 100, "/multiplicity/centrality"}}}); + registry.add("multiplicity", "event multiplicity", {HistType::kTH1F, {{100, 0, 100, "/multiplicity/centrality"}}}); 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 +568,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) { @@ -602,11 +627,10 @@ struct TwoParticleCorrelationsMpi { template bool fillContainerEvent(TTarget target, float multiplicity, CorrelationContainer::CFStep step) { - const float containerMultiplicity = getCorrelationContainerMultiplicity(multiplicity); - if (containerMultiplicity < 0.f) { + if (multiplicity < 0.f) { return false; } - target->fillEvent(containerMultiplicity, step); + target->fillEvent(multiplicity, step); return true; } @@ -1305,7 +1329,7 @@ struct TwoParticleCorrelationsMpi { std::vector discardedPairFills; discardedTriggerFills.reserve(tracks.size()); discardedPairFills.reserve(tracks.size() * tracks.size()); - fillCorrelations(target, tracks, tracks, multiplicity, posZ, magField, 1.0f, &estimate, &discardedTriggerFills, &discardedPairFills, -1.f, false, false); + fillCorrelations(std::move(target), tracks, tracks, multiplicity, posZ, magField, 1.0f, &estimate, &discardedTriggerFills, &discardedPairFills, -1.f, false, false); finalizeEventSeedEstimate(estimate); return estimate.nuncSeeds(); } @@ -1313,8 +1337,7 @@ struct TwoParticleCorrelationsMpi { template void fillCorrelations(TTarget target, TTracks1& tracks1, TTracks2& tracks2, float multiplicity, float posZ, int magField, float eventWeight, EventSeedEstimate* seedEstimate = nullptr, std::vector* pendingTriggerFills = nullptr, std::vector* pendingPairFills = nullptr, double eventSeed = -1.0, bool fillEstimatorAcceptanceQA = true, bool fillLoopQA = true) { - const float containerMultiplicity = getCorrelationContainerMultiplicity(multiplicity); - if (containerMultiplicity < 0.f) { + if (multiplicity < 0.f) { return; } @@ -1392,22 +1415,22 @@ struct TwoParticleCorrelationsMpi { if (cfgUserAxis == EventSeedAxis) { if (pendingTriggerFills) { - pendingTriggerFills->push_back({track1.pt(), containerMultiplicity, posZ, triggerWeight}); + pendingTriggerFills->push_back({track1.pt(), multiplicity, posZ, triggerWeight}); } else { - target->getTriggerHist()->Fill(step, track1.pt(), containerMultiplicity, posZ, eventSeed, triggerWeight); + target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, eventSeed, triggerWeight); } } else if (cfgUserAxis == InvariantMassAxis) { if constexpr (std::experimental::is_detected::value) { - target->getTriggerHist()->Fill(step, track1.pt(), containerMultiplicity, posZ, track1.invMass(), triggerWeight); + target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, track1.invMass(), triggerWeight); } else if constexpr (std::experimental::is_detected::value) { // TParticlePDG *p = pdg->GetParticle(track1.pdgCode()); // target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, p->Mass(), triggerWeight); - target->getTriggerHist()->Fill(step, track1.pt(), containerMultiplicity, posZ, 1.8, triggerWeight); + target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, 1.8, triggerWeight); } else { LOGF(fatal, "Can not fill invariant-mass user axis without invMass column. Disable cfgUserAxis or select another mode."); } } else { - target->getTriggerHist()->Fill(step, track1.pt(), containerMultiplicity, posZ, triggerWeight); + target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, triggerWeight); } const bool triggerHasTemplate = seedEstimate && hasTriggerTemplate(multiplicity, track1.pt()); @@ -1522,20 +1545,20 @@ struct TwoParticleCorrelationsMpi { // last param is the weight if (cfgUserAxis == EventSeedAxis) { if (pendingPairFills) { - pendingPairFills->push_back({deltaEta, track2.pt(), track1.pt(), containerMultiplicity, deltaPhi, posZ, associatedWeight}); + pendingPairFills->push_back({deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, associatedWeight}); } else { - target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), containerMultiplicity, deltaPhi, posZ, eventSeed, associatedWeight); + target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, eventSeed, associatedWeight); } } else if (cfgUserAxis == InvariantMassAxis) { if constexpr (std::experimental::is_detected::value) { - target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), containerMultiplicity, deltaPhi, posZ, track1.invMass(), associatedWeight); + target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, track1.invMass(), associatedWeight); } else if constexpr (std::experimental::is_detected::value) { - target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), containerMultiplicity, deltaPhi, posZ, 1.8, associatedWeight); // p->Mass() + target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, 1.8, associatedWeight); // p->Mass() } else { LOGF(fatal, "Can not fill invariant-mass user axis without invMass column. Disable cfgUserAxis or select another mode."); } } else { - target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), containerMultiplicity, deltaPhi, posZ, associatedWeight); + target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), multiplicity, deltaPhi, posZ, associatedWeight); } } } @@ -1583,11 +1606,6 @@ struct TwoParticleCorrelationsMpi { return eff->GetBinContent(effVars.data()); } - float getCorrelationContainerMultiplicity(float multiplicity) const - { - return multiplicity; - } - template void processSameAODT(TCollision const& collision, TTracks const& tracks, const int* trueNMPI = nullptr) { @@ -1670,22 +1688,24 @@ 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 = getMultiplicity(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) { @@ -1742,13 +1762,19 @@ struct TwoParticleCorrelationsMpi { } } - void processSameDerived(DerivedCollisions::iterator const& collision, soa::Filtered const& tracks) + void processSameDerived(DerivedCollisions::iterator const& collision, DerivedTracks const& tracks) { processSameDerivedT(collision, tracks, tracks); } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processSameDerived, "Process same event on derived data", false); - void processSameDerivedMultSet(soa::Filtered>::iterator const& collision, soa::Filtered const& tracks) + void processSameDerivedCorrected(DerivedCollisionsCorrected::iterator const& collision, DerivedTracks const& tracks) + { + processSameDerivedT(collision, tracks, tracks); + } + PROCESS_SWITCH(TwoParticleCorrelationsMpi, processSameDerivedCorrected, "Process same event on derived data with corrected multiplicity", false); + + void processSameDerivedMultSet(DerivedCollisionsMultSet::iterator const& collision, DerivedTracks const& tracks) { if (!passOutlier(collision)) { return; @@ -1757,6 +1783,15 @@ struct TwoParticleCorrelationsMpi { } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processSameDerivedMultSet, "Process same event on derived data with multiplicity sets", false); + void processSameDerivedMultSetCorrected(DerivedCollisionsMultSetCorrected::iterator const& collision, DerivedTracks 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 +1845,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 +1853,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 +1876,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 +1890,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 +1906,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,83 +1924,23 @@ struct TwoParticleCorrelationsMpi { } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMixedDerived, "Process mixed events on derived data", false); - void processMixedDerivedMultSet(soa::Filtered> const& collisions, DerivedTracks const& tracks) + void processMixedDerivedCorrected(DerivedCollisionsCorrected const& collisions, DerivedTracks const& tracks) { processMixedDerivedT(collisions, tracks); } - PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMixedDerivedMultSet, "Process mixed events on derived data with multiplicity sets", false); + PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMixedDerivedCorrected, "Process mixed events on derived data with corrected multiplicity", false); - int getSpecies(int pdgCode) + void processMixedDerivedMultSet(DerivedCollisionsMultSet const& collisions, DerivedTracks const& tracks) { - switch (pdgCode) { - case 211: // pion - case -211: - return 0; - case 321: // Kaon - case -321: - return 1; - case 2212: // proton - case -2212: - return 2; - default: - break; - } - if (std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), pdgCode) != cfgMcTriggerPDGs->end()) { - return 4; - } - // The efficiency histogram is hardcoded to contain 5 species. Anything special will have the 4th slot. - return 3; + processMixedDerivedT(collisions, tracks); } + PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMixedDerivedMultSet, "Process mixed events on derived data with multiplicity sets", false); - // NOTE SmallGroups includes soa::Filtered always - Preslice perCollision = aod::cftrack::cfCollisionId; - void processMCEfficiency(soa::Filtered::iterator const& mcCollision, aod::CFMcParticles const& mcParticles, soa::SmallGroups const& collisions, aod::CFTracksWithLabel const& tracks) + void processMixedDerivedMultSetCorrected(DerivedCollisionsMultSetCorrected const& collisions, DerivedTracks const& tracks) { - if (cfgVerbosity > 0) { - LOGF(info, "MC collision at vtx-z = %f with %d mc particles and %d reconstructed collisions", mcCollision.posZ(), mcParticles.size(), collisions.size()); - } - - auto multiplicity = mcCollision.multiplicity(); - if (cfgCentBinsForMC > 0) { - if (collisions.size() == 0) { - return; - } - for (const auto& collision : collisions) { - multiplicity = collision.multiplicity(); - } - } - // Primaries - for (const auto& mcParticle : mcParticles) { - if (mcParticle.isPhysicalPrimary() && mcParticle.sign() != 0 && !(std::find(cfgMcTriggerPDGs->begin(), cfgMcTriggerPDGs->end(), mcParticle.pdgCode()) != cfgMcTriggerPDGs->end())) { - same->getTrackHistEfficiency()->Fill(CorrelationContainer::MC, mcParticle.eta(), mcParticle.pt(), getSpecies(mcParticle.pdgCode()), multiplicity, mcCollision.posZ()); - } - } - for (const auto& collision : collisions) { - auto groupedTracks = tracks.sliceBy(perCollision, collision.globalIndex()); - if (cfgVerbosity > 0) { - LOGF(info, " Reconstructed collision at vtx-z = %f", collision.posZ()); - LOGF(info, " which has %d tracks", groupedTracks.size()); - } - - for (const auto& track : groupedTracks) { - if (cfgTrackBitMask > 0 && (track.trackType() & (uint8_t)cfgTrackBitMask) != (uint8_t)cfgTrackBitMask) { - continue; - } - if (track.has_cfMCParticle()) { - const auto& mcParticle = track.cfMCParticle(); - if (mcParticle.isPhysicalPrimary()) { - 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()); - // LOGF(info, "Filled track %d", track.globalIndex()); - } else { - // fake track - same->getTrackHistEfficiency()->Fill(CorrelationContainer::Fake, track.eta(), track.pt(), 0, multiplicity, mcCollision.posZ()); - } - } - } + processMixedDerivedT(collisions, tracks); } - PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMCEfficiency, "MC: Extract efficiencies", false); + PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMixedDerivedMultSetCorrected, "Process mixed events on derived data with corrected multiplicity and multiplicity sets", false); template void processMCSameDerivedT(McCollision const& mcCollision, Particles1 const& mcParticles1, Particles2 const& mcParticles2, soa::SmallGroups const& collisions) @@ -1984,7 +1959,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 {