From fb51a3cbc1151345e5db1a3787db631b5c8d4768 Mon Sep 17 00:00:00 2001 From: Neelkamal Mallick Date: Fri, 29 Aug 2025 13:50:46 +0300 Subject: [PATCH 1/2] [PWGCF] Event shape capabilities Added --- PWGCF/DataModel/CorrelationsDerived.h | 10 ++ PWGCF/TableProducer/filterCorrelations.cxx | 170 ++++++++++++++++++++- PWGCF/Tasks/correlations.cxx | 21 ++- 3 files changed, 199 insertions(+), 2 deletions(-) diff --git a/PWGCF/DataModel/CorrelationsDerived.h b/PWGCF/DataModel/CorrelationsDerived.h index 9be8ac7d58e..0199f2775ac 100644 --- a/PWGCF/DataModel/CorrelationsDerived.h +++ b/PWGCF/DataModel/CorrelationsDerived.h @@ -62,6 +62,16 @@ using CFCollLabel = CFCollLabels::iterator; using CFCollisionsWithLabel = soa::Join; using CFCollisionWithLabel = CFCollisionsWithLabel::iterator; +namespace cfeventshape +{ +DECLARE_SOA_INDEX_COLUMN(CFCollision, cfCollision); //! Index to collision +DECLARE_SOA_COLUMN(Spherocity, spherocity, float); //! Spherocity +} // namespace cfeventshape +DECLARE_SOA_TABLE(CFEventShapes, "AOD", "CFEVENTSHAPE", //! Event shape table + o2::soa::Index<>, + cfeventshape::Spherocity); +using CFEventShape = CFEventShapes::iterator; + namespace cftrack { DECLARE_SOA_INDEX_COLUMN(CFCollision, cfCollision); //! Index to collision diff --git a/PWGCF/TableProducer/filterCorrelations.cxx b/PWGCF/TableProducer/filterCorrelations.cxx index 70a009346c4..1d6b1e281fc 100644 --- a/PWGCF/TableProducer/filterCorrelations.cxx +++ b/PWGCF/TableProducer/filterCorrelations.cxx @@ -574,9 +574,177 @@ struct MultiplicitySelector { PROCESS_SWITCH(MultiplicitySelector, processMCGen, "Select MC particle count as multiplicity", false); }; +struct EventShapeProducer { + Produces OutputEventShape; + + O2_DEFINE_CONFIGURABLE(minTracks, int, 10, "Minimum number of tracks per collision required to estimate spherocity"); + O2_DEFINE_CONFIGURABLE(minPtTrack, float, 0.2, "Minimum pT of tracks required to estimate spherocity"); + O2_DEFINE_CONFIGURABLE(AbsEtaTrack, float, 0.8, "Maximum |eta| of tracks required to estimate spherocity"); + O2_DEFINE_CONFIGURABLE(usePtWeighted, bool, false, "Use pt-weighted spherocity"); + O2_DEFINE_CONFIGURABLE(outputQAHistos, bool, false, "Whether to output QA histograms for spherocity"); + + HistogramRegistry QAhistos{"QAhistos", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry QAhistosCent{"QAhistosCent", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry QAhistosSpherocity{"QAhistosSpherocity", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + + std::vector nx, ny, px, py; + + void init(InitContext const&) + { + AxisSpec axisPt = {100, 0.0, 50.0}; + AxisSpec axisEta = {100, -1.2, 1.2}; + AxisSpec axisPhi = {100, 0.0, 2.0 * M_PI}; + AxisSpec axisNtracks = {100, 0.0, 200.0}; + AxisSpec axisMult = {100, 0.0, 100.0}; + AxisSpec axisSpherocity = {1000, 0.0, 1.0}; + + QAhistos.add("multS0", "",{HistType::kTH2F,{axisMult, axisSpherocity}}); + QAhistos.add("spherocity", "", {HistType::kTH1F, {axisSpherocity}}); + + if(!outputQAHistos) return; + QAhistos.add("pt", "", {HistType::kTH1F, {axisPt}}); + QAhistos.add("eta", "", {HistType::kTH1F, {axisEta}}); + QAhistos.add("phi", "", {HistType::kTH1F, {axisPhi}}); + QAhistos.add("ntracks", "", {HistType::kTH1F, {axisNtracks}}); + QAhistos.add("mult", "", {HistType::kTH1F, {axisMult}}); + QAhistos.add("etaphi", "", {HistType::kTH2F, {axisPhi, axisEta}}); + + QAhistosCent.add("ntracks_0_20", "", {HistType::kTH1F, {axisNtracks}}); + QAhistosCent.add("ntracks_20_40", "", {HistType::kTH1F, {axisNtracks}}); + QAhistosCent.add("ntracks_40_60", "", {HistType::kTH1F, {axisNtracks}}); + QAhistosCent.add("ntracks_60_100", "", {HistType::kTH1F, {axisNtracks}}); + + QAhistosSpherocity.add("spherocity_0_20", "", {HistType::kTH1F, {axisSpherocity}}); + QAhistosSpherocity.add("spherocity_20_40", "", {HistType::kTH1F, {axisSpherocity}}); + QAhistosSpherocity.add("spherocity_40_60", "", {HistType::kTH1F, {axisSpherocity}}); + QAhistosSpherocity.add("spherocity_60_100", "", {HistType::kTH1F, {axisSpherocity}}); + } + + // Spherocity calculation based on the central barrel tracks + template + float CalculateSpherocity(TTracks const& tracks, bool usePtWeighted, bool outputQAHistos) + { + const size_t nTracks = tracks.size(); + if (nTracks < static_cast(minTracks)) return -2.0f; + + if (nx.size() < nTracks) { + nx.resize(nTracks); + ny.resize(nTracks); + px.resize(nTracks); + py.resize(nTracks); + } + + float sumPt = 0.0f; + float retval = 999.0f; + + size_t idx = 0; + for (auto& track : tracks) { + const float cosPhi = std::cos(track.phi()); + const float sinPhi = std::sin(track.phi()); + + nx[idx] = cosPhi; + ny[idx] = sinPhi; + + if (usePtWeighted) { + const float pt = track.pt(); + px[idx] = pt * cosPhi; + py[idx] = pt * sinPhi; + sumPt += pt; + } else { + px[idx] = cosPhi; + py[idx] = sinPhi; + sumPt += 1.0f; + } + if (outputQAHistos) { + QAhistos.fill(HIST("pt"), track.pt()); + QAhistos.fill(HIST("eta"), track.eta()); + QAhistos.fill(HIST("phi"), track.phi()); + QAhistos.fill(HIST("etaphi"), track.phi(), track.eta()); + } + idx++; + } + + if (sumPt == 0.0f) return -1.5f; // no tracks -- avoid division by zero + + // Validation check for non-weighted case + if (!usePtWeighted && std::abs(static_cast(nTracks) - sumPt) > 0.01f) { + LOGF(info, "Spherocity calculation: number of tracks (%zu) does not match sum of pT (%f)", nTracks, sumPt); + return -1.0f; + } + + for (size_t i = 0; i < nTracks; i++) { + float numerator = 0.0f; + const float nyi = ny[i]; + const float nxi = nx[i]; + + for (size_t j = 0; j < nTracks; j++) { + numerator += std::abs(nyi * px[j] - nxi * py[j]); + } + + const float sFull = std::pow(numerator / sumPt, 2); + if (sFull < retval) retval = sFull; + } + + retval = retval * M_PI * M_PI / 4.0f; // normalization factor + + if (retval < 0.0f || retval > 1.0f) { + LOGF(info, "Spherocity value is out of range: %f", retval); + return -0.5f; + } + + return retval; + }//spherocity calculation ends + + Filter cftrackFilter = (aod::cftrack::pt > minPtTrack) && (nabs(aod::cftrack::eta) < AbsEtaTrack); + + using myCollisions = aod::CFCollisions; + using myTracks = soa::Filtered; + + inline int getCentralityBin(float multiplicity) { + if (multiplicity < 20.0f) return 0; + else if (multiplicity < 40.0f) return 1; + else if (multiplicity < 60.0f) return 2; + else if (multiplicity < 100.0f) return 3; + return -1; // outside range + } + + void processSpherocityData(myCollisions::iterator const& coll, myTracks const& tracks) + { + float spherocity = CalculateSpherocity(tracks, usePtWeighted, outputQAHistos); + OutputEventShape(spherocity); + const float multiplicity = coll.multiplicity(); + QAhistos.fill(HIST("multS0"), multiplicity, spherocity); + QAhistos.fill(HIST("spherocity"), spherocity); + + if (!outputQAHistos) return; + + const size_t nTracks = tracks.size(); + const int centralityBin = getCentralityBin(multiplicity); + + QAhistos.fill(HIST("ntracks"), nTracks); + QAhistos.fill(HIST("mult"), multiplicity); + + if (centralityBin == 0) { + QAhistosCent.fill(HIST("ntracks_0_20"), nTracks); + QAhistosSpherocity.fill(HIST("spherocity_0_20"), spherocity); + } else if (centralityBin == 1) { + QAhistosCent.fill(HIST("ntracks_20_40"), nTracks); + QAhistosSpherocity.fill(HIST("spherocity_20_40"), spherocity); + } else if (centralityBin == 2) { + QAhistosCent.fill(HIST("ntracks_40_60"), nTracks); + QAhistosSpherocity.fill(HIST("spherocity_40_60"), spherocity); + } else if (centralityBin == 3) { + QAhistosCent.fill(HIST("ntracks_60_100"), nTracks); + QAhistosSpherocity.fill(HIST("spherocity_60_100"), spherocity); + } + } + PROCESS_SWITCH(EventShapeProducer, processSpherocityData, "Spherocity analysis for Data", false); +}; + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; } diff --git a/PWGCF/Tasks/correlations.cxx b/PWGCF/Tasks/correlations.cxx index 2a6101a5b9a..df41444c104 100644 --- a/PWGCF/Tasks/correlations.cxx +++ b/PWGCF/Tasks/correlations.cxx @@ -109,6 +109,9 @@ struct CorrelationTask { O2_DEFINE_CONFIGURABLE(cfgPtDepMLbkg, std::vector, {}, "pT interval for ML training") O2_DEFINE_CONFIGURABLE(cfgPtCentDepMLbkgSel, std::vector, {}, "Bkg ML selection") + O2_DEFINE_CONFIGURABLE(cfgS0CutLow, float, 0.0, "Lower cut for spherocity"); + O2_DEFINE_CONFIGURABLE(cfgS0CutUp, float, 1.0, "Upper cut for spherocity"); + ConfigurableAxis axisVertex{"axisVertex", {7, -7, 7}, "vertex axis for histograms"}; ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, PIHalf * 3}, "delta phi axis for histograms"}; ConfigurableAxis axisDeltaEta{"axisDeltaEta", {40, -2, 2}, "delta eta axis for histograms"}; @@ -126,7 +129,9 @@ struct CorrelationTask { Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex; // This filter is only applied to AOD Filter collisionVertexTypeFilter = (aod::collision::flags & static_cast(aod::collision::CollisionFlagsRun2::Run2VertexerTracks)) == static_cast(aod::collision::CollisionFlagsRun2::Run2VertexerTracks); - + + Filter collisionSpherocityFilter = (aod::cfeventshape::spherocity > cfgS0CutLow) && (aod::cfeventshape::spherocity < cfgS0CutUp); + // Track filters Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPt) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)); Filter cfTrackFilter = (nabs(aod::cftrack::eta) < cfgCutEta) && (aod::cftrack::pt > cfgCutPt) && ((aod::track::trackType & (uint8_t)cfgTrackBitMask) == (uint8_t)cfgTrackBitMask); @@ -167,6 +172,8 @@ struct CorrelationTask { using DerivedCollisions = soa::Filtered; using DerivedTracks = soa::Filtered; + using DerivedCollisionsWithSpherocity = soa::Filtered>; + void init(o2::framework::InitContext&) { if ((doprocessSame2ProngDerivedML || doprocessSame2Prong2ProngML || doprocessMixed2ProngDerivedML || doprocessMixed2Prong2ProngML) && (cfgPtDepMLbkg->empty() || cfgPtCentDepMLbkgSel->empty())) @@ -859,6 +866,12 @@ struct CorrelationTask { } PROCESS_SWITCH(CorrelationTask, processSameDerivedMultSet, "Process same event on derived data with multiplicity sets", false); + void processSameDerivedWithSpherocity(DerivedCollisionsWithSpherocity::iterator const& collision, soa::Filtered const& tracks) + { + processSameDerivedT(collision, tracks, tracks); + } + PROCESS_SWITCH(CorrelationTask, processSameDerivedWithSpherocity, "Process same event on derived data with spherocity", false); + void processSame2ProngDerived(DerivedCollisions::iterator const& collision, soa::Filtered const& tracks, soa::Filtered const& p2tracks) { processSameDerivedT(collision, p2tracks, tracks); @@ -995,6 +1008,12 @@ struct CorrelationTask { } PROCESS_SWITCH(CorrelationTask, processMixedDerivedMultSet, "Process mixed events on derived data with multiplicity sets", false); + void processMixedDerivedWithSpherocity(DerivedCollisionsWithSpherocity const& collisions, DerivedTracks const& tracks) + { + processMixedDerivedT(collisions, tracks); + } + PROCESS_SWITCH(CorrelationTask, processMixedDerivedWithSpherocity, "Process mixed events on derived data with spherocity", false); + void processMixed2ProngDerived(DerivedCollisions const& collisions, DerivedTracks const& tracks, soa::Filtered const& p2tracks) { processMixedDerivedT(collisions, p2tracks, tracks); From e350fc1f22625c032c634dc37048f7aab5d536d1 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 29 Aug 2025 11:06:01 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- PWGCF/DataModel/CorrelationsDerived.h | 2 +- PWGCF/TableProducer/filterCorrelations.cxx | 68 +++++++++++++--------- PWGCF/Tasks/correlations.cxx | 4 +- 3 files changed, 42 insertions(+), 32 deletions(-) diff --git a/PWGCF/DataModel/CorrelationsDerived.h b/PWGCF/DataModel/CorrelationsDerived.h index 0199f2775ac..1021c2a3c3d 100644 --- a/PWGCF/DataModel/CorrelationsDerived.h +++ b/PWGCF/DataModel/CorrelationsDerived.h @@ -65,7 +65,7 @@ using CFCollisionWithLabel = CFCollisionsWithLabel::iterator; namespace cfeventshape { DECLARE_SOA_INDEX_COLUMN(CFCollision, cfCollision); //! Index to collision -DECLARE_SOA_COLUMN(Spherocity, spherocity, float); //! Spherocity +DECLARE_SOA_COLUMN(Spherocity, spherocity, float); //! Spherocity } // namespace cfeventshape DECLARE_SOA_TABLE(CFEventShapes, "AOD", "CFEVENTSHAPE", //! Event shape table o2::soa::Index<>, diff --git a/PWGCF/TableProducer/filterCorrelations.cxx b/PWGCF/TableProducer/filterCorrelations.cxx index 1d6b1e281fc..1d36a3bfa42 100644 --- a/PWGCF/TableProducer/filterCorrelations.cxx +++ b/PWGCF/TableProducer/filterCorrelations.cxx @@ -582,13 +582,13 @@ struct EventShapeProducer { O2_DEFINE_CONFIGURABLE(AbsEtaTrack, float, 0.8, "Maximum |eta| of tracks required to estimate spherocity"); O2_DEFINE_CONFIGURABLE(usePtWeighted, bool, false, "Use pt-weighted spherocity"); O2_DEFINE_CONFIGURABLE(outputQAHistos, bool, false, "Whether to output QA histograms for spherocity"); - + HistogramRegistry QAhistos{"QAhistos", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; HistogramRegistry QAhistosCent{"QAhistosCent", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; HistogramRegistry QAhistosSpherocity{"QAhistosSpherocity", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; - + std::vector nx, ny, px, py; - + void init(InitContext const&) { AxisSpec axisPt = {100, 0.0, 50.0}; @@ -598,10 +598,11 @@ struct EventShapeProducer { AxisSpec axisMult = {100, 0.0, 100.0}; AxisSpec axisSpherocity = {1000, 0.0, 1.0}; - QAhistos.add("multS0", "",{HistType::kTH2F,{axisMult, axisSpherocity}}); + QAhistos.add("multS0", "", {HistType::kTH2F, {axisMult, axisSpherocity}}); QAhistos.add("spherocity", "", {HistType::kTH1F, {axisSpherocity}}); - if(!outputQAHistos) return; + if (!outputQAHistos) + return; QAhistos.add("pt", "", {HistType::kTH1F, {axisPt}}); QAhistos.add("eta", "", {HistType::kTH1F, {axisEta}}); QAhistos.add("phi", "", {HistType::kTH1F, {axisPhi}}); @@ -613,7 +614,7 @@ struct EventShapeProducer { QAhistosCent.add("ntracks_20_40", "", {HistType::kTH1F, {axisNtracks}}); QAhistosCent.add("ntracks_40_60", "", {HistType::kTH1F, {axisNtracks}}); QAhistosCent.add("ntracks_60_100", "", {HistType::kTH1F, {axisNtracks}}); - + QAhistosSpherocity.add("spherocity_0_20", "", {HistType::kTH1F, {axisSpherocity}}); QAhistosSpherocity.add("spherocity_20_40", "", {HistType::kTH1F, {axisSpherocity}}); QAhistosSpherocity.add("spherocity_40_60", "", {HistType::kTH1F, {axisSpherocity}}); @@ -621,11 +622,12 @@ struct EventShapeProducer { } // Spherocity calculation based on the central barrel tracks - template + template float CalculateSpherocity(TTracks const& tracks, bool usePtWeighted, bool outputQAHistos) - { + { const size_t nTracks = tracks.size(); - if (nTracks < static_cast(minTracks)) return -2.0f; + if (nTracks < static_cast(minTracks)) + return -2.0f; if (nx.size() < nTracks) { nx.resize(nTracks); @@ -641,10 +643,10 @@ struct EventShapeProducer { for (auto& track : tracks) { const float cosPhi = std::cos(track.phi()); const float sinPhi = std::sin(track.phi()); - + nx[idx] = cosPhi; ny[idx] = sinPhi; - + if (usePtWeighted) { const float pt = track.pt(); px[idx] = pt * cosPhi; @@ -664,8 +666,9 @@ struct EventShapeProducer { idx++; } - if (sumPt == 0.0f) return -1.5f; // no tracks -- avoid division by zero - + if (sumPt == 0.0f) + return -1.5f; // no tracks -- avoid division by zero + // Validation check for non-weighted case if (!usePtWeighted && std::abs(static_cast(nTracks) - sumPt) > 0.01f) { LOGF(info, "Spherocity calculation: number of tracks (%zu) does not match sum of pT (%f)", nTracks, sumPt); @@ -676,35 +679,41 @@ struct EventShapeProducer { float numerator = 0.0f; const float nyi = ny[i]; const float nxi = nx[i]; - + for (size_t j = 0; j < nTracks; j++) { numerator += std::abs(nyi * px[j] - nxi * py[j]); } - + const float sFull = std::pow(numerator / sumPt, 2); - if (sFull < retval) retval = sFull; + if (sFull < retval) + retval = sFull; } retval = retval * M_PI * M_PI / 4.0f; // normalization factor - + if (retval < 0.0f || retval > 1.0f) { LOGF(info, "Spherocity value is out of range: %f", retval); return -0.5f; } - + return retval; - }//spherocity calculation ends - + } // spherocity calculation ends + Filter cftrackFilter = (aod::cftrack::pt > minPtTrack) && (nabs(aod::cftrack::eta) < AbsEtaTrack); using myCollisions = aod::CFCollisions; using myTracks = soa::Filtered; - inline int getCentralityBin(float multiplicity) { - if (multiplicity < 20.0f) return 0; - else if (multiplicity < 40.0f) return 1; - else if (multiplicity < 60.0f) return 2; - else if (multiplicity < 100.0f) return 3; + inline int getCentralityBin(float multiplicity) + { + if (multiplicity < 20.0f) + return 0; + else if (multiplicity < 40.0f) + return 1; + else if (multiplicity < 60.0f) + return 2; + else if (multiplicity < 100.0f) + return 3; return -1; // outside range } @@ -715,12 +724,13 @@ struct EventShapeProducer { const float multiplicity = coll.multiplicity(); QAhistos.fill(HIST("multS0"), multiplicity, spherocity); QAhistos.fill(HIST("spherocity"), spherocity); - - if (!outputQAHistos) return; - + + if (!outputQAHistos) + return; + const size_t nTracks = tracks.size(); const int centralityBin = getCentralityBin(multiplicity); - + QAhistos.fill(HIST("ntracks"), nTracks); QAhistos.fill(HIST("mult"), multiplicity); diff --git a/PWGCF/Tasks/correlations.cxx b/PWGCF/Tasks/correlations.cxx index df41444c104..7ec1caec170 100644 --- a/PWGCF/Tasks/correlations.cxx +++ b/PWGCF/Tasks/correlations.cxx @@ -129,9 +129,9 @@ struct CorrelationTask { Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex; // This filter is only applied to AOD Filter collisionVertexTypeFilter = (aod::collision::flags & static_cast(aod::collision::CollisionFlagsRun2::Run2VertexerTracks)) == static_cast(aod::collision::CollisionFlagsRun2::Run2VertexerTracks); - + Filter collisionSpherocityFilter = (aod::cfeventshape::spherocity > cfgS0CutLow) && (aod::cfeventshape::spherocity < cfgS0CutUp); - + // Track filters Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPt) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)); Filter cfTrackFilter = (nabs(aod::cftrack::eta) < cfgCutEta) && (aod::cftrack::pt > cfgCutPt) && ((aod::track::trackType & (uint8_t)cfgTrackBitMask) == (uint8_t)cfgTrackBitMask);