diff --git a/PWGLF/DataModel/LFPhiStrangeCorrelationTables.h b/PWGLF/DataModel/LFPhiStrangeCorrelationTables.h new file mode 100644 index 00000000000..b19b67f9388 --- /dev/null +++ b/PWGLF/DataModel/LFPhiStrangeCorrelationTables.h @@ -0,0 +1,57 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file LFPhiStrangeCorrelationTables.h +/// \brief Data model Phi-Strangeness correlation analysis +/// \author Stefano Cannito (stefano.cannito@cern.ch) + +#ifndef PWGLF_DATAMODEL_LFPHISTRANGECORRELATIONTABLES_H_ +#define PWGLF_DATAMODEL_LFPHISTRANGECORRELATIONTABLES_H_ + +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include + +namespace o2::aod +{ +namespace lf_selection_phi_collision +{ +DECLARE_SOA_COLUMN(PhimesonSel, phimesonSel, bool); +} // namespace lf_selection_phi_collision + +DECLARE_SOA_TABLE(PhimesonSelection, "AOD", "PHIINCOLL", + lf_selection_phi_collision::PhimesonSel); + +namespace lf_selection_phi_candidate +{ +DECLARE_SOA_INDEX_COLUMN(Collision, collision); + +DECLARE_SOA_COLUMN(M, m, float); +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(Y, y, float); +DECLARE_SOA_COLUMN(Phi, phi, float); + +DECLARE_SOA_DYNAMIC_COLUMN(InMassRegion, inMassRegion, + [](float m, float minM, float maxM) -> bool { + return (m >= minM && m <= maxM); + }); +} // namespace lf_selection_phi_candidate + +DECLARE_SOA_TABLE(PhimesonCandidates, "AOD", "PHICANDIDATES", + lf_selection_phi_candidate::CollisionId, + lf_selection_phi_candidate::M, + lf_selection_phi_candidate::Pt, + lf_selection_phi_candidate::Y, + lf_selection_phi_candidate::Phi, + lf_selection_phi_candidate::InMassRegion); +} // namespace o2::aod + +#endif // PWGLF_DATAMODEL_LFPHISTRANGECORRELATIONTABLES_H_ diff --git a/PWGLF/TableProducer/Strangeness/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/CMakeLists.txt index 36914058504..c739a55f6ed 100644 --- a/PWGLF/TableProducer/Strangeness/CMakeLists.txt +++ b/PWGLF/TableProducer/Strangeness/CMakeLists.txt @@ -163,3 +163,8 @@ if(FastJet_FOUND) PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGJECore FastJet::FastJet FastJet::Contrib O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) endif() + +o2physics_add_dpl_workflow(phi-strange-correlator + SOURCES phiStrangeCorrelator.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/PWGLF/TableProducer/Strangeness/phiStrangeCorrelator.cxx b/PWGLF/TableProducer/Strangeness/phiStrangeCorrelator.cxx new file mode 100644 index 00000000000..273b20fe7cf --- /dev/null +++ b/PWGLF/TableProducer/Strangeness/phiStrangeCorrelator.cxx @@ -0,0 +1,462 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file phiStrangeCorrelator.cxx +/// \brief Table producer for Phi-Strangeness correlation analysis +/// \author Stefano Cannito (stefano.cannito@cern.ch) + +#include "PWGLF/DataModel/LFPhiStrangeCorrelationTables.h" +#include "PWGLF/DataModel/mcCentrality.h" +#include "PWGLF/Utils/inelGt.h" + +#include "Common/Core/TableHelper.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/TrackSelectionDefaults.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CCDB/BasicCCDBManager.h" +#include "CommonConstants/PhysicsConstants.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct PhiMesonCandProducer { + // Produce the table with the event selection information + Produces phimesonCandidates; + + HistogramRegistry histos{"phiCandidates", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + + // Configurables for phi's daughter tracks selection + struct : ConfigurableGroup { + Configurable cfgCutCharge{"cfgCutCharge", 0.0f, "Cut on charge"}; + Configurable cfgGlobalWoDCATrack{"cfgGlobalWoDCATrack", true, "Global track selection without DCA"}; + Configurable cfgPVContributor{"cfgPVContributor", true, "PV contributor track selection"}; + Configurable cMinKaonPtcut{"cMinKaonPtcut", 0.15f, "Track minimum pt cut"}; + Configurable etaMax{"etaMax", 0.8f, "eta max"}; + Configurable pTToUseTOF{"pTToUseTOF", 0.5f, "pT above which use TOF"}; + Configurable cMaxDCAzToPVcut{"cMaxDCAzToPVcut", 2.0f, "Track DCAz cut to PV Maximum"}; + Configurable> cMaxDCArToPVPhi{"cMaxDCArToPVPhi", {0.004f, 0.013f, 1.0f}, "Track DCAr cut to PV for Phi"}; + + Configurable cfgIsDCAzParameterized{"cfgIsDCAzParameterized", false, "IsDCAzParameterized"}; + + Configurable nSigmaCutTPCKa{"nSigmaCutTPCKa", 3.0f, "Value of the TPC Nsigma cut for Kaons"}; + Configurable nSigmaCutCombinedKa{"nSigmaCutCombinedKa", 3.0f, "Value of the TPC and TOF Nsigma cut for Kaons"}; + + Configurable minTPCnClsFound{"minTPCnClsFound", 70, "min number of found TPC clusters"}; + } trackConfigs; + + // Configurables on phi selection + struct : ConfigurableGroup { + Configurable maxMPhi{"maxMPhi", 1.5f, "Maximum mass for Phi candidates"}; + Configurable minPhiPt{"minPhiPt", 0.4f, "Minimum pT for Phi candidates"}; + Configurable cfgYAcceptance{"cfgYAcceptance", 0.5f, "Rapidity acceptance"}; + } phiConfigs; + + // Defining the type of the collisions for data and MC + using SelCollisions = soa::Join; + using SimCollisions = soa::Join; + + // Defining the type of the phi's daughter tracks for data and MC + using FullTracks = soa::Join; + using FullMCTracks = soa::Join; + + Partition posTracks = aod::track::signed1Pt > trackConfigs.cfgCutCharge; + Partition negTracks = aod::track::signed1Pt < trackConfigs.cfgCutCharge; + + Partition posMCTracks = aod::track::signed1Pt > trackConfigs.cfgCutCharge; + Partition negMCTracks = aod::track::signed1Pt < trackConfigs.cfgCutCharge; + + // Manual slicing + Preslice trackPerCollision = aod::track::collisionId; + SliceCache cache; + + // Constants + double massKa = o2::constants::physics::MassKPlus; + + void init(InitContext&) + { + AxisSpec massPhiAxis = {200, 0.9f, 1.2f, "#it{M}_{inv} [GeV/#it{c}^{2}]"}; + + histos.add("h1PhiCandidateMass", "Phi candidate invariant mass", kTH1F, {massPhiAxis}); + } + + // Topological track selection + template + bool selectionTrackResonance(const T& track) + { + if (trackConfigs.cfgGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) + return false; + if (trackConfigs.cfgPVContributor && !track.isPVContributor()) + return false; + + if (track.tpcNClsFound() < trackConfigs.minTPCnClsFound) + return false; + + if (track.pt() < trackConfigs.cMinKaonPtcut) + return false; + + if (std::abs(track.dcaXY()) > trackConfigs.cMaxDCArToPVPhi->at(0) + (trackConfigs.cMaxDCArToPVPhi->at(1) / std::pow(track.pt(), trackConfigs.cMaxDCArToPVPhi->at(2)))) + return false; + if (std::abs(track.dcaZ()) > trackConfigs.cMaxDCAzToPVcut) + return false; + + return true; + } + + // PIDKaon track selection + template + bool selectionPIDKaonpTdependent(const T& track) + { + if (track.pt() < trackConfigs.pTToUseTOF && std::abs(track.tpcNSigmaKa()) < trackConfigs.nSigmaCutTPCKa) + return true; + if (track.pt() >= trackConfigs.pTToUseTOF && track.hasTOF() && (std::pow(track.tofNSigmaKa(), 2) + std::pow(track.tpcNSigmaKa(), 2)) < std::pow(trackConfigs.nSigmaCutCombinedKa, 2)) + return true; + + return false; + } + + // Reconstruct the Phi candidate + template + ROOT::Math::PxPyPzMVector recMother(const T& track1, const T& track2, float massdaughter1, float massdaughter2) + { + ROOT::Math::PxPyPzMVector daughter1(track1.px(), track1.py(), track1.pz(), massdaughter1); // set the daughter1 4-momentum + ROOT::Math::PxPyPzMVector daughter2(track2.px(), track2.py(), track2.pz(), massdaughter2); // set the daughter2 4-momentum + ROOT::Math::PxPyPzMVector mother = daughter1 + daughter2; // calculate the mother 4-momentum + + return mother; + } + + void processData(SelCollisions::iterator const& collision, FullTracks const&) + { + auto posThisColl = posTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + auto negThisColl = negTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + + for (const auto& track1 : posThisColl) { + if (!selectionTrackResonance(track1) || !selectionPIDKaonpTdependent(track1)) + continue; + + for (const auto& track2 : negThisColl) { + if (!selectionTrackResonance(track2) || !selectionPIDKaonpTdependent(track2)) + continue; + + ROOT::Math::PxPyPzMVector recPhi = recMother(track1, track2, massKa, massKa); + + if (recPhi.Pt() < phiConfigs.minPhiPt) + continue; + if (recPhi.M() > phiConfigs.maxMPhi) + continue; + if (std::abs(recPhi.Rapidity()) > phiConfigs.cfgYAcceptance) + continue; + + histos.fill(HIST("h1PhiCandidateMass"), recPhi.M()); + + phimesonCandidates(collision.globalIndex(), recPhi.M(), recPhi.Pt(), recPhi.Rapidity(), recPhi.Phi()); + } + } + } + + PROCESS_SWITCH(PhiMesonCandProducer, processData, "Process function to select Phi meson candidates in Data", true); + + void processMCReco(SimCollisions::iterator const& collision, FullMCTracks const&) + { + auto posThisColl = posMCTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + auto negThisColl = negMCTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + + for (const auto& track1 : posThisColl) { + if (!selectionTrackResonance(track1) || !selectionPIDKaonpTdependent(track1)) + continue; + + for (const auto& track2 : negThisColl) { + if (!selectionTrackResonance(track2) || !selectionPIDKaonpTdependent(track2)) + continue; + + ROOT::Math::PxPyPzMVector recPhi = recMother(track1, track2, massKa, massKa); + + if (recPhi.Pt() < phiConfigs.minPhiPt) + continue; + if (recPhi.M() > phiConfigs.maxMPhi) + continue; + if (std::abs(recPhi.Rapidity()) > phiConfigs.cfgYAcceptance) + continue; + + phimesonCandidates(collision.globalIndex(), recPhi.M(), recPhi.Pt(), recPhi.Rapidity(), recPhi.Phi()); + } + } + } + + PROCESS_SWITCH(PhiMesonCandProducer, processMCReco, "Process function to select Phi meson candidates in MCReco", false); + + void processMCGen(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles) + { + for (const auto& mcParticle1 : mcParticles) { + if (!mcParticle1.isPhysicalPrimary() || std::abs(mcParticle1.eta()) > trackConfigs.etaMax) + continue; + + for (const auto& mcParticle2 : mcParticles) { + if (!mcParticle2.isPhysicalPrimary() || std::abs(mcParticle2.eta()) > trackConfigs.etaMax) + continue; + + if (!(mcParticle1.pdgCode() == PDG_t::kKPlus && mcParticle2.pdgCode() == PDG_t::kKMinus) && + !(mcParticle1.pdgCode() == PDG_t::kKMinus && mcParticle2.pdgCode() == PDG_t::kKPlus)) + continue; + + ROOT::Math::PxPyPzMVector genKPair = recMother(mcParticle1, mcParticle2, massKa, massKa); + + if (genKPair.Pt() < phiConfigs.minPhiPt) + continue; + if (genKPair.M() > phiConfigs.maxMPhi) + continue; + if (std::abs(genKPair.Rapidity()) > phiConfigs.cfgYAcceptance) + continue; + + phimesonCandidates(mcCollision.globalIndex(), genKPair.M(), genKPair.Pt(), genKPair.Rapidity(), genKPair.Phi()); + } + } + } + + PROCESS_SWITCH(PhiMesonCandProducer, processMCGen, "Process function to select Phi meson candidates in MCGen", false); +}; + +struct PhiMesonSelCollision { + // Produce the table with the event selection information + Produces phimesonSelection; + + HistogramRegistry histos{"eventSelection", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + + // Configurable for selection type + Configurable selectionType{"selectionType", 1, "Selection type: 0 - default selection only, 1 - default + phi meson selection"}; + + // Configurable for event selection + Configurable cutZVertex{"cutZVertex", 10.0f, "Accepted z-vertex range (cm)"}; + + // Configurable on multiplicity bins + Configurable> binsMult{"binsMult", {0.0, 1.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 70.0, 100.0}, "Multiplicity bin limits"}; + + // Configurable for MC + Configurable cfgiskNoITSROFrameBorder{"cfgiskNoITSROFrameBorder", false, "kNoITSROFrameBorder request on MC collisions"}; + + // Configurables on phi selection + struct : ConfigurableGroup { + Configurable minMPhiSignal{"minMPhiSignal", 1.0095f, "Upper limits on Phi mass for signal extraction"}; + Configurable maxMPhiSignal{"maxMPhiSignal", 1.029f, "Upper limits on Phi mass for signal extraction"}; + } phiConfigs; + + // Defining the type of the collisions for data and MC + using SelCollisions = soa::Join; + using SimCollisions = soa::Join; + using MCCollisions = soa::Join; + + // Defining the type of the phi's daughter tracks for data and MC + using FullTracks = soa::Join; + using FullMCTracks = soa::Join; + + // Necessary service to flag INEL>0 events in GenMC + Service pdgDB; + + void init(InitContext&) + { + // Defining histogram axes + AxisSpec vertexZAxis = {100, -cutZVertex, cutZVertex, "vrtx_{Z} [cm]"}; + AxisSpec multAxis = {120, 0.0f, 120.0f, "centFT0M"}; + AxisSpec binnedmultAxis{(std::vector)binsMult, "centFT0M"}; + + // Booking histograms for event selection QA + // Number of events per selection in Data + histos.add("hEventSelectionData", "hEventSelectionData", kTH1F, {{5, -0.5f, 4.5f}}); + histos.get(HIST("hEventSelectionData"))->GetXaxis()->SetBinLabel(1, "All collisions"); + histos.get(HIST("hEventSelectionData"))->GetXaxis()->SetBinLabel(2, "sel8 cut"); + histos.get(HIST("hEventSelectionData"))->GetXaxis()->SetBinLabel(3, "posZ cut"); + histos.get(HIST("hEventSelectionData"))->GetXaxis()->SetBinLabel(4, "INEL>0 cut"); + histos.get(HIST("hEventSelectionData"))->GetXaxis()->SetBinLabel(5, "With at least a #phi cand"); + + // Number of MC events per selection in MC + histos.add("hEventSelectionMC", "hEventSelectionMC", kTH1F, {{8, -0.5f, 7.5f}}); + histos.get(HIST("hEventSelectionMC"))->GetXaxis()->SetBinLabel(1, "All collisions"); + histos.get(HIST("hEventSelectionMC"))->GetXaxis()->SetBinLabel(2, "kIsTriggerTVX"); + histos.get(HIST("hEventSelectionMC"))->GetXaxis()->SetBinLabel(3, "kNoTimeFrameBorder"); + histos.get(HIST("hEventSelectionMC"))->GetXaxis()->SetBinLabel(4, "kNoITSROFrameBorder"); + histos.get(HIST("hEventSelectionMC"))->GetXaxis()->SetBinLabel(5, "posZ cut"); + histos.get(HIST("hEventSelectionMC"))->GetXaxis()->SetBinLabel(6, "INEL>0 cut"); + histos.get(HIST("hEventSelectionMC"))->GetXaxis()->SetBinLabel(7, "With at least a gen coll"); + histos.get(HIST("hEventSelectionMC"))->GetXaxis()->SetBinLabel(8, "With at least a #phi cand"); + + // Event information + histos.add("hVertexZ", "Vertex Z", kTH1F, {vertexZAxis}); + histos.add("hVertexZWPhi", "Vertex Z with a Phi Candidate", kTH1F, {vertexZAxis}); + histos.add("hMultiplicityPercent", "Multiplicity Percentile", kTH1F, {multAxis}); + histos.add("hMultiplicityPercentWPhi", "Multiplicity Percentile in Events with a Phi Candidate", kTH1F, {multAxis}); + histos.add("h2VertexZvsMult", "Vertex Z vs Multiplicity Percentile", kTH2F, {vertexZAxis, binnedmultAxis}); + histos.add("h2VertexZvsMultWPhi", "Vertex Z vs Multiplicity Percentile with a Phi Candidate", kTH2F, {vertexZAxis, binnedmultAxis}); + + // Phi's daughter tracks information + /*histos.add("hEta", "Eta of Kaon candidates", kTH1F, {{100, -1.0f, 1.0f}}); + histos.add("hNsigmaKaonTPC", "NsigmaKaon TPC distribution vs pt", kTH2F, {{100, 0.0, 5.0, "#it{p} (GeV/#it{c})"}, {100, -10.0f, 10.0f}}); + histos.add("hNsigmaKaonTOF", "NsigmaKaon TOF distribution vs pt", kTH2F, {{100, 0.0, 5.0, "#it{p} (GeV/#it{c})"}, {100, -10.0f, 10.0f}}); + histos.add("h2DauTracksPhiDCAxy", "DCAxy distribution vs pt", kTH2F, {{100, 0.0, 5.0, "#it{p}_{T} (GeV/#it{c})"}, {2000, -0.05, 0.05, "DCA_{xy} (cm)"}}); + histos.add("h2DauTracksPhiDCAz", "DCAz distribution vs pt", kTH2F, {{100, 0.0, 5.0, "#it{p}_{T} (GeV/#it{c})"}, {2000, -0.05, 0.05, "DCA_{xy} (cm)"}});*/ + } + + // Default event selection + template + bool defaultEventSelection(const T& collision) + { + float multPercentile{0.0f}; + + if constexpr (!isMC) { // data event + histos.fill(HIST("hEventSelectionData"), 0); // all collisions + if (!collision.sel8()) + return false; + histos.fill(HIST("hEventSelectionData"), 1); // sel8 collisions + if (std::abs(collision.posZ()) >= cutZVertex) + return false; + histos.fill(HIST("hEventSelectionData"), 2); // vertex-Z selected + if (!collision.isInelGt0()) + return false; + histos.fill(HIST("hEventSelectionData"), 3); // INEL>0 collisions + + multPercentile = collision.centFT0M(); + } else { // MCreco event + static_assert(!std::is_same_v, "Need to set MC_T to MCCollisions for isMC = true"); + + histos.fill(HIST("hEventSelectionMC"), 0); // all collisions + if (!collision.selection_bit(aod::evsel::kIsTriggerTVX)) + return false; + histos.fill(HIST("hEventSelectionMC"), 1); // kIsTriggerTVX collisions + if (!collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) + return false; + histos.fill(HIST("hEventSelectionMC"), 2); // kNoTimeFrameBorder collisions + if (cfgiskNoITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) + return false; + histos.fill(HIST("hEventSelectionMC"), 3); // kNoITSROFrameBorder collisions (by default not requested by the selection) + if (std::abs(collision.posZ()) > cutZVertex) + return false; + histos.fill(HIST("hEventSelectionMC"), 4); // vertex-Z selected + if (!collision.isInelGt0()) + return false; + histos.fill(HIST("hEventSelectionMC"), 5); // INEL>0 collisions + if (!collision.has_mcCollision()) + return false; + histos.fill(HIST("hEventSelectionMC"), 6); // with at least a gen collision + + const auto& mcCollision = collision.template mcCollision_as(); + multPercentile = mcCollision.centFT0M(); + } + + histos.fill(HIST("hVertexZ"), collision.posZ()); + histos.fill(HIST("hMultiplicityPercent"), multPercentile); + histos.fill(HIST("h2VertexZvsMult"), collision.posZ(), multPercentile); + + return true; + } + + // Check if the event has at least one phi candidate + template + bool eventHasPhi(const T1& collision, const T2& phiCandidates) + { + uint16_t nPhi{0}; + + for (const auto& phiCand : phiCandidates) { + + if (phiCand.inMassRegion(phiConfigs.minMPhiSignal, phiConfigs.maxMPhiSignal)) + nPhi++; + + // histos.fill(HIST("hEta"), track1.eta()); + // histos.fill(HIST("hNsigmaKaonTPC"), track1.tpcInnerParam(), track1.tpcNSigmaKa()); + // histos.fill(HIST("hNsigmaKaonTOF"), track1.tpcInnerParam(), track1.tofNSigmaKa()); + // histos.fill(HIST("h2DauTracksPhiDCAxy"), track1.pt(), track1.dcaXY()); + // histos.fill(HIST("h2DauTracksPhiDCAz"), track1.pt(), track1.dcaZ()); + } + + if (nPhi == 0) + return false; + + float multPercentile{0.0f}; + + if constexpr (!isMC) { + histos.fill(HIST("hEventSelectionData"), 4); + + multPercentile = collision.centFT0M(); + } else { + if constexpr (std::is_same_v) { + multPercentile = collision.centFT0M(); + } else { + histos.fill(HIST("hEventSelectionMC"), 7); + + const auto& mcCollision = collision.template mcCollision_as(); + multPercentile = mcCollision.centFT0M(); + } + } + + histos.fill(HIST("hVertexZWPhi"), collision.posZ()); + histos.fill(HIST("hMultiplicityPercentWPhi"), multPercentile); + histos.fill(HIST("h2VertexZvsMultWPhi"), collision.posZ(), multPercentile); + + return true; + } + + void processData(SelCollisions::iterator const& collision, aod::PhimesonCandidates const& phiCandidates) + { + phimesonSelection(defaultEventSelection(collision) && selectionType == 1 ? eventHasPhi(collision, phiCandidates) : true); + } + + PROCESS_SWITCH(PhiMesonSelCollision, processData, "Process function to select events with Phi mesons in Data", true); + + void processMCReco(SimCollisions::iterator const& collision, MCCollisions const&, aod::PhimesonCandidates const& phiCandidates) + { + phimesonSelection(defaultEventSelection(collision) && selectionType == 1 ? eventHasPhi(collision, phiCandidates) : true); + } + + PROCESS_SWITCH(PhiMesonSelCollision, processMCReco, "Process function to select events with Phi mesons in MCReco", false); + + void processMCGen(MCCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles, aod::PhimesonCandidates const& phiCandidates) + { + phimesonSelection(pwglf::isINELgt0mc(mcParticles, pdgDB) && selectionType == 1 ? eventHasPhi(mcCollision, phiCandidates) : true); + } + + PROCESS_SWITCH(PhiMesonSelCollision, processMCGen, "Process function to select events with Phi mesons in MCGen", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGLF/Tasks/Strangeness/CMakeLists.txt b/PWGLF/Tasks/Strangeness/CMakeLists.txt index 333e2600d08..153fa47915c 100644 --- a/PWGLF/Tasks/Strangeness/CMakeLists.txt +++ b/PWGLF/Tasks/Strangeness/CMakeLists.txt @@ -89,6 +89,11 @@ o2physics_add_dpl_workflow(phik0shortanalysis PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(phi-strange-correlation + SOURCES phiStrangeCorrelation.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(lambdapolarization SOURCES lambdapolarization.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx b/PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx new file mode 100644 index 00000000000..bfa459f6502 --- /dev/null +++ b/PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx @@ -0,0 +1,708 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file phiStrangeCorrelation.cxx +/// \brief Analysis task for phi-strangeness rapidity correlations analysis +/// \author Stefano Cannito (stefano.cannito@cern.ch) + +#include "PWGLF/DataModel/LFPhiStrangeCorrelationTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/mcCentrality.h" +#include "PWGLF/Utils/inelGt.h" + +#include "Common/Core/TableHelper.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/TrackSelectionDefaults.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CCDB/BasicCCDBManager.h" +#include "CommonConstants/PhysicsConstants.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +static constexpr std::array phiMassRegionLabels{"Signal", "Sideband"}; + +enum ParticleOfInterest { + Phi = 0, + K0S, + PionTPC, + PionTPCTOF +}; + +static constexpr std::array particleOfInterestLabels{"Phi", "K0S", "PionTPC", "PionTPCTOF"}; + +/* +#define LIST_OF_PARTICLES_OF_INTEREST \ + X(Phi) \ + X(K0S) \ + X(PionTPC) \ + X(PionTPCTOF) + +enum ParticleOfInterest { +#define X(name) name, + LIST_OF_PARTICLES_OF_INTEREST +#undef X +}; + +static constexpr auto particleOfInterestLabels = std::to_array({ +#define X(name) #name, + LIST_OF_PARTICLES_OF_INTEREST +#undef X +}); +*/ + +struct BoundEfficiencyMap { + using CoordsTuple = std::tuple; + + const TH3* effMap; + CoordsTuple coords; + + BoundEfficiencyMap(const std::shared_ptr& effMap, float x, float y, float z) : effMap(effMap.get()), coords(x, y, z) {} + BoundEfficiencyMap(const std::shared_ptr& effMap, const CoordsTuple& coords) : effMap(effMap.get()), coords(coords) {} + + float getBinEfficiency() const + { + if (!effMap) { + return 1.0f; + } + + const auto& [x, y, z] = coords; + return effMap->GetBinContent(effMap->FindFixBin(x, y, z)); + } + + float interpolateEfficiency() const + { + if (!effMap) { + return 1.0f; + } + + const auto& [x, y, z] = coords; + return effMap->Interpolate(x, y, z); + } +}; + +/* +struct AnalysisRegion { + std::string suffix; + float minMass; + float maxMass; +}; +*/ + +struct PhiStrangenessCorrelation { + HistogramRegistry histos{"phiStrangenessCorrelation", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + + // Configurable for selection type + Configurable selectionType{"selectionType", 1, "Selection type: 0 - default selection only, 1 - default + phi meson selection"}; + + // Configurable for analysis mode + Configurable analysisMode{"analysisMode", 1, "Analysis mode: 0 - old method with online normalization, 1 - new method with correlations"}; + + // Configurable for event selection + Configurable cutZVertex{"cutZVertex", 10.0f, "Accepted z-vertex range (cm)"}; // TO BE REMOVED + + // Configurable on multiplicity bins + Configurable> binsMult{"binsMult", {0.0, 1.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 70.0, 100.0}, "Multiplicity bin limits"}; + + // Configurables for tracks selection + struct : ConfigurableGroup { + Configurable cfgGlobalWoDCATrack{"cfgGlobalWoDCATrack", true, "Global track selection without DCA"}; + Configurable cfgPVContributor{"cfgPVContributor", true, "PV contributor track selection"}; + Configurable cMinKaonPtcut{"cMinKaonPtcut", 0.15f, "Track minimum pt cut"}; + Configurable etaMax{"etaMax", 0.8f, "eta max"}; + Configurable pTToUseTOF{"pTToUseTOF", 0.5f, "pT above which use TOF"}; + Configurable cMaxDCAzToPVcut{"cMaxDCAzToPVcut", 2.0f, "Track DCAz cut to PV Maximum"}; + Configurable> cMaxDCArToPVPhi{"cMaxDCArToPVPhi", {0.004f, 0.013f, 1.0f}, "Track DCAr cut to PV for Phi"}; + + Configurable cfgIsTOFChecked{"cfgIsTOFChecked", true, "Is TOF checked in PID for pions"}; + Configurable> cMaxDCArToPVPion{"cMaxDCArToPVPion", {0.004f, 0.013f, 1.0f}, "Track DCAr cut to PV for Pions"}; + Configurable cfgIsDCAzParameterized{"cfgIsDCAzParameterized", false, "IsDCAzParameterized"}; + Configurable> cMaxDCAzToPVPion{"cMaxDCAzToPVPion", {0.004f, 0.013f, 1.0f}, "Track DCAz cut to PV for Pions"}; + + Configurable nSigmaCutTPCKa{"nSigmaCutTPCKa", 2.0f, "Value of the TPC Nsigma cut for Kaons"}; + Configurable nSigmaCutCombinedKa{"nSigmaCutCombinedKa", 2.0f, "Value of the TPC and TOF Nsigma cut for Kaons"}; + + Configurable nSigmaCutTPCPrimPion{"nSigmaCutTPCPrimPion", 2.0f, "Value of the TPC Nsigma cut for primary Pions"}; + Configurable nSigmaCutTPCSecPion{"nSigmaCutTPCSecPion", 4.0f, "Value of the TPC Nsigma cut for secondary Pions"}; + Configurable nSigmaCutCombinedPi{"nSigmaCutCombinedPi", 2.0f, "Value of the TPC and TOF Nsigma cut for Pions"}; + Configurable cMinPionPtcut{"cMinPionPtcut", 0.2f, "Track minimum pt cut"}; + + Configurable minTPCnClsFound{"minTPCnClsFound", 70, "min number of found TPC clusters"}; + Configurable minNCrossedRowsTPC{"minNCrossedRowsTPC", 70, "min number of TPC crossed rows"}; + Configurable maxChi2TPC{"maxChi2TPC", 4.0f, "max chi2 per cluster TPC"}; + Configurable minITSnCls{"minITSnCls", 4, "min number of ITS clusters"}; + Configurable maxChi2ITS{"maxChi2ITS", 36.0f, "max chi2 per cluster ITS"}; + } trackConfigs; + + // Configurables on phi selection + struct : ConfigurableGroup { + Configurable minPhiPt{"minPhiPt", 0.4f, "Minimum pT for Phi candidates"}; + Configurable> rangeMPhiSignal{"rangeMPhiSignal", {1.0095f, 1.029f}, "Phi mass range for signal extraction"}; + Configurable> rangeMPhiSideband{"rangeMPhiSideband", {1.1f, 1.2f}, "Phi mass range for sideband extraction"}; + } phiConfigs; + + // Configurables on phi pT bins + Configurable> binspTPhi{"binspTPhi", {0.4, 0.8, 1.4, 2.0, 2.8, 4.0, 6.0, 10.0}, "pT bin limits for Phi"}; + + // Configurables for V0 selection + struct : ConfigurableGroup { + Configurable v0SettingCosPA{"v0SettingCosPA", 0.98f, "V0 CosPA"}; + Configurable v0SettingRadius{"v0SettingRadius", 0.5f, "v0radius"}; + Configurable v0SettingDCAV0Dau{"v0SettingDCAV0Dau", 1.0f, "DCA V0 Daughters"}; + Configurable v0SettingDCAPosToPV{"v0SettingDCAPosToPV", 0.1f, "DCA Pos To PV"}; + Configurable v0SettingDCANegToPV{"v0SettingDCANegToPV", 0.1f, "DCA Neg To PV"}; + Configurable v0SettingMinPt{"v0SettingMinPt", 0.1f, "V0 min pt"}; + + Configurable cfgFurtherV0Selection{"cfgFurtherV0Selection", false, "Further V0 selection"}; + Configurable ctauK0s{"ctauK0s", 20.0f, "C tau K0s(cm)"}; + Configurable paramArmenterosCut{"paramArmenterosCut", 0.2f, "parameter Armenteros Cut"}; + Configurable v0rejK0s{"v0rejK0s", 0.005f, "V0 rej K0s"}; + + Configurable lowMK0S{"lowMK0S", 0.48f, "Lower limit on K0Short mass"}; + Configurable upMK0S{"upMK0S", 0.52f, "Upper limit on K0Short mass"}; + } v0Configs; + + // Configurable on K0S pT bins + Configurable> binspTK0S{"binspTK0S", {0.1, 0.5, 0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 4.0, 6.0}, "pT bin limits for K0S"}; + + // Configurable on pion pT bins + Configurable> binspTPi{"binspTPi", {0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.5, 2.0, 3.0}, "pT bin limits for pions"}; + + // Configurables for delta y selection + struct : ConfigurableGroup { + Configurable nBinsY{"nBinsY", 20, "Number of bins in y axis"}; + Configurable nBinsDeltaY{"nBinsDeltaY", 20, "Number of bins in deltay axis"}; + Configurable cfgYAcceptance{"cfgYAcceptance", 0.5f, "Rapidity acceptance"}; + Configurable> cfgDeltaYAcceptanceBins{"cfgDeltaYAcceptanceBins", {0.5f}, "Rapidity acceptance bins"}; + } yConfigs; + + // Configurables to apply efficiency online and how to + Configurable applyEfficiency{"applyEfficiency", false, "Use efficiency for filling histograms"}; + Configurable useEffInterpolation{"useEffInterpolation", false, "If true, interpolates efficiency map, else uses bin center"}; + + // Configurable for event mixing + Configurable cfgNoMixedEvents{"cfgNoMixedEvents", 5, "Number of mixed events per event"}; + + // Configurables for CCDB + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository to use"}; + Configurable ccdbEfficiencyPath{"ccdbEfficiencyPath", "Users/s/scannito/Efficiencies", "Correction path to file"}; + + // Constants + double massPi = o2::constants::physics::MassPiPlus; + double massK0S = o2::constants::physics::MassK0Short; + double massLambda = o2::constants::physics::MassLambda0; + + // Filter on phi selected collisions + Filter collisionFilter = aod::lf_selection_phi_collision::phimesonSel == true; + + // Defining filters on V0s (cannot filter on dynamic columns) + Filter v0PreFilter = (nabs(aod::v0data::dcapostopv) > v0Configs.v0SettingDCAPosToPV && nabs(aod::v0data::dcanegtopv) > v0Configs.v0SettingDCANegToPV && aod::v0data::dcaV0daughters < v0Configs.v0SettingDCAV0Dau); + + // Defining the type of the collisions for data and MC + using SelCollisions = soa::Filtered>; + using SimCollisions = soa::Join; + using MCCollisions = soa::Filtered>; + + // Defining the type of the V0s and corresponding daughter tracks for data and MC + using FullV0s = soa::Filtered; + using FullMCV0s = soa::Join; + + using V0DauTracks = soa::Join; + using V0DauMCTracks = soa::Join; + + // Defining the type of the tracks for data and MC + using FullTracks = soa::Join; + using FullMCTracks = soa::Join; + + // using FilteredTracks = soa::Filtered; + // using FilteredMCTracks = soa::Filtered; + + // Preslice for manual slicing + struct : PresliceGroup { + PresliceUnsorted collPerMCCollision = aod::mccollisionlabel::mcCollisionId; + Preslice phiCandPerCollision = aod::lf_selection_phi_candidate::collisionId; + Preslice v0PerCollision = aod::v0::collisionId; + Preslice trackPerCollision = aod::track::collisionId; + // Preslice mcPartPerMCCollision = aod::mcparticle::mcCollisionId; + } preslices; + + // Necessary service to retrieve efficiency maps from CCDB + Service ccdb; + + // Efficiency maps + /*std::shared_ptr effMapPhi = nullptr; + std::shared_ptr effMapK0S = nullptr; + std::shared_ptr effMapPionTPC = nullptr; + std::shared_ptr effMapPionTPCTOF = nullptr;*/ + + std::array, 4> effMaps{}; + + void init(InitContext&) + { + AxisSpec vertexZAxis = {100, -cutZVertex, cutZVertex, "vrtx_{Z} [cm]"}; // TO BE REMOVED + AxisSpec yAxis = {yConfigs.nBinsY, -yConfigs.cfgYAcceptance, yConfigs.cfgYAcceptance, "#it{y}"}; + AxisSpec deltayAxis = {yConfigs.nBinsDeltaY, -1.0f, 1.0f, "#Delta#it{y}"}; + AxisSpec deltaphiAxis = {72, -o2::constants::math::PIHalf, o2::constants::math::PIHalf * 3, "#Delta#varphi"}; + AxisSpec multAxis = {120, 0.0f, 120.0f, "centFT0M"}; + AxisSpec binnedmultAxis{(std::vector)binsMult, "centFT0M"}; + AxisSpec massPhiAxis = {200, 0.9f, 1.2f, "#it{M}_{inv} [GeV/#it{c}^{2}]"}; + AxisSpec pTPhiAxis = {120, 0.0f, 12.0f, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec binnedpTPhiAxis{(std::vector)binspTPhi, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec pTK0SAxis = {100, 0.0f, 10.0f, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec binnedpTK0SAxis{(std::vector)binspTK0S, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec pTPiAxis = {50, 0.0f, 5.0f, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec binnedpTPiAxis{(std::vector)binspTPi, "#it{p}_{T} (GeV/#it{c})"}; + + histos.add("phi/h3PhiData", "Invariant mass of Phi in Data", kTH3F, {binnedmultAxis, binnedpTPhiAxis, massPhiAxis}); + for (const auto& label : phiMassRegionLabels) { + histos.add(fmt::format("phiK0S/h5PhiK0SData2PartCorr{}", label).c_str(), "Deltay vs deltaphi for Phi and K0Short in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTK0SAxis, deltayAxis, deltaphiAxis}); + histos.add(fmt::format("phiPi/h5PhiPiData2PartCorr{}", label).c_str(), "Deltay vs deltaphi for Phi and Pion in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, deltaphiAxis}); + } + + // histos.add("phiK0S/h5PhiK0SDataNewProc", "2D Invariant mass of Phi and K0Short in Data", kTHnSparseF, {deltayAxis, binnedmultAxis, binnedpTK0SAxis, massK0SAxis, massPhiAxis}); + // histos.add("phiPi/h5PhiPiTPCDataNewProc", "Phi Invariant mass vs Pion nSigma TPC in Data", kTHnSparseF, {deltayAxis, binnedmultAxis, binnedpTPiAxis, nSigmaPiAxis, massPhiAxis}); + // histos.add("phiPi/h5PhiPiTOFDataNewProc", "Phi Invariant mass vs Pion nSigma TOF in Data", kTHnSparseF, {deltayAxis, binnedmultAxis, binnedpTPiAxis, nSigmaPiAxis, massPhiAxis}); + + histos.add("event/hRecoMCMultiplicityPercent", "RecoMC Multiplicity Percentile", kTH1F, {binnedmultAxis}); + histos.add("event/h2RecoMCVertexZvsMult", "RecoMC Vertex Z vs Multiplicity Percentile", kTH2F, {vertexZAxis, binnedmultAxis}); + histos.add("event/hSplitVertexZ", "Split in z-vtx", kTH1F, {{100, -5.0f, 5.0f}}); + histos.add("event/hGenMCMultiplicityPercent", "Generated MC Multiplicity Percentile", kTH1F, {binnedmultAxis}); + histos.add("event/hGenMCAssocRecoMultiplicityPercent", "Generated MC associated Multiplicity Percentile", kTH1F, {binnedmultAxis}); + histos.add("event/h2GenMCAssocRecoVertexZvsMult", "Generated MC associated reco Vertex Z vs multiplicity", kTH2F, {vertexZAxis, binnedmultAxis}); + + histos.add("phi/h4PhiMCReco", "Phi in MC Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTPhiAxis, yAxis}); + histos.add("phi/h3PhiMCGen", "Phi in MC Gen", kTH3F, {binnedmultAxis, binnedpTPhiAxis, yAxis}); + histos.add("phi/h4PhiMCGenAssocReco", "Phi in MC Gen Assoc Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTPhiAxis, yAxis}); + + histos.add("k0s/h4K0SMCReco", "K0S in MC Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTK0SAxis, yAxis}); + histos.add("k0s/h3K0SMCGen", "K0S in MC Gen", kTH3F, {binnedmultAxis, binnedpTK0SAxis, yAxis}); + histos.add("k0s/h4K0SMCGenAssocReco", "K0S in MC Gen Assoc Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTK0SAxis, yAxis}); + + histos.add("pi/h4PiMCReco", "Pion in MC Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTPiAxis, yAxis}); + histos.add("pi/h3PiMCGen", "Pion in MC Gen", kTH3F, {binnedmultAxis, binnedpTPiAxis, yAxis}); + histos.add("pi/h4PiMCGenAssocReco", "Pion in MC Gen Assoc Reco", kTHnSparseF, {vertexZAxis, binnedmultAxis, binnedpTPiAxis, yAxis}); + + histos.add("pi/h2RecMCDCAxyPrimPi", "Dcaxy distribution vs pt for Primary Pions", kTH2F, {binnedpTPiAxis, {2000, -0.05, 0.05, "DCA_{xy} (cm)"}}); + histos.add("pi/h2RecMCDCAxySecWeakDecayPi", "Dcaz distribution vs pt for Secondary Pions from Weak Decay", kTH2F, {binnedpTPiAxis, {2000, -0.05, 0.05, "DCA_{xy} (cm)"}}); + histos.add("pi/h2RecMCDCAxySecMaterialPi", "Dcaxy distribution vs pt for Secondary Pions from Material", kTH2F, {binnedpTPiAxis, {2000, -0.05, 0.05, "DCA_{xy} (cm)"}}); + + // Load efficiency maps from CCDB + if (applyEfficiency) { + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + // getEfficiencyMapsFromCCDB(); + + for (int i = 0; i < 4; ++i) { + loadEfficiencyMapFromCCDB(static_cast(i)); + } + } + } + + void loadEfficiencyMapFromCCDB(ParticleOfInterest poi) + { + effMaps[poi] = std::shared_ptr(ccdb->get(fmt::format("{}/h3EffMap{}", ccdbEfficiencyPath.value, particleOfInterestLabels[poi]))); + if (!effMaps[poi]) + LOG(fatal) << "Could not load efficiency map for " << particleOfInterestLabels[poi] << "!"; + LOG(info) << "Efficiency map for " << particleOfInterestLabels[poi] << " loaded from CCDB"; + } + + /* + void getEfficiencyMapsFromCCDB() + { + } + */ + + // Compute weight based on efficiencies + template + float computeWeight(const BoundEffMaps&... boundEffMaps) + { + if (!applyEfficiency) + return 1.0f; + + float totalEfficiency = ((useEffInterpolation ? boundEffMaps.interpolateEfficiency() : boundEffMaps.getBinEfficiency()) * ...); + + return totalEfficiency <= 0.0f ? 1.0f : 1.0f / totalEfficiency; + } + + // Single track selection for strangeness sector + template + bool selectionTrackStrangeness(const T& track) + { + if (!track.hasTPC()) + return false; + if (track.tpcNClsFound() < trackConfigs.minTPCnClsFound) + return false; + if (track.tpcNClsCrossedRows() < trackConfigs.minNCrossedRowsTPC) + return false; + if (track.tpcChi2NCl() > trackConfigs.maxChi2TPC) + return false; + + if (std::abs(track.eta()) > trackConfigs.etaMax) + return false; + return true; + } + + // V0 selection + template + bool selectionV0(const T1& v0, const T2& collision) + { + const auto& posDaughterTrack = v0.template posTrack_as(); + const auto& negDaughterTrack = v0.template negTrack_as(); + + if (!selectionTrackStrangeness(posDaughterTrack) || !selectionTrackStrangeness(negDaughterTrack)) + return false; + + if constexpr (!isMC) { + if (std::abs(posDaughterTrack.tpcNSigmaPi()) > trackConfigs.nSigmaCutTPCSecPion) + return false; + if (std::abs(negDaughterTrack.tpcNSigmaPi()) > trackConfigs.nSigmaCutTPCSecPion) + return false; + } + + if (v0.v0cosPA() < v0Configs.v0SettingCosPA) + return false; + if (v0.v0radius() < v0Configs.v0SettingRadius) + return false; + if (v0.pt() < v0Configs.v0SettingMinPt) + return false; + + if (v0Configs.cfgFurtherV0Selection) { + if (v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * massK0S > v0Configs.ctauK0s) + return false; + if (v0.qtarm() < (v0Configs.paramArmenterosCut * std::abs(v0.alpha()))) + return false; + if (std::abs(v0.mLambda() - massLambda) < v0Configs.v0rejK0s) + return false; + } + + if (std::abs(v0.yK0Short()) > yConfigs.cfgYAcceptance) + return false; + + return true; + } + + // Topological selection for pions + template + bool selectionPion(const T& track) + { + if (!track.isGlobalTrackWoDCA()) + return false; + + if (track.itsNCls() < trackConfigs.minITSnCls) + return false; + if (track.tpcNClsFound() < trackConfigs.minTPCnClsFound) + return false; + + if (track.pt() < trackConfigs.cMinPionPtcut) + return false; + + if (std::abs(track.dcaXY()) > trackConfigs.cMaxDCArToPVPion->at(0) + (trackConfigs.cMaxDCArToPVPion->at(1) / std::pow(track.pt(), trackConfigs.cMaxDCArToPVPion->at(2)))) + return false; + if (trackConfigs.cfgIsDCAzParameterized) { + if (std::abs(track.dcaZ()) > trackConfigs.cMaxDCAzToPVPion->at(0) + (trackConfigs.cMaxDCAzToPVPion->at(1) / std::pow(track.pt(), trackConfigs.cMaxDCAzToPVPion->at(2)))) + return false; + } else { + if (std::abs(track.dcaZ()) > trackConfigs.cMaxDCAzToPVcut) + return false; + } + + if (trackConfigs.cfgIsTOFChecked && track.pt() >= trackConfigs.pTToUseTOF && !track.hasTOF()) + return false; + + if (analysisMode == 1) { + if (track.pt() < trackConfigs.pTToUseTOF && std::abs(track.tpcNSigmaPi()) >= trackConfigs.nSigmaCutTPCPrimPion) + return false; + if (trackConfigs.cfgIsTOFChecked && track.pt() >= trackConfigs.pTToUseTOF && (std::pow(track.tofNSigmaPi(), 2) + std::pow(track.tpcNSigmaPi(), 2)) >= std::pow(trackConfigs.nSigmaCutCombinedPi, 2)) + return false; + } + + if (std::abs(track.rapidity(massPi)) > yConfigs.cfgYAcceptance) + return false; + + return true; + } + + /* + void processPhiK0SPionDeltayDeltaphiData2D(SelCollisions::iterator const& collision, aod::PhimesonCandidates const& phiCandidates, FullTracks const& fullTracks, FullV0s const& V0s, V0DauTracks const&) + { + float multiplicity = collision.centFT0M(); + + std::vector analysisRegions = { + {"Signal", phiConfigs.rangeMPhiSignal.first, phiConfigs.rangeMPhiSignal.second}, + {"Sideband", phiConfigs.rangeMPhiSideband.first, phiConfigs.rangeMPhiSideband.second}}; + + // Loop over all positive tracks + for (const auto& phiCand : phiCandidates) { + float weightPhi = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y())); + + histos.fill(HIST("phi/h3PhiData"), multiplicity, phiCand.pt(), phiCand.m(), weightPhi); + + for (const auto& region : analysisRegions) { + if (!phiCand.inMassRegion(region.minMass, region.maxMass)) + continue; + + // V0 already reconstructed by the builder + for (const auto& v0 : V0s) { + // Cut on V0 dynamic columns + if (!selectionV0(v0, collision)) + continue; + + float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + BoundEfficiencyMap(effMapK0S, multiplicity, v0.pt(), v0.yK0Short())); + + histos.fill(HIST("phiK0S/h5PhiK0SData2PartCorr"), multiplicity, phiCand.pt(), v0.pt(), phiCand.y() - v0.yK0Short(), phiCand.phi() - v0.phi(), weightPhiK0S); + } + + // Loop over all primary pion candidates + for (const auto& track : fullTracks) { + if (!selectionPion(track)) + continue; + + float weightPhiPion = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + track.pt() < trackConfigs.pTToUseTOF ? BoundEfficiencyMap(effMapPionTPC, multiplicity, track.pt(), track.rapidity(massPi)) : BoundEfficiencyMap(effMapPionTPCTOF, multiplicity, track.pt(), track.rapidity(massPi))); + + histos.fill(HIST("phiPi/h5PhiPiData2PartCorr"), multiplicity, phiCand.pt(), track.pt(), phiCand.y() - track.rapidity(massPi), phiCand.phi() - track.phi(), weightPhiPion); + } + } + } + } + + PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiK0SPionDeltayDeltaphiData2D, "Process function for Phi-K0S and Phi-Pion Deltay and Deltaphi 2D Correlations in Data", true); + */ + + void processPhiK0SPionDeltayDeltaphiData2D(SelCollisions::iterator const& collision, aod::PhimesonCandidates const& phiCandidates, FullTracks const& fullTracks, FullV0s const& V0s, V0DauTracks const&) + { + float multiplicity = collision.centFT0M(); + + const std::array, 2> phiMassRegions = {phiConfigs.rangeMPhiSignal, phiConfigs.rangeMPhiSideband}; + + // Loop over all positive tracks + for (const auto& phiCand : phiCandidates) { + float weightPhi = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y())); + + // float weightPhi = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y())); + + histos.fill(HIST("phi/h3PhiData"), multiplicity, phiCand.pt(), phiCand.m(), weightPhi); + + static_for<0, phiMassRegionLabels.size() - 1>([&](auto i_idx) { + constexpr unsigned int i = i_idx.value; + + const auto& [minMass, maxMass] = phiMassRegions[i]; + if (!phiCand.inMassRegion(minMass, maxMass)) + return; + + // V0 already reconstructed by the builder + for (const auto& v0 : V0s) { + // Cut on V0 dynamic columns + if (!selectionV0(v0, collision)) + continue; + + float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()), + BoundEfficiencyMap(effMaps[K0S], multiplicity, v0.pt(), v0.yK0Short())); + + /*float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + BoundEfficiencyMap(effMapK0S, multiplicity, v0.pt(), v0.yK0Short()));*/ + + histos.fill(HIST("phiK0S/h5PhiK0SData2PartCorr") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), v0.pt(), phiCand.y() - v0.yK0Short(), phiCand.phi() - v0.phi(), weightPhiK0S); + } + + // Loop over all primary pion candidates + for (const auto& track : fullTracks) { + if (!selectionPion(track)) + continue; + + auto Pion = track.pt() < trackConfigs.pTToUseTOF ? PionTPC : PionTPCTOF; + + float weightPhiPion = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()), + BoundEfficiencyMap(effMaps[Pion], multiplicity, track.pt(), track.rapidity(massPi))); + + /*auto effMapPion = track.pt() < trackConfigs.pTToUseTOF ? effMapPionTPC : effMapPionTPCTOF; + + float weightPhiPion = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()), + BoundEfficiencyMap(effMapPion, multiplicity, track.pt(), track.rapidity(massPi)));*/ + + histos.fill(HIST("phiPi/h5PhiPiData2PartCorr") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), track.pt(), phiCand.y() - track.rapidity(massPi), phiCand.phi() - track.phi(), weightPhiPion); + } + }); + } + } + + PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiK0SPionDeltayDeltaphiData2D, "Process function for Phi-K0S and Phi-Pion Deltay and Deltaphi 2D Correlations in Data", true); + + void processParticleEfficiency(MCCollisions::iterator const& mcCollision, SimCollisions const& collisions, FullMCTracks const& fullMCTracks, FullMCV0s const& V0s, V0DauMCTracks const&, aod::McParticles const& mcParticles, aod::PhimesonCandidates const& phiCandidates) + { + uint16_t numberAssocColls{0}; + std::vector zVtxs; + + const auto collsThisMCColl = collisions.sliceBy(preslices.collPerMCCollision, mcCollision.globalIndex()); + + for (const auto& collision : collsThisMCColl) { + histos.fill(HIST("event/hRecoMCMultiplicityPercent"), mcCollision.centFT0M()); + histos.fill(HIST("event/h2RecoMCVertexZvsMult"), collision.posZ(), mcCollision.centFT0M()); + + zVtxs.push_back(collision.posZ()); + + if (selectionType == 0) { + const auto phiCandidatesThisColl = phiCandidates.sliceBy(preslices.phiCandPerCollision, collision.globalIndex()); + for (const auto& phiCand : phiCandidatesThisColl) { + histos.fill(HIST("phi/h4PhiMCReco"), collision.posZ(), mcCollision.centFT0M(), phiCand.pt(), phiCand.y()); + } + } + + const auto v0sThisColl = V0s.sliceBy(preslices.v0PerCollision, collision.globalIndex()); + const auto fullMCTracksThisColl = fullMCTracks.sliceBy(preslices.trackPerCollision, collision.globalIndex()); + + for (const auto& v0 : v0sThisColl) { + if (!selectionV0(v0, collision)) + continue; + + if (!v0.has_mcParticle()) + continue; + + const auto& v0McParticle = mcParticles.rawIteratorAt(v0.mcParticleId()); + if (std::abs(v0McParticle.pdgCode()) != PDG_t::kK0Short || !v0McParticle.isPhysicalPrimary()) + continue; + + histos.fill(HIST("k0s/h4K0SMCReco"), collision.posZ(), mcCollision.centFT0M(), v0McParticle.pt(), v0McParticle.y()); + } + + for (const auto& track : fullMCTracksThisColl) { + if (!selectionPion(track)) + continue; + + if (!track.has_mcParticle()) + continue; + + const auto& trackMcParticle = mcParticles.rawIteratorAt(track.mcParticleId()); + if (std::abs(trackMcParticle.pdgCode()) != PDG_t::kPiPlus) + continue; + + if (trackMcParticle.isPhysicalPrimary()) { + histos.fill(HIST("pi/h2RecMCDCAxyPrimPi"), track.pt(), track.dcaXY()); + } else { + if (trackMcParticle.getProcess() == TMCProcess::kPDecay) { // Selection of secondary pions from weak decay + histos.fill(HIST("pi/h2RecMCDCAxySecWeakDecayPi"), track.pt(), track.dcaXY()); + } else { // Selection of secondary pions from material interactions + histos.fill(HIST("pi/h2RecMCDCAxySecMaterialPi"), track.pt(), track.dcaXY()); + } + continue; + } + + histos.fill(HIST("pi/h4PiMCReco"), collision.posZ(), mcCollision.centFT0M(), trackMcParticle.pt(), trackMcParticle.y()); + } + + numberAssocColls++; + } + + histos.fill(HIST("event/hGenMCMultiplicityPercent"), mcCollision.centFT0M()); + + const bool hasAssoc = (numberAssocColls > 0); + const float zVtxRef = hasAssoc ? zVtxs[0] : 0.0f; + + //////TOBECHANGED////// + if (hasAssoc) { + if (zVtxs.size() > 1) { + for (size_t i = 1; i < zVtxs.size(); ++i) { + histos.fill(HIST("event/hSplitVertexZ"), zVtxs[i] - zVtxRef); + } + } + + histos.fill(HIST("event/hGenMCAssocRecoMultiplicityPercent"), mcCollision.centFT0M()); + histos.fill(HIST("event/h2GenMCAssocRecoVertexZvsMult"), zVtxRef, mcCollision.centFT0M()); + } + /////////////////////// + + auto inYAcceptance = [&](const auto& mcParticle) { + return std::abs(mcParticle.y()) <= yConfigs.cfgYAcceptance; + }; + + auto fillGenHistos = [&](auto h3Key, auto h4Key, const auto& mcParticle) { + histos.fill(h3Key, mcCollision.centFT0M(), mcParticle.pt(), mcParticle.y()); + if (hasAssoc) + histos.fill(h4Key, zVtxRef, mcCollision.centFT0M(), mcParticle.pt(), mcParticle.y()); + }; + + for (const auto& mcParticle : mcParticles /*| std::views::filter(inYAcceptance)*/) { + if (!inYAcceptance(mcParticle)) + continue; + + switch (std::abs(mcParticle.pdgCode())) { + case o2::constants::physics::Pdg::kPhi: + if (selectionType == 0 && mcParticle.isPhysicalPrimary() && mcParticle.pt() >= phiConfigs.minPhiPt) + fillGenHistos(HIST("phi/h3PhiMCGen"), HIST("phi/h4PhiMCGenAssocReco"), mcParticle); + break; + case PDG_t::kK0Short: + if (mcParticle.isPhysicalPrimary() && mcParticle.pt() >= v0Configs.v0SettingMinPt) + fillGenHistos(HIST("k0s/h3K0SMCGen"), HIST("k0s/h4K0SMCGenAssocReco"), mcParticle); + break; + case PDG_t::kPiPlus: + if (mcParticle.isPhysicalPrimary() && mcParticle.pt() >= trackConfigs.cMinPionPtcut) + fillGenHistos(HIST("pi/h3PiMCGen"), HIST("pi/h4PiMCGenAssocReco"), mcParticle); + break; + default: + break; + } + } + } + + PROCESS_SWITCH(PhiStrangenessCorrelation, processParticleEfficiency, "Process function for Efficiency Computation for Particles of Interest", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx b/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx index cacee6ced5d..2166d02ca28 100644 --- a/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx +++ b/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx @@ -3483,6 +3483,5 @@ struct Phik0shortanalysis { WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; }