From da7c51b98d41c6a5f9daa86b1ee9665f3eb2a09c Mon Sep 17 00:00:00 2001 From: PaolaVT Date: Thu, 6 Aug 2026 10:15:31 -0600 Subject: [PATCH 01/10] replace Framework/Constants.h with CommonConstants/MathConstants.h for fattenicity task --- PWGMM/UE/Tasks/flattenicityTask.cxx | 321 ++++++++++++++++++++++++++++ 1 file changed, 321 insertions(+) create mode 100644 PWGMM/UE/Tasks/flattenicityTask.cxx diff --git a/PWGMM/UE/Tasks/flattenicityTask.cxx b/PWGMM/UE/Tasks/flattenicityTask.cxx new file mode 100644 index 00000000000..c1cb6ea100e --- /dev/null +++ b/PWGMM/UE/Tasks/flattenicityTask.cxx @@ -0,0 +1,321 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file flattenicityTask.cxx +/// \brief Flattenicity analysis task for UE studies +/// \author Eisha Rani +/// \since August 2026 + +#include "Common/Core/TrackSelection.h" +#include "Common/Core/TrackSelectionDefaults.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct FlattenicityTask { + + // --- Flattenicity constants --- + static constexpr int NCH_A = 96; + static constexpr int NCH_C = 112; + static constexpr int NCELL = NCH_A + NCH_C; + static constexpr int NPHISECTORS = 8; + static constexpr int NETA_A = NCH_A / NPHISECTORS; + static constexpr int NETA_C = NCH_C / NPHISECTORS; + + // FT0 acceptance + static constexpr float FT0A_ETA_MIN = 3.5f; + static constexpr float FT0A_ETA_MAX = 4.9f; + static constexpr float FT0C_ETA_MIN = -3.3f; + static constexpr float FT0C_ETA_MAX = -2.1f; + + // --- Event selection constants --- + static constexpr float VERTEX_CUT = 10.0f; + static constexpr float FLAT_MIN = 0.0f; + static constexpr float INEL_ETA_CUT = 1.0f; + static constexpr float MIDRAP_ETA_CUT = 0.8f; + + // --- Configurables --- + Configurable cfgTrkEtaCut{"cfgTrkEtaCut", 0.8f, "Eta range for tracks"}; + Configurable cfgTrkLowPtCut{"cfgTrkLowPtCut", 0.15f, "Minimum pT"}; + + Configurable isRun3{"isRun3", true, "is Run3 dataset"}; + Configurable timeEvsel{"timeEvsel", true, "TPC Time frame boundary cut"}; + Configurable piluprejection{"piluprejection", true, "Pileup rejection"}; + Configurable goodzvertex{"goodzvertex", true, "Good Z vertex"}; + + // --- Track selection --- + TrackSelection mySelectionPrim; + + // --- Histograms --- + HistogramRegistry registry; + + // --- Init --- + void init(InitContext const&) override + { + // Initialize track selection + mySelectionPrim = myTrackSelectionPrim(); + + // Define histograms + AxisSpec flatBins = {40, 0.0, 1.0, "#rho"}; + AxisSpec nchBins = {100, -0.5, 99.5, "N_{ch}"}; + + registry.add("hFlattenicityTruth", "Truth flattenicity; 1-#rho; Events", + HistType::kTH1D, {flatBins}); + registry.add("hFlattenicityReco", "Reco flattenicity; 1-#rho; Events", + HistType::kTH1D, {flatBins}); + registry.add("hFlattenicityCorrelation", "Truth vs Reco; 1-#rho_{truth}; 1-#rho_{reco}", + HistType::kTH2D, {flatBins, flatBins}); + registry.add("hNch", "Reco Nch distribution; N_{ch}; Events", + HistType::kTH1D, {nchBins}); + registry.add("hNchTruth", "Truth Nch distribution; N_{ch}; Events", + HistType::kTH1D, {nchBins}); + } + + // --- Track selection function --- + TrackSelection myTrackSelectionPrim() + { + TrackSelection selectedTracks; + selectedTracks.SetPtRange(0.1f, 1e10f); + selectedTracks.SetEtaRange(-0.8f, 0.8f); + selectedTracks.SetRequireITSRefit(true); + selectedTracks.SetRequireTPCRefit(true); + selectedTracks.SetMinNCrossedRowsTPC(70); + selectedTracks.SetMinNCrossedRowsOverFindableClustersTPC(0.4f); + selectedTracks.SetMaxChi2PerClusterTPC(4.0f); + selectedTracks.SetRequireHitsInITSLayers(1, {0, 1}); + selectedTracks.SetMaxChi2PerClusterITS(36.0f); + selectedTracks.SetMaxDcaXYPtDep([](float pt) { return 0.0105f + 0.0350f / std::pow(pt, 1.1f); }); + selectedTracks.SetMaxDcaZ(2.0f); + return selectedTracks; + } + + // --- Helper: Get cell ID for a particle in FT0 acceptance --- + int getCellId(float eta, float phi) + { + // Check if in FT0-A acceptance + if (eta > FT0A_ETA_MIN && eta < FT0A_ETA_MAX) { + int phiBin = static_cast(std::floor(phi / (2.0f * o2::constants::math::PI / static_cast(NPHISECTORS)))); + phiBin = std::clamp(phiBin, 0, NPHISECTORS - 1); + int etaBin = static_cast(std::floor((eta - FT0A_ETA_MIN) / ((FT0A_ETA_MAX - FT0A_ETA_MIN) / static_cast(NETA_A)))); + etaBin = std::clamp(etaBin, 0, NETA_A - 1); + return etaBin * NPHISECTORS + phiBin; + } + + // Check if in FT0-C acceptance + if (eta > FT0C_ETA_MIN && eta < FT0C_ETA_MAX) { + int phiBin = static_cast(std::floor(phi / (2.0f * o2::constants::math::PI / static_cast(NPHISECTORS)))); + phiBin = std::clamp(phiBin, 0, NPHISECTORS - 1); + int etaBin = static_cast(std::floor((eta - FT0C_ETA_MIN) / ((FT0C_ETA_MAX - FT0C_ETA_MIN) / static_cast(NETA_C)))); + etaBin = std::clamp(etaBin, 0, NETA_C - 1); + return NCH_A + etaBin * NPHISECTORS + phiBin; + } + + return -1; // Not in FT0 acceptance + } + + // --- Flattenicity calculation --- + float calculateFlattenicity(const std::vector& counts) + { + if (counts.size() != static_cast(NCELL)) { + return -1.0f; + } + + float total = 0.0f; + for (const auto& c : counts) { + total += c; + } + if (total <= 0.0f) { + return -1.0f; + } + + float mean = total / static_cast(NCELL); + if (mean <= 0.0f) { + return -1.0f; + } + + float sumSq = 0.0f; + for (const auto& c : counts) { + sumSq += (c - mean) * (c - mean); + } + + float rho = std::sqrt(sumSq / (static_cast(NCELL) * static_cast(NCELL))) / mean; + return 1.0f - rho; + } + + // --- Process Data --- + void processData(aod::Collision const& collision, + soa::Filtered const& tracks, + aod::FT0s const& ft0s) + { + // Event selection (Paola/Jesus) + if (!collision.sel8()) { + return; + } + if (std::abs(collision.posZ()) >= VERTEX_CUT) { + return; + } + + // Track loop for Nch + int nch = 0; + for (const auto& track : tracks) { + if (!mySelectionPrim.IsSelected(track)) { + continue; + } + nch++; + } + registry.fill(HIST("hNch"), nch); + + // FT0 flattenicity + auto ft0 = collision.ft0(); + if (ft0.hasAmplitudeA() && ft0.hasAmplitudeC()) { + auto ampA = ft0.amplitudeA(); + auto ampC = ft0.amplitudeC(); + + std::vector counts(NCELL, 0.0f); + for (int i = 0; i < static_cast(ampA.size()) && i < NCH_A; ++i) { + counts[i] = ampA[i]; + } + for (int i = 0; i < static_cast(ampC.size()) && i < NCH_C; ++i) { + counts[NCH_A + i] = ampC[i]; + } + + float flat = calculateFlattenicity(counts); + if (flat >= FLAT_MIN) { + registry.fill(HIST("hFlattenicityReco"), flat); + } + } + } + PROCESS_SWITCH(FlattenicityTask, processData, "Process data", true); + + // --- Process MC --- + void processMC(aod::McCollision const& mcCollision, + aod::McParticles const& particles, + soa::SmallGroups> const& collisions, + aod::FT0s const& ft0s, + aod::BCs const& /*bcs*/) + { + // ---- Truth-level processing ---- + bool inel = false; + int nchTruth = 0; + std::vector truthCounts(NCELL, 0.0f); + + for (const auto& particle : particles) { + // Check if physical primary + if (!particle.isPhysicalPrimary()) { + continue; + } + + // Check if charged + if (std::abs(particle.pdgCode()) == 0) { + continue; + } + + // Check pT > 0 + if (particle.pt() <= 0.0f) { + continue; + } + + // INEL>0 check: primary charged with |eta| < INEL_ETA_CUT + if (std::abs(particle.eta()) < INEL_ETA_CUT) { + inel = true; + } + + // Nch at midrapidity: |eta| < MIDRAP_ETA_CUT, pT > cfgTrkLowPtCut + if (std::abs(particle.eta()) < MIDRAP_ETA_CUT && particle.pt() > cfgTrkLowPtCut) { + nchTruth++; + } + + // Flattenicity: particles in FT0 acceptance + int cellId = getCellId(particle.eta(), particle.phi()); + if (cellId >= 0 && cellId < NCELL) { + truthCounts[cellId] += 1.0f; + } + } + + // Apply truth-level event selection (Paola/Jesus) + if (!inel) { + return; + } + if (std::abs(mcCollision.posZ()) >= VERTEX_CUT) { + return; + } + + // Fill truth multiplicity + registry.fill(HIST("hNchTruth"), nchTruth); + registry.fill(HIST("hNch"), nchTruth); + + // Calculate truth flattenicity + float truthFlat = calculateFlattenicity(truthCounts); + if (truthFlat >= FLAT_MIN) { + registry.fill(HIST("hFlattenicityTruth"), truthFlat); + } + + // ---- Reconstructed-level processing for matched collisions ---- + for (const auto& collision : collisions) { + // Apply reconstruction-level event selection + if (!collision.sel8()) { + continue; + } + if (std::abs(collision.posZ()) >= VERTEX_CUT) { + continue; + } + + // Get FT0 flattenicity for this collision + auto ft0 = collision.ft0(); + if (!ft0.hasAmplitudeA() || !ft0.hasAmplitudeC()) { + continue; + } + + auto ampA = ft0.amplitudeA(); + auto ampC = ft0.amplitudeC(); + + std::vector recoCounts(NCELL, 0.0f); + for (int i = 0; i < static_cast(ampA.size()) && i < NCH_A; ++i) { + recoCounts[i] = ampA[i]; + } + for (int i = 0; i < static_cast(ampC.size()) && i < NCH_C; ++i) { + recoCounts[NCH_A + i] = ampC[i]; + } + + float recoFlat = calculateFlattenicity(recoCounts); + if (recoFlat >= FLAT_MIN && truthFlat >= FLAT_MIN) { + registry.fill(HIST("hFlattenicityReco"), recoFlat); + registry.fill(HIST("hFlattenicityCorrelation"), truthFlat, recoFlat); + } + } + } + PROCESS_SWITCH(FlattenicityTask, processMC, "Process MC", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{}; + workflow.push_back(adaptAnalysisTask(cfgc)); + return workflow; +} From 2216d5e52911eb4ce415cb2e6d61b9c459a05199 Mon Sep 17 00:00:00 2001 From: PaolaVT Date: Thu, 6 Aug 2026 10:24:15 -0600 Subject: [PATCH 02/10] Add flattenicityTask task to CMakeLists --- PWGMM/UE/Tasks/CMakeLists.txt | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/PWGMM/UE/Tasks/CMakeLists.txt b/PWGMM/UE/Tasks/CMakeLists.txt index f9ee32c6e10..27d74e095c1 100644 --- a/PWGMM/UE/Tasks/CMakeLists.txt +++ b/PWGMM/UE/Tasks/CMakeLists.txt @@ -22,4 +22,9 @@ o2physics_add_dpl_workflow(ue-zdc-analysis o2physics_add_dpl_workflow(dedx-analysis SOURCES dedxAnalysis.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore - COMPONENT_NAME Analysis) \ No newline at end of file + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(flattenicity-task + SOURCES flattenicityTask.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) From e193c75d94eb4d0d23ab06ad32eeb8dc47a6d3c6 Mon Sep 17 00:00:00 2001 From: PaolaVT Date: Fri, 7 Aug 2026 00:22:30 -0600 Subject: [PATCH 03/10] Fix some errors for flattenicityTask --- PWGMM/UE/Tasks/flattenicityTask.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/PWGMM/UE/Tasks/flattenicityTask.cxx b/PWGMM/UE/Tasks/flattenicityTask.cxx index c1cb6ea100e..0569f1b201e 100644 --- a/PWGMM/UE/Tasks/flattenicityTask.cxx +++ b/PWGMM/UE/Tasks/flattenicityTask.cxx @@ -1,3 +1,8 @@ +/// \file flattenicityTask.cxx +/// \brief Flattenicity analysis task for UE studies +/// \author Eisha Rani +/// \since August 2026 + // Copyright 2019-2020 CERN and copyright holders of ALICE O2. // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. // All rights not expressly granted are reserved. @@ -9,11 +14,6 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file flattenicityTask.cxx -/// \brief Flattenicity analysis task for UE studies -/// \author Eisha Rani -/// \since August 2026 - #include "Common/Core/TrackSelection.h" #include "Common/Core/TrackSelectionDefaults.h" #include "Common/DataModel/EventSelection.h" From 33bec9e73094d7aa8ea87239619d651e9b5bc553 Mon Sep 17 00:00:00 2001 From: PaolaVT Date: Fri, 7 Aug 2026 00:42:15 -0600 Subject: [PATCH 04/10] Fix copyright for flattenicityTask --- PWGMM/UE/Tasks/flattenicityTask.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/PWGMM/UE/Tasks/flattenicityTask.cxx b/PWGMM/UE/Tasks/flattenicityTask.cxx index 0569f1b201e..c1cb6ea100e 100644 --- a/PWGMM/UE/Tasks/flattenicityTask.cxx +++ b/PWGMM/UE/Tasks/flattenicityTask.cxx @@ -1,8 +1,3 @@ -/// \file flattenicityTask.cxx -/// \brief Flattenicity analysis task for UE studies -/// \author Eisha Rani -/// \since August 2026 - // Copyright 2019-2020 CERN and copyright holders of ALICE O2. // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. // All rights not expressly granted are reserved. @@ -14,6 +9,11 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file flattenicityTask.cxx +/// \brief Flattenicity analysis task for UE studies +/// \author Eisha Rani +/// \since August 2026 + #include "Common/Core/TrackSelection.h" #include "Common/Core/TrackSelectionDefaults.h" #include "Common/DataModel/EventSelection.h" From cd97f05d819cde18779749eb843a8f4c0e62b7f3 Mon Sep 17 00:00:00 2001 From: PaolaVT Date: Fri, 7 Aug 2026 13:37:44 -0600 Subject: [PATCH 05/10] Fix errors for flattenicityTask --- PWGMM/UE/Tasks/flattenicityTask.cxx | 572 ++++++++++++++++------------ 1 file changed, 338 insertions(+), 234 deletions(-) diff --git a/PWGMM/UE/Tasks/flattenicityTask.cxx b/PWGMM/UE/Tasks/flattenicityTask.cxx index c1cb6ea100e..77f9dd02236 100644 --- a/PWGMM/UE/Tasks/flattenicityTask.cxx +++ b/PWGMM/UE/Tasks/flattenicityTask.cxx @@ -14,303 +14,407 @@ /// \author Eisha Rani /// \since August 2026 -#include "Common/Core/TrackSelection.h" -#include "Common/Core/TrackSelectionDefaults.h" +#include "PWGDQ/DataModel/ReducedInfoTables.h" + +#include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include +#include "CCDB/BasicCCDBManager.h" +#include "CommonConstants/PhysicsConstants.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" + +#include +#include +#include +#include +#include +#include + +#include #include using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; +using namespace o2::soa; +using namespace o2::constants::physics; -struct FlattenicityTask { +using FullTracks = soa::Join; - // --- Flattenicity constants --- - static constexpr int NCH_A = 96; - static constexpr int NCH_C = 112; - static constexpr int NCELL = NCH_A + NCH_C; - static constexpr int NPHISECTORS = 8; - static constexpr int NETA_A = NCH_A / NPHISECTORS; - static constexpr int NETA_C = NCH_C / NPHISECTORS; - - // FT0 acceptance - static constexpr float FT0A_ETA_MIN = 3.5f; - static constexpr float FT0A_ETA_MAX = 4.9f; - static constexpr float FT0C_ETA_MIN = -3.3f; - static constexpr float FT0C_ETA_MAX = -2.1f; - - // --- Event selection constants --- - static constexpr float VERTEX_CUT = 10.0f; - static constexpr float FLAT_MIN = 0.0f; - static constexpr float INEL_ETA_CUT = 1.0f; - static constexpr float MIDRAP_ETA_CUT = 0.8f; - - // --- Configurables --- - Configurable cfgTrkEtaCut{"cfgTrkEtaCut", 0.8f, "Eta range for tracks"}; - Configurable cfgTrkLowPtCut{"cfgTrkLowPtCut", 0.15f, "Minimum pT"}; - - Configurable isRun3{"isRun3", true, "is Run3 dataset"}; - Configurable timeEvsel{"timeEvsel", true, "TPC Time frame boundary cut"}; - Configurable piluprejection{"piluprejection", true, "Pileup rejection"}; - Configurable goodzvertex{"goodzvertex", true, "Good Z vertex"}; - - // --- Track selection --- - TrackSelection mySelectionPrim; - - // --- Histograms --- - HistogramRegistry registry; - - // --- Init --- - void init(InitContext const&) override - { - // Initialize track selection - mySelectionPrim = myTrackSelectionPrim(); - - // Define histograms - AxisSpec flatBins = {40, 0.0, 1.0, "#rho"}; - AxisSpec nchBins = {100, -0.5, 99.5, "N_{ch}"}; - - registry.add("hFlattenicityTruth", "Truth flattenicity; 1-#rho; Events", - HistType::kTH1D, {flatBins}); - registry.add("hFlattenicityReco", "Reco flattenicity; 1-#rho; Events", - HistType::kTH1D, {flatBins}); - registry.add("hFlattenicityCorrelation", "Truth vs Reco; 1-#rho_{truth}; 1-#rho_{reco}", - HistType::kTH2D, {flatBins, flatBins}); - registry.add("hNch", "Reco Nch distribution; N_{ch}; Events", - HistType::kTH1D, {nchBins}); - registry.add("hNchTruth", "Truth Nch distribution; N_{ch}; Events", - HistType::kTH1D, {nchBins}); - } +struct FlattenicityTask { - // --- Track selection function --- - TrackSelection myTrackSelectionPrim() + // ============================================ + // FT0 Constants + // ============================================ + static constexpr int N_PHI_SECTORS = 8; + static constexpr int N_ETA_A = 12; + static constexpr int N_ETA_C = 14; + static constexpr int N_CH_A = 96; + static constexpr int N_CH_C = 112; + static constexpr int N_CELL = N_CH_A + N_CH_C; // 208 + + static constexpr float FT0A_ETA_MIN = 3.5; + static constexpr float FT0A_ETA_MAX = 4.9; + static constexpr float FT0C_ETA_MIN = -3.3; + static constexpr float FT0C_ETA_MAX = -2.1; + + static constexpr int PHYSICAL_PRIMARY_BIT = 0x4; + + // ============================================ + // Histogram Definitions + // ============================================ + HistogramRegistry histos{ + "histos", + { + // Event-level + {"hEvents", "Event selection;;Counts", {HistType::kTH1F, {{5, 0, 5}}}}, + + // dNch/deta + {"hNch_INEL", "Nch distribution (INEL>0);N_{ch};Entries", {HistType::kTH1F, {{100, -0.5, 99.5}}}}, + {"hNch_FT0", "Nch distribution (INEL>0 & FT0);N_{ch};Entries", {HistType::kTH1F, {{100, -0.5, 99.5}}}}, + + // Flattenicity + {"hFlattenicity", "Flattenicity distribution;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicity_vs_Nch", "Flattenicity vs Nch;N_{ch};1-#rho", {HistType::kTH2F, {{50, -0.5, 99.5}, {50, 0.0, 1.0}}}}, + + // FT0 cell occupancy + {"hCellOccupancy", "FT0 cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{N_CELL, 0, N_CELL}}}}, + {"hCellOccupancyFT0A", "FT0-A cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{N_CH_A, 0, N_CH_A}}}}, + {"hCellOccupancyFT0C", "FT0-C cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{N_CH_C, 0, N_CH_C}}}}, + + // Multiplicity classes + {"hFlattenicity_0_10", "Flattenicity 0-10%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicity_10_20", "Flattenicity 10-20%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicity_20_30", "Flattenicity 20-30%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicity_30_40", "Flattenicity 30-40%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicity_40_50", "Flattenicity 40-50%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicity_50_60", "Flattenicity 50-60%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicity_60_70", "Flattenicity 60-70%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicity_70_80", "Flattenicity 70-80%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicity_80_90", "Flattenicity 80-90%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicity_90_100", "Flattenicity 90-100%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + }}; + + // ============================================ + // Configurables + // ============================================ + Configurable cfgPtMin{"cfgPtMin", 0.1, "Minimum pT for tracks"}; + Configurable cfgEtaMax{"cfgEtaMax", 0.8, "Maximum |eta| for tracks"}; + Configurable cfgVzMax{"cfgVzMax", 10.0, "Maximum |vz| for collisions"}; + Configurable cfgNCrossedRowsTPC{"cfgNCrossedRowsTPC", 70, "Minimum TPC crossed rows"}; + Configurable cfgChi2PerClusterTPC{"cfgChi2PerClusterTPC", 4.0, "Maximum TPC chi2 per cluster"}; + Configurable cfgChi2PerClusterITS{"cfgChi2PerClusterITS", 36.0, "Maximum ITS chi2 per cluster"}; + Configurable cfgDCAZ{"cfgDCAZ", 0.1, "Maximum DCA z"}; + Configurable cfgRequireGoldenChi2{"cfgRequireGoldenChi2", true, "Require golden chi2"}; + + // ============================================ + // Particle charge function + // ============================================ + int getCharge(int pdgCode) { - TrackSelection selectedTracks; - selectedTracks.SetPtRange(0.1f, 1e10f); - selectedTracks.SetEtaRange(-0.8f, 0.8f); - selectedTracks.SetRequireITSRefit(true); - selectedTracks.SetRequireTPCRefit(true); - selectedTracks.SetMinNCrossedRowsTPC(70); - selectedTracks.SetMinNCrossedRowsOverFindableClustersTPC(0.4f); - selectedTracks.SetMaxChi2PerClusterTPC(4.0f); - selectedTracks.SetRequireHitsInITSLayers(1, {0, 1}); - selectedTracks.SetMaxChi2PerClusterITS(36.0f); - selectedTracks.SetMaxDcaXYPtDep([](float pt) { return 0.0105f + 0.0350f / std::pow(pt, 1.1f); }); - selectedTracks.SetMaxDcaZ(2.0f); - return selectedTracks; + switch (std::abs(pdgCode)) { + case 211: // pion + case 321: // kaon + case 2212: // proton + return 1; + case 11: // electron + case 13: // muon + return -1; + default: + return 0; + } } - // --- Helper: Get cell ID for a particle in FT0 acceptance --- - int getCellId(float eta, float phi) + // ============================================ + // Flattenicity calculation + // ============================================ + float computeFlattenicity(const std::array& counts) { - // Check if in FT0-A acceptance - if (eta > FT0A_ETA_MIN && eta < FT0A_ETA_MAX) { - int phiBin = static_cast(std::floor(phi / (2.0f * o2::constants::math::PI / static_cast(NPHISECTORS)))); - phiBin = std::clamp(phiBin, 0, NPHISECTORS - 1); - int etaBin = static_cast(std::floor((eta - FT0A_ETA_MIN) / ((FT0A_ETA_MAX - FT0A_ETA_MIN) / static_cast(NETA_A)))); - etaBin = std::clamp(etaBin, 0, NETA_A - 1); - return etaBin * NPHISECTORS + phiBin; + float total = 0.0; + for (int i = 0; i < N_CELL; i++) { + total += counts[i]; } - // Check if in FT0-C acceptance - if (eta > FT0C_ETA_MIN && eta < FT0C_ETA_MAX) { - int phiBin = static_cast(std::floor(phi / (2.0f * o2::constants::math::PI / static_cast(NPHISECTORS)))); - phiBin = std::clamp(phiBin, 0, NPHISECTORS - 1); - int etaBin = static_cast(std::floor((eta - FT0C_ETA_MIN) / ((FT0C_ETA_MAX - FT0C_ETA_MIN) / static_cast(NETA_C)))); - etaBin = std::clamp(etaBin, 0, NETA_C - 1); - return NCH_A + etaBin * NPHISECTORS + phiBin; + if (total <= 0) + return -1.0; + + float mean = total / N_CELL; + if (mean <= 0) + return -1.0; + + float sumSq = 0.0; + for (int i = 0; i < N_CELL; i++) { + sumSq += (counts[i] - mean) * (counts[i] - mean); } - return -1; // Not in FT0 acceptance + float rho = std::sqrt(sumSq) / (N_CELL * mean); + return rho; } - // --- Flattenicity calculation --- - float calculateFlattenicity(const std::vector& counts) + // ============================================ + // Assign particle to FT0 cell + // ============================================ + int assignToFT0Cell(float eta, float phi, bool& isFT0A) { - if (counts.size() != static_cast(NCELL)) { - return -1.0f; + // Check if in FT0 acceptance + bool inFT0A = (eta > FT0A_ETA_MIN && eta < FT0A_ETA_MAX); + bool inFT0C = (eta > FT0C_ETA_MIN && eta < FT0C_ETA_MAX); + + if (!inFT0A && !inFT0C) + return -1; + + isFT0A = inFT0A; + + // Phi bin + int phiBin = static_cast(std::floor(phi / (2 * TMath::Pi() / N_PHI_SECTORS))); + phiBin = std::max(0, std::min(phiBin, N_PHI_SECTORS - 1)); + + int cellId = -1; + + if (inFT0A) { + // FT0-A: cells 0-95 + float etaWidth = (FT0A_ETA_MAX - FT0A_ETA_MIN) / N_ETA_A; + int etaBin = static_cast(std::floor((eta - FT0A_ETA_MIN) / etaWidth)); + etaBin = std::max(0, std::min(etaBin, N_ETA_A - 1)); + cellId = etaBin * N_PHI_SECTORS + phiBin; + } else if (inFT0C) { + // FT0-C: cells 96-207 + float etaWidth = (FT0C_ETA_MAX - FT0C_ETA_MIN) / N_ETA_C; + int etaBin = static_cast(std::floor((eta - FT0C_ETA_MIN) / etaWidth)); + etaBin = std::max(0, std::min(etaBin, N_ETA_C - 1)); + cellId = N_CH_A + etaBin * N_PHI_SECTORS + phiBin; } - float total = 0.0f; - for (const auto& c : counts) { - total += c; - } - if (total <= 0.0f) { - return -1.0f; - } + return cellId; + } - float mean = total / static_cast(NCELL); - if (mean <= 0.0f) { - return -1.0f; - } + // ============================================ + // Track selection (Paola/Jesus criteria) + // ============================================ + template + bool isSelectedTrack(const T& track) + { + // pT selection + if (track.pt() < cfgPtMin) + return false; - float sumSq = 0.0f; - for (const auto& c : counts) { - sumSq += (c - mean) * (c - mean); - } + // Eta selection + if (std::abs(track.eta()) > cfgEtaMax) + return false; - float rho = std::sqrt(sumSq / (static_cast(NCELL) * static_cast(NCELL))) / mean; - return 1.0f - rho; - } + // TPC crossed rows + if (track.tpcNClsCrossedRows() < cfgNCrossedRowsTPC) + return false; - // --- Process Data --- - void processData(aod::Collision const& collision, - soa::Filtered const& tracks, - aod::FT0s const& ft0s) - { - // Event selection (Paola/Jesus) - if (!collision.sel8()) { - return; - } - if (std::abs(collision.posZ()) >= VERTEX_CUT) { - return; - } + // TPC chi2 per cluster + if (track.tpcChi2NCl() > cfgChi2PerClusterTPC) + return false; - // Track loop for Nch - int nch = 0; - for (const auto& track : tracks) { - if (!mySelectionPrim.IsSelected(track)) { - continue; - } - nch++; - } - registry.fill(HIST("hNch"), nch); + // ITS chi2 per cluster + if (track.itsChi2NCl() > cfgChi2PerClusterITS) + return false; - // FT0 flattenicity - auto ft0 = collision.ft0(); - if (ft0.hasAmplitudeA() && ft0.hasAmplitudeC()) { - auto ampA = ft0.amplitudeA(); - auto ampC = ft0.amplitudeC(); + // DCA z + if (std::abs(track.dcaZ()) > cfgDCAZ) + return false; - std::vector counts(NCELL, 0.0f); - for (int i = 0; i < static_cast(ampA.size()) && i < NCH_A; ++i) { - counts[i] = ampA[i]; - } - for (int i = 0; i < static_cast(ampC.size()) && i < NCH_C; ++i) { - counts[NCH_A + i] = ampC[i]; - } + // Golden chi2 (global track) + if (cfgRequireGoldenChi2 && !track.isGlobalTrack()) + return false; - float flat = calculateFlattenicity(counts); - if (flat >= FLAT_MIN) { - registry.fill(HIST("hFlattenicityReco"), flat); - } - } + return true; } - PROCESS_SWITCH(FlattenicityTask, processData, "Process data", true); - - // --- Process MC --- - void processMC(aod::McCollision const& mcCollision, - aod::McParticles const& particles, - soa::SmallGroups> const& collisions, - aod::FT0s const& ft0s, - aod::BCs const& /*bcs*/) + + // ============================================ + // Process MC collisions + // ============================================ + void processMC( + aod::McCollision const& /* mcCollision */, + aod::McParticles const& mcParticles) { - // ---- Truth-level processing ---- - bool inel = false; - int nchTruth = 0; - std::vector truthCounts(NCELL, 0.0f); - - for (const auto& particle : particles) { - // Check if physical primary - if (!particle.isPhysicalPrimary()) { + // Initialize counters + std::array truthCounts; + truthCounts.fill(0.0); + + int nch_INEL = 0; + int nch_FT0 = 0; + bool hasFT0A = false; + bool hasFT0C = false; + + // Loop over MC particles + for (const auto& particle : mcParticles) { + // Check if primary + if (!(particle.flags() & PHYSICAL_PRIMARY_BIT)) continue; - } // Check if charged - if (std::abs(particle.pdgCode()) == 0) { + int charge = getCharge(particle.pdgCode()); + if (charge == 0) continue; - } - // Check pT > 0 - if (particle.pt() <= 0.0f) { + // pT > 0.1 + if (particle.pt() < cfgPtMin) continue; - } - // INEL>0 check: primary charged with |eta| < INEL_ETA_CUT - if (std::abs(particle.eta()) < INEL_ETA_CUT) { - inel = true; + // INEL>0: |eta| < 1 + if (std::abs(particle.eta()) < 1.0) { + nch_INEL++; } - // Nch at midrapidity: |eta| < MIDRAP_ETA_CUT, pT > cfgTrkLowPtCut - if (std::abs(particle.eta()) < MIDRAP_ETA_CUT && particle.pt() > cfgTrkLowPtCut) { - nchTruth++; + // dNch/deta: |eta| < 0.8 + if (std::abs(particle.eta()) < cfgEtaMax) { + nch_FT0++; } - // Flattenicity: particles in FT0 acceptance - int cellId = getCellId(particle.eta(), particle.phi()); - if (cellId >= 0 && cellId < NCELL) { - truthCounts[cellId] += 1.0f; + // FT0 acceptance + bool inFT0A = (particle.eta() > FT0A_ETA_MIN && particle.eta() < FT0A_ETA_MAX); + bool inFT0C = (particle.eta() > FT0C_ETA_MIN && particle.eta() < FT0C_ETA_MAX); + + if (inFT0A) + hasFT0A = true; + if (inFT0C) + hasFT0C = true; + + // Assign to FT0 cell + bool isFT0A = false; + int cellId = assignToFT0Cell(particle.eta(), particle.phi(), isFT0A); + + if (cellId >= 0 && cellId < N_CELL) { + truthCounts[cellId] += 1.0; + histos.fill(HIST("hCellOccupancy"), cellId); + if (isFT0A) { + histos.fill(HIST("hCellOccupancyFT0A"), cellId); + } else { + histos.fill(HIST("hCellOccupancyFT0C"), cellId - N_CH_A); + } } } - // Apply truth-level event selection (Paola/Jesus) - if (!inel) { - return; - } - if (std::abs(mcCollision.posZ()) >= VERTEX_CUT) { - return; - } + // Event selection + bool isINEL = (nch_INEL > 0); + bool isFT0 = (hasFT0A && hasFT0C); - // Fill truth multiplicity - registry.fill(HIST("hNchTruth"), nchTruth); - registry.fill(HIST("hNch"), nchTruth); + histos.fill(HIST("hEvents"), 0); // All events - // Calculate truth flattenicity - float truthFlat = calculateFlattenicity(truthCounts); - if (truthFlat >= FLAT_MIN) { - registry.fill(HIST("hFlattenicityTruth"), truthFlat); + if (isINEL) { + histos.fill(HIST("hEvents"), 1); // INEL>0 + histos.fill(HIST("hNch_INEL"), nch_FT0); } - // ---- Reconstructed-level processing for matched collisions ---- - for (const auto& collision : collisions) { - // Apply reconstruction-level event selection - if (!collision.sel8()) { - continue; + if (isINEL && isFT0) { + histos.fill(HIST("hEvents"), 2); // INEL>0 & FT0 + histos.fill(HIST("hNch_FT0"), nch_FT0); + + // Compute flattenicity + float rho = computeFlattenicity(truthCounts); + if (rho > 0) { + float flattenicity = 1.0 - rho; + histos.fill(HIST("hFlattenicity"), flattenicity); + histos.fill(HIST("hFlattenicity_vs_Nch"), nch_FT0, flattenicity); + + // Multiplicity classes (based on Nch) + if (nch_FT0 < 5) { + histos.fill(HIST("hFlattenicity_0_10"), flattenicity); + } else if (nch_FT0 < 8) { + histos.fill(HIST("hFlattenicity_10_20"), flattenicity); + } else if (nch_FT0 < 11) { + histos.fill(HIST("hFlattenicity_20_30"), flattenicity); + } else if (nch_FT0 < 14) { + histos.fill(HIST("hFlattenicity_30_40"), flattenicity); + } else if (nch_FT0 < 17) { + histos.fill(HIST("hFlattenicity_40_50"), flattenicity); + } else if (nch_FT0 < 20) { + histos.fill(HIST("hFlattenicity_50_60"), flattenicity); + } else if (nch_FT0 < 24) { + histos.fill(HIST("hFlattenicity_60_70"), flattenicity); + } else if (nch_FT0 < 28) { + histos.fill(HIST("hFlattenicity_70_80"), flattenicity); + } else if (nch_FT0 < 33) { + histos.fill(HIST("hFlattenicity_80_90"), flattenicity); + } else { + histos.fill(HIST("hFlattenicity_90_100"), flattenicity); + } } - if (std::abs(collision.posZ()) >= VERTEX_CUT) { - continue; + } + } + + PROCESS_SWITCH(FlattenicityTask, processMC, "Process MC events", true); + + // ============================================ + // Process data collisions + // ============================================ + void processData( + aod::Collision const& collision, + aod::FT0s const& ft0s, + FullTracks const& tracks) + { + // Event selection: |vz| < 10 cm + if (std::abs(collision.posZ()) > cfgVzMax) + return; + + // Find FT0 matching this collision's BC + auto ft0 = ft0s.begin(); + bool foundFT0 = false; + for (const auto& f : ft0s) { + if (f.bcId() == collision.bcId()) { + ft0 = f; + foundFT0 = true; + break; } + } + if (!foundFT0) + return; - // Get FT0 flattenicity for this collision - auto ft0 = collision.ft0(); - if (!ft0.hasAmplitudeA() || !ft0.hasAmplitudeC()) { + // Track selection and counting + // int nTracks = 0; + std::array recoCounts; + recoCounts.fill(0.0); + + for (const auto& track : tracks) { + if (!isSelectedTrack(track)) continue; + // nTracks++; + + // Assign to FT0 cell using track extrapolation + bool isFT0A = false; + int cellId = assignToFT0Cell(track.eta(), track.phi(), isFT0A); + if (cellId >= 0 && cellId < N_CELL) { + recoCounts[cellId] += 1.0; } + } - auto ampA = ft0.amplitudeA(); - auto ampC = ft0.amplitudeC(); + histos.fill(HIST("hEvents"), 3); // Data events - std::vector recoCounts(NCELL, 0.0f); - for (int i = 0; i < static_cast(ampA.size()) && i < NCH_A; ++i) { - recoCounts[i] = ampA[i]; - } - for (int i = 0; i < static_cast(ampC.size()) && i < NCH_C; ++i) { - recoCounts[NCH_A + i] = ampC[i]; - } + // Compute flattenicity from FT0 amplitudes + // The FT0 channels are stored as arrays of (channel, amplitude) pairs + // The channelA() and channelC() return vectors of pairs + std::array ft0Counts; + ft0Counts.fill(0.0); - float recoFlat = calculateFlattenicity(recoCounts); - if (recoFlat >= FLAT_MIN && truthFlat >= FLAT_MIN) { - registry.fill(HIST("hFlattenicityReco"), recoFlat); - registry.fill(HIST("hFlattenicityCorrelation"), truthFlat, recoFlat); - } + // FT0-A channels (0-95) + for (int i = 0; i < N_CH_A; i++) { + ft0Counts[i] = ft0.channelA()[i]; + } + + // FT0-C channels (96-207) + for (int i = 0; i < N_CH_C; i++) { + ft0Counts[N_CH_A + i] = ft0.channelC()[i]; + } + + float rho = computeFlattenicity(ft0Counts); + if (rho > 0) { + histos.fill(HIST("hFlattenicity"), 1.0 - rho); } } - PROCESS_SWITCH(FlattenicityTask, processMC, "Process MC", true); + + PROCESS_SWITCH(FlattenicityTask, processData, "Process data events", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From be418d248ceaff23269da00078296bbd7cc7f87a Mon Sep 17 00:00:00 2001 From: PaolaVT Date: Fri, 7 Aug 2026 14:24:08 -0600 Subject: [PATCH 06/10] Fix linker errors for flattenicityTask --- PWGMM/UE/Tasks/flattenicityTask.cxx | 147 +++++++++++++++------------- 1 file changed, 80 insertions(+), 67 deletions(-) diff --git a/PWGMM/UE/Tasks/flattenicityTask.cxx b/PWGMM/UE/Tasks/flattenicityTask.cxx index 77f9dd02236..7602b86bbf0 100644 --- a/PWGMM/UE/Tasks/flattenicityTask.cxx +++ b/PWGMM/UE/Tasks/flattenicityTask.cxx @@ -22,13 +22,13 @@ #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "CCDB/BasicCCDBManager.h" -#include "CommonConstants/PhysicsConstants.h" -#include "Framework/ASoAHelpers.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" -#include "ReconstructionDataFormats/Track.h" +#include +#include +#include +#include +#include +#include +#include #include #include @@ -53,19 +53,32 @@ struct FlattenicityTask { // ============================================ // FT0 Constants // ============================================ - static constexpr int N_PHI_SECTORS = 8; - static constexpr int N_ETA_A = 12; - static constexpr int N_ETA_C = 14; - static constexpr int N_CH_A = 96; - static constexpr int N_CH_C = 112; - static constexpr int N_CELL = N_CH_A + N_CH_C; // 208 + static constexpr int NPhiSectors = 8; + static constexpr int NEtaA = 12; + static constexpr int NEtaC = 14; + static constexpr int NchA = 96; + static constexpr int NchC = 112; + static constexpr int NCell = NchA + NchC; // 208 - static constexpr float FT0A_ETA_MIN = 3.5; - static constexpr float FT0A_ETA_MAX = 4.9; - static constexpr float FT0C_ETA_MIN = -3.3; - static constexpr float FT0C_ETA_MAX = -2.1; + static constexpr float FT0AEtaMin = 3.5; + static constexpr float FT0AEtaMax = 4.9; + static constexpr float FT0CEtaMin = -3.3; + static constexpr float FT0CEtaMax = -2.1; - static constexpr int PHYSICAL_PRIMARY_BIT = 0x4; + static constexpr int NPhysicalPrimaryBit = 0x4; + + // ============================================ + // Multiplicity class boundaries (Nch cuts defining the 0-10% ... 90-100% classes) + // ============================================ + static constexpr int NchBound10 = 5; + static constexpr int NchBound20 = 8; + static constexpr int NchBound30 = 11; + static constexpr int NchBound40 = 14; + static constexpr int NchBound50 = 17; + static constexpr int NchBound60 = 20; + static constexpr int NchBound70 = 24; + static constexpr int NchBound80 = 28; + static constexpr int NchBound90 = 33; // ============================================ // Histogram Definitions @@ -85,9 +98,9 @@ struct FlattenicityTask { {"hFlattenicity_vs_Nch", "Flattenicity vs Nch;N_{ch};1-#rho", {HistType::kTH2F, {{50, -0.5, 99.5}, {50, 0.0, 1.0}}}}, // FT0 cell occupancy - {"hCellOccupancy", "FT0 cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{N_CELL, 0, N_CELL}}}}, - {"hCellOccupancyFT0A", "FT0-A cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{N_CH_A, 0, N_CH_A}}}}, - {"hCellOccupancyFT0C", "FT0-C cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{N_CH_C, 0, N_CH_C}}}}, + {"hCellOccupancy", "FT0 cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{NCell, 0, NCell}}}}, + {"hCellOccupancyFT0A", "FT0-A cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{NchA, 0, NchA}}}}, + {"hCellOccupancyFT0C", "FT0-C cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{NchC, 0, NchC}}}}, // Multiplicity classes {"hFlattenicity_0_10", "Flattenicity 0-10%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, @@ -135,26 +148,26 @@ struct FlattenicityTask { // ============================================ // Flattenicity calculation // ============================================ - float computeFlattenicity(const std::array& counts) + float computeFlattenicity(const std::array& counts) { float total = 0.0; - for (int i = 0; i < N_CELL; i++) { + for (int i = 0; i < NCell; i++) { total += counts[i]; } if (total <= 0) return -1.0; - float mean = total / N_CELL; + float mean = total / NCell; if (mean <= 0) return -1.0; float sumSq = 0.0; - for (int i = 0; i < N_CELL; i++) { + for (int i = 0; i < NCell; i++) { sumSq += (counts[i] - mean) * (counts[i] - mean); } - float rho = std::sqrt(sumSq) / (N_CELL * mean); + float rho = std::sqrt(sumSq) / (NCell * mean); return rho; } @@ -164,8 +177,8 @@ struct FlattenicityTask { int assignToFT0Cell(float eta, float phi, bool& isFT0A) { // Check if in FT0 acceptance - bool inFT0A = (eta > FT0A_ETA_MIN && eta < FT0A_ETA_MAX); - bool inFT0C = (eta > FT0C_ETA_MIN && eta < FT0C_ETA_MAX); + bool inFT0A = (eta > FT0AEtaMin && eta < FT0AEtaMax); + bool inFT0C = (eta > FT0CEtaMin && eta < FT0CEtaMax); if (!inFT0A && !inFT0C) return -1; @@ -173,23 +186,23 @@ struct FlattenicityTask { isFT0A = inFT0A; // Phi bin - int phiBin = static_cast(std::floor(phi / (2 * TMath::Pi() / N_PHI_SECTORS))); - phiBin = std::max(0, std::min(phiBin, N_PHI_SECTORS - 1)); + int phiBin = static_cast(std::floor(phi / (o2::constants::math::TwoPI / NPhiSectors))); + phiBin = std::max(0, std::min(phiBin, NPhiSectors - 1)); int cellId = -1; if (inFT0A) { // FT0-A: cells 0-95 - float etaWidth = (FT0A_ETA_MAX - FT0A_ETA_MIN) / N_ETA_A; - int etaBin = static_cast(std::floor((eta - FT0A_ETA_MIN) / etaWidth)); - etaBin = std::max(0, std::min(etaBin, N_ETA_A - 1)); - cellId = etaBin * N_PHI_SECTORS + phiBin; + float etaWidth = (FT0AEtaMax - FT0AEtaMin) / NEtaA; + int etaBin = static_cast(std::floor((eta - FT0AEtaMin) / etaWidth)); + etaBin = std::max(0, std::min(etaBin, NEtaA - 1)); + cellId = etaBin * NPhiSectors + phiBin; } else if (inFT0C) { // FT0-C: cells 96-207 - float etaWidth = (FT0C_ETA_MAX - FT0C_ETA_MIN) / N_ETA_C; - int etaBin = static_cast(std::floor((eta - FT0C_ETA_MIN) / etaWidth)); - etaBin = std::max(0, std::min(etaBin, N_ETA_C - 1)); - cellId = N_CH_A + etaBin * N_PHI_SECTORS + phiBin; + float etaWidth = (FT0CEtaMax - FT0CEtaMin) / NEtaC; + int etaBin = static_cast(std::floor((eta - FT0CEtaMin) / etaWidth)); + etaBin = std::max(0, std::min(etaBin, NEtaC - 1)); + cellId = NchA + etaBin * NPhiSectors + phiBin; } return cellId; @@ -240,18 +253,18 @@ struct FlattenicityTask { aod::McParticles const& mcParticles) { // Initialize counters - std::array truthCounts; + std::array truthCounts; truthCounts.fill(0.0); - int nch_INEL = 0; - int nch_FT0 = 0; + int nchINEL = 0; + int nchFT0 = 0; bool hasFT0A = false; bool hasFT0C = false; // Loop over MC particles for (const auto& particle : mcParticles) { // Check if primary - if (!(particle.flags() & PHYSICAL_PRIMARY_BIT)) + if (!(particle.flags() & NPhysicalPrimaryBit)) continue; // Check if charged @@ -265,17 +278,17 @@ struct FlattenicityTask { // INEL>0: |eta| < 1 if (std::abs(particle.eta()) < 1.0) { - nch_INEL++; + nchINEL++; } // dNch/deta: |eta| < 0.8 if (std::abs(particle.eta()) < cfgEtaMax) { - nch_FT0++; + nchFT0++; } // FT0 acceptance - bool inFT0A = (particle.eta() > FT0A_ETA_MIN && particle.eta() < FT0A_ETA_MAX); - bool inFT0C = (particle.eta() > FT0C_ETA_MIN && particle.eta() < FT0C_ETA_MAX); + bool inFT0A = (particle.eta() > FT0AEtaMin && particle.eta() < FT0AEtaMax); + bool inFT0C = (particle.eta() > FT0CEtaMin && particle.eta() < FT0CEtaMax); if (inFT0A) hasFT0A = true; @@ -286,57 +299,57 @@ struct FlattenicityTask { bool isFT0A = false; int cellId = assignToFT0Cell(particle.eta(), particle.phi(), isFT0A); - if (cellId >= 0 && cellId < N_CELL) { + if (cellId >= 0 && cellId < NCell) { truthCounts[cellId] += 1.0; histos.fill(HIST("hCellOccupancy"), cellId); if (isFT0A) { histos.fill(HIST("hCellOccupancyFT0A"), cellId); } else { - histos.fill(HIST("hCellOccupancyFT0C"), cellId - N_CH_A); + histos.fill(HIST("hCellOccupancyFT0C"), cellId - NchA); } } } // Event selection - bool isINEL = (nch_INEL > 0); + bool isINEL = (nchINEL > 0); bool isFT0 = (hasFT0A && hasFT0C); histos.fill(HIST("hEvents"), 0); // All events if (isINEL) { histos.fill(HIST("hEvents"), 1); // INEL>0 - histos.fill(HIST("hNch_INEL"), nch_FT0); + histos.fill(HIST("hNch_INEL"), nchFT0); } if (isINEL && isFT0) { histos.fill(HIST("hEvents"), 2); // INEL>0 & FT0 - histos.fill(HIST("hNch_FT0"), nch_FT0); + histos.fill(HIST("hNch_FT0"), nchFT0); // Compute flattenicity float rho = computeFlattenicity(truthCounts); if (rho > 0) { float flattenicity = 1.0 - rho; histos.fill(HIST("hFlattenicity"), flattenicity); - histos.fill(HIST("hFlattenicity_vs_Nch"), nch_FT0, flattenicity); + histos.fill(HIST("hFlattenicity_vs_Nch"), nchFT0, flattenicity); // Multiplicity classes (based on Nch) - if (nch_FT0 < 5) { + if (nchFT0 < NchBound10) { histos.fill(HIST("hFlattenicity_0_10"), flattenicity); - } else if (nch_FT0 < 8) { + } else if (nchFT0 < NchBound20) { histos.fill(HIST("hFlattenicity_10_20"), flattenicity); - } else if (nch_FT0 < 11) { + } else if (nchFT0 < NchBound30) { histos.fill(HIST("hFlattenicity_20_30"), flattenicity); - } else if (nch_FT0 < 14) { + } else if (nchFT0 < NchBound40) { histos.fill(HIST("hFlattenicity_30_40"), flattenicity); - } else if (nch_FT0 < 17) { + } else if (nchFT0 < NchBound50) { histos.fill(HIST("hFlattenicity_40_50"), flattenicity); - } else if (nch_FT0 < 20) { + } else if (nchFT0 < NchBound60) { histos.fill(HIST("hFlattenicity_50_60"), flattenicity); - } else if (nch_FT0 < 24) { + } else if (nchFT0 < NchBound70) { histos.fill(HIST("hFlattenicity_60_70"), flattenicity); - } else if (nch_FT0 < 28) { + } else if (nchFT0 < NchBound80) { histos.fill(HIST("hFlattenicity_70_80"), flattenicity); - } else if (nch_FT0 < 33) { + } else if (nchFT0 < NchBound90) { histos.fill(HIST("hFlattenicity_80_90"), flattenicity); } else { histos.fill(HIST("hFlattenicity_90_100"), flattenicity); @@ -374,7 +387,7 @@ struct FlattenicityTask { // Track selection and counting // int nTracks = 0; - std::array recoCounts; + std::array recoCounts; recoCounts.fill(0.0); for (const auto& track : tracks) { @@ -385,7 +398,7 @@ struct FlattenicityTask { // Assign to FT0 cell using track extrapolation bool isFT0A = false; int cellId = assignToFT0Cell(track.eta(), track.phi(), isFT0A); - if (cellId >= 0 && cellId < N_CELL) { + if (cellId >= 0 && cellId < NCell) { recoCounts[cellId] += 1.0; } } @@ -395,17 +408,17 @@ struct FlattenicityTask { // Compute flattenicity from FT0 amplitudes // The FT0 channels are stored as arrays of (channel, amplitude) pairs // The channelA() and channelC() return vectors of pairs - std::array ft0Counts; + std::array ft0Counts; ft0Counts.fill(0.0); // FT0-A channels (0-95) - for (int i = 0; i < N_CH_A; i++) { + for (int i = 0; i < NchA; i++) { ft0Counts[i] = ft0.channelA()[i]; } // FT0-C channels (96-207) - for (int i = 0; i < N_CH_C; i++) { - ft0Counts[N_CH_A + i] = ft0.channelC()[i]; + for (int i = 0; i < NchC; i++) { + ft0Counts[NchA + i] = ft0.channelC()[i]; } float rho = computeFlattenicity(ft0Counts); From cc48bf186e624c943e48ff74da880019b37cdb47 Mon Sep 17 00:00:00 2001 From: PaolaVT Date: Fri, 7 Aug 2026 14:34:45 -0600 Subject: [PATCH 07/10] Fix Megalinker errors for flattenicityTask --- PWGMM/UE/Tasks/flattenicityTask.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGMM/UE/Tasks/flattenicityTask.cxx b/PWGMM/UE/Tasks/flattenicityTask.cxx index 7602b86bbf0..309f95a93b1 100644 --- a/PWGMM/UE/Tasks/flattenicityTask.cxx +++ b/PWGMM/UE/Tasks/flattenicityTask.cxx @@ -37,6 +37,7 @@ #include #include +#include #include #include From da883cdf05360b9606922fa905be5cfc8a09a59c Mon Sep 17 00:00:00 2001 From: PaolaVT Date: Sat, 8 Aug 2026 10:22:00 -0600 Subject: [PATCH 08/10] Fix code-check errors for flattenicityTask --- PWGMM/UE/Tasks/flattenicityTask.cxx | 102 ++++++++++++++-------------- 1 file changed, 51 insertions(+), 51 deletions(-) diff --git a/PWGMM/UE/Tasks/flattenicityTask.cxx b/PWGMM/UE/Tasks/flattenicityTask.cxx index 309f95a93b1..07c5a3c8712 100644 --- a/PWGMM/UE/Tasks/flattenicityTask.cxx +++ b/PWGMM/UE/Tasks/flattenicityTask.cxx @@ -14,30 +14,20 @@ /// \author Eisha Rani /// \since August 2026 -#include "PWGDQ/DataModel/ReducedInfoTables.h" - -#include "Common/DataModel/Centrality.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" -#include -#include -#include +#include #include #include +#include +#include #include -#include -#include -#include -#include -#include -#include -#include +#include #include +#include +#include #include #include @@ -50,7 +40,6 @@ using namespace o2::constants::physics; using FullTracks = soa::Join; struct FlattenicityTask { - // ============================================ // FT0 Constants // ============================================ @@ -69,7 +58,7 @@ struct FlattenicityTask { static constexpr int NPhysicalPrimaryBit = 0x4; // ============================================ - // Multiplicity class boundaries (Nch cuts defining the 0-10% ... 90-100% classes) + // Multiplicity class boundaries // ============================================ static constexpr int NchBound10 = 5; static constexpr int NchBound20 = 8; @@ -134,12 +123,12 @@ struct FlattenicityTask { int getCharge(int pdgCode) { switch (std::abs(pdgCode)) { - case 211: // pion - case 321: // kaon - case 2212: // proton + case PDG_t::kPiPlus: // 211 + case PDG_t::kKPlus: // 321 + case PDG_t::kProton: // 2212 return 1; - case 11: // electron - case 13: // muon + case PDG_t::kElectron: // 11 + case PDG_t::kMuonPlus: // 13 return -1; default: return 0; @@ -156,12 +145,14 @@ struct FlattenicityTask { total += counts[i]; } - if (total <= 0) + if (total <= 0) { return -1.0; + } float mean = total / NCell; - if (mean <= 0) + if (mean <= 0) { return -1.0; + } float sumSq = 0.0; for (int i = 0; i < NCell; i++) { @@ -181,8 +172,9 @@ struct FlattenicityTask { bool inFT0A = (eta > FT0AEtaMin && eta < FT0AEtaMax); bool inFT0C = (eta > FT0CEtaMin && eta < FT0CEtaMax); - if (!inFT0A && !inFT0C) + if (!inFT0A && !inFT0C) { return -1; + } isFT0A = inFT0A; @@ -216,32 +208,39 @@ struct FlattenicityTask { bool isSelectedTrack(const T& track) { // pT selection - if (track.pt() < cfgPtMin) + if (track.pt() < cfgPtMin) { return false; + } // Eta selection - if (std::abs(track.eta()) > cfgEtaMax) + if (std::abs(track.eta()) > cfgEtaMax) { return false; + } // TPC crossed rows - if (track.tpcNClsCrossedRows() < cfgNCrossedRowsTPC) + if (track.tpcNClsCrossedRows() < cfgNCrossedRowsTPC) { return false; + } // TPC chi2 per cluster - if (track.tpcChi2NCl() > cfgChi2PerClusterTPC) + if (track.tpcChi2NCl() > cfgChi2PerClusterTPC) { return false; + } // ITS chi2 per cluster - if (track.itsChi2NCl() > cfgChi2PerClusterITS) + if (track.itsChi2NCl() > cfgChi2PerClusterITS) { return false; + } // DCA z - if (std::abs(track.dcaZ()) > cfgDCAZ) + if (std::abs(track.dcaZ()) > cfgDCAZ) { return false; + } // Golden chi2 (global track) - if (cfgRequireGoldenChi2 && !track.isGlobalTrack()) + if (cfgRequireGoldenChi2 && !track.isGlobalTrack()) { return false; + } return true; } @@ -254,8 +253,7 @@ struct FlattenicityTask { aod::McParticles const& mcParticles) { // Initialize counters - std::array truthCounts; - truthCounts.fill(0.0); + std::array truthCounts{}; int nchINEL = 0; int nchFT0 = 0; @@ -265,17 +263,20 @@ struct FlattenicityTask { // Loop over MC particles for (const auto& particle : mcParticles) { // Check if primary - if (!(particle.flags() & NPhysicalPrimaryBit)) + if (!(particle.flags() & NPhysicalPrimaryBit)) { continue; + } // Check if charged int charge = getCharge(particle.pdgCode()); - if (charge == 0) + if (charge == 0) { continue; + } // pT > 0.1 - if (particle.pt() < cfgPtMin) + if (particle.pt() < cfgPtMin) { continue; + } // INEL>0: |eta| < 1 if (std::abs(particle.eta()) < 1.0) { @@ -291,10 +292,12 @@ struct FlattenicityTask { bool inFT0A = (particle.eta() > FT0AEtaMin && particle.eta() < FT0AEtaMax); bool inFT0C = (particle.eta() > FT0CEtaMin && particle.eta() < FT0CEtaMax); - if (inFT0A) + if (inFT0A) { hasFT0A = true; - if (inFT0C) + } + if (inFT0C) { hasFT0C = true; + } // Assign to FT0 cell bool isFT0A = false; @@ -370,8 +373,9 @@ struct FlattenicityTask { FullTracks const& tracks) { // Event selection: |vz| < 10 cm - if (std::abs(collision.posZ()) > cfgVzMax) + if (std::abs(collision.posZ()) > cfgVzMax) { return; + } // Find FT0 matching this collision's BC auto ft0 = ft0s.begin(); @@ -383,18 +387,17 @@ struct FlattenicityTask { break; } } - if (!foundFT0) + if (!foundFT0) { return; + } // Track selection and counting - // int nTracks = 0; - std::array recoCounts; - recoCounts.fill(0.0); + std::array recoCounts{}; for (const auto& track : tracks) { - if (!isSelectedTrack(track)) + if (!isSelectedTrack(track)) { continue; - // nTracks++; + } // Assign to FT0 cell using track extrapolation bool isFT0A = false; @@ -407,10 +410,7 @@ struct FlattenicityTask { histos.fill(HIST("hEvents"), 3); // Data events // Compute flattenicity from FT0 amplitudes - // The FT0 channels are stored as arrays of (channel, amplitude) pairs - // The channelA() and channelC() return vectors of pairs - std::array ft0Counts; - ft0Counts.fill(0.0); + std::array ft0Counts{}; // FT0-A channels (0-95) for (int i = 0; i < NchA; i++) { From fb1402ec3bdebbb741ae1c2a7a7b9a243ec75b0d Mon Sep 17 00:00:00 2001 From: PaolaVT Date: Sat, 8 Aug 2026 10:45:26 -0600 Subject: [PATCH 09/10] Fix code-check errors for flattenicityTask --- PWGMM/UE/Tasks/flattenicityTask.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/PWGMM/UE/Tasks/flattenicityTask.cxx b/PWGMM/UE/Tasks/flattenicityTask.cxx index 07c5a3c8712..e6a720c7d7b 100644 --- a/PWGMM/UE/Tasks/flattenicityTask.cxx +++ b/PWGMM/UE/Tasks/flattenicityTask.cxx @@ -21,6 +21,7 @@ #include #include #include +#include #include #include @@ -263,7 +264,7 @@ struct FlattenicityTask { // Loop over MC particles for (const auto& particle : mcParticles) { // Check if primary - if (!(particle.flags() & NPhysicalPrimaryBit)) { + if ((particle.flags() & NPhysicalPrimaryBit) == 0) { continue; } From 811df3b785b7fae39605fedc48f358be8e12f480 Mon Sep 17 00:00:00 2001 From: PaolaVT Date: Wed, 12 Aug 2026 19:40:31 -0600 Subject: [PATCH 10/10] Add single flattenicity (FT0-only) definition and update percentiles --- PWGMM/UE/Tasks/flattenicityTask.cxx | 281 +++++++++++++++++----------- 1 file changed, 168 insertions(+), 113 deletions(-) diff --git a/PWGMM/UE/Tasks/flattenicityTask.cxx b/PWGMM/UE/Tasks/flattenicityTask.cxx index e6a720c7d7b..d4507cd76b0 100644 --- a/PWGMM/UE/Tasks/flattenicityTask.cxx +++ b/PWGMM/UE/Tasks/flattenicityTask.cxx @@ -10,7 +10,10 @@ // or submit itself to any jurisdiction. /// \file flattenicityTask.cxx -/// \brief Flattenicity analysis task for UE studies +/// \brief Flattenicity analysis task for UE studies, with two flattenicity +/// definitions: particle-based (charged particles/tracks mapped into +/// FT0 cells) and FT0-detector-amplitude-based (real FT0 channel +/// signals), for Hyperloop /// \author Eisha Rani /// \since August 2026 @@ -29,8 +32,6 @@ #include #include #include -#include -#include using namespace o2; using namespace o2::framework; @@ -58,19 +59,6 @@ struct FlattenicityTask { static constexpr int NPhysicalPrimaryBit = 0x4; - // ============================================ - // Multiplicity class boundaries - // ============================================ - static constexpr int NchBound10 = 5; - static constexpr int NchBound20 = 8; - static constexpr int NchBound30 = 11; - static constexpr int NchBound40 = 14; - static constexpr int NchBound50 = 17; - static constexpr int NchBound60 = 20; - static constexpr int NchBound70 = 24; - static constexpr int NchBound80 = 28; - static constexpr int NchBound90 = 33; - // ============================================ // Histogram Definitions // ============================================ @@ -84,32 +72,44 @@ struct FlattenicityTask { {"hNch_INEL", "Nch distribution (INEL>0);N_{ch};Entries", {HistType::kTH1F, {{100, -0.5, 99.5}}}}, {"hNch_FT0", "Nch distribution (INEL>0 & FT0);N_{ch};Entries", {HistType::kTH1F, {{100, -0.5, 99.5}}}}, - // Flattenicity - {"hFlattenicity", "Flattenicity distribution;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, - {"hFlattenicity_vs_Nch", "Flattenicity vs Nch;N_{ch};1-#rho", {HistType::kTH2F, {{50, -0.5, 99.5}, {50, 0.0, 1.0}}}}, + // Definition 1: flattenicity from charged particles/tracks mapped into FT0 cells + {"hFlattenicityParticles", "Flattenicity from charged particles in FT0 acceptance;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityParticles_vs_Nch", "Flattenicity (particles) vs Nch;N_{ch};1-#rho", {HistType::kTH2F, {{50, -0.5, 99.5}, {50, 0.0, 1.0}}}}, - // FT0 cell occupancy + // Definition 2: flattenicity from FT0 detector amplitudes only + {"hFlattenicityFT0", "Flattenicity from FT0 detector amplitudes;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + + // FT0 cell occupancy (from particle-level mapping) {"hCellOccupancy", "FT0 cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{NCell, 0, NCell}}}}, {"hCellOccupancyFT0A", "FT0-A cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{NchA, 0, NchA}}}}, {"hCellOccupancyFT0C", "FT0-C cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{NchC, 0, NchC}}}}, - // Multiplicity classes - {"hFlattenicity_0_10", "Flattenicity 0-10%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, - {"hFlattenicity_10_20", "Flattenicity 10-20%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, - {"hFlattenicity_20_30", "Flattenicity 20-30%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, - {"hFlattenicity_30_40", "Flattenicity 30-40%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, - {"hFlattenicity_40_50", "Flattenicity 40-50%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, - {"hFlattenicity_50_60", "Flattenicity 50-60%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, - {"hFlattenicity_60_70", "Flattenicity 60-70%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, - {"hFlattenicity_70_80", "Flattenicity 70-80%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, - {"hFlattenicity_80_90", "Flattenicity 80-90%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, - {"hFlattenicity_90_100", "Flattenicity 90-100%;1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + // Percentile classes (I-VIII: 0-1,1-5,5-10,10-20,20-30,30-40,40-50,50-100%) + // -- particle-based definition + {"hFlattenicityParticles_0_1", "Flattenicity (particles) class I (0-1%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityParticles_1_5", "Flattenicity (particles) class II (1-5%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityParticles_5_10", "Flattenicity (particles) class III (5-10%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityParticles_10_20", "Flattenicity (particles) class IV (10-20%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityParticles_20_30", "Flattenicity (particles) class V (20-30%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityParticles_30_40", "Flattenicity (particles) class VI (30-40%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityParticles_40_50", "Flattenicity (particles) class VII (40-50%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityParticles_50_100", "Flattenicity (particles) class VIII (50-100%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + + // -- FT0-amplitude-based definition + {"hFlattenicityFT0_0_1", "Flattenicity (FT0) class I (0-1%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityFT0_1_5", "Flattenicity (FT0) class II (1-5%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityFT0_5_10", "Flattenicity (FT0) class III (5-10%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityFT0_10_20", "Flattenicity (FT0) class IV (10-20%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityFT0_20_30", "Flattenicity (FT0) class V (20-30%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityFT0_30_40", "Flattenicity (FT0) class VI (30-40%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityFT0_40_50", "Flattenicity (FT0) class VII (40-50%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, + {"hFlattenicityFT0_50_100", "Flattenicity (FT0) class VIII (50-100%);1-#rho;Entries", {HistType::kTH1F, {{50, 0.0, 1.0}}}}, }}; // ============================================ // Configurables // ============================================ - Configurable cfgPtMin{"cfgPtMin", 0.1, "Minimum pT for tracks"}; + Configurable cfgPtMin{"cfgPtMin", 0.0, "Minimum pT for tracks/particles (GeV/c)"}; Configurable cfgEtaMax{"cfgEtaMax", 0.8, "Maximum |eta| for tracks"}; Configurable cfgVzMax{"cfgVzMax", 10.0, "Maximum |vz| for collisions"}; Configurable cfgNCrossedRowsTPC{"cfgNCrossedRowsTPC", 70, "Minimum TPC crossed rows"}; @@ -118,6 +118,32 @@ struct FlattenicityTask { Configurable cfgDCAZ{"cfgDCAZ", 0.1, "Maximum DCA z"}; Configurable cfgRequireGoldenChi2{"cfgRequireGoldenChi2", true, "Require golden chi2"}; + // Percentile-class boundaries (in 1-rho), one full set of 7 edges per + // definition (classes I-VIII: 0-1,1-5,5-10,10-20,20-30,30-40,40-50,50-100%). + // Kept as individually named Configurables for Hyperloop -- these are NOT + // computed on the fly; they must be set from an actual percentile + // calibration of the corresponding 1-rho distribution (particle-based and + // FT0-amplitude-based distributions are different quantities and will not + // share the same boundary values). The defaults below are placeholders. + + // -- particle-based definition boundaries + Configurable cfgParticlesCut0To1{"cfgParticlesCut0To1", 0.90, "Particles 1-rho lower edge for 0-1%"}; + Configurable cfgParticlesCut1To5{"cfgParticlesCut1To5", 0.80, "Particles 1-rho lower edge for 1-5%"}; + Configurable cfgParticlesCut5To10{"cfgParticlesCut5To10", 0.75, "Particles 1-rho lower edge for 5-10%"}; + Configurable cfgParticlesCut10To20{"cfgParticlesCut10To20", 0.65, "Particles 1-rho lower edge for 10-20%"}; + Configurable cfgParticlesCut20To30{"cfgParticlesCut20To30", 0.55, "Particles 1-rho lower edge for 20-30%"}; + Configurable cfgParticlesCut30To40{"cfgParticlesCut30To40", 0.45, "Particles 1-rho lower edge for 30-40%"}; + Configurable cfgParticlesCut40To50{"cfgParticlesCut40To50", 0.35, "Particles 1-rho lower edge for 40-50%"}; + + // -- FT0-amplitude-based definition boundaries + Configurable cfgFT0Cut0To1{"cfgFT0Cut0To1", 0.904, "FT0 1-rho lower edge for 0-1%"}; + Configurable cfgFT0Cut1To5{"cfgFT0Cut1To5", 0.888, "FT0 1-rho lower edge for 1-5%"}; + Configurable cfgFT0Cut5To10{"cfgFT0Cut5To10", 0.840, "FT0 1-rho lower edge for 5-10%"}; + Configurable cfgFT0Cut10To20{"cfgFT0Cut10To20", 0.780, "FT0 1-rho lower edge for 10-20%"}; + Configurable cfgFT0Cut20To30{"cfgFT0Cut20To30", 0.720, "FT0 1-rho lower edge for 20-30%"}; + Configurable cfgFT0Cut30To40{"cfgFT0Cut30To40", 0.660, "FT0 1-rho lower edge for 30-40%"}; + Configurable cfgFT0Cut40To50{"cfgFT0Cut40To50", 0.600, "FT0 1-rho lower edge for 40-50%"}; + // ============================================ // Particle charge function // ============================================ @@ -137,69 +163,119 @@ struct FlattenicityTask { } // ============================================ - // Flattenicity calculation + // Generic ALICE flattenicity formula, shared by both definitions // ============================================ - float computeFlattenicity(const std::array& counts) + float computeRhoFromCells(const std::array& counts) { - float total = 0.0; - for (int i = 0; i < NCell; i++) { + float total = 0.0f; + for (int i = 0; i < NCell; ++i) { total += counts[i]; } - - if (total <= 0) { - return -1.0; + if (total <= 0.0f) { + return -1.0f; } float mean = total / NCell; - if (mean <= 0) { - return -1.0; + if (mean <= 0.0f) { + return -1.0f; } - float sumSq = 0.0; - for (int i = 0; i < NCell; i++) { + float sumSq = 0.0f; + for (int i = 0; i < NCell; ++i) { sumSq += (counts[i] - mean) * (counts[i] - mean); } - float rho = std::sqrt(sumSq) / (NCell * mean); - return rho; + return std::sqrt(sumSq) / (NCell * mean); + } + + // Definition 1: charged-particle/track flattenicity in FT0 acceptance. + float computeFlattenicityParticles(const std::array& particleCounts) + { + float rho = computeRhoFromCells(particleCounts); + return (rho >= 0.0f) ? (1.0f - rho) : -1.0f; + } + + // Definition 2: FT0-detector-only flattenicity from amplitudes. + float computeFlattenicityFT0(const std::array& ft0Counts) + { + float rho = computeRhoFromCells(ft0Counts); + return (rho >= 0.0f) ? (1.0f - rho) : -1.0f; } // ============================================ - // Assign particle to FT0 cell + // Assign particle/track to FT0 cell // ============================================ int assignToFT0Cell(float eta, float phi, bool& isFT0A) { - // Check if in FT0 acceptance bool inFT0A = (eta > FT0AEtaMin && eta < FT0AEtaMax); bool inFT0C = (eta > FT0CEtaMin && eta < FT0CEtaMax); - if (!inFT0A && !inFT0C) { return -1; } isFT0A = inFT0A; - // Phi bin int phiBin = static_cast(std::floor(phi / (o2::constants::math::TwoPI / NPhiSectors))); phiBin = std::max(0, std::min(phiBin, NPhiSectors - 1)); - int cellId = -1; - if (inFT0A) { // FT0-A: cells 0-95 float etaWidth = (FT0AEtaMax - FT0AEtaMin) / NEtaA; int etaBin = static_cast(std::floor((eta - FT0AEtaMin) / etaWidth)); etaBin = std::max(0, std::min(etaBin, NEtaA - 1)); - cellId = etaBin * NPhiSectors + phiBin; - } else if (inFT0C) { - // FT0-C: cells 96-207 - float etaWidth = (FT0CEtaMax - FT0CEtaMin) / NEtaC; - int etaBin = static_cast(std::floor((eta - FT0CEtaMin) / etaWidth)); - etaBin = std::max(0, std::min(etaBin, NEtaC - 1)); - cellId = NchA + etaBin * NPhiSectors + phiBin; + return etaBin * NPhiSectors + phiBin; + } + + // FT0-C: cells 96-207 + float etaWidth = (FT0CEtaMax - FT0CEtaMin) / NEtaC; + int etaBin = static_cast(std::floor((eta - FT0CEtaMin) / etaWidth)); + etaBin = std::max(0, std::min(etaBin, NEtaC - 1)); + return NchA + etaBin * NPhiSectors + phiBin; + } + + // ============================================ + // Percentile-class filling helpers + // ============================================ + void fillParticlesPercentileHistograms(float flat) + { + if (flat >= cfgParticlesCut0To1) { + histos.fill(HIST("hFlattenicityParticles_0_1"), flat); + } else if (flat >= cfgParticlesCut1To5) { + histos.fill(HIST("hFlattenicityParticles_1_5"), flat); + } else if (flat >= cfgParticlesCut5To10) { + histos.fill(HIST("hFlattenicityParticles_5_10"), flat); + } else if (flat >= cfgParticlesCut10To20) { + histos.fill(HIST("hFlattenicityParticles_10_20"), flat); + } else if (flat >= cfgParticlesCut20To30) { + histos.fill(HIST("hFlattenicityParticles_20_30"), flat); + } else if (flat >= cfgParticlesCut30To40) { + histos.fill(HIST("hFlattenicityParticles_30_40"), flat); + } else if (flat >= cfgParticlesCut40To50) { + histos.fill(HIST("hFlattenicityParticles_40_50"), flat); + } else { + histos.fill(HIST("hFlattenicityParticles_50_100"), flat); } + } - return cellId; + void fillFT0PercentileHistograms(float flat) + { + if (flat >= cfgFT0Cut0To1) { + histos.fill(HIST("hFlattenicityFT0_0_1"), flat); + } else if (flat >= cfgFT0Cut1To5) { + histos.fill(HIST("hFlattenicityFT0_1_5"), flat); + } else if (flat >= cfgFT0Cut5To10) { + histos.fill(HIST("hFlattenicityFT0_5_10"), flat); + } else if (flat >= cfgFT0Cut10To20) { + histos.fill(HIST("hFlattenicityFT0_10_20"), flat); + } else if (flat >= cfgFT0Cut20To30) { + histos.fill(HIST("hFlattenicityFT0_20_30"), flat); + } else if (flat >= cfgFT0Cut30To40) { + histos.fill(HIST("hFlattenicityFT0_30_40"), flat); + } else if (flat >= cfgFT0Cut40To50) { + histos.fill(HIST("hFlattenicityFT0_40_50"), flat); + } else { + histos.fill(HIST("hFlattenicityFT0_50_100"), flat); + } } // ============================================ @@ -208,8 +284,8 @@ struct FlattenicityTask { template bool isSelectedTrack(const T& track) { - // pT selection - if (track.pt() < cfgPtMin) { + // pT > 0 + if (track.pt() <= 0.0f || track.pt() < cfgPtMin) { return false; } @@ -248,13 +324,14 @@ struct FlattenicityTask { // ============================================ // Process MC collisions + // (only the particle-based definition is available here: there is no + // FT0 detector object at MC-truth level in this process signature) // ============================================ void processMC( aod::McCollision const& /* mcCollision */, aod::McParticles const& mcParticles) { - // Initialize counters - std::array truthCounts{}; + std::array particleCounts{}; int nchINEL = 0; int nchFT0 = 0; @@ -274,8 +351,8 @@ struct FlattenicityTask { continue; } - // pT > 0.1 - if (particle.pt() < cfgPtMin) { + // pT > 0 + if (particle.pt() <= 0.0f || particle.pt() < cfgPtMin) { continue; } @@ -305,7 +382,7 @@ struct FlattenicityTask { int cellId = assignToFT0Cell(particle.eta(), particle.phi(), isFT0A); if (cellId >= 0 && cellId < NCell) { - truthCounts[cellId] += 1.0; + particleCounts[cellId] += 1.0f; histos.fill(HIST("hCellOccupancy"), cellId); if (isFT0A) { histos.fill(HIST("hCellOccupancyFT0A"), cellId); @@ -330,35 +407,11 @@ struct FlattenicityTask { histos.fill(HIST("hEvents"), 2); // INEL>0 & FT0 histos.fill(HIST("hNch_FT0"), nchFT0); - // Compute flattenicity - float rho = computeFlattenicity(truthCounts); - if (rho > 0) { - float flattenicity = 1.0 - rho; - histos.fill(HIST("hFlattenicity"), flattenicity); - histos.fill(HIST("hFlattenicity_vs_Nch"), nchFT0, flattenicity); - - // Multiplicity classes (based on Nch) - if (nchFT0 < NchBound10) { - histos.fill(HIST("hFlattenicity_0_10"), flattenicity); - } else if (nchFT0 < NchBound20) { - histos.fill(HIST("hFlattenicity_10_20"), flattenicity); - } else if (nchFT0 < NchBound30) { - histos.fill(HIST("hFlattenicity_20_30"), flattenicity); - } else if (nchFT0 < NchBound40) { - histos.fill(HIST("hFlattenicity_30_40"), flattenicity); - } else if (nchFT0 < NchBound50) { - histos.fill(HIST("hFlattenicity_40_50"), flattenicity); - } else if (nchFT0 < NchBound60) { - histos.fill(HIST("hFlattenicity_50_60"), flattenicity); - } else if (nchFT0 < NchBound70) { - histos.fill(HIST("hFlattenicity_60_70"), flattenicity); - } else if (nchFT0 < NchBound80) { - histos.fill(HIST("hFlattenicity_70_80"), flattenicity); - } else if (nchFT0 < NchBound90) { - histos.fill(HIST("hFlattenicity_80_90"), flattenicity); - } else { - histos.fill(HIST("hFlattenicity_90_100"), flattenicity); - } + float flatParticles = computeFlattenicityParticles(particleCounts); + if (flatParticles >= 0.0f) { + histos.fill(HIST("hFlattenicityParticles"), flatParticles); + histos.fill(HIST("hFlattenicityParticles_vs_Nch"), nchFT0, flatParticles); + fillParticlesPercentileHistograms(flatParticles); } } } @@ -367,6 +420,8 @@ struct FlattenicityTask { // ============================================ // Process data collisions + // (both definitions are available here: particle-based from selected + // reconstructed tracks, and FT0-amplitude-based from the real FT0 signal) // ============================================ void processData( aod::Collision const& collision, @@ -392,40 +447,40 @@ struct FlattenicityTask { return; } - // Track selection and counting - std::array recoCounts{}; - + // Definition 1: selected charged tracks mapped into FT0 cells + std::array recoParticleCounts{}; for (const auto& track : tracks) { if (!isSelectedTrack(track)) { continue; } - - // Assign to FT0 cell using track extrapolation bool isFT0A = false; int cellId = assignToFT0Cell(track.eta(), track.phi(), isFT0A); if (cellId >= 0 && cellId < NCell) { - recoCounts[cellId] += 1.0; + recoParticleCounts[cellId] += 1.0f; } } - histos.fill(HIST("hEvents"), 3); // Data events - - // Compute flattenicity from FT0 amplitudes + // Definition 2: FT0 detector amplitudes only std::array ft0Counts{}; - - // FT0-A channels (0-95) - for (int i = 0; i < NchA; i++) { + for (int i = 0; i < NchA; ++i) { ft0Counts[i] = ft0.channelA()[i]; } - - // FT0-C channels (96-207) - for (int i = 0; i < NchC; i++) { + for (int i = 0; i < NchC; ++i) { ft0Counts[NchA + i] = ft0.channelC()[i]; } - float rho = computeFlattenicity(ft0Counts); - if (rho > 0) { - histos.fill(HIST("hFlattenicity"), 1.0 - rho); + histos.fill(HIST("hEvents"), 3); // Data events + + float flatParticles = computeFlattenicityParticles(recoParticleCounts); + if (flatParticles >= 0.0f) { + histos.fill(HIST("hFlattenicityParticles"), flatParticles); + fillParticlesPercentileHistograms(flatParticles); + } + + float flatFT0 = computeFlattenicityFT0(ft0Counts); + if (flatFT0 >= 0.0f) { + histos.fill(HIST("hFlattenicityFT0"), flatFT0); + fillFT0PercentileHistograms(flatFT0); } }