From 86859349e6c619701acac4f779073e901dd5d6f6 Mon Sep 17 00:00:00 2001 From: Prottay Das Date: Wed, 26 Nov 2025 15:02:29 +0100 Subject: [PATCH] QA code for eff calculation --- PWGLF/Tasks/Strangeness/CMakeLists.txt | 5 + PWGLF/Tasks/Strangeness/lambdak0seff.cxx | 458 +++++++++++++++++++++++ 2 files changed, 463 insertions(+) create mode 100644 PWGLF/Tasks/Strangeness/lambdak0seff.cxx diff --git a/PWGLF/Tasks/Strangeness/CMakeLists.txt b/PWGLF/Tasks/Strangeness/CMakeLists.txt index 2d902a8e597..333e2600d08 100644 --- a/PWGLF/Tasks/Strangeness/CMakeLists.txt +++ b/PWGLF/Tasks/Strangeness/CMakeLists.txt @@ -99,6 +99,11 @@ o2physics_add_dpl_workflow(lambdapolsp PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(lambdak0seff + SOURCES lambdak0seff.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(task-lambda-spin-corr SOURCES taskLambdaSpinCorr.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGLF/Tasks/Strangeness/lambdak0seff.cxx b/PWGLF/Tasks/Strangeness/lambdak0seff.cxx new file mode 100644 index 00000000000..2a3864f2ede --- /dev/null +++ b/PWGLF/Tasks/Strangeness/lambdak0seff.cxx @@ -0,0 +1,458 @@ +// 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. +// Fast Lambda k0s eff QA task for correlation analysis +// prottay.das@cern.ch + +#include "PWGLF/DataModel/LFStrangenessPIDTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/SPCalibrationTables.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/FT0Corrected.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CCDB/BasicCCDBManager.h" +#include "CommonConstants/PhysicsConstants.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsParameters/GRPObject.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include "Framework/StepTHn.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" + +#include "Math/GenVector/Boost.h" +#include "Math/Vector3D.h" +#include "Math/Vector4D.h" +#include "TF1.h" +#include "TRandom3.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include // <<< CHANGED: for dedup sets +#include +#include +#include // <<< CHANGED: for seenMap +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using std::array; + +struct lambdak0seff { + + struct : ConfigurableGroup { + Configurable cfgURL{"cfgURL", "http://alice-ccdb.cern.ch", "Address of the CCDB to browse"}; + Configurable nolaterthan{"ccdb-no-later-than", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "Latest acceptable timestamp of creation for the object"}; + } cfgCcdbParam; + + int mRunNumber; + Service ccdb; + Service pdg; + o2::ccdb::CcdbApi ccdbApi; + TH1D* hwgtAL; + // fill output + struct : ConfigurableGroup { + Configurable additionalEvSel{"additionalEvSel", false, "additionalEvSel"}; + Configurable additionalEvSel2{"additionalEvSel2", false, "additionalEvSel2"}; + Configurable additionalEvSel3{"additionalEvSel3", false, "additionalEvSel3"}; + } evselGrp; + // events + Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; + Configurable cfgCutCentralityMax{"cfgCutCentralityMax", 50.0f, "Accepted maximum Centrality"}; + Configurable cfgCutCentralityMin{"cfgCutCentralityMin", 30.0f, "Accepted minimum Centrality"}; + // proton track cut + Configurable cfgCutPT{"cfgCutPT", 0.15, "PT cut on daughter track"}; + Configurable cfgCutEta{"cfgCutEta", 0.8, "Eta cut on daughter track"}; + Configurable cfgCutDCAxy{"cfgCutDCAxy", 0.1f, "DCAxy range for tracks"}; + Configurable cfgCutDCAz{"cfgCutDCAz", 0.1f, "DCAz range for tracks"}; + Configurable cfgITScluster{"cfgITScluster", 5, "Number of ITS cluster"}; + Configurable cfgTPCcluster{"cfgTPCcluster", 70, "Number of TPC cluster"}; + Configurable isPVContributor{"isPVContributor", true, "is PV contributor"}; + // Configs for V0 + Configurable ConfV0PtMin{"ConfV0PtMin", 0.f, "Minimum transverse momentum of V0"}; + Configurable ConfV0Rap{"ConfV0Rap", 0.8f, "Rapidity range of V0"}; + Configurable ConfV0DCADaughMax{"ConfV0DCADaughMax", 0.2f, "Maximum DCA between the V0 daughters"}; + Configurable ConfV0CPAMin{"ConfV0CPAMin", 0.9998f, "Minimum CPA of V0"}; + Configurable ConfV0TranRadV0Min{"ConfV0TranRadV0Min", 1.5f, "Minimum transverse radius"}; + Configurable ConfV0TranRadV0Max{"ConfV0TranRadV0Max", 100.f, "Maximum transverse radius"}; + Configurable cMaxV0DCA{"cMaxV0DCA", 1.2, "Maximum V0 DCA to PV"}; + Configurable cMinV0DCAPr{"cMinV0DCAPr", 0.05, "Minimum V0 daughters DCA to PV for Pr"}; + Configurable cMinV0DCAPi{"cMinV0DCAPi", 0.05, "Minimum V0 daughters DCA to PV for Pi"}; + Configurable cMaxV0LifeTime{"cMaxV0LifeTime", 20, "Maximum V0 life time"}; + Configurable analyzeLambda{"analyzeLambda", true, "flag for lambda analysis"}; + Configurable analyzeK0s{"analyzeK0s", false, "flag for K0s analysis"}; + Configurable qtArmenterosMinForK0{"qtArmenterosMinForK0", 0.2, "Armenterous cut for K0s"}; + // config for V0 daughters + Configurable ConfDaughEta{"ConfDaughEta", 0.8f, "V0 Daugh sel: max eta"}; + Configurable cfgDaughPrPt{"cfgDaughPrPt", 0.4, "minimum daughter proton pt"}; + Configurable cfgDaughPiPt{"cfgDaughPiPt", 0.2, "minimum daughter pion pt"}; + Configurable rcrfc{"rcrfc", 0.8f, "Ratio of CR to FC"}; + Configurable ConfDaughTPCnclsMin{"ConfDaughTPCnclsMin", 50.f, "V0 Daugh sel: Min. nCls TPC"}; + Configurable ConfDaughPIDCuts{"ConfDaughPIDCuts", 3, "PID selections for Lambda daughters"}; + + struct : ConfigurableGroup { + Configurable IMNbins{"IMNbins", 100, "Number of bins in invariant mass"}; + Configurable lbinIM{"lbinIM", 1.0, "lower bin value in IM histograms"}; + Configurable hbinIM{"hbinIM", 1.2, "higher bin value in IM histograms"}; + } binGrp; + struct : ConfigurableGroup { + ConfigurableAxis configcentAxis{"configcentAxis", {VARIABLE_WIDTH, 0.0, 10.0, 40.0, 80.0, 150, 300}, "Cent FT0C"}; + ConfigurableAxis configthnAxispT{"configthnAxisPt", {VARIABLE_WIDTH, 0.2, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.5, 8.0, 10.0, 100.0}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis configetaAxis{"configetaAxis", {VARIABLE_WIDTH, -0.8, -0.4, -0.2, 0, 0.2, 0.4, 0.8}, "Eta"}; + ConfigurableAxis configvzAxis{"configvzAxis", {VARIABLE_WIDTH, -10, -5, -0.0, 5, 10}, "Vz"}; + } axisGrp; + + SliceCache cache; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + void init(o2::framework::InitContext&) + { + AxisSpec thnAxisInvMass{binGrp.IMNbins, binGrp.lbinIM, binGrp.hbinIM, "#it{M} (GeV/#it{c}^{2})"}; + + std::vector runaxes2 = {thnAxisInvMass, axisGrp.configthnAxispT, axisGrp.configetaAxis, axisGrp.configvzAxis, axisGrp.configcentAxis}; + + histos.add("hCentrality", "Centrality distribution", kTH1F, {{axisGrp.configcentAxis}}); + histos.add("hSparseGenLambda", "hSparseGenLambda", HistType::kTHnSparseF, runaxes2, true); + histos.add("hSparseGenAntiLambda", "hSparseGenAntiLambda", HistType::kTHnSparseF, runaxes2, true); + histos.add("hSparseRecLambda", "hSparseRecLambda", HistType::kTHnSparseF, runaxes2, true); + histos.add("hSparseRecAntiLambda", "hSparseRecAntiLambda", HistType::kTHnSparseF, runaxes2, true); + histos.add("hSparseGenK0s", "hSparseGenK0s", HistType::kTHnSparseF, runaxes2, true); + histos.add("hSparseRecK0s", "hSparseRecK0s", HistType::kTHnSparseF, runaxes2, true); + + ccdb->setURL(cfgCcdbParam.cfgURL); + ccdbApi.init("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + LOGF(info, "Getting alignment offsets from the CCDB..."); + } + + template + bool selectionTrack(const T& candidate) + { + if (!(candidate.isGlobalTrack() && candidate.isPVContributor() && candidate.itsNCls() > cfgITScluster && candidate.tpcNClsFound() > cfgTPCcluster && candidate.itsNClsInnerBarrel() >= 1)) { + return false; + } + return true; + } + + template + bool SelectionV0(Collision const& collision, V0 const& candidate) + { + if (TMath::Abs(candidate.dcav0topv()) > cMaxV0DCA) { + return false; + } + const float pT = candidate.pt(); + const float tranRad = candidate.v0radius(); + const float dcaDaughv0 = TMath::Abs(candidate.dcaV0daughters()); + const float cpav0 = candidate.v0cosPA(); + + float CtauLambda = candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * massLambda; + float CtauK0s = candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * massK0s; + + if (pT < ConfV0PtMin) { + return false; + } + if (dcaDaughv0 > ConfV0DCADaughMax) { + return false; + } + if (cpav0 < ConfV0CPAMin) { + return false; + } + if (tranRad < ConfV0TranRadV0Min) { + return false; + } + if (tranRad > ConfV0TranRadV0Max) { + return false; + } + if (analyzeLambda && TMath::Abs(CtauLambda) > cMaxV0LifeTime) { + return false; + } + if (analyzeK0s && TMath::Abs(CtauK0s) > cMaxV0LifeTime) { + return false; + } + if (analyzeLambda && TMath::Abs(candidate.yLambda()) > ConfV0Rap) { + return false; + } + if (analyzeK0s && TMath::Abs(candidate.yK0Short()) > ConfV0Rap) { + return false; + } + return true; + } + + template + bool isSelectedV0Daughter(V0 const& candidate, T const& track, int pid, int pid2) + { + const auto tpcNClsF = track.tpcNClsFound(); + if (track.tpcNClsCrossedRows() < cfgTPCcluster) { + return false; + } + + if (tpcNClsF < ConfDaughTPCnclsMin) { + return false; + } + if (track.tpcCrossedRowsOverFindableCls() < rcrfc) { + return false; + } + + if (analyzeLambda && pid == 0 && TMath::Abs(track.tpcNSigmaPr()) > ConfDaughPIDCuts) { + return false; + } + if (analyzeLambda && pid == 1 && TMath::Abs(track.tpcNSigmaPi()) > ConfDaughPIDCuts) { + return false; + } + if (analyzeK0s && TMath::Abs(track.tpcNSigmaPi()) > ConfDaughPIDCuts) { + return false; + } + if (pid == 0 && (candidate.positivept() < cfgDaughPrPt || candidate.negativept() < cfgDaughPiPt)) { + return false; // doesn´t pass lambda pT sels + } + if (pid == 1 && (candidate.positivept() < cfgDaughPiPt || candidate.negativept() < cfgDaughPrPt)) { + return false; // doesn´t pass antilambda pT sels + } + if (std::abs(candidate.positiveeta()) > ConfDaughEta || std::abs(candidate.negativeeta()) > ConfDaughEta) { + return false; + } + + if (analyzeLambda && pid2 == 0 && (TMath::Abs(candidate.dcapostopv()) < cMinV0DCAPr || TMath::Abs(candidate.dcanegtopv()) < cMinV0DCAPi)) { + return false; + } + if (analyzeLambda && pid2 == 1 && (TMath::Abs(candidate.dcapostopv()) < cMinV0DCAPi || TMath::Abs(candidate.dcanegtopv()) < cMinV0DCAPr)) { + return false; + } + if (analyzeK0s && (TMath::Abs(candidate.dcapostopv()) < cMinV0DCAPi || TMath::Abs(candidate.dcanegtopv()) < cMinV0DCAPi)) { + return false; + } + if (analyzeK0s && (candidate.qtarm() / (std::abs(candidate.alpha()))) < 0.2) { + return false; + } + + return true; + } + + bool shouldReject(bool LambdaTag, bool aLambdaTag, + const ROOT::Math::PxPyPzMVector& Lambdadummy, + const ROOT::Math::PxPyPzMVector& AntiLambdadummy) + { + const double minMass = 1.105; + const double maxMass = 1.125; + return (LambdaTag && aLambdaTag && + (Lambdadummy.M() > minMass && Lambdadummy.M() < maxMass) && + (AntiLambdadummy.M() > minMass && AntiLambdadummy.M() < maxMass)); + } + + ROOT::Math::PxPyPzMVector Lambda, AntiLambda, Lambdadummy, AntiLambdadummy, Proton, Pion, AntiProton, AntiPion, K0sdummy, K0s; + double massLambda = o2::constants::physics::MassLambda; + double massK0s = o2::constants::physics::MassK0Short; + double massPr = o2::constants::physics::MassProton; + double massPi = o2::constants::physics::MassPionCharged; + + Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; + Filter centralityFilter = (nabs(aod::cent::centFT0C) < cfgCutCentralityMax && nabs(aod::cent::centFT0C) > cfgCutCentralityMin); + Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && nabs(aod::track::pt) > cfgCutPT); + + using EventCandidatesMC = soa::Filtered>; + using AllTrackCandidates = soa::Filtered>; + using ResoV0s = aod::V0Datas; + + using TrackMCTrueTable = aod::McParticles; + ROOT::Math::PxPyPzMVector lambdadummymc, antiLambdadummymc, kshortdummymc, protonmc, pionmc, antiProtonmc, antiPionmc; + + void processMC(EventCandidatesMC::iterator const& collision, AllTrackCandidates const& /*tracks*/, TrackMCTrueTable const& GenParticles, ResoV0s const& V0s) + { + if (!collision.sel8()) { + return; + } + double centrality = -999.; + centrality = collision.centFT0C(); + // centrality = collision.multiplicity(); + double vz = collision.posZ(); + + if (evselGrp.additionalEvSel && (!collision.selection_bit(aod::evsel::kNoSameBunchPileup) || !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))) { + return; + } + if (evselGrp.additionalEvSel2 && (!collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || !collision.selection_bit(aod::evsel::kNoITSROFrameBorder))) { + return; + } + if (evselGrp.additionalEvSel3 && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { + return; + } + + histos.fill(HIST("hCentrality"), centrality); + + for (const auto& v0 : V0s) { + + auto postrack = v0.template posTrack_as(); + auto negtrack = v0.template negTrack_as(); + + if (analyzeLambda && analyzeK0s) + continue; + if (!analyzeLambda && !analyzeK0s) + continue; + + int LambdaTag = 0; + int aLambdaTag = 0; + int K0sTag = 0; + + const auto signpos = postrack.sign(); + const auto signneg = negtrack.sign(); + + if (signpos < 0 || signneg > 0) { + continue; + } + if (analyzeLambda) { + if (isSelectedV0Daughter(v0, postrack, 0, 0) && isSelectedV0Daughter(v0, negtrack, 1, 0)) { + LambdaTag = 1; + } + if (isSelectedV0Daughter(v0, negtrack, 0, 1) && isSelectedV0Daughter(v0, postrack, 1, 1)) { + aLambdaTag = 1; + } + } + if (analyzeK0s) { + if (isSelectedV0Daughter(v0, postrack, 0, 0) && isSelectedV0Daughter(v0, negtrack, 1, 0)) { + K0sTag = 1; + } + } + + if (analyzeLambda && (!LambdaTag && !aLambdaTag)) + continue; + if (analyzeK0s && (!K0sTag)) + continue; + + if (!SelectionV0(collision, v0)) { + continue; + } + + if (analyzeLambda) { + if (LambdaTag) { + Proton = ROOT::Math::PxPyPzMVector(v0.pxpos(), v0.pypos(), v0.pzpos(), massPr); + AntiPion = ROOT::Math::PxPyPzMVector(v0.pxneg(), v0.pyneg(), v0.pzneg(), massPi); + Lambdadummy = Proton + AntiPion; + } + if (aLambdaTag) { + AntiProton = ROOT::Math::PxPyPzMVector(v0.pxneg(), v0.pyneg(), v0.pzneg(), massPr); + Pion = ROOT::Math::PxPyPzMVector(v0.pxpos(), v0.pypos(), v0.pzpos(), massPi); + AntiLambdadummy = AntiProton + Pion; + } + + if (shouldReject(LambdaTag, aLambdaTag, Lambdadummy, AntiLambdadummy)) { + continue; + } + } + + if (analyzeK0s) { + if (K0sTag) { + Pion = ROOT::Math::PxPyPzMVector(v0.pxpos(), v0.pypos(), v0.pzpos(), massPi); + AntiPion = ROOT::Math::PxPyPzMVector(v0.pxneg(), v0.pyneg(), v0.pzneg(), massPi); + K0sdummy = Pion + AntiPion; + } + } + + if (TMath::Abs(v0.eta()) > 0.8) + continue; + + if (LambdaTag) { + Lambda = Proton + AntiPion; + histos.fill(HIST("hSparseRecLambda"), v0.mLambda(), v0.pt(), v0.eta(), vz, centrality); + } + if (aLambdaTag) { + AntiLambda = AntiProton + Pion; + histos.fill(HIST("hSparseRecAntiLambda"), v0.mAntiLambda(), v0.pt(), v0.eta(), vz, centrality); + } + if (K0sTag) { + histos.fill(HIST("hSparseRecK0s"), v0.mK0Short(), v0.pt(), v0.eta(), vz, centrality); + } + } + + for (const auto& mcParticle : GenParticles) { + if (analyzeLambda && std::abs(mcParticle.pdgCode()) != PDG_t::kLambda0) { + continue; + } + if (analyzeK0s && std::abs(mcParticle.pdgCode()) != PDG_t::kK0Short) { + continue; + } + if (std::abs(mcParticle.y()) > ConfV0Rap) { + continue; + } + auto pdg1 = mcParticle.pdgCode(); + auto kDaughters = mcParticle.daughters_as(); + int daughsize = 2; + if (kDaughters.size() != daughsize) { + continue; + } + for (const auto& kCurrentDaughter : kDaughters) { + + if (std::abs(kCurrentDaughter.pdgCode()) != PDG_t::kProton && std::abs(kCurrentDaughter.pdgCode()) != PDG_t::kPiPlus) { + continue; + } + if (kCurrentDaughter.pdgCode() == PDG_t::kProton) { + protonmc = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassProton); + } + if (kCurrentDaughter.pdgCode() == PDG_t::kPiMinus) { + antiPionmc = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassPionCharged); + } + + if (kCurrentDaughter.pdgCode() == PDG_t::kProtonBar) { + antiProtonmc = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassProton); + } + if (kCurrentDaughter.pdgCode() == PDG_t::kPiPlus) { + pionmc = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassPionCharged); + } + } + if (pdg1 == PDG_t::kLambda0) { + lambdadummymc = protonmc + antiPionmc; + histos.fill(HIST("hSparseGenLambda"), lambdadummymc.M(), lambdadummymc.Pt(), lambdadummymc.Eta(), vz, centrality); + } + + if (pdg1 == PDG_t::kLambda0Bar) { + antiLambdadummymc = antiProtonmc + pionmc; + histos.fill(HIST("hSparseGenAntiLambda"), antiLambdadummymc.M(), antiLambdadummymc.Pt(), lambdadummymc.Eta(), vz, centrality); + } + if (pdg1 == PDG_t::kK0Short) { + kshortdummymc = antiPionmc + pionmc; + histos.fill(HIST("hSparseGenK0s"), kshortdummymc.M(), kshortdummymc.Pt(), kshortdummymc.Eta(), vz, centrality); + } + } + } + PROCESS_SWITCH(lambdak0seff, processMC, "Process MC", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"lambdak0seff"})}; +}