diff --git a/PWGCF/FemtoUniverse/Core/FemtoUniverseCollisionSelection.h b/PWGCF/FemtoUniverse/Core/FemtoUniverseCollisionSelection.h index 40edf1bd87d..81aa05a86b0 100644 --- a/PWGCF/FemtoUniverse/Core/FemtoUniverseCollisionSelection.h +++ b/PWGCF/FemtoUniverse/Core/FemtoUniverseCollisionSelection.h @@ -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 @@ -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 @@ -162,7 +160,48 @@ class FemtoUniverseCollisionSelection template 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: diff --git a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx index 2d5a07f46e7..55e62e3a566 100644 --- a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx +++ b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx @@ -110,6 +110,8 @@ struct femtoUniverseProducerTask { Configurable 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 ConfDoSpher{"ConfDoSpher", false, "Calculate sphericity. If false sphericity will take value of 2."}; + /// Event cuts FemtoUniverseCollisionSelection colCuts; Configurable ConfEvtUseTPCmult{"ConfEvtUseTPCmult", false, "Use multiplicity based on the number of tracks with TPC information"}; @@ -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) { @@ -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 @@ -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) { @@ -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 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) { @@ -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 @@ -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()); @@ -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 const& tracks, - o2::aod::V0Datas const& fullV0s) + void processTrackV0(aod::FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + soa::Filtered const& tracks, + o2::aod::V0Datas const& fullV0s) { // get magnetic field for run getMagneticFieldTesla(col.bc_as()); @@ -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 const& tracks, - aod::McCollisions const& mcCollisions, - aod::McParticles const& mcParticles, - soa::Join const& fullV0s) + void processFullMC(aod::FemtoFullCollisionMC const& col, + aod::BCsWithTimestamps const&, + soa::Join const& tracks, + aod::McCollisions const& mcCollisions, + aod::McParticles const& mcParticles, + soa::Join const& fullV0s) { // get magnetic field for run getMagneticFieldTesla(col.bc_as()); @@ -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 const& tracks, - aod::McCollisions const& mcCollisions, - aod::McParticles const&) + void processTrackMC(aod::FemtoFullCollisionMC const& col, + aod::BCsWithTimestamps const&, + soa::Join const& tracks, + aod::McCollisions const& mcCollisions, + aod::McParticles const&) { // get magnetic field for run getMagneticFieldTesla(col.bc_as()); @@ -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()); @@ -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 const& tracks, - soa::Join const& candidates) + void processTrackD0mesonData(aod::FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + soa::Filtered const& tracks, + soa::Join const& candidates) { // get magnetic field for run getMagneticFieldTesla(col.bc_as()); @@ -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> const& collisions, - aod::McParticles const& mcParticles, - aod::BCsWithTimestamps const&) + void processTrackMCTruth(aod::McCollision const& mcCol, + soa::SmallGroups> const& collisions, + aod::McParticles const& mcParticles, + aod::BCsWithTimestamps const&) { // magnetic field for run not needed for mc truth // fill the tables @@ -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()); @@ -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()); diff --git a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackMultKtExtended.cxx b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackMultKtExtended.cxx index 508ce7dbdf2..777266a111b 100644 --- a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackMultKtExtended.cxx +++ b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackMultKtExtended.cxx @@ -135,8 +135,11 @@ struct femtoUniversePairTaskTrackTrackMultKtExtended { /// Event part Configurable ConfV0MLow{"ConfV0MLow", 0.0, "Lower limit for V0M multiplicity"}; Configurable ConfV0MHigh{"ConfV0MHigh", 25000.0, "Upper limit for V0M multiplicity"}; + Configurable ConfSphericityCutMin{"ConfSphericityCutMin", 0, "Min. sphericity"}; + Configurable 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 @@ -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 { @@ -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); @@ -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);