Skip to content

[PWGDQ] Implementation of Transverse Spherocity in PWGDQ #18299

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
akavumpa wants to merge 3 commits into
base: master
Choose a base branch
from
Draft

[PWGDQ] Implementation of Transverse Spherocity in PWGDQ #18299

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f19a6ee
Committing older uncom changes
Aug 31, 2026
e2c8d48
Changes after PR type 12 sept 2026
Sep 12, 2026
05232cb
Implementation of Spherocity in PWGDQ
Oct 14, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file added PWGDQ.zip
View file
Open in desktop
Binary file not shown.
29 changes: 29 additions & 0 deletions PWGDQ/Core/HistogramsLibrary.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,6 +121,26 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
if (subGroupStr.Contains("randomplane")) {
hm->AddHistogram(histClass, "random plane angle", "random Psi", false, 100, -TMath::Pi() / 2, TMath::Pi() / 2, VarManager::kRandomPsi2);
}
if (subGroupStr.Contains("spherocity")) {
hm->AddHistogram(histClass, "Spherocity", "Transverse Spherocity", false, 100, 0.0, 1.0, VarManager::kSpherocity);
hm->AddHistogram(histClass, "SpherocityPtWeighted", "pT-weighted Spherocity", false, 100, 0.0, 1.0, VarManager::kSpherocityPtWeighted); // Spherocity vs event properties
hm->AddHistogram(histClass, "Spherocity_VtxZ", "Spherocity vs Vtx Z", false, 100, 0.0, 1.0, VarManager::kSpherocity, 60, -15.0, 15.0, VarManager::kVtxZ);
hm->AddHistogram(histClass, "MultITSTPC_VtxZ", "MultITSTPC vs Vtx Z", false, 200, 0.0, 200.0, VarManager::kMultNTracksITSTPC, 60, -15.0, 15.0, VarManager::kVtxZ);
hm->AddHistogram(histClass, "Spherocity_VtxNContrib", "Spherocity vs Vtx n contributors", false, 100, 0.0, 1.0, VarManager::kSpherocity, 100, 0.0, 100.0, VarManager::kVtxNcontrib);
hm->AddHistogram(histClass, "MultITSTPC_VtxNContrib", "MultITSTPC vs Vtx n contributors", false, 200, 0.0, 200.0, VarManager::kMultNTracksITSTPC, 100, 0.0, 100.0, VarManager::kVtxNcontrib);
hm->AddHistogram(histClass, "Spherocity_MultITS", "Spherocity vs ITS only Multiplicity", false, 100, 0.0, 1.0, VarManager::kSpherocity, 250, 0.0, 250.0, VarManager::kMultNTracksITSOnly);
hm->AddHistogram(histClass, "Spherocity_MultTPC", "Spherocity vs TPC only Multiplicity", false, 100, 0.0, 1.0, VarManager::kSpherocity, 250, 0.0, 5000.0, VarManager::kMultTPC);
hm->AddHistogram(histClass, "Spherocity_MultITSTPC", "Spherocity vs ITS-TPC Multiplicity", false, 100, 0.0, 1.0, VarManager::kSpherocity, 250, 0.0, 250.0, VarManager::kMultNTracksITSTPC);
hm->AddHistogram(histClass, "Spherocity_CentFT0C", "Spherocity vs Centrality FT0C", false, 100, 0.0, 1.0, VarManager::kSpherocity, 20, 0.0, 100.0, VarManager::kCentFT0C);
hm->AddHistogram(histClass, "SpherocityPtWeighted_CentFT0C", "pT-weighted Spherocity vs Centrality FT0C", false, 100, 0.0, 1.0, VarManager::kSpherocityPtWeighted, 20, 0.0, 100.0, VarManager::kCentFT0C);
hm->AddHistogram(histClass, "Spherocity_CentFT0C_3D", "Spherocity vs Centrality vs Mult", false, 50, 0.0, 1.0, VarManager::kSpherocity, 20, 0.0, 100.0, VarManager::kCentFT0C, 100, 0.0, 50000.0, VarManager::kMultTPC);

if (subGroupStr.Contains("mySpheroCent")) {
hm->AddHistogram(histClass, "Spherocity_CentFT0C", "Spherocity vs Centrality FT0C", false, 100, 0.0, 1.0, VarManager::kSpherocity, 20, 0.0, 100.0, VarManager::kCentFT0C);
hm->AddHistogram(histClass, "SpherocityPtWeighted_CentFT0C", "pT-weighted Spherocity vs Centrality FT0C", false, 100, 0.0, 1.0, VarManager::kSpherocityPtWeighted, 20, 0.0, 100.0, VarManager::kCentFT0C);
hm->AddHistogram(histClass, "Spherocity_CentFT0C_3D", "Spherocity vs Centrality vs Mult", false, 50, 0.0, 1.0, VarManager::kSpherocity, 20, 0.0, 100.0, VarManager::kCentFT0C, 100, 0.0, 50000.0, VarManager::kMultTPC);
}
}
if (subGroupStr.Contains("mult")) {
if (subGroupStr.Contains("pp")) {
hm->AddHistogram(histClass, "MultTPC", "MultTPC", false, 250, 0.0, 500.0, VarManager::kMultTPC);
Expand Down Expand Up @@ -1632,6 +1652,15 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
hm->AddHistogram(histClass, "Mass_Pt", "", false, 750, 0.0, 15.0, VarManager::kMass, 120, 0.0, 30.0, VarManager::kPt);
hm->AddHistogram(histClass, "Mass_Rapidity", "", false, 750, 0.0, 15.0, VarManager::kMass, 150, 2.5, 4.0, VarManager::kRap);
hm->AddHistogram(histClass, "Mass_Phi", "", false, 750, 0.0, 15.0, VarManager::kMass, 180, constants::math::PI, 2 * constants::math::PI, VarManager::kPhi);
hm->AddHistogram(histClass, "Mass", "", false, 750, 0.0, 5.0, VarManager::kMass);
if (subGroupStr.Contains("dimuon-spherocity")) {
hm->AddHistogram(histClass, "Spherocity", "Transverse Spherocity", false, 100, 0.0, 1.0, VarManager::kSpherocity);
hm->AddHistogram(histClass, "SpherocityPtWeighted", "pT-weighted Spherocity", false, 100, 0.0, 1.0, VarManager::kSpherocityPtWeighted);
hm->AddHistogram(histClass, "Mass_Spherocity", "", false, 750, 0.0, 5.0, VarManager::kMass, 100, 0.0, 1.0, VarManager::kSpherocity);
hm->AddHistogram(histClass, "Mass_SpherocityPtWeighted", "", false, 750, 0.0, 5.0, VarManager::kMass, 100, 0.0, 1.0, VarManager::kSpherocityPtWeighted);
hm->AddHistogram(histClass, "Mass_Pt_Spherocity", "", false, 750, 0.0, 5.0, VarManager::kMass, 120, 0.0, 30.0, VarManager::kPt, 100, 0.0, 1.0, VarManager::kSpherocity);
hm->AddHistogram(histClass, "Mass_Pt_SpherocityPtWeighted", "", false, 200, 0.0, 5.0, VarManager::kMass, 30, 0.0, 30.0, VarManager::kPt, 50, 0.0, 1.0, VarManager::kSpherocityPtWeighted);
}
if (subGroupStr.Contains("dimuon-tag-and-probe")) {
hm->AddHistogram(histClass, "Mass_PtProbe", "mass vs probe pT", false, 750, 0.0, 15.0, VarManager::kMass, 120, 0.0, 30.0, VarManager::kPt2); // Warning: in the tag-and-probe task t2 is always the probe
hm->AddHistogram(histClass, "Mass_EtaProbe", "mass vs probe eta", false, 750, 0.0, 15.0, VarManager::kMass, 120, 0.0, 30.0, VarManager::kEta2); // Warning: in the tag-and-probe task t2 is always the probe
Expand Down
6 changes: 6 additions & 0 deletions PWGDQ/Core/VarManager.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -816,6 +816,10 @@ void VarManager::SetDefaultVarNames()
fgVariableUnits[kMultDimuonsME] = "";
fgVariableNames[kCentFT0C] = "Centrality FT0C";
fgVariableUnits[kCentFT0C] = "%";
fgVariableNames[kSpherocity] = "Spherocity";
fgVariableUnits[kSpherocity] = "";
fgVariableNames[kSpherocityPtWeighted] = "SpherocityPtWeighted";
fgVariableUnits[kSpherocityPtWeighted] = "";
fgVariableNames[kMCEventGeneratorId] = "MC Generator ID";
fgVariableNames[kMCEventSubGeneratorId] = "MC SubGenerator ID";
fgVariableNames[kMCVtxX] = "MC Vtx X";
Expand Down Expand Up @@ -2096,6 +2100,8 @@ void VarManager::SetDefaultVarNames()
fgVarNamesMap["kVtxChi2"] = kVtxChi2;
fgVarNamesMap["kCentVZERO"] = kCentVZERO;
fgVarNamesMap["kCentFT0C"] = kCentFT0C;
fgVarNamesMap["kSpherocity"] = kSpherocity;
fgVarNamesMap["kSpherocityPtWeighted"] = kSpherocityPtWeighted;
fgVarNamesMap["kCentFT0A"] = kCentFT0A;
fgVarNamesMap["kCentFT0M"] = kCentFT0M;
fgVarNamesMap["kMultTPC"] = kMultTPC;
Expand Down
82 changes: 82 additions & 0 deletions PWGDQ/Core/VarManager.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,7 @@ class VarManager : public TObject
RapidityGapFilter = BIT(21),
Fit = BIT(22),
ReducedFit = BIT(23),
ReducedEventSpherocity = BIT(24),
Track = BIT(0),
TrackCov = BIT(1),
TrackExtra = BIT(2),
Expand Down Expand Up @@ -283,6 +284,8 @@ class VarManager : public TObject
kMultNTracksITSTPC,
kMultNTracksPVeta1,
kMultNTracksPVetaHalf,
kSpherocity, // transverse spherocity
kSpherocityPtWeighted, // pT-weighted transverse spherocity
kTrackOccupancyInTimeRange,
kFT0COccupancyInTimeRange,
kNoCollInTimeRangeStandard,
Expand Down Expand Up @@ -1438,6 +1441,10 @@ class VarManager : public TObject
{
return RecoDecay::constrainAngle(psi1 - psi2, -o2::constants::math::PI / harmonic, harmonic);
}

template <typename T>
static float CalculateSpherocity(T const& tracks, float ptMin = 0.15, float ptMax = 10.0, float etaMin = -0.8, float etaMax = 0.8, int minMult = 10, bool usePtWeight = false);

template <typename T, typename T1>
static o2::dataformats::VertexBase RecalculatePrimaryVertex(T const& track0, T const& track1, const T1& collision);

Expand Down Expand Up @@ -2316,6 +2323,11 @@ void VarManager::FillEvent(T const& event, float* values)
}
}

if constexpr ((fillMap & ObjTypes::ReducedEventSpherocity) > 0) {
values[kSpherocity] = event.spherocity();
values[kSpherocityPtWeighted] = event.spherocityPtWeighted();
}

if constexpr ((fillMap & ReducedEventVtxCov) > 0) {
values[kVtxCovXX] = event.covXX();
values[kVtxCovXY] = event.covXY();
Expand Down Expand Up @@ -7996,5 +8008,75 @@ void VarManager::FillTripletVertexingALICE3(C const& collision, T const& t1, T c
values[kVertexingTauxyzProjected] = values[kVertexingLxyzProjected] * v123.M() / (v123.P());
}
}
template <typename T>
float VarManager::CalculateSpherocity(T const& tracks, float ptMin, float ptMax, float etaMin, float etaMax, int minMult, bool usePtWeight)
{
std::vector<float> pxNorm;
std::vector<float> pyNorm;
std::vector<float> ptValues;
float sumPt = 0.0;
int nTracks = 0;

// Collect valid tracks
for (const auto& track : tracks) {
float pt = track.pt();
float eta = track.eta();

if (pt < ptMin || pt > ptMax)
continue;
if (eta < etaMin || eta > etaMax)
continue;

float px = track.px();
float py = track.py();

if (pt <= 0.f) {
continue;
}
pxNorm.push_back(px / pt);
pyNorm.push_back(py / pt);
ptValues.push_back(pt);
sumPt += pt;
nTracks++;
}

// Check minimum multiplicity
if (nTracks < minMult) {
return -1.0f;
}

// Minimize over azimuthal angle, hence nSteps = 360 for 1 degree steps
float spherocity = 5.0f;
const int nSteps = 360;

for (int iStep = 0; iStep < nSteps; iStep++) {
float phi = 2.0f * M_PI * static_cast<float>(iStep) / static_cast<float>(nSteps);
float nx = std::cos(phi);
float ny = std::sin(phi);

float numerator = 0.0f;
float denominator = 0.0f;

for (int iTrk = 0; iTrk < nTracks; iTrk++) {
float crossProduct = std::abs(pxNorm[iTrk] * ny - pyNorm[iTrk] * nx); // Cross product magnitude: |p_T x n| = |px*ny - py*nx|

if (usePtWeight) {
numerator += ptValues[iTrk] * crossProduct;
denominator += ptValues[iTrk];
} else {
numerator += crossProduct;
denominator += 1.0f;
}
}

float sph = std::pow(numerator / denominator, 2);
if (sph < spherocity) {
spherocity = sph;
}
}

spherocity *= (M_PI * M_PI / 4.0f);
return spherocity;
}

#endif // PWGDQ_CORE_VARMANAGER_H_
6 changes: 6 additions & 0 deletions PWGDQ/DataModel/ReducedInfoTables.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,8 @@ DECLARE_SOA_COLUMN(NTPCoccupMeanTimeShortA, nTPCoccupMeanTimeShortA, float);
DECLARE_SOA_COLUMN(NTPCoccupMeanTimeShortC, nTPCoccupMeanTimeShortC, float); //! TPC pileup mean time on C side (short time range)
DECLARE_SOA_COLUMN(NTPCoccupMedianTimeShortA, nTPCoccupMedianTimeShortA, float); //! TPC pileup median time on A side (short time range)
DECLARE_SOA_COLUMN(NTPCoccupMedianTimeShortC, nTPCoccupMedianTimeShortC, float); //! TPC pileup median time on C side (short time range)
DECLARE_SOA_COLUMN(Spherocity, spherocity, float); //! transverse spherocity
DECLARE_SOA_COLUMN(SpherocityPtWeighted, spherocityPtWeighted, float); //! pT-weighted transverse spherocity
DECLARE_SOA_COLUMN(DCAzBimodalityCoefficient, dcazBimodalityCoefficient, float); //! Bimodality coefficient of the DCAz distribution of the tracks in the event
DECLARE_SOA_COLUMN(DCAzBimodalityCoefficientBinned, dcazBimodalityCoefficientBinned, float); //! Bimodality coefficient of the DCAz distribution of the tracks in the event, binned
DECLARE_SOA_COLUMN(DCAzBimodalityCoefficientBinnedTrimmed1, dcazBimodalityCoefficientBinnedTrimmed1, float); //! Bimodality coefficient of the DCAz distribution of the tracks in the event, binned and trimmed 1
Expand Down Expand Up @@ -179,6 +181,10 @@ DECLARE_SOA_TABLE_VERSIONED(ReducedEventsExtended_001, "AOD", "REEXTENDED", 1, /

using ReducedEventsExtended = ReducedEventsExtended_001;

DECLARE_SOA_TABLE(ReducedEventsSpherocity, "AOD", "RESPHEROCITY",
reducedevent::Spherocity,
reducedevent::SpherocityPtWeighted);

DECLARE_SOA_TABLE(ReducedEventsMultPV_000, "AOD", "REMULTPV", //! Multiplicity information for primary vertex
mult::MultNTracksHasITS, mult::MultNTracksHasTPC, mult::MultNTracksHasTOF, mult::MultNTracksHasTRD,
mult::MultNTracksITSOnly, mult::MultNTracksTPCOnly, mult::MultNTracksITSTPC,
Expand Down
72 changes: 68 additions & 4 deletions PWGDQ/TableProducer/tableMaker_withAssoc.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -222,6 +222,7 @@ struct TableMaker {
// Q-vector related tables, to be filled only if the user selects the corresponding option
Produces<ReducedEventsQvectorCentr> eventQvectorCentr;
Produces<ReducedEventsQvectorCentrExtra> eventQvectorCentrExtra;
Produces<ReducedEventsSpherocity> eventSpherocity;
} outTables;

OutputObj<THashList> fOutputList{"output"}; //! the histogram manager output list
Expand Down Expand Up @@ -333,6 +334,17 @@ struct TableMaker {
Configurable<bool> fExcludeShort{"cfgTPCExcludeShort", true, "Exclude short term from long term occupancy (micro-seconds)"};
} fConfigVariousOptions;

// Spherocity configurables
struct : ConfigurableGroup {
Configurable<bool> fConfigSpherocityCompute{"cfgSpherocityCompute", false, "Compute and store event spherocity (off by default)"};
Configurable<float> fConfigSpherocityPtMin{"cfgSpherocityPtMin", 0.15f, "Min pT for spherocity calculation"};
Configurable<float> fConfigSpherocityPtMax{"cfgSpherocityPtMax", 10.0f, "Max pT for spherocity calculation"};
Configurable<float> fConfigSpherocityEtaMin{"cfgSpherocityEtaMin", -0.8f, "Min eta for spherocity calculation"};
Configurable<float> fConfigSpherocityEtaMax{"cfgSpherocityEtaMax", 0.8f, "Max eta for spherocity calculation"};
Configurable<int> fConfigSpherocityMinMult{"cfgSpherocityMinMult", 5, "Min multiplicity for spherocity calculation"};
Configurable<bool> fConfigSpherocityUsePtWeight{"cfgSpherocityUsePtWeight", false, "Use pT weighting for spherocity"};
} fConfigSpherocity;

Service<o2::ccdb::BasicCCDBManager> fCCDB{};
o2::ccdb::CcdbApi fCCDBApi;

Expand Down Expand Up @@ -1033,10 +1045,11 @@ struct TableMaker {
}
}

template <uint32_t TEventFillMap, typename TEvents, typename TBCs,
typename TZdcs, typename TFt0s, typename TFv0as, typename TFdds>
template <uint32_t TEventFillMap, uint32_t TTrackFillMap, typename TEvents, typename TBCs,
typename TZdcs, typename TFt0s, typename TFv0as, typename TFdds, typename TTrackAssoc, typename TTracks>
void skimCollisions(TEvents const& collisions, TBCs const& bcs, TZdcs const& /*zdcs*/,
TFt0s const& ft0s, TFv0as const& fv0as, TFdds const& fdds)
TFt0s const& ft0s, TFv0as const& fv0as, TFdds const& fdds,
TTrackAssoc const& trackAssocs, TTracks const& tracksBarrel)
{
// Skim collisions
// NOTE: So far, collisions are filtered based on the user specified analysis cuts AND the filterPP or Zorro event filter.
Expand Down Expand Up @@ -1222,6 +1235,55 @@ struct TableMaker {
}
(dynamic_cast<TH2D*>(fStatsList->At(kStatsEvent)))->Fill(3.0, static_cast<float>(o2::aod::evsel::kNsel));

// Calculate spherocity for this collision
float spherocityValue = -1.0f;
float spherocityPtWeighted = -1.0f;

if (fConfigSpherocity.fConfigSpherocityCompute) {
if constexpr (static_cast<bool>(TTrackFillMap)) {
auto groupedTrackIndices = trackAssocs.sliceBy(trackIndicesPerCollision, collision.globalIndex());

if (static_cast<int>(groupedTrackIndices.size()) >= fConfigSpherocity.fConfigSpherocityMinMult) {
std::vector<typename std::decay_t<decltype(tracksBarrel.begin())>> trackVec;
trackVec.reserve(groupedTrackIndices.size());

int validTracks = 0;
for (auto const& assoc : groupedTrackIndices) {
if (assoc.trackId() < tracksBarrel.size()) {
auto track = tracksBarrel.rawIteratorAt(assoc.trackId());
trackVec.push_back(track);
validTracks++;
}
}

if (validTracks >= fConfigSpherocity.fConfigSpherocityMinMult) {
// Calculate unweighted spherocity
spherocityValue = VarManager::CalculateSpherocity(
trackVec,
fConfigSpherocity.fConfigSpherocityPtMin,
fConfigSpherocity.fConfigSpherocityPtMax,
fConfigSpherocity.fConfigSpherocityEtaMin,
fConfigSpherocity.fConfigSpherocityEtaMax,
fConfigSpherocity.fConfigSpherocityMinMult,
false);

// Calculate pT-weighted spherocity
spherocityPtWeighted = VarManager::CalculateSpherocity(
trackVec,
fConfigSpherocity.fConfigSpherocityPtMin,
fConfigSpherocity.fConfigSpherocityPtMax,
fConfigSpherocity.fConfigSpherocityEtaMin,
fConfigSpherocity.fConfigSpherocityEtaMax,
fConfigSpherocity.fConfigSpherocityMinMult,
true);
}
}
}
}

VarManager::fgValues[VarManager::kSpherocity] = spherocityValue;
VarManager::fgValues[VarManager::kSpherocityPtWeighted] = spherocityPtWeighted;

fHistMan->FillHistClass("Event_AfterCuts", dqtablemaker_helpers::varValues());

// create the event tables
Expand All @@ -1248,6 +1310,7 @@ struct TableMaker {
multTPC, multFV0A, multFV0C, multFT0A, multFT0C, multFDDA, multFDDC, multZNA, multZNC, multTracklets, multTracksPV, centFT0C, centFT0A, centFT0M);
outTables.eventVtxCov(collision.covXX(), collision.covXY(), collision.covXZ(), collision.covYY(), collision.covYZ(), collision.covZZ(), collision.chi2());
outTables.eventInfo(collision.globalIndex());
outTables.eventSpherocity(spherocityValue, spherocityPtWeighted);

if constexpr ((TEventFillMap & VarManager::ObjTypes::CollisionQvectCentr) > 0) {
if (fConfigQvector.fConfigQvectCalibAvailable) {
Expand Down Expand Up @@ -1914,8 +1977,9 @@ struct TableMaker {
outTables.event.reserve(collisions.size());
outTables.eventExtended.reserve(collisions.size());
outTables.eventVtxCov.reserve(collisions.size());
outTables.eventSpherocity.reserve(collisions.size());

skimCollisions<TEventFillMap>(collisions, bcs, zdcs, ft0s, fv0as, fdds);
skimCollisions<TEventFillMap, TTrackFillMap>(collisions, bcs, zdcs, ft0s, fv0as, fdds, trackAssocs, tracksBarrel);
if (fCollIndexMap.empty()) {
return;
}
Expand Down
Loading
Loading