diff --git a/PWGMM/UE/Tasks/flattenicityTask.cxx b/PWGMM/UE/Tasks/flattenicityTask.cxx index afd21b7a580..f279984cb9f 100644 --- a/PWGMM/UE/Tasks/flattenicityTask.cxx +++ b/PWGMM/UE/Tasks/flattenicityTask.cxx @@ -14,6 +14,7 @@ /// \author Eisha Rani /// \since August 2026 +#include "Common/DataModel/Centrality.h" #include "Common/DataModel/TrackSelectionTables.h" #include @@ -41,6 +42,8 @@ using namespace o2::soa; using namespace o2::constants::physics; using FullTracks = soa::Join; +using CollisionsWithCent = soa::Join; +using CollisionsWithCentAndMcLabel = soa::Join; struct FlattenicityTask { // ============================================ @@ -60,21 +63,22 @@ 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; + // Centrality percentile class boundaries (used in fillMultiplicityClass) + static constexpr float CentBound1 = 1.0f; + static constexpr float CentBound5 = 5.0f; + static constexpr float CentBound10 = 10.0f; + static constexpr float CentBound20 = 20.0f; + static constexpr float CentBound30 = 30.0f; + static constexpr float CentBound40 = 40.0f; + static constexpr float CentBound50 = 50.0f; + static constexpr float CentBound60 = 60.0f; + static constexpr float CentBound70 = 70.0f; + static constexpr float CentBound80 = 80.0f; + static constexpr float CentBound90 = 90.0f; + static constexpr float CentBound95 = 95.0f; // ============================================ - // Histogram Definitions + // Histogram Definitions - 100 BINS (bin width = 0.01) // ============================================ HistogramRegistry histos{ "histos", @@ -86,26 +90,73 @@ 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}}}}, + // ============================================ + // Flattenicity Histograms - 100 BINS! + // Bin width = 0.01, allowing 0-1% and 1-5% classes + // ============================================ + {"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}}}}, // FT0 cell occupancy {"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}}}}, + // ============================================ + // 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}}}}, + {"hFlattenicityFT0_0_1", "Flattenicity (FT0) class 0-1%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 1-5% + {"hFlattenicityParticles_1_5", "Flattenicity (particles) class 1-5%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0_1_5", "Flattenicity (FT0) class 1-5%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 5-10% + {"hFlattenicityParticles_5_10", "Flattenicity (particles) class 5-10%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0_5_10", "Flattenicity (FT0) class 5-10%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 10-20% + {"hFlattenicityParticles_10_20", "Flattenicity (particles) class 10-20%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0_10_20", "Flattenicity (FT0) class 10-20%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 20-30% + {"hFlattenicityParticles_20_30", "Flattenicity (particles) class 20-30%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0_20_30", "Flattenicity (FT0) class 20-30%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 30-40% + {"hFlattenicityParticles_30_40", "Flattenicity (particles) class 30-40%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0_30_40", "Flattenicity (FT0) class 30-40%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 40-50% + {"hFlattenicityParticles_40_50", "Flattenicity (particles) class 40-50%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0_40_50", "Flattenicity (FT0) class 40-50%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 50-60% + {"hFlattenicityParticles_50_60", "Flattenicity (particles) class 50-60%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0_50_60", "Flattenicity (FT0) class 50-60%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 60-70% + {"hFlattenicityParticles_60_70", "Flattenicity (particles) class 60-70%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0_60_70", "Flattenicity (FT0) class 60-70%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 70-80% + {"hFlattenicityParticles_70_80", "Flattenicity (particles) class 70-80%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0_70_80", "Flattenicity (FT0) class 70-80%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 80-90% + {"hFlattenicityParticles_80_90", "Flattenicity (particles) class 80-90%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0_80_90", "Flattenicity (FT0) class 80-90%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 90-95% + {"hFlattenicityParticles_90_95", "Flattenicity (particles) class 90-95%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"hFlattenicityFT0_90_95", "Flattenicity (FT0) class 90-95%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + // 95-100% + {"hFlattenicityParticles_95_100", "Flattenicity (particles) class 95-100%;1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, + {"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. + // ============================================ + {"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}}}}, + {"hFlattenicityGen_MatchedToReco", "Gen-level flattenicity (only events with a reco match);1-#rho;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}}, }}; // ============================================ @@ -249,7 +300,73 @@ struct FlattenicityTask { } // ============================================ - // Process MC collisions + // Helper function to fill multiplicity class histograms + // ============================================ + void fillMultiplicityClass(float flattenicity, float centrality, bool isParticle) + { + // Particle histograms + if (isParticle) { + if (centrality < CentBound1) { + histos.fill(HIST("hFlattenicityParticles_0_1"), flattenicity); + } else if (centrality < CentBound5) { + histos.fill(HIST("hFlattenicityParticles_1_5"), flattenicity); + } else if (centrality < CentBound10) { + histos.fill(HIST("hFlattenicityParticles_5_10"), flattenicity); + } else if (centrality < CentBound20) { + histos.fill(HIST("hFlattenicityParticles_10_20"), flattenicity); + } else if (centrality < CentBound30) { + histos.fill(HIST("hFlattenicityParticles_20_30"), flattenicity); + } else if (centrality < CentBound40) { + histos.fill(HIST("hFlattenicityParticles_30_40"), flattenicity); + } else if (centrality < CentBound50) { + histos.fill(HIST("hFlattenicityParticles_40_50"), flattenicity); + } else if (centrality < CentBound60) { + histos.fill(HIST("hFlattenicityParticles_50_60"), flattenicity); + } else if (centrality < CentBound70) { + histos.fill(HIST("hFlattenicityParticles_60_70"), flattenicity); + } else if (centrality < CentBound80) { + histos.fill(HIST("hFlattenicityParticles_70_80"), flattenicity); + } else if (centrality < CentBound90) { + histos.fill(HIST("hFlattenicityParticles_80_90"), flattenicity); + } else if (centrality < CentBound95) { + histos.fill(HIST("hFlattenicityParticles_90_95"), flattenicity); + } else { + histos.fill(HIST("hFlattenicityParticles_95_100"), flattenicity); + } + } else { + // FT0 histograms + if (centrality < CentBound1) { + histos.fill(HIST("hFlattenicityFT0_0_1"), flattenicity); + } else if (centrality < CentBound5) { + histos.fill(HIST("hFlattenicityFT0_1_5"), flattenicity); + } else if (centrality < CentBound10) { + histos.fill(HIST("hFlattenicityFT0_5_10"), flattenicity); + } else if (centrality < CentBound20) { + histos.fill(HIST("hFlattenicityFT0_10_20"), flattenicity); + } else if (centrality < CentBound30) { + histos.fill(HIST("hFlattenicityFT0_20_30"), flattenicity); + } else if (centrality < CentBound40) { + histos.fill(HIST("hFlattenicityFT0_30_40"), flattenicity); + } else if (centrality < CentBound50) { + histos.fill(HIST("hFlattenicityFT0_40_50"), flattenicity); + } else if (centrality < CentBound60) { + histos.fill(HIST("hFlattenicityFT0_50_60"), flattenicity); + } else if (centrality < CentBound70) { + histos.fill(HIST("hFlattenicityFT0_60_70"), flattenicity); + } else if (centrality < CentBound80) { + histos.fill(HIST("hFlattenicityFT0_70_80"), flattenicity); + } else if (centrality < CentBound90) { + histos.fill(HIST("hFlattenicityFT0_80_90"), flattenicity); + } else if (centrality < CentBound95) { + histos.fill(HIST("hFlattenicityFT0_90_95"), flattenicity); + } else { + histos.fill(HIST("hFlattenicityFT0_95_100"), flattenicity); + } + } + } + + // ============================================ + // Process MC collisions (generator level only) // ============================================ void processMC( aod::McCollisions const& /* mcCollisions */, @@ -341,31 +458,8 @@ struct FlattenicityTask { 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); - } + histos.fill(HIST("hFlattenicityParticles"), flattenicity); + histos.fill(HIST("hFlattenicityParticles_vs_Nch"), nchFT0, flattenicity); } } } @@ -373,10 +467,10 @@ struct FlattenicityTask { PROCESS_SWITCH(FlattenicityTask, processMC, "Process MC events", false); // ============================================ - // Process data collisions + // Process data collisions (reconstruction level only) // ============================================ void processData( - aod::Collision const& collision, + CollisionsWithCent::iterator const& collision, aod::FT0s const& ft0s, FullTracks const& tracks) { @@ -385,6 +479,9 @@ struct FlattenicityTask { return; } + // Get centrality percentile (0-100) + float centrality = collision.centFT0M(); + // Find FT0 matching this collision's BC auto ft0 = ft0s.begin(); bool foundFT0 = false; @@ -417,28 +514,153 @@ struct FlattenicityTask { histos.fill(HIST("hEvents"), 3); // Data events - // FT0-A channels (0-95) + // Compute particle flattenicity + float rhoParticles = computeFlattenicity(recoCounts); + if (rhoParticles > 0) { + float flattenicity = 1.0 - rhoParticles; + histos.fill(HIST("hFlattenicityParticles"), flattenicity); + fillMultiplicityClass(flattenicity, centrality, true); + } + + // Compute flattenicity from FT0 amplitudes std::array ft0Counts{}; - 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]; + + // FT0-A channels (0-95) + 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]; + } } } + // FT0-C channels (96-207) - 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]; + 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]; + } } } - float rho = computeFlattenicity(ft0Counts); - if (rho > 0) { - histos.fill(HIST("hFlattenicity"), 1.0 - rho); + + float rhoFT0 = computeFlattenicity(ft0Counts); + if (rhoFT0 > 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 + // ============================================ + Preslice perMCCol = aod::mcparticle::mcCollisionId; + SliceCache cache; + + void processGenRecCorrelation( + CollisionsWithCentAndMcLabel::iterator const& collision, + aod::McCollisions const&, + aod::McParticles const& mcParticles, + aod::FT0s const& ft0s, + FullTracks const& tracks) + { + if (std::abs(collision.posZ()) > cfgVzMax) { + return; + } + 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) { + continue; + } + if (getCharge(particle.pdgCode()) == 0) { + continue; + } + if (particle.pt() < cfgPtMin) { + continue; + } + bool isFT0A = false; + int cellId = assignToFT0Cell(particle.eta(), particle.phi(), isFT0A); + if (cellId >= 0 && cellId < NCell) { + truthCounts[cellId] += 1.0; + } + } + float rhoGen = computeFlattenicity(truthCounts); + 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) { + if (f.bcId() == collision.bcId()) { + ft0 = f; + foundFT0 = true; + break; + } + } + + std::array recoCounts{}; + for (const auto& track : tracks) { + if (!isSelectedTrack(track)) { + continue; + } + bool isFT0A = false; + int cellId = assignToFT0Cell(track.eta(), track.phi(), isFT0A); + if (cellId >= 0 && cellId < NCell) { + recoCounts[cellId] += 1.0; + } + } + float rhoRecParticles = computeFlattenicity(recoCounts); + + // Fill gen-level histogram restricted to events with a reco match + histos.fill(HIST("hFlattenicityGen_MatchedToReco"), flatGen); + + if (rhoRecParticles > 0) { + float flatRecParticles = 1.0 - rhoRecParticles; + histos.fill(HIST("hFlatParticles_Gen_vs_Rec"), flatGen, flatRecParticles); + } + + if (foundFT0) { + std::array ft0Counts{}; + 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]; + } + } + } + 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]; + } + } + } + float rhoRecFT0 = computeFlattenicity(ft0Counts); + if (rhoRecFT0 > 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); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)