From da7c51b98d41c6a5f9daa86b1ee9665f3eb2a09c Mon Sep 17 00:00:00 2001 From: PaolaVT Date: Thu, 6 Aug 2026 10:15:31 -0600 Subject: [PATCH 1/7] 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 2/7] 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 3/7] 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 4/7] 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 5/7] 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 6/7] 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 7/7] 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