|
| 1 | +//============================================================================= |
| 2 | +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. |
| 3 | +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. |
| 4 | +// All rights not expressly granted are reserved. |
| 5 | +// |
| 6 | +// This software is distributed under the terms of the GNU General Public |
| 7 | +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". |
| 8 | +// |
| 9 | +// In applying this license CERN does not waive the privileges and immunities |
| 10 | +// granted to it by virtue of its status as an Intergovernmental Organization |
| 11 | +// or submit itself to any jurisdiction. |
| 12 | +/// |
| 13 | +/// \author Eishah Rani |
| 14 | +/// \since August 5, 2026 |
| 15 | +/// \file flattenicityTask.cxx |
| 16 | +/// \brief Generator-Level vs Reconstruction-Level Flattenicity |
| 17 | +// Based on: |
| 18 | +// PWGMM/UE/Tasks/uecharged.cxx |
| 19 | +//============================================================================= |
| 20 | + |
| 21 | +#include "Common/Core/TrackSelection.h" |
| 22 | +#include "Common/Core/TrackSelectionDefaults.h" |
| 23 | +#include "Common/DataModel/EventSelection.h" |
| 24 | +#include "Common/DataModel/Multiplicity.h" |
| 25 | +#include "Common/DataModel/TrackSelectionTables.h" |
| 26 | + |
| 27 | +#include "Framework/ASoAHelpers.h" |
| 28 | +#include "Framework/AnalysisDataModel.h" |
| 29 | +#include "Framework/AnalysisTask.h" |
| 30 | +#include "Framework/Configurable.h" |
| 31 | +#include "Framework/Constants.h" |
| 32 | +#include "Framework/HistogramRegistry.h" |
| 33 | +#include "Framework/InitContext.h" |
| 34 | +#include "Framework/O2DatabasePDGPlugin.h" |
| 35 | +#include "Framework/runDataProcessing.h" |
| 36 | +#include "ReconstructionDataFormats/Track.h" |
| 37 | + |
| 38 | +#include <cmath> |
| 39 | +#include <vector> |
| 40 | + |
| 41 | +using namespace o2; |
| 42 | +using namespace o2::framework; |
| 43 | +using namespace o2::framework::expressions; |
| 44 | + |
| 45 | +struct FlattenicityTask { |
| 46 | + |
| 47 | + // --- Flattenicity constants --- |
| 48 | + static constexpr int NCH_A = 96; |
| 49 | + static constexpr int NCH_C = 112; |
| 50 | + static constexpr int NCELL = NCH_A + NCH_C; |
| 51 | + static constexpr int NPHISECTORS = 8; |
| 52 | + static constexpr int NETA_A = NCH_A / NPHISECTORS; |
| 53 | + static constexpr int NETA_C = NCH_C / NPHISECTORS; |
| 54 | + |
| 55 | + // FT0 acceptance |
| 56 | + static constexpr float FT0A_ETA_MIN = 3.5f; |
| 57 | + static constexpr float FT0A_ETA_MAX = 4.9f; |
| 58 | + static constexpr float FT0C_ETA_MIN = -3.3f; |
| 59 | + static constexpr float FT0C_ETA_MAX = -2.1f; |
| 60 | + |
| 61 | + // --- Event selection constants --- |
| 62 | + static constexpr float VERTEX_CUT = 10.0f; |
| 63 | + static constexpr float FLAT_MIN = 0.0f; |
| 64 | + static constexpr float INEL_ETA_CUT = 1.0f; |
| 65 | + static constexpr float MIDRAP_ETA_CUT = 0.8f; |
| 66 | + |
| 67 | + // --- Configurables --- |
| 68 | + Configurable<float> cfgTrkEtaCut{"cfgTrkEtaCut", 0.8f, "Eta range for tracks"}; |
| 69 | + Configurable<float> cfgTrkLowPtCut{"cfgTrkLowPtCut", 0.15f, "Minimum pT"}; |
| 70 | + |
| 71 | + Configurable<bool> isRun3{"isRun3", true, "is Run3 dataset"}; |
| 72 | + Configurable<bool> timeEvsel{"timeEvsel", true, "TPC Time frame boundary cut"}; |
| 73 | + Configurable<bool> piluprejection{"piluprejection", true, "Pileup rejection"}; |
| 74 | + Configurable<bool> goodzvertex{"goodzvertex", true, "Good Z vertex"}; |
| 75 | + |
| 76 | + // --- Track selection --- |
| 77 | + TrackSelection mySelectionPrim; |
| 78 | + |
| 79 | + // --- Histograms --- |
| 80 | + HistogramRegistry registry; |
| 81 | + |
| 82 | + // --- Init --- |
| 83 | + void init(InitContext const&) override |
| 84 | + { |
| 85 | + // Initialize track selection |
| 86 | + mySelectionPrim = myTrackSelectionPrim(); |
| 87 | + |
| 88 | + // Define histograms |
| 89 | + AxisSpec flatBins = {40, 0.0, 1.0, "#rho"}; |
| 90 | + AxisSpec nchBins = {100, -0.5, 99.5, "N_{ch}"}; |
| 91 | + |
| 92 | + registry.add("hFlattenicityTruth", "Truth flattenicity; 1-#rho; Events", |
| 93 | + HistType::kTH1D, {flatBins}); |
| 94 | + registry.add("hFlattenicityReco", "Reco flattenicity; 1-#rho; Events", |
| 95 | + HistType::kTH1D, {flatBins}); |
| 96 | + registry.add("hFlattenicityCorrelation", "Truth vs Reco; 1-#rho_{truth}; 1-#rho_{reco}", |
| 97 | + HistType::kTH2D, {flatBins, flatBins}); |
| 98 | + registry.add("hNch", "Reco Nch distribution; N_{ch}; Events", |
| 99 | + HistType::kTH1D, {nchBins}); |
| 100 | + registry.add("hNchTruth", "Truth Nch distribution; N_{ch}; Events", |
| 101 | + HistType::kTH1D, {nchBins}); |
| 102 | + } |
| 103 | + |
| 104 | + // --- Track selection function --- |
| 105 | + TrackSelection myTrackSelectionPrim() |
| 106 | + { |
| 107 | + TrackSelection selectedTracks; |
| 108 | + selectedTracks.SetPtRange(0.1f, 1e10f); |
| 109 | + selectedTracks.SetEtaRange(-0.8f, 0.8f); |
| 110 | + selectedTracks.SetRequireITSRefit(true); |
| 111 | + selectedTracks.SetRequireTPCRefit(true); |
| 112 | + selectedTracks.SetMinNCrossedRowsTPC(70); |
| 113 | + selectedTracks.SetMinNCrossedRowsOverFindableClustersTPC(0.4f); |
| 114 | + selectedTracks.SetMaxChi2PerClusterTPC(4.0f); |
| 115 | + selectedTracks.SetRequireHitsInITSLayers(1, {0, 1}); |
| 116 | + selectedTracks.SetMaxChi2PerClusterITS(36.0f); |
| 117 | + selectedTracks.SetMaxDcaXYPtDep([](float pt) { return 0.0105f + 0.0350f / std::pow(pt, 1.1f); }); |
| 118 | + selectedTracks.SetMaxDcaZ(2.0f); |
| 119 | + return selectedTracks; |
| 120 | + } |
| 121 | + |
| 122 | + // --- Helper: Get cell ID for a particle in FT0 acceptance --- |
| 123 | + int getCellId(float eta, float phi) |
| 124 | + { |
| 125 | + // Check if in FT0-A acceptance |
| 126 | + if (eta > FT0A_ETA_MIN && eta < FT0A_ETA_MAX) { |
| 127 | + int phiBin = static_cast<int>(std::floor(phi / (2.0f * o2::constants::math::PI / static_cast<float>(NPHISECTORS)))); |
| 128 | + phiBin = std::clamp(phiBin, 0, NPHISECTORS - 1); |
| 129 | + int etaBin = static_cast<int>(std::floor((eta - FT0A_ETA_MIN) / ((FT0A_ETA_MAX - FT0A_ETA_MIN) / static_cast<float>(NETA_A)))); |
| 130 | + etaBin = std::clamp(etaBin, 0, NETA_A - 1); |
| 131 | + return etaBin * NPHISECTORS + phiBin; |
| 132 | + } |
| 133 | + |
| 134 | + // Check if in FT0-C acceptance |
| 135 | + if (eta > FT0C_ETA_MIN && eta < FT0C_ETA_MAX) { |
| 136 | + int phiBin = static_cast<int>(std::floor(phi / (2.0f * o2::constants::math::PI / static_cast<float>(NPHISECTORS)))); |
| 137 | + phiBin = std::clamp(phiBin, 0, NPHISECTORS - 1); |
| 138 | + int etaBin = static_cast<int>(std::floor((eta - FT0C_ETA_MIN) / ((FT0C_ETA_MAX - FT0C_ETA_MIN) / static_cast<float>(NETA_C)))); |
| 139 | + etaBin = std::clamp(etaBin, 0, NETA_C - 1); |
| 140 | + return NCH_A + etaBin * NPHISECTORS + phiBin; |
| 141 | + } |
| 142 | + |
| 143 | + return -1; // Not in FT0 acceptance |
| 144 | + } |
| 145 | + |
| 146 | + // --- Flattenicity calculation --- |
| 147 | + float calculateFlattenicity(const std::vector<float>& counts) |
| 148 | + { |
| 149 | + if (counts.size() != static_cast<size_t>(NCELL)) { |
| 150 | + return -1.0f; |
| 151 | + } |
| 152 | + |
| 153 | + float total = 0.0f; |
| 154 | + for (const auto& c : counts) { |
| 155 | + total += c; |
| 156 | + } |
| 157 | + if (total <= 0.0f) { |
| 158 | + return -1.0f; |
| 159 | + } |
| 160 | + |
| 161 | + float mean = total / static_cast<float>(NCELL); |
| 162 | + if (mean <= 0.0f) { |
| 163 | + return -1.0f; |
| 164 | + } |
| 165 | + |
| 166 | + float sumSq = 0.0f; |
| 167 | + for (const auto& c : counts) { |
| 168 | + sumSq += (c - mean) * (c - mean); |
| 169 | + } |
| 170 | + |
| 171 | + float rho = std::sqrt(sumSq / (static_cast<float>(NCELL) * static_cast<float>(NCELL))) / mean; |
| 172 | + return 1.0f - rho; |
| 173 | + } |
| 174 | + |
| 175 | + // --- Process Data --- |
| 176 | + void processData(aod::Collision const& collision, |
| 177 | + soa::Filtered<aod::Tracks> const& tracks, |
| 178 | + aod::FT0s const& ft0s) |
| 179 | + { |
| 180 | + // Event selection (Paola/Jesus) |
| 181 | + if (!collision.sel8()) { |
| 182 | + return; |
| 183 | + } |
| 184 | + if (std::abs(collision.posZ()) >= VERTEX_CUT) { |
| 185 | + return; |
| 186 | + } |
| 187 | + |
| 188 | + // Track loop for Nch |
| 189 | + int nch = 0; |
| 190 | + for (const auto& track : tracks) { |
| 191 | + if (!mySelectionPrim.IsSelected(track)) { |
| 192 | + continue; |
| 193 | + } |
| 194 | + nch++; |
| 195 | + } |
| 196 | + registry.fill(HIST("hNch"), nch); |
| 197 | + |
| 198 | + // FT0 flattenicity |
| 199 | + auto ft0 = collision.ft0(); |
| 200 | + if (ft0.hasAmplitudeA() && ft0.hasAmplitudeC()) { |
| 201 | + auto ampA = ft0.amplitudeA(); |
| 202 | + auto ampC = ft0.amplitudeC(); |
| 203 | + |
| 204 | + std::vector<float> counts(NCELL, 0.0f); |
| 205 | + for (int i = 0; i < static_cast<int>(ampA.size()) && i < NCH_A; ++i) { |
| 206 | + counts[i] = ampA[i]; |
| 207 | + } |
| 208 | + for (int i = 0; i < static_cast<int>(ampC.size()) && i < NCH_C; ++i) { |
| 209 | + counts[NCH_A + i] = ampC[i]; |
| 210 | + } |
| 211 | + |
| 212 | + float flat = calculateFlattenicity(counts); |
| 213 | + if (flat >= FLAT_MIN) { |
| 214 | + registry.fill(HIST("hFlattenicityReco"), flat); |
| 215 | + } |
| 216 | + } |
| 217 | + } |
| 218 | + PROCESS_SWITCH(FlattenicityTask, processData, "Process data", true); |
| 219 | + |
| 220 | + // --- Process MC --- |
| 221 | + void processMC(aod::McCollision const& mcCollision, |
| 222 | + aod::McParticles const& particles, |
| 223 | + soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions>> const& collisions, |
| 224 | + aod::FT0s const& ft0s, |
| 225 | + aod::BCs const& /*bcs*/) |
| 226 | + { |
| 227 | + // ---- Truth-level processing ---- |
| 228 | + bool inel = false; |
| 229 | + int nchTruth = 0; |
| 230 | + std::vector<float> truthCounts(NCELL, 0.0f); |
| 231 | + |
| 232 | + for (const auto& particle : particles) { |
| 233 | + // Check if physical primary |
| 234 | + if (!particle.isPhysicalPrimary()) { |
| 235 | + continue; |
| 236 | + } |
| 237 | + |
| 238 | + // Check if charged |
| 239 | + if (std::abs(particle.pdgCode()) == 0) { |
| 240 | + continue; |
| 241 | + } |
| 242 | + |
| 243 | + // Check pT > 0 |
| 244 | + if (particle.pt() <= 0.0f) { |
| 245 | + continue; |
| 246 | + } |
| 247 | + |
| 248 | + // INEL>0 check: primary charged with |eta| < INEL_ETA_CUT |
| 249 | + if (std::abs(particle.eta()) < INEL_ETA_CUT) { |
| 250 | + inel = true; |
| 251 | + } |
| 252 | + |
| 253 | + // Nch at midrapidity: |eta| < MIDRAP_ETA_CUT, pT > cfgTrkLowPtCut |
| 254 | + if (std::abs(particle.eta()) < MIDRAP_ETA_CUT && particle.pt() > cfgTrkLowPtCut) { |
| 255 | + nchTruth++; |
| 256 | + } |
| 257 | + |
| 258 | + // Flattenicity: particles in FT0 acceptance |
| 259 | + int cellId = getCellId(particle.eta(), particle.phi()); |
| 260 | + if (cellId >= 0 && cellId < NCELL) { |
| 261 | + truthCounts[cellId] += 1.0f; |
| 262 | + } |
| 263 | + } |
| 264 | + |
| 265 | + // Apply truth-level event selection (Paola/Jesus) |
| 266 | + if (!inel) { |
| 267 | + return; |
| 268 | + } |
| 269 | + if (std::abs(mcCollision.posZ()) >= VERTEX_CUT) { |
| 270 | + return; |
| 271 | + } |
| 272 | + |
| 273 | + // Fill truth multiplicity |
| 274 | + registry.fill(HIST("hNchTruth"), nchTruth); |
| 275 | + registry.fill(HIST("hNch"), nchTruth); |
| 276 | + |
| 277 | + // Calculate truth flattenicity |
| 278 | + float truthFlat = calculateFlattenicity(truthCounts); |
| 279 | + if (truthFlat >= FLAT_MIN) { |
| 280 | + registry.fill(HIST("hFlattenicityTruth"), truthFlat); |
| 281 | + } |
| 282 | + |
| 283 | + // ---- Reconstructed-level processing for matched collisions ---- |
| 284 | + for (const auto& collision : collisions) { |
| 285 | + // Apply reconstruction-level event selection |
| 286 | + if (!collision.sel8()) { |
| 287 | + continue; |
| 288 | + } |
| 289 | + if (std::abs(collision.posZ()) >= VERTEX_CUT) { |
| 290 | + continue; |
| 291 | + } |
| 292 | + |
| 293 | + // Get FT0 flattenicity for this collision |
| 294 | + auto ft0 = collision.ft0(); |
| 295 | + if (!ft0.hasAmplitudeA() || !ft0.hasAmplitudeC()) { |
| 296 | + continue; |
| 297 | + } |
| 298 | + |
| 299 | + auto ampA = ft0.amplitudeA(); |
| 300 | + auto ampC = ft0.amplitudeC(); |
| 301 | + |
| 302 | + std::vector<float> recoCounts(NCELL, 0.0f); |
| 303 | + for (int i = 0; i < static_cast<int>(ampA.size()) && i < NCH_A; ++i) { |
| 304 | + recoCounts[i] = ampA[i]; |
| 305 | + } |
| 306 | + for (int i = 0; i < static_cast<int>(ampC.size()) && i < NCH_C; ++i) { |
| 307 | + recoCounts[NCH_A + i] = ampC[i]; |
| 308 | + } |
| 309 | + |
| 310 | + float recoFlat = calculateFlattenicity(recoCounts); |
| 311 | + if (recoFlat >= FLAT_MIN && truthFlat >= FLAT_MIN) { |
| 312 | + registry.fill(HIST("hFlattenicityReco"), recoFlat); |
| 313 | + registry.fill(HIST("hFlattenicityCorrelation"), truthFlat, recoFlat); |
| 314 | + } |
| 315 | + } |
| 316 | + } |
| 317 | + PROCESS_SWITCH(FlattenicityTask, processMC, "Process MC", true); |
| 318 | +}; |
| 319 | + |
| 320 | +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) |
| 321 | +{ |
| 322 | + WorkflowSpec workflow{}; |
| 323 | + workflow.push_back(adaptAnalysisTask<FlattenicityTask>(cfgc)); |
| 324 | + return workflow; |
| 325 | +} |
0 commit comments