diff --git a/PWGCF/Flow/Tasks/flowEventPlane.cxx b/PWGCF/Flow/Tasks/flowEventPlane.cxx index 8287f3df49b..ea895bc8232 100644 --- a/PWGCF/Flow/Tasks/flowEventPlane.cxx +++ b/PWGCF/Flow/Tasks/flowEventPlane.cxx @@ -13,6 +13,8 @@ /// \brief Flow calculation using event plane. /// \author Yash Patley +#include "PWGLF/DataModel/LFStrangenessTables.h" + #include "Common/Core/RecoDecay.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/CollisionAssociationTables.h" @@ -56,16 +58,18 @@ DECLARE_SOA_TABLE(ColSPExt, "AOD", "COLSPEXT", o2::soa::Index<>, colspext::Xc, colspext::Yc); -namespace trackidext +namespace tracksid { +DECLARE_SOA_COLUMN(IsCharged, isCharged, bool); DECLARE_SOA_COLUMN(IsPion, isPion, bool); DECLARE_SOA_COLUMN(IsKaon, isKaon, bool); DECLARE_SOA_COLUMN(IsProton, isProton, bool); -} // namespace trackidext -DECLARE_SOA_TABLE(TrackIdExt, "AOD", "TRACKIDEXT", o2::soa::Index<>, - trackidext::IsPion, - trackidext::IsKaon, - trackidext::IsProton); +} // namespace tracksid +DECLARE_SOA_TABLE(TracksId, "AOD", "TRACKSID", o2::soa::Index<>, + tracksid::IsCharged, + tracksid::IsPion, + tracksid::IsKaon, + tracksid::IsProton); } // namespace o2::aod enum GainClibCorr { @@ -102,9 +106,15 @@ enum ParticleType { kNPart }; +enum V0Type { + kLambda = 0, + kAntiLambda +}; + struct SpectatorPlaneTableProducer { // Table producer Produces colSPExtTable; + Produces tracksIdTable; // Configurables // Collisions @@ -147,6 +157,26 @@ struct SpectatorPlaneTableProducer { Configurable cCcdbUrl{"cCcdbUrl", "http://ccdb-test.cern.ch:8080", "url of ccdb"}; Configurable cCcdbPath{"cCcdbPath", "Users/y/ypatley/DFOO", "Path for ccdb-object"}; + // Tracks + Configurable cTrackMinPt{"cTrackMinPt", 0.1, "p_{T} minimum"}; + Configurable cTrackMaxPt{"cTrackMaxPt", 10.0, "p_{T} maximum"}; + Configurable cNEtaBins{"cNEtaBins", 7, "# of eta bins"}; + Configurable cTrackEtaCut{"cTrackEtaCut", 0.8, "Pseudorapidity cut"}; + Configurable cTrackGlobal{"cTrackGlobal", true, "Global Track"}; + Configurable cTrackDcaXYCut{"cTrackDcaXYCut", 0.1, "DcaXY Cut"}; + Configurable cTrackDcaZCut{"cTrackDcaZCut", 1., "DcaXY Cut"}; + + // Track PID + Configurable cTpcElRejCutMin{"cTpcElRejCutMin", -3., "Electron Rejection Cut Minimum"}; + Configurable cTpcElRejCutMax{"cTpcElRejCutMax", 5., "Electron Rejection Cut Maximum"}; + Configurable cTpcNSigmaCut{"cTpcNSigmaCut", 2, "TPC NSigma Cut"}; + Configurable cTpcRejCut{"cTpcRejCut", 3, "TPC Rej Cut"}; + Configurable cTofNSigmaCut{"cTofNSigmaCut", 2, "TOF NSigma Cut"}; + Configurable cTofRejCut{"cTofRejCut", 3, "TOF Rej Cut"}; + Configurable cPionPtCut{"cPionPtCut", 0.6, "Pion TPC pT cutoff"}; + Configurable cKaonPtCut{"cKaonPtCut", 0.6, "Kaon TPC pT cutoff"}; + Configurable cProtonPtCut{"cProtonPtCut", 1.1, "Proton TPC pT cutoff"}; + // Initialize CCDB Service Service ccdbService; @@ -260,6 +290,9 @@ struct SpectatorPlaneTableProducer { histos.add("Checks/hYaYc", "Y^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent}); histos.add("Checks/hXaYc", "X^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent}); histos.add("Checks/hYaXc", "Y^{A}_{1}X^{C}_{1}", kTProfile, {axisCent}); + + // Directed flow QXY vector + histos.add("DF/hQaQc", "X^{A}_{1}X^{C}_{1} + Y^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent}); } template @@ -568,10 +601,73 @@ struct SpectatorPlaneTableProducer { return true; } + // Track Selection + template + bool selectTrack(T const& track) + { + if (track.pt() <= cTrackMinPt || track.pt() >= cTrackMaxPt || std::abs(track.eta()) >= cTrackEtaCut) { + return false; + } + + if (cTrackGlobal && !track.isGlobalTrackWoDCA()) { + return false; + } + + if (std::abs(track.dcaXY()) >= cTrackDcaXYCut || std::abs(track.dcaZ()) >= cTrackDcaZCut) { + return false; + } + + return true; + } + + template + bool checkTrackPid(float const& ptCut, float const& trackPt, std::vector const& vTpcNsig, std::vector const& vTofNsig, bool const& tofFlag) + { + bool retFlag = false; + if (tofFlag) { + if (vTofNsig[part1] < cTofNSigmaCut && vTofNsig[part2] > cTofRejCut && vTofNsig[part3] > cTofRejCut && vTpcNsig[part1] < cTpcNSigmaCut) { + retFlag = true; + } + } else { + if (trackPt < ptCut && vTpcNsig[part1] < cTpcNSigmaCut && vTpcNsig[part2] > cTpcRejCut && vTpcNsig[part3] > cTpcRejCut) { + retFlag = true; + } + } + return retFlag; + } + + template + bool identifyTrack(T const& track) + { + // Electron rejection + if (std::abs(track.tpcNSigmaPi()) > cTpcRejCut && std::abs(track.tpcNSigmaKa()) > cTpcRejCut && std::abs(track.tpcNSigmaPr()) > cTpcRejCut && track.tpcNSigmaEl() > cTpcElRejCutMin && track.tpcNSigmaEl() < cTpcElRejCutMax) { + return false; + } + + // Pion, Kaon, Proton Identification + std::vector vPtCut = {cPionPtCut, cKaonPtCut, cProtonPtCut}; + std::vector vTpcNsig = {std::abs(track.tpcNSigmaPi()), std::abs(track.tpcNSigmaKa()), std::abs(track.tpcNSigmaPr())}; + std::vector vTofNsig = {std::abs(track.tofNSigmaPi()), std::abs(track.tofNSigmaKa()), std::abs(track.tofNSigmaPr())}; + bool retFlag = false; + + if (partType == kPi && checkTrackPid(vPtCut[kPi], track.pt(), vTpcNsig, vTofNsig, track.hasTOF())) { + retFlag = true; + } else if (partType == kKa && checkTrackPid(vPtCut[kKa], track.pt(), vTpcNsig, vTofNsig, track.hasTOF())) { + retFlag = true; + } else if (partType == kPr && checkTrackPid(vPtCut[kPr], track.pt(), vTpcNsig, vTofNsig, track.hasTOF())) { + retFlag = true; + } else { + return false; + } + + return retFlag; + } + using BCsRun3 = soa::Join; using CollisionsRun3 = soa::Join; + using Tracks = soa::Join; - void process(CollisionsRun3::iterator const& collision, BCsRun3 const&, aod::Zdcs const&) + void processSpectatorPlane(CollisionsRun3::iterator const& collision, BCsRun3 const&, aod::Zdcs const&) { // Analyze collision and get Spectator Plane Vector std::array vSP = {0., 0., 0., 0.}; @@ -593,40 +689,69 @@ struct SpectatorPlaneTableProducer { histos.fill(HIST("Checks/hYaYc"), cent, (vSP[kYa] * vSP[kYc])); histos.fill(HIST("Checks/hXaYc"), cent, (vSP[kXa] * vSP[kYc])); histos.fill(HIST("Checks/hYaXc"), cent, (vSP[kYa] * vSP[kXc])); + + // Directed flow QXY vector + float qac = (vSP[kXa] * vSP[kXc]) + (vSP[kYa] * vSP[kYc]); + histos.fill(HIST("DF/hQaQc"), cent, qac); } // Fill table colSPExtTable(colSPExtFlag, vSP[kXa], vSP[kYa], vSP[kXc], vSP[kYc]); } -}; -struct IdHadronFlow { - // Table producer - Produces trackIdExtTable; + PROCESS_SWITCH(SpectatorPlaneTableProducer, processSpectatorPlane, "Spectator Plane Process", true); + + void processIdHadrons(Tracks const& tracks) + { + for (auto const& track : tracks) { + bool chargedFlag = selectTrack(track); + bool pionFlag = identifyTrack(track); + bool kaonFlag = identifyTrack(track); + bool protonFlag = identifyTrack(track); + tracksIdTable(chargedFlag, pionFlag, kaonFlag, protonFlag); + } + } + + PROCESS_SWITCH(SpectatorPlaneTableProducer, processIdHadrons, "Hadron Id Process", false); +}; +struct FlowEventPlane { // Tracks - Configurable cTrackMinPt{"cTrackMinPt", 0.1, "p_{T} minimum"}; - Configurable cTrackMaxPt{"cTrackMaxPt", 10.0, "p_{T} maximum"}; Configurable cNEtaBins{"cNEtaBins", 7, "# of eta bins"}; - Configurable cTrackEtaCut{"cTrackEtaCut", 0.8, "Pseudorapidity cut"}; - Configurable cTrackGlobal{"cTrackGlobal", true, "Global Track"}; - Configurable cTrackDcaXYCut{"cTrackDcaXYCut", 0.1, "DcaXY Cut"}; - Configurable cTrackDcaZCut{"cTrackDcaZCut", 1., "DcaXY Cut"}; - // Track PID - Configurable cTpcNSigmaCut{"cTpcNSigmaCut", 2, "TPC NSigma Cut"}; - Configurable cTpcRejCut{"cTpcRejCut", 3, "TPC Rej Cut"}; - Configurable cTofNSigmaCut{"cTofNSigmaCut", 2, "TOF NSigma Cut"}; - Configurable cTofRejCut{"cTofRejCut", 3, "TOF Rej Cut"}; - Configurable cPionPtCut{"cPionPtCut", 0.6, "Pion TPC pT cutoff"}; - Configurable cKaonPtCut{"cKaonPtCut", 0.6, "Kaon TPC pT cutoff"}; - Configurable cProtonPtCut{"cProtonPtCut", 1.1, "Proton TPC pT cutoff"}; + // Resonance + Configurable cNRapBins{"cNRapBins", 5, "# of y bins"}; + Configurable cNInvMassBins{"cNInvMassBins", 500, "# of m bins"}; + Configurable cResRapCut{"cResRapCut", 0.5, "Resonance rapidity cut"}; + + // V0 + // Tracks + Configurable cTrackPtCut{"cTrackPtCut", 0.1, "p_{T} minimum"}; + Configurable cTrackEtaCut{"cTrackEtaCut", 0.8, "Pseudorapidity cut"}; + Configurable cTpcNsigmaCut{"cTpcNsigmaCut", 3.0, "TPC NSigma Selection Cut"}; + + // V0s + Configurable cV0TypeSelection{"cV0TypeSelection", 1, "V0 Type Selection"}; + Configurable cMinDcaProtonToPV{"cMinDcaProtonToPV", 0.02, "Minimum Proton DCAr to PV"}; + Configurable cMinDcaPionToPV{"cMinDcaPionToPV", 0.06, "Minimum Pion DCAr to PV"}; + Configurable cDcaV0Dau{"cDcaV0Dau", 1., "DCA between V0 daughters"}; + Configurable cDcaV0ToPv{"cDcaV0ToPv", 0.1, "DCA V0 to PV"}; + Configurable cMinV0Radius{"cMinV0Radius", 0.5, "Minimum V0 radius from PV"}; + Configurable cMaxV0Radius{"cMaxV0Radius", 999.0, "Maximum V0 radius from PV"}; + Configurable cV0CTau{"cV0CTau", 30.0, "Decay length cut"}; + Configurable cV0CosPA{"cV0CosPA", 0.995, "V0 CosPA to PV"}; + Configurable cK0SMassRej{"cK0SMassRej", 0.01, "Reject K0Short Candidates"}; + Configurable cLambdaMassWin{"cLambdaMassWin", 0.007, "Lambda Mass Window"}; + Configurable cMinV0Pt{"cMinV0Pt", 0.5, "Min v0 pT"}; + Configurable cMaxV0Pt{"cMaxV0Pt", 4.5, "Max v0 pT"}; + Configurable cV0RapCut{"cV0RapCut", 0.5, "V0 rap cut"}; // Histogram registry: an object to hold your histograms HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; // Global objects float cent = 0.; + std::array vSP = {0., 0., 0., 0.}; void init(InitContext const&) { @@ -643,232 +768,220 @@ struct IdHadronFlow { const AxisSpec axisTrackdEdx{360, 20, 200, "#frac{dE}{dx}"}; const AxisSpec axisTrackNSigma{161, -4.025, 4.025, {"n#sigma"}}; - // Create histograms - // Track QA - histos.add("TrackQA/hPtDcaXY", "DCA_{XY} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaXY}); - histos.add("TrackQA/hPtDcaZ", "DCA_{Z} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaZ}); - histos.add("TrackQA/hTrackTPCdEdX", "hTrackTPCdEdX", kTH2F, {axisTrackPt, axisTrackdEdx}); - - // Charged particle directed flow - histos.add("DF/hQaQc", "X^{A}_{1}X^{C}_{1} + Y^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent}); - histos.add("DF/hAQu", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta}); - histos.add("DF/hCQu", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta}); - histos.add("DF/hAQuPos", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta}); - histos.add("DF/hCQuPos", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta}); - histos.add("DF/hAQuNeg", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta}); - histos.add("DF/hCQuNeg", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta}); - - // Identified particle - histos.add("PartId/Pion/hdEdX", "PartId/Pion/hdEdX", kTH2F, {axisTrackPt, axisTrackdEdx}); - histos.add("PartId/Pion/hTPCNSigma", "PartId/Pion/hTPCNSigma", kTH2F, {axisTrackPt, axisTrackNSigma}); - histos.add("PartId/Pion/hTOFNSigma", "PartId/Pion/hTOFNSigma", kTH2F, {axisTrackPt, axisTrackNSigma}); - histos.add("PartId/Pion/hAQuPos", "PartId/Pion/hAQuPos", kTProfile2D, {axisCent, axisTrackEta}); - histos.add("PartId/Pion/hAQuNeg", "PartId/Pion/hAQuNeg", kTProfile2D, {axisCent, axisTrackEta}); - histos.add("PartId/Pion/hCQuPos", "PartId/Pion/hCQuPos", kTProfile2D, {axisCent, axisTrackEta}); - histos.add("PartId/Pion/hCQuNeg", "PartId/Pion/hCQuNeg", kTProfile2D, {axisCent, axisTrackEta}); - histos.addClone("PartId/Pion/", "PartId/Kaon/"); - histos.addClone("PartId/Pion/", "PartId/Proton/"); - } - - // Track Selection - template - bool selectTrack(T const& track) - { - if (track.pt() <= cTrackMinPt || track.pt() >= cTrackMaxPt || std::abs(track.eta()) >= cTrackEtaCut) { - return false; - } + const AxisSpec axisTrackRap{cNRapBins, -0.5, 0.5, "y"}; + const AxisSpec axisInvMass{cNInvMassBins, 0.87, 1.12, "M_{KK} (GeV/#it{c}^{2}"}; + const AxisSpec axisMomPID(80, 0, 4, "p_{T} (GeV/#it{c})"); + const AxisSpec axisNsigma(401, -10.025, 10.025, {"n#sigma"}); + const AxisSpec axisdEdx(360, 20, 200, "#frac{dE}{dx}"); + const AxisSpec axisTrackTofSignal(240, 0, 1.2, "#beta"); + const AxisSpec axisRadius(2000, 0, 200, "r(cm)"); + const AxisSpec axisCosPA(300, 0.97, 1.0, "cos(#theta_{PA})"); + const AxisSpec axisDcaV0PV(1000, 0., 10., "dca (cm)"); + const AxisSpec axisDcaProngPV(5000, -50., 50., "dca (cm)"); + const AxisSpec axisDcaDau(75, 0., 1.5, "Daug DCA (#sigma)"); + const AxisSpec axisCTau(2000, 0, 200, "c#tau (cm)"); + const AxisSpec axisAlpha(40, -1, 1, "#alpha"); + const AxisSpec axisQtarm(40, 0, 0.4, "q_{T}"); + const AxisSpec axisLambdaPt(50, 0, 10, "p_{T} (GeV/#it{c})"); + const AxisSpec axisLambdaInvMass(140, 1.08, 1.15, "M_{p#pi} (GeV/#it{c}^{2})"); - if (cTrackGlobal && !track.isGlobalTrackWoDCA()) { - return false; + // Create histograms + // Charged particles + if (doprocessChargedFlow) { + histos.add("TrackQA/hPtDcaXY", "DCA_{XY} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaXY}); + histos.add("TrackQA/hPtDcaZ", "DCA_{Z} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaZ}); + histos.add("TrackQA/hTrackTPCdEdX", "hTrackTPCdEdX", kTH2F, {axisMomPID, axisdEdx}); + histos.add("DF/hQaQc", "X^{A}_{1}X^{C}_{1} + Y^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent}); + histos.add("DF/hAQu", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta}); + histos.add("DF/hCQu", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta}); + histos.add("DF/hAQuPos", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta}); + histos.add("DF/hCQuPos", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta}); + histos.add("DF/hAQuNeg", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta}); + histos.add("DF/hCQuNeg", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta}); } - if (std::abs(track.dcaXY()) > cTrackDcaXYCut || std::abs(track.dcaZ()) > cTrackDcaZCut) { - return false; + // Identified hadrons + if (doprocessIdFlow) { + histos.add("PartId/Pion/hdEdX", "dE/dx vs pT", kTH2F, {axisMomPID, axisTrackdEdx}); + histos.add("PartId/Pion/hTOFSignal", "#beta_{TOF} vs p_{T}", kTH2F, {axisMomPID, axisTrackTofSignal}); + histos.add("PartId/Pion/hTPCNSigma", "n#sigma_{TPC} vs p_{T}", kTH2F, {axisMomPID, axisTrackNSigma}); + histos.add("PartId/Pion/hTOFNSigma", "n#sigma_{TOF} vs p_{T}", kTH2F, {axisMomPID, axisTrackNSigma}); + histos.add("PartId/Pion/hAQuPos", "PartId/Pion/hAQuPos", kTProfile2D, {axisCent, axisTrackEta}); + histos.add("PartId/Pion/hAQuNeg", "PartId/Pion/hAQuNeg", kTProfile2D, {axisCent, axisTrackEta}); + histos.add("PartId/Pion/hCQuPos", "PartId/Pion/hCQuPos", kTProfile2D, {axisCent, axisTrackEta}); + histos.add("PartId/Pion/hCQuNeg", "PartId/Pion/hCQuNeg", kTProfile2D, {axisCent, axisTrackEta}); + histos.addClone("PartId/Pion/", "PartId/Kaon/"); + histos.addClone("PartId/Pion/", "PartId/Proton/"); } - return true; - } - - template - bool checkTrackPid(float const& ptCut, float const& trackPt, std::vector const& vTpcNsig, std::vector const& vTofNsig, bool const& tofFlag) - { - bool retFlag = false; - if (tofFlag) { - if (vTofNsig[part1] < cTofNSigmaCut && vTofNsig[part2] > cTofRejCut && vTofNsig[part3] > cTofRejCut && vTpcNsig[part1] < cTpcNSigmaCut) { - retFlag = true; - } - } else { - if (trackPt < ptCut && vTpcNsig[part1] < cTpcNSigmaCut && vTpcNsig[part2] > cTpcRejCut && vTpcNsig[part3] > cTpcRejCut) { - retFlag = true; - } - } - return retFlag; - } - - template - bool identifyTrack(T const& track) - { - std::vector vPtCut = {cPionPtCut, cKaonPtCut, cProtonPtCut}; - std::vector vTpcNsig = {std::abs(track.tpcNSigmaPi()), std::abs(track.tpcNSigmaKa()), std::abs(track.tpcNSigmaPr())}; - std::vector vTofNsig = {std::abs(track.tofNSigmaPi()), std::abs(track.tofNSigmaKa()), std::abs(track.tofNSigmaPr())}; - bool retFlag = false; - - if (partType == kPi && checkTrackPid(vPtCut[kPi], track.pt(), vTpcNsig, vTofNsig, track.hasTOF())) { - retFlag = true; - } else if (partType == kKa && checkTrackPid(vPtCut[kKa], track.pt(), vTpcNsig, vTofNsig, track.hasTOF())) { - retFlag = true; - } else if (partType == kPr && checkTrackPid(vPtCut[kPr], track.pt(), vTpcNsig, vTofNsig, track.hasTOF())) { - retFlag = true; - } else { - return false; + // Resonance + if (doprocessResoFlow) { + histos.add("Reso/Phi/hSigCentPtInvMass", "hUSCentPtInvMass", kTH3F, {axisCent, axisTrackPt, axisInvMass}); + histos.add("Reso/Phi/hBkgCentPtInvMass", "hLSCentPtInvMass", kTH3F, {axisCent, axisTrackPt, axisInvMass}); + histos.add("Reso/Phi/Sig/hPhiQuA", "hPhiQuA", kTProfile3D, {axisCent, axisTrackRap, axisInvMass}); + histos.add("Reso/Phi/Sig/hPhiQuC", "hPhiQuC", kTProfile3D, {axisCent, axisTrackRap, axisInvMass}); + histos.add("Reso/Phi/Bkg/hPhiQuA", "hPhiQuA", kTProfile3D, {axisCent, axisTrackRap, axisInvMass}); + histos.add("Reso/Phi/Bkg/hPhiQuC", "hPhiQuC", kTProfile3D, {axisCent, axisTrackRap, axisInvMass}); } - return retFlag; - } - - template - void getIdHadronFlow(float const& cent, T const& track, float const& v1a, float const& v1c) - { - static constexpr std::string_view SubDir[] = {"Pion/", "Kaon/", "Proton/"}; - float tpcNsigma = 0., tofNsigma = 0.; - if (part == kPi) { - tpcNsigma = track.tpcNSigmaPi(); - tofNsigma = track.tofNSigmaPi(); - } else if (part == kKa) { - tpcNsigma = track.tpcNSigmaKa(); - tofNsigma = track.tofNSigmaKa(); - } else if (part == kPr) { - tpcNsigma = track.tpcNSigmaPr(); - tofNsigma = track.tofNSigmaPr(); - } else { - return; - } - histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hdEdX"), track.pt(), track.tpcSignal()); - histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hTPCNSigma"), track.pt(), tpcNsigma); - if (track.hasTOF()) { - histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hTOFNSigma"), track.pt(), tofNsigma); - } - if (track.sign() > 0) { - histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hAQuPos"), cent, track.eta(), v1a); - histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hCQuPos"), cent, track.eta(), v1c); - } else { - histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hAQuNeg"), cent, track.eta(), v1a); - histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hCQuNeg"), cent, track.eta(), v1c); + // Lambda + if (doprocessLambdaFlow) { + histos.add("Lambda/QA/hQtVsAlpha", "Armentros-Podolanski Plot", kTH2F, {axisAlpha, axisQtarm}); + histos.add("Lambda/QA/hDcaV0Dau", "DCA between V0 daughters", kTH1F, {axisDcaDau}); + histos.add("Lambda/QA/hDcaPosToPv", "DCA positive prong to PV", kTH1F, {axisDcaProngPV}); + histos.add("Lambda/QA/hDcaNegToPv", "DCA negative prong to PV", kTH1F, {axisDcaProngPV}); + histos.add("Lambda/QA/hDcaV0ToPv", "DCA V0 to PV", kTH1F, {axisDcaV0PV}); + histos.add("Lambda/QA/hCosPa", "cos(#theta_{PA})", kTH1F, {axisCosPA}); + histos.add("Lambda/QA/hRxy", "V_{0} Decay Radius in XY plane", kTH1F, {axisRadius}); + histos.add("Lambda/QA/hCTau", "V_{0} c#tau", kTH1F, {axisCTau}); + histos.add("Lambda/QA/hPosdEdXVsP", "TPC Signal Pos-Prong", kTH2F, {axisMomPID, axisdEdx}); + histos.add("Lambda/QA/hNegdEdXVsP", "TPC Signal Neg-Prong", kTH2F, {axisMomPID, axisdEdx}); + histos.add("Lambda/QA/hPosNsigPrVsP", "TPC n#sigma Pos Prong", kTH2F, {axisMomPID, axisNsigma}); + histos.add("Lambda/QA/hNegNsigPrVsP", "TPC n#sigma Neg Prong", kTH2F, {axisMomPID, axisNsigma}); + histos.add("Lambda/QA/hPosNsigPiVsP", "TPC n#sigma Pos Prong", kTH2F, {axisMomPID, axisNsigma}); + histos.add("Lambda/QA/hNegNsigPiVsP", "TPC n#sigma Neg Prong", kTH2F, {axisMomPID, axisNsigma}); + histos.add("Lambda/hInvMassVsPt", "hInvMassVsPt", kTH3F, {axisCent, axisLambdaInvMass, axisLambdaPt}); + histos.add("Lambda/Flow/hQuA", "hPhiQuA", kTProfile3D, {axisCent, axisTrackRap, axisLambdaInvMass}); + histos.add("Lambda/Flow/hQuC", "hPhiQuC", kTProfile3D, {axisCent, axisTrackRap, axisLambdaInvMass}); + histos.addClone("Lambda/", "AntiLambda/"); } } - template - void fillTrackHist(T const& track) - { - histos.fill(HIST("TrackQA/hPtDcaZ"), track.pt(), track.dcaZ()); - histos.fill(HIST("TrackQA/hPtDcaXY"), track.pt(), track.dcaXY()); - } - - using CollisionsRun3 = soa::Join; - using Tracks = soa::Join; - - void processDummy(CollisionsRun3::iterator const&) {} - - PROCESS_SWITCH(IdHadronFlow, processDummy, "Dummy process", true); - - void processIdHadronFlow(CollisionsRun3::iterator const& collision, Tracks const& tracks) + template + bool selCollision(C const& collision, std::array& v) { // Check collision if (!collision.selColFlag()) { - return; + return false; } // Set centrality cent = collision.centFT0C(); // Flow vectors - std::array vSP = {0., 0., 0., 0.}; - vSP[kXa] = collision.xa(); - vSP[kYa] = collision.ya(); - vSP[kXc] = collision.xc(); - vSP[kYc] = collision.yc(); + v[kXa] = collision.xa(); + v[kYa] = collision.ya(); + v[kXc] = collision.xc(); + v[kYc] = collision.yc(); - // Directed flow QXY vector - float qac = (vSP[kXa] * vSP[kXc]) + (vSP[kYa] * vSP[kYc]); - histos.fill(HIST("DF/hQaQc"), cent, qac); + return true; + } - // Loop over tracks + template + void getIdFlow(T const& tracks) + { float ux = 0., uy = 0., v1a = 0., v1c = 0.; + float tpcNsigma = 0., tofNsigma = 0.; for (auto const& track : tracks) { - // Select track - if (!selectTrack(track)) { - trackIdExtTable(false, false, false); - continue; + static constexpr std::string_view SubDir[] = {"Pion/", "Kaon/", "Proton/"}; + if (part == kPi) { + tpcNsigma = track.tpcNSigmaPi(); + tofNsigma = track.tofNSigmaPi(); + } else if (part == kKa) { + tpcNsigma = track.tpcNSigmaKa(); + tofNsigma = track.tofNSigmaKa(); + } else if (part == kPr) { + tpcNsigma = track.tpcNSigmaPr(); + tofNsigma = track.tofNSigmaPr(); + } else { + return; + } + histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hdEdX"), track.pt(), track.tpcSignal()); + histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hTPCNSigma"), track.pt(), tpcNsigma); + if (track.hasTOF()) { + histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hTOFSignal"), track.pt(), track.beta()); + histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hTOFNSigma"), track.pt(), tofNsigma); } - // Fill track QA - fillTrackHist(track); - - // Get directed flow ux = std::cos(track.phi()); uy = std::sin(track.phi()); v1a = ux * vSP[kXa] + uy * vSP[kYa]; v1c = ux * vSP[kXc] + uy * vSP[kYc]; - // Charged particle directed flow - histos.fill(HIST("DF/hAQu"), cent, track.eta(), v1a); - histos.fill(HIST("DF/hCQu"), cent, track.eta(), v1c); if (track.sign() > 0) { - histos.fill(HIST("DF/hAQuPos"), cent, track.eta(), v1a); - histos.fill(HIST("DF/hCQuPos"), cent, track.eta(), v1c); + histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hAQuPos"), cent, track.eta(), v1a); + histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hCQuPos"), cent, track.eta(), v1c); } else { - histos.fill(HIST("DF/hAQuNeg"), cent, track.eta(), v1a); - histos.fill(HIST("DF/hCQuNeg"), cent, track.eta(), v1c); - } - - // Identified directed flow - if (identifyTrack(track)) { - trackIdExtTable(true, false, false); - getIdHadronFlow(cent, track, v1a, v1c); - } else if (identifyTrack(track)) { - trackIdExtTable(false, true, false); - getIdHadronFlow(cent, track, v1a, v1c); - } else if (identifyTrack(track)) { - trackIdExtTable(false, false, true); - getIdHadronFlow(cent, track, v1a, v1c); + histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hAQuNeg"), cent, track.eta(), v1a); + histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hCQuNeg"), cent, track.eta(), v1c); } } } - PROCESS_SWITCH(IdHadronFlow, processIdHadronFlow, "Identified hadron flow process", false); -}; -struct FlowEventPlane { - // Resonance - Configurable cNRapBins{"cNRapBins", 5, "# of y bins"}; - Configurable cNInvMassBins{"cNInvMassBins", 500, "# of m bins"}; - Configurable cResRapCut{"cResRapCut", 0.5, "Resonance rapidity cut"}; + template + bool selV0DauTracks(V const& v0, T const& postrack, T const& negtrack) + { + // Kinematic selection + if (postrack.pt() <= cTrackPtCut || negtrack.pt() <= cTrackPtCut || std::abs(postrack.eta()) >= cTrackEtaCut || std::abs(negtrack.eta()) >= cTrackEtaCut) { + return false; + } - // Histogram registry: an object to hold your histograms - HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + // Apply DCA Selection on Daughter Tracks Based on Lambda/AntiLambda daughters + float dcaProton = 0., dcaPion = 0.; + if (part == kLambda) { + dcaProton = std::abs(v0.dcapostopv()); + dcaPion = std::abs(v0.dcanegtopv()); + } else if (part == kAntiLambda) { + dcaPion = std::abs(v0.dcapostopv()); + dcaProton = std::abs(v0.dcanegtopv()); + } else { + return false; + } - // Global objects - float cent = 0.; + if (dcaProton <= cMinDcaProtonToPV || dcaPion <= cMinDcaPionToPV) { + return false; + } - void init(InitContext const&) - { - // Define axes - const AxisSpec axisCent{100, 0., 100, "FT0C%"}; + // Daughter track PID + float tpcNSigmaPr = 0., tpcNSigmaPi = 0.; + + switch (part) { + // postrack = Proton, negtrack = Pion + case kLambda: + tpcNSigmaPr = postrack.tpcNSigmaPr(); + tpcNSigmaPi = negtrack.tpcNSigmaPi(); + break; + + // negtrack = Proton, postrack = Pion + case kAntiLambda: + tpcNSigmaPr = negtrack.tpcNSigmaPr(); + tpcNSigmaPi = postrack.tpcNSigmaPi(); + break; + } - const AxisSpec axisXYac{600, -6, 6, "Q^{t}Q^{p}"}; - const AxisSpec axisV1{400, -4, 4, "v_{1}"}; + if (std::abs(tpcNSigmaPr) >= cTpcNsigmaCut || std::abs(tpcNSigmaPi) >= cTpcNsigmaCut) { + return false; + } - const AxisSpec axisTrackPt{100, 0., 10., "p_{T} (GeV/#it{c})"}; - const AxisSpec axisTrackRap{cNRapBins, -0.5, 0.5, "y"}; - const AxisSpec axisInvMass{cNInvMassBins, 0.87, 1.12, "M_{KK} (GeV/#it{c}^{2}"}; + return true; + } - // Create histograms - // Resonance - histos.add("Reso/Phi/hSigCentPtInvMass", "hUSCentPtInvMass", kTH3F, {axisCent, axisTrackPt, axisInvMass}); - histos.add("Reso/Phi/hBkgCentPtInvMass", "hLSCentPtInvMass", kTH3F, {axisCent, axisTrackPt, axisInvMass}); - histos.add("Reso/Phi/Sig/hPhiQuA", "hPhiQuA", kTProfile3D, {axisCent, axisTrackRap, axisInvMass}); - histos.add("Reso/Phi/Sig/hPhiQuC", "hPhiQuC", kTProfile3D, {axisCent, axisTrackRap, axisInvMass}); - histos.add("Reso/Phi/Bkg/hPhiQuA", "hPhiQuA", kTProfile3D, {axisCent, axisTrackRap, axisInvMass}); - histos.add("Reso/Phi/Bkg/hPhiQuC", "hPhiQuC", kTProfile3D, {axisCent, axisTrackRap, axisInvMass}); + template + void fillV0QAHist(C const& col, V const& v0, T const&) + { + static constexpr std::string_view SubDir[] = {"Lambda/QA/", "AntiLambda/QA/"}; + + // daugthers + auto postrack = v0.template posTrack_as(); + auto negtrack = v0.template negTrack_as(); + + // ctau + float ctau = v0.distovertotmom(col.posX(), col.posY(), col.posZ()) * MassLambda0; + + histos.fill(HIST(SubDir[part]) + HIST("hQtVsAlpha"), v0.alpha(), v0.qtarm()); + histos.fill(HIST(SubDir[part]) + HIST("hDcaV0Dau"), v0.dcaV0daughters()); + histos.fill(HIST(SubDir[part]) + HIST("hDcaPosToPv"), v0.dcapostopv()); + histos.fill(HIST(SubDir[part]) + HIST("hDcaNegToPv"), v0.dcanegtopv()); + histos.fill(HIST(SubDir[part]) + HIST("hDcaV0ToPv"), v0.dcav0topv()); + histos.fill(HIST(SubDir[part]) + HIST("hCosPa"), v0.v0cosPA()); + histos.fill(HIST(SubDir[part]) + HIST("hRxy"), v0.v0radius()); + histos.fill(HIST(SubDir[part]) + HIST("hCTau"), ctau); + histos.fill(HIST(SubDir[part]) + HIST("hPosdEdXVsP"), postrack.tpcInnerParam(), postrack.tpcSignal()); + histos.fill(HIST(SubDir[part]) + HIST("hNegdEdXVsP"), negtrack.tpcInnerParam(), negtrack.tpcSignal()); + histos.fill(HIST(SubDir[part]) + HIST("hPosNsigPrVsP"), postrack.tpcInnerParam(), postrack.tpcNSigmaPr()); + histos.fill(HIST(SubDir[part]) + HIST("hNegNsigPrVsP"), negtrack.tpcInnerParam(), negtrack.tpcNSigmaPr()); + histos.fill(HIST(SubDir[part]) + HIST("hPosNsigPiVsP"), postrack.tpcInnerParam(), postrack.tpcNSigmaPi()); + histos.fill(HIST(SubDir[part]) + HIST("hNegNsigPiVsP"), negtrack.tpcInnerParam(), negtrack.tpcNSigmaPi()); } template @@ -927,31 +1040,77 @@ struct FlowEventPlane { } using CollisionsRun3 = soa::Join; - using Tracks = soa::Join; + using Tracks = soa::Join; + using TracksV0s = soa::Join; + // Partitions SliceCache cache; - Partition kaonTrackPartition = (aod::trackidext::isKaon == true); + Partition chargedTrackPartition = (aod::tracksid::isCharged == true); + Partition pionTrackPartition = (aod::tracksid::isCharged == true) && (aod::tracksid::isPion == true); + Partition kaonTrackPartition = (aod::tracksid::isCharged == true) && (aod::tracksid::isKaon == true); + Partition protonTrackPartition = (aod::tracksid::isCharged == true) && (aod::tracksid::isProton == true); void processDummy(CollisionsRun3::iterator const&) {} PROCESS_SWITCH(FlowEventPlane, processDummy, "Dummy process", true); - void processResoFlow(CollisionsRun3::iterator const& collision, Tracks const&) + void processChargedFlow(CollisionsRun3::iterator const& collision, Tracks const&) { // Check collision - if (!collision.selColFlag()) { + if (!selCollision(collision, vSP)) { return; } + // Loop over tracks + auto chargedTracks = chargedTrackPartition->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + float ux = 0., uy = 0., v1a = 0., v1c = 0.; + for (auto const& track : chargedTracks) { + // Fill track QA + histos.fill(HIST("TrackQA/hPtDcaZ"), track.pt(), track.dcaZ()); + histos.fill(HIST("TrackQA/hPtDcaXY"), track.pt(), track.dcaXY()); + histos.fill(HIST("TrackQA/hTrackTPCdEdX"), track.pt(), track.tpcSignal()); - // Set centrality - cent = collision.centFT0C(); + // Get directed flow + ux = std::cos(track.phi()); + uy = std::sin(track.phi()); + v1a = ux * vSP[kXa] + uy * vSP[kYa]; + v1c = ux * vSP[kXc] + uy * vSP[kYc]; - // Flow vectors - std::array vSP = {0., 0., 0., 0.}; - vSP[kXa] = collision.xa(); - vSP[kYa] = collision.ya(); - vSP[kXc] = collision.xc(); - vSP[kYc] = collision.yc(); + // Charged particle directed flow + histos.fill(HIST("DF/hAQu"), cent, track.eta(), v1a); + histos.fill(HIST("DF/hCQu"), cent, track.eta(), v1c); + if (track.sign() > 0) { + histos.fill(HIST("DF/hAQuPos"), cent, track.eta(), v1a); + histos.fill(HIST("DF/hCQuPos"), cent, track.eta(), v1c); + } else { + histos.fill(HIST("DF/hAQuNeg"), cent, track.eta(), v1a); + histos.fill(HIST("DF/hCQuNeg"), cent, track.eta(), v1c); + } + } + } + + PROCESS_SWITCH(FlowEventPlane, processChargedFlow, "Charged particle flow process", false); + + void processIdFlow(CollisionsRun3::iterator const& collision, Tracks const&) + { + if (!selCollision(collision, vSP)) { + return; + } + // Loop over tracks + auto pionTracks = pionTrackPartition->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + auto kaonTracks = kaonTrackPartition->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + auto protonTracks = protonTrackPartition->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + getIdFlow(pionTracks); + getIdFlow(kaonTracks); + getIdFlow(protonTracks); + } + + PROCESS_SWITCH(FlowEventPlane, processIdFlow, "Identified particle flow process", false); + + void processResoFlow(CollisionsRun3::iterator const& collision, Tracks const&) + { + if (!selCollision(collision, vSP)) { + return; + } // Track partitions auto kaonTracks = kaonTrackPartition->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); @@ -960,12 +1119,83 @@ struct FlowEventPlane { getResoFlow(kaonTracks, kaonTracks, vSP); } PROCESS_SWITCH(FlowEventPlane, processResoFlow, "Resonance flow process", false); + + void processLambdaFlow(CollisionsRun3::iterator const& collision, aod::V0Datas const& V0s, TracksV0s const& tracks) + { + if (!selCollision(collision, vSP)) { + return; + } + + // Loop over v0s + for (auto const& v0 : V0s) { + // V0 kinematic selection + if (v0.pt() <= cMinV0Pt || v0.pt() >= cMaxV0Pt || std::abs(v0.yLambda()) >= cV0RapCut) { + continue; + } + + // Topological selections + float ctau = v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * MassLambda0; + if (v0.dcaV0daughters() >= cDcaV0Dau || v0.dcav0topv() >= cDcaV0ToPv || + v0.v0radius() <= cMinV0Radius || v0.v0radius() >= cMaxV0Radius || + v0.v0cosPA() <= cV0CosPA || ctau >= cV0CTau || v0.v0Type() != cV0TypeSelection) { + continue; + } + + // Ks Mass Rejection + if (std::abs(v0.mK0Short() - MassK0Short) <= cK0SMassRej) { + continue; + } + + // Initialize daughter tracks + auto postrack = v0.template posTrack_as(); + auto negtrack = v0.template negTrack_as(); + + // Initialize selection flags + bool lambdaFlag = false, antiLambdaFlag = false; + + // Get v0 track as lambda + if ((std::abs(v0.mLambda() - MassLambda0) < cLambdaMassWin) && (selV0DauTracks(v0, postrack, negtrack))) { + lambdaFlag = true; + } + + // Get v0 track as anti-lambda + if ((std::abs(v0.mAntiLambda() - MassLambda0) < cLambdaMassWin) && (selV0DauTracks(v0, postrack, negtrack))) { + antiLambdaFlag = true; + } + + // Lambda/Anti-Lambda selection checks + if (!lambdaFlag && !antiLambdaFlag) { // neither Lambda nor Anti-Lambda + continue; + } else if (lambdaFlag && antiLambdaFlag) { // check if the track is identified as lambda and anti-lambda both (DISCARD THIS TRACK) + continue; + } + + // We have a Lambda/Anti-Lambda + // Directed flow + float ux = std::cos(v0.phi()); + float uy = std::sin(v0.phi()); + float v1a = ux * vSP[kXa] + uy * vSP[kYa]; + float v1c = ux * vSP[kXc] + uy * vSP[kYc]; + + if (lambdaFlag) { + fillV0QAHist(collision, v0, tracks); + histos.fill(HIST("Lambda/hInvMassVsPt"), cent, v0.mLambda(), v0.pt()); + histos.fill(HIST("Lambda/Flow/hQuA"), cent, v0.yLambda(), v0.mLambda(), v1a); + histos.fill(HIST("Lambda/Flow/hQuC"), cent, v0.yLambda(), v0.mLambda(), v1c); + } else if (antiLambdaFlag) { + fillV0QAHist(collision, v0, tracks); + histos.fill(HIST("AntiLambda/hInvMassVsPt"), cent, v0.mAntiLambda(), v0.pt()); + histos.fill(HIST("AntiLambda/Flow/hQuA"), cent, v0.yLambda(), v0.mAntiLambda(), v1a); + histos.fill(HIST("AntiLambda/Flow/hQuC"), cent, v0.yLambda(), v0.mAntiLambda(), v1c); + } + } + } + PROCESS_SWITCH(FlowEventPlane, processLambdaFlow, "Lambda flow process", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc)}; }