diff --git a/PWGMM/UE/Tasks/flattenicityTask.cxx b/PWGMM/UE/Tasks/flattenicityTask.cxx index f279984cb9f..875cef702a9 100644 --- a/PWGMM/UE/Tasks/flattenicityTask.cxx +++ b/PWGMM/UE/Tasks/flattenicityTask.cxx @@ -15,6 +15,7 @@ /// \since August 2026 #include "Common/DataModel/Centrality.h" +#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/TrackSelectionTables.h" #include @@ -30,8 +31,6 @@ #include #include #include -#include -#include #include #include @@ -47,14 +46,23 @@ using CollisionsWithCentAndMcLabel = soa::Join= ChannelsPerSector * i_sec && i_ch <= (ChannelsPerSector - 1) + ChannelsPerSector * i_sec) { + return i_sec; + } + } + return -1; + } + + int getT0CSector(int i_ch) + { + for (int i_sec = 0; i_sec < NSectorsC; ++i_sec) { + if (i_ch >= ChannelsPerSector * i_sec && i_ch <= (ChannelsPerSector - 1) + ChannelsPerSector * i_sec) { + return i_sec; + } + } + return -1; + } + // ============================================ // Histogram Definitions - 100 BINS (bin width = 0.01) // ============================================ @@ -96,7 +127,9 @@ struct FlattenicityTask { // ============================================ {"hFlattenicityParticles", "Flattenicity from charged particles;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, {"hFlattenicityParticles_vs_Nch", "Flattenicity (particles) vs Nch;N_{ch};1-#rho", {HistType::kTH2F, {{50, -0.5, 99.5}, {100, 0.0, 1.0}}}}, - {"hFlattenicityFT0", "Flattenicity from FT0 detector amplitudes;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0", "Flattenicity from FT0 detector amplitudes (avg of FT0-A and FT0-C);1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0A", "Flattenicity from FT0-A amplitudes (24 sectors);1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0C", "Flattenicity from FT0-C amplitudes (28 sectors);1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, // FT0 cell occupancy {"hCellOccupancy", "FT0 cell occupancy;Cell ID;Entries", {HistType::kTH1F, {{NCell, 0, NCell}}}}, @@ -106,7 +139,6 @@ struct FlattenicityTask { // ============================================ // MULTIPLICITY CLASSES USING PERCENTILES // All classes have 100 bins (bin width = 0.01) - // Following Antonio's publication style // ============================================ // 0-1% {"hFlattenicityParticles_0_1", "Flattenicity (particles) class 0-1%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, @@ -149,10 +181,7 @@ struct FlattenicityTask { {"hFlattenicityFT0_95_100", "Flattenicity (FT0) class 95-100%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, // ============================================ - // GEN vs RECO CORRELATION - // Antonio's original request: compare flattenicity at - // generator (MC truth) level vs reconstruction level, - // for the same collision. + // GEN vs RECO CORRELATION (3 plots) // ============================================ {"hFlatParticles_Gen_vs_Rec", "Flattenicity (particles): gen vs reco;1-#rho (gen);1-#rho (reco)", {HistType::kTH2F, {{100, 0.0, 1.0}, {100, 0.0, 1.0}}}}, {"hFlatFT0_Gen_vs_Rec", "Flattenicity: gen (particles) vs reco (FT0 amplitudes);1-#rho (gen);1-#rho (reco FT0)", {HistType::kTH2F, {{100, 0.0, 1.0}, {100, 0.0, 1.0}}}}, @@ -191,67 +220,63 @@ struct FlattenicityTask { // ============================================ // Flattenicity calculation + // Templated on array size so it can be reused for the 208-cell + // particle-level grid as well as the 24-cell (FT0-A) and + // 28-cell (FT0-C) detector-amplitude arrays. // ============================================ - float computeFlattenicity(const std::array& counts) + template + float computeFlattenicity(const std::array& counts) { float total = 0.0; - for (int i = 0; i < NCell; i++) { + for (std::size_t i = 0; i < N; i++) { total += counts[i]; } - if (total <= 0) { return -1.0; } - - float mean = total / NCell; + float mean = total / N; if (mean <= 0) { return -1.0; } - float sumSq = 0.0; - for (int i = 0; i < NCell; i++) { + for (std::size_t i = 0; i < N; i++) { sumSq += (counts[i] - mean) * (counts[i] - mean); } - - float rho = std::sqrt(sumSq) / (NCell * mean); + float rho = std::sqrt(sumSq) / (N * mean); return rho; } // ============================================ - // Assign particle to FT0 cell + // Assign particle to FT0 cell (Updated to match multFilter.cxx) // ============================================ 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 + // Map phi to sector (0-7) 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 + // FT0-A: NEtaA eta bins x NPhiSectors phi bins = NchA cells (0..NchA-1) 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; + int sector = static_cast(std::floor((eta - FT0AEtaMin) / etaWidth)); + sector = std::max(0, std::min(sector, NEtaA - 1)); + cellId = sector * NPhiSectors + phiBin; } else if (inFT0C) { - // FT0-C: cells 96-207 + // FT0-C: NEtaC eta bins x NPhiSectors phi bins = NchC cells, + // offset by NchA so FT0-A and FT0-C cell IDs never overlap. 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; + int sector = static_cast(std::floor((eta - FT0CEtaMin) / etaWidth)); + sector = std::max(0, std::min(sector, NEtaC - 1)); + cellId = NchA + sector * NPhiSectors + phiBin; } - return cellId; } @@ -261,41 +286,20 @@ struct FlattenicityTask { template 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; } @@ -304,7 +308,6 @@ struct FlattenicityTask { // ============================================ void fillMultiplicityClass(float flattenicity, float centrality, bool isParticle) { - // Particle histograms if (isParticle) { if (centrality < CentBound1) { histos.fill(HIST("hFlattenicityParticles_0_1"), flattenicity); @@ -334,7 +337,6 @@ struct FlattenicityTask { histos.fill(HIST("hFlattenicityParticles_95_100"), flattenicity); } } else { - // FT0 histograms if (centrality < CentBound1) { histos.fill(HIST("hFlattenicityFT0_0_1"), flattenicity); } else if (centrality < CentBound5) { @@ -368,66 +370,40 @@ struct FlattenicityTask { // ============================================ // Process MC collisions (generator level only) // ============================================ - void processMC( - aod::McCollisions const& /* mcCollisions */, - aod::McParticles const& mcParticles) + void processMC(aod::McCollisions const& /* mcCollisions */, aod::McParticles const& mcParticles) { if (mcParticles.size() == 0) { LOG(warning) << "No MC particles found in this data frame"; return; } - - // Initialize counters std::array truthCounts{}; - 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() & NPhysicalPrimaryBit) == 0) { + if ((particle.flags() & NPhysicalPrimaryBit) == 0) continue; - } - - // Check if charged - int charge = getCharge(particle.pdgCode()); - if (charge == 0) { + if (getCharge(particle.pdgCode()) == 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) { + if (std::abs(particle.eta()) < 1.0) nchINEL++; - } - - // dNch/deta: |eta| < 0.8 - if (std::abs(particle.eta()) < cfgEtaMax) { + if (std::abs(particle.eta()) < cfgEtaMax) nchFT0++; - } - // FT0 acceptance 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; int cellId = assignToFT0Cell(particle.eta(), particle.phi(), isFT0A); - if (cellId >= 0 && cellId < NCell) { truthCounts[cellId] += 1.0; histos.fill(HIST("hCellOccupancy"), cellId); @@ -439,22 +415,16 @@ struct FlattenicityTask { } } - // Event selection bool isINEL = (nchINEL > 0); bool isFT0 = (hasFT0A && hasFT0C); - - histos.fill(HIST("hEvents"), 0); // All events - + histos.fill(HIST("hEvents"), 0); if (isINEL) { - histos.fill(HIST("hEvents"), 1); // INEL>0 + histos.fill(HIST("hEvents"), 1); histos.fill(HIST("hNch_INEL"), nchFT0); } - if (isINEL && isFT0) { - histos.fill(HIST("hEvents"), 2); // INEL>0 & FT0 + histos.fill(HIST("hEvents"), 2); histos.fill(HIST("hNch_FT0"), nchFT0); - - // Compute flattenicity float rho = computeFlattenicity(truthCounts); if (rho > 0) { float flattenicity = 1.0 - rho; @@ -463,26 +433,17 @@ struct FlattenicityTask { } } } - - PROCESS_SWITCH(FlattenicityTask, processMC, "Process MC events", false); + PROCESS_SWITCH(FlattenicityTask, processMC, "Process MC events", true); // ============================================ // Process data collisions (reconstruction level only) // ============================================ - void processData( - CollisionsWithCent::iterator const& collision, - aod::FT0s const& ft0s, - FullTracks const& tracks) + void processData(CollisionsWithCent::iterator const& collision, aod::FT0s const& ft0s, FullTracks const& tracks) { - // Event selection: |vz| < 10 cm - if (std::abs(collision.posZ()) > cfgVzMax) { + if (std::abs(collision.posZ()) > cfgVzMax) return; - } - - // Get centrality percentile (0-100) float centrality = collision.centFT0M(); - // Find FT0 matching this collision's BC auto ft0 = ft0s.begin(); bool foundFT0 = false; for (const auto& f : ft0s) { @@ -492,29 +453,21 @@ struct FlattenicityTask { break; } } - if (!foundFT0) { + if (!foundFT0) return; - } - // Track selection and counting std::array recoCounts{}; - for (const auto& track : tracks) { - if (!isSelectedTrack(track)) { + 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; } } + histos.fill(HIST("hEvents"), 3); - histos.fill(HIST("hEvents"), 3); // Data events - - // Compute particle flattenicity float rhoParticles = computeFlattenicity(recoCounts); if (rhoParticles > 0) { float flattenicity = 1.0 - rhoParticles; @@ -522,41 +475,52 @@ struct FlattenicityTask { fillMultiplicityClass(flattenicity, centrality, true); } - // Compute flattenicity from FT0 amplitudes - std::array ft0Counts{}; - - // FT0-A channels (0-95) + // ============================================ + // FT0 detector-amplitude flattenicity. + // Matches multFilter.cxx (Antonio): FT0-A (24 sectors) and + // FT0-C (28 sectors) are each turned into their own flattenicity + // value, then averaged. No further phi subdivision is applied + // here, unlike the particle-level cell grid above. + // ============================================ + std::array ft0CountsA{}; + std::array ft0CountsC{}; if (ft0.amplitudeA().size() > 0) { for (std::size_t i = 0; i < ft0.amplitudeA().size(); i++) { uint8_t channel = ft0.channelA()[i]; - if (channel < NchA) { - ft0Counts[channel] = ft0.amplitudeA()[i]; + int sector = getT0ASector(channel); + if (sector >= 0 && sector < NSectorsA) { + ft0CountsA[sector] += ft0.amplitudeA()[i]; } } } - - // FT0-C channels (96-207) if (ft0.amplitudeC().size() > 0) { for (std::size_t i = 0; i < ft0.amplitudeC().size(); i++) { uint8_t channel = ft0.channelC()[i]; - if (channel < NchC) { - ft0Counts[NchA + channel] = ft0.amplitudeC()[i]; + int sector = getT0CSector(channel); + if (sector >= 0 && sector < NSectorsC) { + ft0CountsC[sector] += ft0.amplitudeC()[i]; } } } - - float rhoFT0 = computeFlattenicity(ft0Counts); - if (rhoFT0 > 0) { + float rhoFT0A = computeFlattenicity(ft0CountsA); + float rhoFT0C = computeFlattenicity(ft0CountsC); + if (rhoFT0A > 0) { + histos.fill(HIST("hFlattenicityFT0A"), 1.0 - rhoFT0A); + } + if (rhoFT0C > 0) { + histos.fill(HIST("hFlattenicityFT0C"), 1.0 - rhoFT0C); + } + if (rhoFT0A > 0 && rhoFT0C > 0) { + float rhoFT0 = (rhoFT0A + rhoFT0C) / 2.0; float flattenicity = 1.0 - rhoFT0; histos.fill(HIST("hFlattenicityFT0"), flattenicity); fillMultiplicityClass(flattenicity, centrality, false); } } - PROCESS_SWITCH(FlattenicityTask, processData, "Process data events", true); // ============================================ - // Process GEN vs RECO correlation + // Process GEN vs RECO correlation (3 additional plots) // ============================================ Preslice perMCCol = aod::mcparticle::mcCollisionId; SliceCache cache; @@ -568,28 +532,22 @@ struct FlattenicityTask { aod::FT0s const& ft0s, FullTracks const& tracks) { - if (std::abs(collision.posZ()) > cfgVzMax) { + if (std::abs(collision.posZ()) > cfgVzMax) return; - } - if (!collision.has_mcCollision()) { + if (!collision.has_mcCollision()) return; - } - // ---- Generator level: same truth-flattenicity logic as processMC ---- const auto& mcCollision = collision.mcCollision(); const auto& particlesInThisCollision = mcParticles.sliceBy(perMCCol, mcCollision.globalIndex()); std::array truthCounts{}; for (const auto& particle : particlesInThisCollision) { - if ((particle.flags() & NPhysicalPrimaryBit) == 0) { + if ((particle.flags() & NPhysicalPrimaryBit) == 0) continue; - } - if (getCharge(particle.pdgCode()) == 0) { + if (getCharge(particle.pdgCode()) == 0) continue; - } - if (particle.pt() < cfgPtMin) { + if (particle.pt() < cfgPtMin) continue; - } bool isFT0A = false; int cellId = assignToFT0Cell(particle.eta(), particle.phi(), isFT0A); if (cellId >= 0 && cellId < NCell) { @@ -597,12 +555,10 @@ struct FlattenicityTask { } } float rhoGen = computeFlattenicity(truthCounts); - if (rhoGen <= 0) { + if (rhoGen <= 0) return; - } float flatGen = 1.0 - rhoGen; - // ---- Reconstruction level: same reco logic as processData ---- auto ft0 = ft0s.begin(); bool foundFT0 = false; for (const auto& f : ft0s) { @@ -615,9 +571,8 @@ struct FlattenicityTask { std::array recoCounts{}; for (const auto& track : tracks) { - if (!isSelectedTrack(track)) { + if (!isSelectedTrack(track)) continue; - } bool isFT0A = false; int cellId = assignToFT0Cell(track.eta(), track.phi(), isFT0A); if (cellId >= 0 && cellId < NCell) { @@ -625,8 +580,6 @@ struct FlattenicityTask { } } float rhoRecParticles = computeFlattenicity(recoCounts); - - // Fill gen-level histogram restricted to events with a reco match histos.fill(HIST("hFlattenicityGen_MatchedToReco"), flatGen); if (rhoRecParticles > 0) { @@ -635,32 +588,37 @@ struct FlattenicityTask { } if (foundFT0) { - std::array ft0Counts{}; + // Same FT0-A/FT0-C split-then-average fix as processData above. + std::array ft0CountsA{}; + std::array ft0CountsC{}; if (ft0.amplitudeA().size() > 0) { for (std::size_t i = 0; i < ft0.amplitudeA().size(); i++) { uint8_t channel = ft0.channelA()[i]; - if (channel < NchA) { - ft0Counts[channel] = ft0.amplitudeA()[i]; + int sector = getT0ASector(channel); + if (sector >= 0 && sector < NSectorsA) { + ft0CountsA[sector] += ft0.amplitudeA()[i]; } } } if (ft0.amplitudeC().size() > 0) { for (std::size_t i = 0; i < ft0.amplitudeC().size(); i++) { uint8_t channel = ft0.channelC()[i]; - if (channel < NchC) { - ft0Counts[NchA + channel] = ft0.amplitudeC()[i]; + int sector = getT0CSector(channel); + if (sector >= 0 && sector < NSectorsC) { + ft0CountsC[sector] += ft0.amplitudeC()[i]; } } } - float rhoRecFT0 = computeFlattenicity(ft0Counts); - if (rhoRecFT0 > 0) { + float rhoRecFT0A = computeFlattenicity(ft0CountsA); + float rhoRecFT0C = computeFlattenicity(ft0CountsC); + if (rhoRecFT0A > 0 && rhoRecFT0C > 0) { + float rhoRecFT0 = (rhoRecFT0A + rhoRecFT0C) / 2.0; float flatRecFT0 = 1.0 - rhoRecFT0; histos.fill(HIST("hFlatFT0_Gen_vs_Rec"), flatGen, flatRecFT0); } } } - - PROCESS_SWITCH(FlattenicityTask, processGenRecCorrelation, "Process gen-vs-reco correlation (needs MC file)", false); + PROCESS_SWITCH(FlattenicityTask, processGenRecCorrelation, "Process gen-vs-reco correlation (needs MC file)", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)