Skip to content

PWGCF: Implement sphericity calculations and add sphericity cut #5080

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

Merged
zchochul merged 1 commit into from
Mar 11, 2024
Merged

PWGCF: Implement sphericity calculations and add sphericity cut #5080

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
47 changes: 43 additions & 4 deletions PWGCF/FemtoUniverse/Core/FemtoUniverseCollisionSelection.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,7 @@ class FemtoUniverseCollisionSelection
mHistogramRegistry->add("Event/MultNTracksPV", "; vMultNTracksPV; Entries", kTH1F, {{120, 0, 120}});
mHistogramRegistry->add("Event/MultNTracklets", "; vMultNTrackslets; Entries", kTH1F, {{300, 0, 300}});
mHistogramRegistry->add("Event/MultTPC", "; vMultTPC; Entries", kTH1I, {{600, 0, 600}});
mHistogramRegistry->add("Event/Sphericity", "; Sphericity; Entries", kTH1I, {{200, 0, 3}});
}

/// Print some debug information
Expand Down Expand Up @@ -150,10 +151,7 @@ class FemtoUniverseCollisionSelection
}
}

/// \todo to be implemented!
/// Compute the sphericity of an event
/// Important here is that the filter on tracks does not interfere here!
/// In Run 2 we used here global tracks within |eta| < 0.8
/// \tparam T1 type of the collision
/// \tparam T2 type of the tracks
/// \param col Collision
Expand All @@ -162,7 +160,48 @@ class FemtoUniverseCollisionSelection
template <typename T1, typename T2>
float computeSphericity(T1 const& col, T2 const& tracks)
{
return 2.f;
double S00 = 0;
double S11 = 0;
double S10 = 0;
double sumPt = 0;
int partNumber = 0;
double spher = 0;

for (auto& p : tracks) {
double phi = p.phi();
double pT = p.pt();
double px = pT * TMath::Cos(phi);
double py = pT * TMath::Sin(phi);

S00 = S00 + px * px / pT;
S11 = S11 + py * py / pT;
S10 = S10 + px * py / pT;
sumPt = sumPt + pT;
partNumber++;
}

if (sumPt != 0) {
S00 = S00 / sumPt;
S11 = S11 / sumPt;
S10 = S10 / sumPt;

double lambda1 = (S00 + S11 + TMath::Sqrt((S00 + S11) * (S00 + S11) - 4.0 * (S00 * S11 - S10 * S10))) / 2.0;
double lambda2 = (S00 + S11 - TMath::Sqrt((S00 + S11) * (S00 + S11) - 4.0 * (S00 * S11 - S10 * S10))) / 2.0;

if ((lambda1 + lambda2) != 0 && partNumber > 2) {
spher = 2 * lambda2 / (lambda1 + lambda2);
} else {
spher = 2;
}
} else {
spher = 2;
}

if (mHistogramRegistry) {
mHistogramRegistry->fill(HIST("Event/Sphericity"), spher);
}

return spher;
}

private:
Expand Down
126 changes: 71 additions & 55 deletions PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,8 @@ struct femtoUniverseProducerTask {

Configurable<bool> ConfIsForceGRP{"ConfIsForceGRP", false, "Set true if the magnetic field configuration is not available in the usual CCDB directory (e.g. for Run 2 converted data or unanchorad Monte Carlo)"};

Configurable<bool> ConfDoSpher{"ConfDoSpher", false, "Calculate sphericity. If false sphericity will take value of 2."};

/// Event cuts
FemtoUniverseCollisionSelection colCuts;
Configurable<bool> ConfEvtUseTPCmult{"ConfEvtUseTPCmult", false, "Use multiplicity based on the number of tracks with TPC information"};
Expand Down Expand Up @@ -504,7 +506,6 @@ struct femtoUniverseProducerTask {
void fillCollisions(CollisionType const& col, TrackType const& tracks)
{
const auto vtxZ = col.posZ();
const auto spher = colCuts.computeSphericity(col, tracks);
int mult = 0;
int multNtr = 0;
if (ConfIsRun3) {
Expand All @@ -526,13 +527,21 @@ struct femtoUniverseProducerTask {
// particle candidates for such collisions
if (!colCuts.isSelected(col)) {
if (ConfIsTrigger) {
outputCollision(vtxZ, mult, multNtr, spher, mMagField);
if (ConfDoSpher) {
outputCollision(vtxZ, mult, multNtr, colCuts.computeSphericity(col, tracks), mMagField);
} else {
outputCollision(vtxZ, mult, multNtr, 2, mMagField);
}
}
return;
}

colCuts.fillQA(col);
outputCollision(vtxZ, mult, multNtr, spher, mMagField);
if (ConfDoSpher) {
outputCollision(vtxZ, mult, multNtr, colCuts.computeSphericity(col, tracks), mMagField);
} else {
outputCollision(vtxZ, mult, multNtr, 2, mMagField);
}
}

template <typename CollisionType, typename TrackType>
Expand All @@ -553,7 +562,6 @@ struct femtoUniverseProducerTask {
void fillCollisionsCentRun2(CollisionType const& col, TrackType const& tracks)
{
const auto vtxZ = col.posZ();
const auto spher = colCuts.computeSphericity(col, tracks);
int cent = 0;
int multNtr = 0;
if (!ConfIsRun3) {
Expand All @@ -568,20 +576,28 @@ struct femtoUniverseProducerTask {
// particle candidates for such collisions
if (!colCuts.isSelected(col)) {
if (ConfIsTrigger) {
outputCollision(vtxZ, cent, multNtr, spher, mMagField); //////
if (ConfDoSpher) {
outputCollision(vtxZ, cent, multNtr, colCuts.computeSphericity(col, tracks), mMagField);
} else {
outputCollision(vtxZ, cent, multNtr, 2, mMagField);
}
}

return;
}

// colCuts.fillQA(col); //for now, TODO: create a configurable so in the FemroUniverseCollisionSelection.h there is an option to plot QA just for the posZ
outputCollision(vtxZ, cent, multNtr, spher, mMagField);
if (ConfDoSpher) {
outputCollision(vtxZ, cent, multNtr, colCuts.computeSphericity(col, tracks), mMagField);
} else {
outputCollision(vtxZ, cent, multNtr, 2, mMagField);
}
}

template <bool isMC, typename CollisionType, typename TrackType>
void fillCollisionsCentRun3(CollisionType const& col, TrackType const& tracks)
{
const auto vtxZ = col.posZ();
const auto spher = colCuts.computeSphericity(col, tracks);
int cent = 0;
int multNtr = 0;
if (ConfIsRun3) {
Expand All @@ -596,13 +612,22 @@ struct femtoUniverseProducerTask {
// particle candidates for such collisions
if (!colCuts.isSelectedRun3(col)) {
if (ConfIsTrigger) {
outputCollision(vtxZ, cent, multNtr, spher, mMagField); //////
}
if (ConfDoSpher) {
outputCollision(vtxZ, cent, multNtr, colCuts.computeSphericity(col, tracks), mMagField);
} else {
outputCollision(vtxZ, cent, multNtr, 2, mMagField);
}
} //////

return;
}

// colCuts.fillQA(col); //for now, TODO: create a configurable so in the FemroUniverseCollisionSelection.h there is an option to plot QA just for the posZ
outputCollision(vtxZ, cent, multNtr, spher, mMagField);
if (ConfDoSpher) {
outputCollision(vtxZ, cent, multNtr, colCuts.computeSphericity(col, tracks), mMagField);
} else {
outputCollision(vtxZ, cent, multNtr, 2, mMagField);
}
}

template <bool isMC, typename TrackType>
Expand Down Expand Up @@ -1039,11 +1064,10 @@ struct femtoUniverseProducerTask {
}
}

void
processFullData(aod::FemtoFullCollision const& col,
aod::BCsWithTimestamps const&,
aod::FemtoFullTracks const& tracks,
o2::aod::V0Datas const& fullV0s)
void processFullData(aod::FemtoFullCollision const& col,
aod::BCsWithTimestamps const&,
aod::FemtoFullTracks const& tracks,
o2::aod::V0Datas const& fullV0s)
{
// get magnetic field for run
getMagneticFieldTesla(col.bc_as<aod::BCsWithTimestamps>());
Expand All @@ -1052,11 +1076,10 @@ struct femtoUniverseProducerTask {
}
PROCESS_SWITCH(femtoUniverseProducerTask, processFullData, "Provide experimental data", false);

void
processTrackV0(aod::FemtoFullCollision const& col,
aod::BCsWithTimestamps const&,
soa::Filtered<aod::FemtoFullTracks> const& tracks,
o2::aod::V0Datas const& fullV0s)
void processTrackV0(aod::FemtoFullCollision const& col,
aod::BCsWithTimestamps const&,
soa::Filtered<aod::FemtoFullTracks> const& tracks,
o2::aod::V0Datas const& fullV0s)
{
// get magnetic field for run
getMagneticFieldTesla(col.bc_as<aod::BCsWithTimestamps>());
Expand All @@ -1065,13 +1088,12 @@ struct femtoUniverseProducerTask {
}
PROCESS_SWITCH(femtoUniverseProducerTask, processTrackV0, "Provide experimental data for track v0", false);

void
processFullMC(aod::FemtoFullCollisionMC const& col,
aod::BCsWithTimestamps const&,
soa::Join<aod::FemtoFullTracks, aod::McTrackLabels> const& tracks,
aod::McCollisions const& mcCollisions,
aod::McParticles const& mcParticles,
soa::Join<o2::aod::V0Datas, aod::McV0Labels> const& fullV0s)
void processFullMC(aod::FemtoFullCollisionMC const& col,
aod::BCsWithTimestamps const&,
soa::Join<aod::FemtoFullTracks, aod::McTrackLabels> const& tracks,
aod::McCollisions const& mcCollisions,
aod::McParticles const& mcParticles,
soa::Join<o2::aod::V0Datas, aod::McV0Labels> const& fullV0s)
{
// get magnetic field for run
getMagneticFieldTesla(col.bc_as<aod::BCsWithTimestamps>());
Expand All @@ -1080,12 +1102,11 @@ struct femtoUniverseProducerTask {
}
PROCESS_SWITCH(femtoUniverseProducerTask, processFullMC, "Provide MC data (tracks, V0, Phi)", false);

void
processTrackMC(aod::FemtoFullCollisionMC const& col,
aod::BCsWithTimestamps const&,
soa::Join<aod::FemtoFullTracks, aod::McTrackLabels> const& tracks,
aod::McCollisions const& mcCollisions,
aod::McParticles const&)
void processTrackMC(aod::FemtoFullCollisionMC const& col,
aod::BCsWithTimestamps const&,
soa::Join<aod::FemtoFullTracks, aod::McTrackLabels> const& tracks,
aod::McCollisions const& mcCollisions,
aod::McParticles const&)
{
// get magnetic field for run
getMagneticFieldTesla(col.bc_as<aod::BCsWithTimestamps>());
Expand All @@ -1095,10 +1116,9 @@ struct femtoUniverseProducerTask {
}
PROCESS_SWITCH(femtoUniverseProducerTask, processTrackMC, "Provide MC data for track analysis", false);

void
processTrackData(aod::FemtoFullCollision const& col,
aod::BCsWithTimestamps const&,
aod::FemtoFullTracks const& tracks)
void processTrackData(aod::FemtoFullCollision const& col,
aod::BCsWithTimestamps const&,
aod::FemtoFullTracks const& tracks)
{
// get magnetic field for run
getMagneticFieldTesla(col.bc_as<aod::BCsWithTimestamps>());
Expand All @@ -1124,11 +1144,10 @@ struct femtoUniverseProducerTask {
PROCESS_SWITCH(femtoUniverseProducerTask, processTrackPhiData,
"Provide experimental data for track phi", false);

void
processTrackD0mesonData(aod::FemtoFullCollision const& col,
aod::BCsWithTimestamps const&,
soa::Filtered<aod::FemtoFullTracks> const& tracks,
soa::Join<aod::HfCand2Prong, aod::HfSelD0> const& candidates)
void processTrackD0mesonData(aod::FemtoFullCollision const& col,
aod::BCsWithTimestamps const&,
soa::Filtered<aod::FemtoFullTracks> const& tracks,
soa::Join<aod::HfCand2Prong, aod::HfSelD0> const& candidates)
{
// get magnetic field for run
getMagneticFieldTesla(col.bc_as<aod::BCsWithTimestamps>());
Expand All @@ -1140,11 +1159,10 @@ struct femtoUniverseProducerTask {
PROCESS_SWITCH(femtoUniverseProducerTask, processTrackD0mesonData,
"Provide experimental data for track D0 meson", false);

void
processTrackMCTruth(aod::McCollision const& mcCol,
soa::SmallGroups<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::McCollisionLabels>> const& collisions,
aod::McParticles const& mcParticles,
aod::BCsWithTimestamps const&)
void processTrackMCTruth(aod::McCollision const& mcCol,
soa::SmallGroups<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::McCollisionLabels>> const& collisions,
aod::McParticles const& mcParticles,
aod::BCsWithTimestamps const&)
{
// magnetic field for run not needed for mc truth
// fill the tables
Expand All @@ -1153,10 +1171,9 @@ struct femtoUniverseProducerTask {
}
PROCESS_SWITCH(femtoUniverseProducerTask, processTrackMCTruth, "Provide MC data for MC truth track analysis", false);

void
processTrackCentRun2Data(aod::FemtoFullCollisionCentRun2 const& col,
aod::BCsWithTimestamps const&,
aod::FemtoFullTracks const& tracks)
void processTrackCentRun2Data(aod::FemtoFullCollisionCentRun2 const& col,
aod::BCsWithTimestamps const&,
aod::FemtoFullTracks const& tracks)
{
// get magnetic field for run
getMagneticFieldTesla(col.bc_as<aod::BCsWithTimestamps>());
Expand All @@ -1166,10 +1183,9 @@ struct femtoUniverseProducerTask {
}
PROCESS_SWITCH(femtoUniverseProducerTask, processTrackCentRun2Data, "Provide experimental data for Run 2 with centrality for track track", false);

void
processTrackCentRun3Data(aod::FemtoFullCollisionCentRun3 const& col,
aod::BCsWithTimestamps const&,
aod::FemtoFullTracks const& tracks)
void processTrackCentRun3Data(aod::FemtoFullCollisionCentRun3 const& col,
aod::BCsWithTimestamps const&,
aod::FemtoFullTracks const& tracks)
{
// get magnetic field for run
getMagneticFieldTesla(col.bc_as<aod::BCsWithTimestamps>());
Expand Down
8 changes: 8 additions & 0 deletions PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackMultKtExtended.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,8 +135,11 @@ struct femtoUniversePairTaskTrackTrackMultKtExtended {
/// Event part
Configurable<float> ConfV0MLow{"ConfV0MLow", 0.0, "Lower limit for V0M multiplicity"};
Configurable<float> ConfV0MHigh{"ConfV0MHigh", 25000.0, "Upper limit for V0M multiplicity"};
Configurable<float> ConfSphericityCutMin{"ConfSphericityCutMin", 0, "Min. sphericity"};
Configurable<float> ConfSphericityCutMax{"ConfSphericityCutMax", 3, "Max. sphericity"};

Filter collV0Mfilter = ((o2::aod::femtouniversecollision::multV0M > ConfV0MLow) && (o2::aod::femtouniversecollision::multV0M < ConfV0MHigh));
Filter colSpherfilter = ((o2::aod::femtouniversecollision::sphericity > ConfSphericityCutMin) && (o2::aod::femtouniversecollision::sphericity < ConfSphericityCutMax));
// Filter trackAdditionalfilter = (nabs(aod::femtouniverseparticle::eta) < twotracksconfigs.ConfEtaMax); // example filtering on configurable

/// Particle part
Expand Down Expand Up @@ -214,6 +217,8 @@ struct femtoUniversePairTaskTrackTrackMultKtExtended {
HistogramRegistry SameMultRegistryMM{"SameMultRegistryMM", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
HistogramRegistry MixedMultRegistryMM{"MixedMultRegistryMM", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};

HistogramRegistry sphericityRegistry{"SphericityHisto", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};

// PID for protons
bool IsProtonNSigma(float mom, float nsigmaTPCPr, float nsigmaTOFPr) // previous version from: https://github.com/alisw/AliPhysics/blob/master/PWGCF/FEMTOSCOPY/AliFemtoUser/AliFemtoMJTrackCut.cxx
{
Expand Down Expand Up @@ -343,6 +348,8 @@ struct femtoUniversePairTaskTrackTrackMultKtExtended {
void init(InitContext&)
{
eventHisto.init(&qaRegistry);
sphericityRegistry.add("sphericity", ";Sphericity;Entries", kTH1F, {{150, 0.0, 3, "Sphericity"}});

trackHistoPartOne.init(&qaRegistry, ConfTempFitVarpTBins, ConfTempFitVarBins, twotracksconfigs.ConfIsMC, trackonefilter.ConfPDGCodePartOne, true);

trackHistoPartTwo.init(&qaRegistry, ConfTempFitVarpTBins, ConfTempFitVarBins, twotracksconfigs.ConfIsMC, tracktwofilter.ConfPDGCodePartTwo, true);
Expand Down Expand Up @@ -528,6 +535,7 @@ struct femtoUniversePairTaskTrackTrackMultKtExtended {
FilteredFemtoFullParticles& parts)
{
fillCollision(col);
sphericityRegistry.fill(HIST("sphericity"), col.sphericity());

auto thegroupPartsOne = partsOne->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
auto thegroupPartsTwo = partsTwo->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
Expand Down