diff --git a/PWGCF/Flow/Tasks/CMakeLists.txt b/PWGCF/Flow/Tasks/CMakeLists.txt index a7cfbe3c95f..0f5f97f64fb 100644 --- a/PWGCF/Flow/Tasks/CMakeLists.txt +++ b/PWGCF/Flow/Tasks/CMakeLists.txt @@ -113,3 +113,9 @@ o2physics_add_dpl_workflow(pid-flow-pt-corr SOURCES pidFlowPtCorr.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::GFWCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(flow-gpol-ese + SOURCES flowDirectedFlowTask.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + diff --git a/PWGCF/Flow/Tasks/flowDirectedFlowTask.cxx b/PWGCF/Flow/Tasks/flowDirectedFlowTask.cxx new file mode 100644 index 00000000000..28f32417010 --- /dev/null +++ b/PWGCF/Flow/Tasks/flowDirectedFlowTask.cxx @@ -0,0 +1,594 @@ +// 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. + +/// \author Junlee Kim (jikim1290@gmail.com) +/// \file flowDirectedFlowTask.cxx +/// \brief Task for q1-dependent directed flow and global polarization +/// \since May 2026 + +#include "PWGLF/DataModel/LFStrangenessPIDTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/SPCalibrationTables.h" + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/CCDB/RCTSelectionFlags.h" +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.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 +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::aod::rctsel; +using namespace o2::constants::physics; + +struct flowDirectedFlowTask { + Configurable cfgUrl{"cfgUrl", "http://alice-ccdb.cern.ch", "CCDB URL"}; + Configurable nolaterthan{"nolaterthan", -1, "Latest acceptable CCDB timestamp; -1 = current"}; + + Configurable cfgPVSel{"cfgPVSel", true, "flag for the primary vertex selection"}; + Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; + Configurable cfgAddEvtSelPileup{"cfgAddEvtSelPileup", true, "flag for the pilup rejection"}; + Configurable cfgMaxOccupancy{"cfgMaxOccupancy", 999999, "maximum occupancy of tracks in neighbouring collisions in a given time range"}; + Configurable cfgMinOccupancy{"cfgMinOccupancy", 0, "maximum occupancy of tracks in neighbouring collisions in a given time range"}; + + Configurable cfgq1mode{"cfgq1mode", 0, "0: |QC-QA|, 1: |QA|, 2: |QC|, 3: 0.5*(|QA|+|QC|)"}; + + Configurable cfgUseAccCorr{"cfgUseAccCorr", true, "Use acceptance correction profile"}; + Configurable cfgAccPathL{"cfgAccPathL", "Users/p/prottay/My/Object/acccorrL", "CCDB path to Lambda acceptance correction"}; + Configurable cfgAccPathAL{"cfgAccPathAL", "Users/p/prottay/My/Object/acccorrAL", "CCDB path to AntiLambda acceptance correction"}; + + Configurable cfgV0PtMin{"cfgV0PtMin", 0.f, "Minimum V0 pT"}; + Configurable cfgV0Rap{"cfgV0Rap", 0.8f, "V0 rapidity range"}; + Configurable cfgV0DCADaughMax{"cfgV0DCADaughMax", 0.2, "Maximum DCA between V0 daughters"}; + Configurable cfgV0CPAMin{"cfgV0CPAMin", 0.9998, "Minimum V0 CPA"}; + Configurable cfgV0TranRadV0Min{"cfgV0TranRadV0Min", 1.5f, "Minimum V0 transverse radius"}; + Configurable cfgV0TranRadV0Max{"cfgV0TranRadV0Max", 100.f, "Maximum V0 transverse radius"}; + Configurable cfgMaxV0DCA{"cfgMaxV0DCA", 1.2, "Maximum V0 DCA to PV"}; + Configurable cfgMinV0DCAPr{"cfgMinV0DCAPr", 0.05, "Minimum proton DCA to PV"}; + Configurable cfgMinV0DCAPi{"cfgMinV0DCAPi", 0.05, "Minimum pion DCA to PV"}; + Configurable cfgMaxV0LifeTime{"cfgMaxV0LifeTime", 20, "Maximum V0 lifetime"}; + + Configurable cfgCutPT{"cfgCutPT", 0.15, "PT cut on daughter track"}; + Configurable cfgCutEta{"cfgCutEta", 0.8, "Eta cut on daughter track"}; + + Configurable cfgDaughEta{"cfgDaughEta", 0.8f, "Daughter eta range"}; + Configurable cfgDaughPrPt{"cfgDaughPrPt", 0.4f, "Minimum daughter proton pT"}; + Configurable cfgDaughPiPt{"cfgDaughPiPt", 0.2f, "Minimum daughter pion pT"}; + Configurable cfgDaughTPCnclsMin{"cfgDaughTPCnclsMin", 50.f, "Minimum TPC found clusters"}; + Configurable cfgTPCcluster{"cfgTPCcluster", 70, "Minimum TPC crossed rows"}; + Configurable cfgRcrfc{"cfgRcrfc", 0.8f, "TPC crossed rows / findable clusters"}; + Configurable cfgDaughPIDCuts{"cfgDaughPIDCuts", 3.f, "TPC PID nsigma cut"}; + + Configurable cfgTrackPtMin{"cfgTrackPtMin", 0.2f, "Minimum track pT for v1"}; + Configurable cfgTrackPtMax{"cfgTrackPtMax", 5.0f, "Maximum track pT for v1"}; + Configurable cfgTrackEtaMax{"cfgTrackEtaMax", 0.8f, "Maximum |eta| for v1"}; + Configurable cfgMinTPCClustersTrack{"cfgMinTPCClustersTrack", 70, "Minimum TPC clusters for v1 tracks"}; + + struct : ConfigurableGroup { + Configurable cfgRequireRCTFlagChecker{"cfgRequireRCTFlagChecker", true, "Check event quality in run condition table"}; + Configurable cfgEvtRCTFlagCheckerLabel{"cfgEvtRCTFlagCheckerLabel", "CBT_hadronPID", "Evt sel: RCT flag checker label"}; + Configurable cfgEvtRCTFlagCheckerZDCCheck{"cfgEvtRCTFlagCheckerZDCCheck", true, "Evt sel: RCT flag checker ZDC check"}; + Configurable cfgEvtRCTFlagCheckerLimitAcceptAsBad{"cfgEvtRCTFlagCheckerLimitAcceptAsBad", false, "Evt sel: RCT flag checker treat Limited Acceptance As Bad"}; + } rctCut; + + float massLambda = o2::constants::physics::MassLambda0; + float massPr = o2::constants::physics::MassProton; + float massPi = o2::constants::physics::MassPionCharged; + + ConfigurableAxis centAxis{"centAxis", {VARIABLE_WIDTH, 0., 5., 10., 20., 30., 40., 50., 80.}, "Centrality (%)"}; + ConfigurableAxis ptAxis{"ptAxis", {VARIABLE_WIDTH, 0.2, 0.5, 1., 1.5, 2., 2.5, 3., 4., 5., 6.5, 8., 10.}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis etaAxis{"etaAxis", {VARIABLE_WIDTH, -0.8, -0.4, -0.2, 0., 0.2, 0.4, 0.8}, "#eta"}; + ConfigurableAxis massAxis{"massAxis", {100, 1.0, 1.2}, "#it{M} (GeV/#it{c}^{2})"}; + ConfigurableAxis polAxis{"polAxis", {200, -1.0, 1.0}, "polarization observable"}; + ConfigurableAxis spAxis{"spAxis", {400, -10.0, 10.0}, "SP observable"}; + ConfigurableAxis q1Axis{"q1Axis", {VARIABLE_WIDTH, 0., 0.5, 1., 1.5, 2., 3., 5., 10.}, "q_{1}^{ZDC}"}; + ConfigurableAxis qAxis{"qAxis", {200, -10., 10.}, "Q"}; + ConfigurableAxis resAxis{"resAxis", {200, -1., 1.}, "resolution / correlation"}; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + Service ccdb; + TProfile2D* accprofileL = nullptr; + TProfile2D* accprofileAL = nullptr; + int currentRunNumber = -999; + int lastRunNumber = -999; + + Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; + Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && nabs(aod::track::pt) > cfgCutPT); + + RCTFlagsChecker rctChecker; + + using EventCandidates = soa::Filtered>; + using AllTrackCandidates = soa::Filtered>; + using ResoV0s = aod::V0Datas; + using BCsRun3 = soa::Join; + + void init(InitContext&) + { + ccdb->setURL(cfgUrl.value); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + + rctChecker.init(rctCut.cfgEvtRCTFlagCheckerLabel, rctCut.cfgEvtRCTFlagCheckerZDCCheck, rctCut.cfgEvtRCTFlagCheckerLimitAcceptAsBad); + + histos.add("hCentrality", "Centrality", kTH1F, {{centAxis}}); + histos.add("hQ1ZDCvscent", "q1 vs centrality", kTH2F, {{centAxis}, {q1Axis}}); + histos.add("hQxAvscent", "Qx A vs centrality", kTH2F, {{centAxis}, {qAxis}}); + histos.add("hQyAvscent", "Qy A vs centrality", kTH2F, {{centAxis}, {qAxis}}); + histos.add("hQxCvscent", "Qx C vs centrality", kTH2F, {{centAxis}, {qAxis}}); + histos.add("hQyCvscent", "Qy C vs centrality", kTH2F, {{centAxis}, {qAxis}}); + + histos.add("hpResCosAC", "cos(#Psi_{A}-#Psi_{C}) vs centrality", kTH2F, {{centAxis}, {resAxis}}); + histos.add("hpResDotAC", "Q_{A}#upoint Q_{C} vs centrality", kTH2F, {{centAxis}, {resAxis}}); + histos.add("hpQxAQxC", "QxA QxC", kTH2F, {{centAxis}, {resAxis}}); + histos.add("hpQyAQyC", "QyA QyC", kTH2F, {{centAxis}, {resAxis}}); + histos.add("hpQxAQyC", "QxA QyC", kTH2F, {{centAxis}, {resAxis}}); + histos.add("hpQxCQyA", "QxC QyA", kTH2F, {{centAxis}, {resAxis}}); + + std::vector axesV1SP = {centAxis, ptAxis, etaAxis, q1Axis, spAxis}; + std::vector axesV1EP = {centAxis, ptAxis, etaAxis, q1Axis, polAxis}; + + histos.add("hV1SPFullQ1", "charged v1 SP numerator u#upoint(QC-QA) vs q1", HistType::kTHnSparseF, axesV1SP, true); + histos.add("hV1SPAQ1", "charged v1 SP numerator u#upointQA vs q1", HistType::kTHnSparseF, axesV1SP, true); + histos.add("hV1SPCQ1", "charged v1 SP numerator u#upointQC vs q1", HistType::kTHnSparseF, axesV1SP, true); + histos.add("hV1EPFullQ1", "charged v1 EP QA cos(phi-PsiZDC) vs q1", HistType::kTHnSparseF, axesV1EP, true); + histos.add("hV1EPAQ1", "charged v1 EP QA cos(phi-PsiA) vs q1", HistType::kTHnSparseF, axesV1EP, true); + histos.add("hV1EPCQ1", "charged v1 EP QA cos(phi-PsiC) vs q1", HistType::kTHnSparseF, axesV1EP, true); + + std::vector axesPolSPQ1 = {massAxis, ptAxis, spAxis, centAxis, q1Axis}; + + histos.add("hSparseLambdaPolSPQ1", "Lambda SP polarization numerator vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaPolSPwgtQ1", "Lambda SP polarization numerator / acceptance vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaPolSPAQ1", "Lambda SP polarization numerator A vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaPolSPCQ1", "Lambda SP polarization numerator C vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaPolSPAwgtQ1", "Lambda SP polarization numerator A / acceptance vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaPolSPCwgtQ1", "Lambda SP polarization numerator C / acceptance vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + + histos.add("hSparseAntiLambdaPolSPQ1", "AntiLambda SP polarization numerator vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaPolSPwgtQ1", "AntiLambda SP polarization numerator / acceptance vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaPolSPAQ1", "AntiLambda SP polarization numerator A vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaPolSPCQ1", "AntiLambda SP polarization numerator C vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaPolSPAwgtQ1", "AntiLambda SP polarization numerator A / acceptance vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaPolSPCwgtQ1", "AntiLambda SP polarization numerator C / acceptance vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + + histos.add("hSparseLambdaPolEPQ1", "Lambda EP QA sin(phi*-PsiZDC) vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaPolEPAQ1", "Lambda EP QA sin(phi*-PsiA) vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaPolEPCQ1", "Lambda EP QA sin(phi*-PsiC) vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaPolEPQ1", "AntiLambda EP QA sin(phi*-PsiZDC) vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaPolEPAQ1", "AntiLambda EP QA sin(phi*-PsiA) vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaPolEPCQ1", "AntiLambda EP QA sin(phi*-PsiC) vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + + histos.add("hSparseLambdaCosPsiQ1", "Lambda cos(PsiZDC) vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaSinPsiQ1", "Lambda sin(PsiZDC) vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaCosPsiQ1", "AntiLambda cos(PsiZDC) vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaSinPsiQ1", "AntiLambda sin(PsiZDC) vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + + histos.add("hSparseLambdaCorrSinPhiStarQ1", "Lambda sin(phi*) acceptance correction vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaCorrCosPhiStarQ1", "Lambda cos(phi*) acceptance correction vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaCorrSinThetaStarQ1", "Lambda sin(theta*) acceptance correction vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaAvgUxQ1", "Lambda vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseLambdaAvgUyQ1", "Lambda vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + + histos.add("hSparseAntiLambdaCorrSinPhiStarQ1", "AntiLambda sin(phi*) acceptance correction vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaCorrCosPhiStarQ1", "AntiLambda cos(phi*) acceptance correction vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaCorrSinThetaStarQ1", "AntiLambda sin(theta*) acceptance correction vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaAvgUxQ1", "AntiLambda vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + histos.add("hSparseAntiLambdaAvgUyQ1", "AntiLambda vs q1", HistType::kTHnSparseF, axesPolSPQ1, true); + } + + float getPhiInRange(float phi) + { + while (phi <= -o2::constants::math::PI) { + phi += o2::constants::math::TwoPI; + } + while (phi > o2::constants::math::PI) { + phi -= o2::constants::math::TwoPI; + } + return phi; + } + + template + bool eventSelected(TCollision collision) + { + if (!collision.sel8()) { + return 0; + } + if (!collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { + return 0; + } + if (!collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { + return 0; + } + if (cfgPVSel && std::abs(collision.posZ()) > cfgCutVertex) { + return 0; + } + if (cfgAddEvtSelPileup && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { + return 0; + } + if (collision.trackOccupancyInTimeRange() > cfgMaxOccupancy || collision.trackOccupancyInTimeRange() < cfgMinOccupancy) { + return 0; + } + if (rctCut.cfgRequireRCTFlagChecker && !rctChecker(collision)) { + return 0; + } + + return 1; + } + + template + bool getZDCQ1(Collision const& collision, float& qxA, float& qyA, float& qxC, float& qyC, + float& psiA, float& psiC, float& psiFull, float& q1) + { + const float rawQxA = collision.qxZDCA(); + const float rawQyA = collision.qyZDCA(); + const float rawQxC = collision.qxZDCC(); + const float rawQyC = collision.qyZDCC(); + + psiA = collision.psiZDCA(); + psiC = collision.psiZDCC(); + + const float magA = std::sqrt(rawQxA * rawQxA + rawQyA * rawQyA); + const float magC = std::sqrt(rawQxC * rawQxC + rawQyC * rawQyC); + qxA = magA * std::cos(psiA); + qyA = magA * std::sin(psiA); + qxC = magC * std::cos(psiC); + qyC = magC * std::sin(psiC); + + const float qxFull = qxC - qxA; + const float qyFull = qyC - qyA; + psiFull = std::atan2(qyFull, qxFull); + + const float q1A = std::sqrt(qxA * qxA + qyA * qyA); + const float q1C = std::sqrt(qxC * qxC + qyC * qyC); + const float q1Full = std::sqrt(qxFull * qxFull + qyFull * qyFull); + const float q1Mean = 0.5 * (q1A + q1C); + + if (cfgq1mode == 1) { + q1 = q1A; + } else if (cfgq1mode == 2) { + q1 = q1C; + } else if (cfgq1mode == 3) { + q1 = q1Mean; + } else { + q1 = q1Full; + } + + return std::isfinite(qxA) && std::isfinite(qyA) && std::isfinite(qxC) && std::isfinite(qyC) && std::isfinite(psiA) && std::isfinite(psiC) && std::isfinite(psiFull) && std::isfinite(q1); + } + + template + bool selectV1Track(Track const& track) + { + if (!track.isGlobalTrack()) { + return false; + } + if (track.pt() < cfgTrackPtMin || track.pt() > cfgTrackPtMax) { + return false; + } + if (std::abs(track.eta()) > cfgTrackEtaMax) { + return false; + } + if (track.tpcNClsFound() < cfgMinTPCClustersTrack) { + return false; + } + return true; + } + + template + bool isCompatible(V0 const& v0, int pid) + { + // pid = 0: Lambda, pid = 1: AntiLambda + auto posTrack = v0.template posTrack_as(); + auto negTrack = v0.template negTrack_as(); + + if (pid == 0 && (v0.positivept() < cfgDaughPrPt || v0.negativept() < cfgDaughPiPt)) { + return false; + } + if (pid == 1 && (v0.positivept() < cfgDaughPiPt || v0.negativept() < cfgDaughPrPt)) { + return false; + } + if (std::abs(v0.positiveeta()) > cfgDaughEta || std::abs(v0.negativeeta()) > cfgDaughEta) { + return false; + } + + if (posTrack.tpcNClsCrossedRows() < cfgTPCcluster || negTrack.tpcNClsCrossedRows() < cfgTPCcluster) { + return false; + } + if (posTrack.tpcNClsFound() < cfgDaughTPCnclsMin || negTrack.tpcNClsFound() < cfgDaughTPCnclsMin) { + return false; + } + if (posTrack.tpcCrossedRowsOverFindableCls() < cfgRcrfc || negTrack.tpcCrossedRowsOverFindableCls() < cfgRcrfc) { + return false; + } + + if (pid == 0 && (std::abs(posTrack.tpcNSigmaPr()) > cfgDaughPIDCuts || std::abs(negTrack.tpcNSigmaPi()) > cfgDaughPIDCuts)) { + return false; + } + if (pid == 1 && (std::abs(posTrack.tpcNSigmaPi()) > cfgDaughPIDCuts || std::abs(negTrack.tpcNSigmaPr()) > cfgDaughPIDCuts)) { + return false; + } + + if (pid == 0 && (std::abs(v0.dcapostopv()) < cfgMinV0DCAPr || std::abs(v0.dcanegtopv()) < cfgMinV0DCAPi)) { + return false; + } + if (pid == 1 && (std::abs(v0.dcapostopv()) < cfgMinV0DCAPi || std::abs(v0.dcanegtopv()) < cfgMinV0DCAPr)) { + return false; + } + return true; + } + + template + bool selectV0(Collision const& collision, V0 const& v0) + { + const float ctauLambda = v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * massLambda; + if (std::abs(v0.dcav0topv()) > cfgMaxV0DCA) { + return false; + } + if (v0.pt() < cfgV0PtMin) { + return false; + } + if (std::abs(v0.dcaV0daughters()) > cfgV0DCADaughMax) { + return false; + } + if (v0.v0cosPA() < cfgV0CPAMin) { + return false; + } + if (v0.v0radius() < cfgV0TranRadV0Min || v0.v0radius() > cfgV0TranRadV0Max) { + return false; + } + if (std::abs(ctauLambda) > cfgMaxV0LifeTime) { + return false; + } + if (std::abs(v0.yLambda()) > cfgV0Rap) { + return false; + } + if (std::abs(v0.eta()) > 0.8) { + return false; + } + return true; + } + + float acceptanceWeight(bool isLambda, float eta, float pt) + { + if (!cfgUseAccCorr) { + return 1.0; + } + TProfile2D* prof = isLambda ? accprofileL : accprofileAL; + if (!prof) { + return 1.0; + } + const int binx = prof->GetXaxis()->FindBin(eta); + const int biny = prof->GetYaxis()->FindBin(pt); + const float acc = prof->GetBinContent(binx, biny); + if (acc <= 0.0 || !std::isfinite(acc)) { + return 1.0; + } + return 1.0 / acc; + } + + template + void fillPolarization(bool isLambda, V0 const& v0, float centrality, float q1, float psiA, float psiC, float psiFull, float qxA, float qyA, float qxC, float qyC) + { + ROOT::Math::PxPyPzMVector daughter; + ROOT::Math::PxPyPzMVector pion; + float mass = 0.0; + + if (isLambda) { + daughter = ROOT::Math::PxPyPzMVector(v0.pxpos(), v0.pypos(), v0.pzpos(), massPr); + pion = ROOT::Math::PxPyPzMVector(v0.pxneg(), v0.pyneg(), v0.pzneg(), massPi); + mass = v0.mLambda(); + } else { + daughter = ROOT::Math::PxPyPzMVector(v0.pxneg(), v0.pyneg(), v0.pzneg(), massPr); + pion = ROOT::Math::PxPyPzMVector(v0.pxpos(), v0.pypos(), v0.pzpos(), massPi); + mass = v0.mAntiLambda(); + } + + const auto hyperon = daughter + pion; + ROOT::Math::Boost boostToRest(hyperon.BoostToCM()); + const auto daughterStar = boostToRest(daughter); + + float phiStar = std::atan2(daughterStar.Py(), daughterStar.Px()); + const float ux = std::cos(phiStar); + const float uy = std::sin(phiStar); + + const float cosThetaStar = daughterStar.Pz() / daughterStar.P(); + const float sinThetaStar = std::sqrt(std::max(0.0, 1.0 - cosThetaStar * cosThetaStar)); + const float sinPhiStar = std::sin(getPhiInRange(phiStar)); + const float cosPhiStar = std::cos(getPhiInRange(phiStar)); + + float polEP_A = std::sin(getPhiInRange(phiStar - psiA)); + float polEP_C = std::sin(getPhiInRange(phiStar - psiC)); + float polEP = std::sin(getPhiInRange(phiStar - psiFull)); + + const float qxFull = qxC - qxA; + const float qyFull = qyC - qyA; + + float polSP = uy * qxFull - ux * qyFull; + float polSP_A = uy * qxA - ux * qyA; + float polSP_C = uy * qxC - ux * qyC; + + float accDen = acceptanceWeight(isLambda, v0.eta(), v0.pt()); + float wgt = 1.0; // efficiency + + polSP /= accDen; + polSP_A /= accDen; + polSP_C /= accDen; + + polEP /= accDen; + polEP_A /= accDen; + polEP_C /= accDen; + + const float cosPsi = std::cos(getPhiInRange(psiFull)); + const float sinPsi = std::sin(getPhiInRange(psiFull)); + + if (isLambda) { + histos.fill(HIST("hSparseLambdaPolSPQ1"), mass, v0.pt(), polSP, centrality, q1, wgt); + histos.fill(HIST("hSparseLambdaPolSPAQ1"), mass, v0.pt(), polSP_A, centrality, q1, wgt); + histos.fill(HIST("hSparseLambdaPolSPCQ1"), mass, v0.pt(), polSP_C, centrality, q1, wgt); + + histos.fill(HIST("hSparseLambdaPolEPQ1"), mass, v0.pt(), polEP, centrality, q1, wgt); + histos.fill(HIST("hSparseLambdaPolEPAQ1"), mass, v0.pt(), polEP_A, centrality, q1, wgt); + histos.fill(HIST("hSparseLambdaPolEPCQ1"), mass, v0.pt(), polEP_C, centrality, q1, wgt); + + histos.fill(HIST("hSparseLambdaCosPsiQ1"), mass, v0.pt(), cosPsi, centrality, q1, wgt); + histos.fill(HIST("hSparseLambdaSinPsiQ1"), mass, v0.pt(), sinPsi, centrality, q1, wgt); + histos.fill(HIST("hSparseLambdaCorrSinPhiStarQ1"), mass, v0.pt(), sinPhiStar, centrality, q1, wgt); + histos.fill(HIST("hSparseLambdaCorrCosPhiStarQ1"), mass, v0.pt(), cosPhiStar, centrality, q1, wgt); + histos.fill(HIST("hSparseLambdaCorrSinThetaStarQ1"), mass, v0.pt(), sinThetaStar, centrality, q1, wgt); + histos.fill(HIST("hSparseLambdaAvgUxQ1"), mass, v0.pt(), ux, centrality, q1, wgt); + histos.fill(HIST("hSparseLambdaAvgUyQ1"), mass, v0.pt(), uy, centrality, q1, wgt); + } else { + histos.fill(HIST("hSparseAntiLambdaPolSPQ1"), mass, v0.pt(), polSP, centrality, q1, wgt); + histos.fill(HIST("hSparseAntiLambdaPolSPAQ1"), mass, v0.pt(), polSP_A, centrality, q1, wgt); + histos.fill(HIST("hSparseAntiLambdaPolSPCQ1"), mass, v0.pt(), polSP_C, centrality, q1, wgt); + + histos.fill(HIST("hSparseAntiLambdaPolEPQ1"), mass, v0.pt(), polEP, centrality, q1, wgt); + histos.fill(HIST("hSparseAntiLambdaPolEPAQ1"), mass, v0.pt(), polEP_A, centrality, q1, wgt); + histos.fill(HIST("hSparseAntiLambdaPolEPCQ1"), mass, v0.pt(), polEP_C, centrality, q1, wgt); + + histos.fill(HIST("hSparseAntiLambdaCosPsiQ1"), mass, v0.pt(), cosPsi, centrality, q1, wgt); + histos.fill(HIST("hSparseAntiLambdaSinPsiQ1"), mass, v0.pt(), sinPsi, centrality, q1, wgt); + histos.fill(HIST("hSparseAntiLambdaCorrSinPhiStarQ1"), mass, v0.pt(), sinPhiStar, centrality, q1, wgt); + histos.fill(HIST("hSparseAntiLambdaCorrCosPhiStarQ1"), mass, v0.pt(), cosPhiStar, centrality, q1, wgt); + histos.fill(HIST("hSparseAntiLambdaCorrSinThetaStarQ1"), mass, v0.pt(), sinThetaStar, centrality, q1, wgt); + histos.fill(HIST("hSparseAntiLambdaAvgUxQ1"), mass, v0.pt(), ux, centrality, q1, wgt); + histos.fill(HIST("hSparseAntiLambdaAvgUyQ1"), mass, v0.pt(), uy, centrality, q1, wgt); + } + } + + void processData(EventCandidates::iterator const& collision, AllTrackCandidates const& tracks, ResoV0s const& v0s, BCsRun3 const&) + { + if (!eventSelected(collision)) { + return; + } + const float centrality = collision.centFT0C(); + + if (cfgUseAccCorr) { + auto bc = collision.bc_as(); + currentRunNumber = bc.runNumber(); + if (currentRunNumber != lastRunNumber) { + accprofileL = ccdb->getForTimeStamp(cfgAccPathL.value, bc.timestamp()); + accprofileAL = ccdb->getForTimeStamp(cfgAccPathAL.value, bc.timestamp()); + lastRunNumber = currentRunNumber; + } + } + + float qxA = 0.0, qyA = 0.0, qxC = 0.0, qyC = 0.0; + float psiA = 0.0, psiC = 0.0, psiFull = 0.0, q1 = 0.0; + if (!getZDCQ1(collision, qxA, qyA, qxC, qyC, psiA, psiC, psiFull, q1)) { + return; + } + + histos.fill(HIST("hCentrality"), centrality); + histos.fill(HIST("hQ1ZDCvscent"), centrality, q1); + histos.fill(HIST("hQxAvscent"), centrality, qxA); + histos.fill(HIST("hQyAvscent"), centrality, qyA); + histos.fill(HIST("hQxCvscent"), centrality, qxC); + histos.fill(HIST("hQyCvscent"), centrality, qyC); + + float magA = std::sqrt(qxA * qxA + qyA * qyA); + float magC = std::sqrt(qxC * qxC + qyC * qyC); + float qxFull = qxC - qxA; + float qyFull = qyC - qyA; + float dotAC = qxA * qxC + qyA * qyC; + float resDot = dotAC / (magA * magC); + + histos.fill(HIST("hpResCosAC"), centrality, std::cos(getPhiInRange(psiA - psiC))); + histos.fill(HIST("hpResDotAC"), centrality, resDot); + histos.fill(HIST("hpQxAQxC"), centrality, qxA * qxC); + histos.fill(HIST("hpQyAQyC"), centrality, qyA * qyC); + histos.fill(HIST("hpQxAQyC"), centrality, qxA * qyC); + histos.fill(HIST("hpQxCQyA"), centrality, qxC * qyA); + + for (const auto& track : tracks) { + if (!selectV1Track(track)) { + continue; + } + + float phi = track.phi(); + float ux = std::cos(phi); + float uy = std::sin(phi); + + float v1SPFull = ux * qxFull + uy * qyFull; + float v1SPA = ux * qxA + uy * qyA; + float v1SPC = ux * qxC + uy * qyC; + + histos.fill(HIST("hV1SPFullQ1"), centrality, track.pt(), track.eta(), q1, v1SPFull); + histos.fill(HIST("hV1SPAQ1"), centrality, track.pt(), track.eta(), q1, v1SPA); + histos.fill(HIST("hV1SPCQ1"), centrality, track.pt(), track.eta(), q1, v1SPC); + + histos.fill(HIST("hV1EPFullQ1"), centrality, track.pt(), track.eta(), q1, std::cos(getPhiInRange(phi - psiFull))); + histos.fill(HIST("hV1EPAQ1"), centrality, track.pt(), track.eta(), q1, std::cos(getPhiInRange(phi - psiA))); + histos.fill(HIST("hV1EPCQ1"), centrality, track.pt(), track.eta(), q1, std::cos(getPhiInRange(phi - psiC))); + } + + for (const auto& v0 : v0s) { + bool lambdaTag = isCompatible(v0, 0); + bool antiLambdaTag = isCompatible(v0, 1); + if (!lambdaTag && !antiLambdaTag) { + continue; + } + if (!selectV0(collision, v0)) { + continue; + } + if (lambdaTag) { + fillPolarization(true, v0, centrality, q1, psiA, psiC, psiFull, qxA, qyA, qxC, qyC); + } + if (antiLambdaTag) { + fillPolarization(false, v0, centrality, q1, psiA, psiC, psiFull, qxA, qyA, qxC, qyC); + } + } + } + PROCESS_SWITCH(flowDirectedFlowTask, processData, "Process for data", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}