From e74bb352f7a95da000691b539b59491086a2da10 Mon Sep 17 00:00:00 2001 From: Francesco Mazzaschi Date: Thu, 13 Aug 2026 13:48:06 +0200 Subject: [PATCH 1/2] snapshot, double omega tree creator added --- PWGLF/DataModel/LFDoubleCascTables.h | 72 -- PWGLF/DataModel/LFDoubleOmegaTables.h | 112 +++ .../TableProducer/Strangeness/CMakeLists.txt | 4 +- .../Strangeness/doubleCascTreeCreator.cxx | 352 -------- .../Strangeness/doubleOmegaTreeCreator.cxx | 760 ++++++++++++++++++ 5 files changed, 874 insertions(+), 426 deletions(-) delete mode 100644 PWGLF/DataModel/LFDoubleCascTables.h create mode 100644 PWGLF/DataModel/LFDoubleOmegaTables.h delete mode 100644 PWGLF/TableProducer/Strangeness/doubleCascTreeCreator.cxx create mode 100644 PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx diff --git a/PWGLF/DataModel/LFDoubleCascTables.h b/PWGLF/DataModel/LFDoubleCascTables.h deleted file mode 100644 index 03f77e8f34c..00000000000 --- a/PWGLF/DataModel/LFDoubleCascTables.h +++ /dev/null @@ -1,72 +0,0 @@ -// 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. - -#ifndef PWGLF_DATAMODEL_LFDOUBLECASCTABLES_H_ -#define PWGLF_DATAMODEL_LFDOUBLECASCTABLES_H_ - -#include - -namespace o2::aod -{ - -namespace DoubleCascTables -{ -DECLARE_SOA_COLUMN(PtCasc1, ptCasc1, float); // signed pt of the cascade -DECLARE_SOA_COLUMN(EtaCasc1, etaCasc1, float); -DECLARE_SOA_COLUMN(PhiCasc1, phiCasc1, float); -DECLARE_SOA_COLUMN(CascDecLength1, cascDecLength1, float); -DECLARE_SOA_COLUMN(OmegaMassCasc1, omegaMassCasc1, float); -DECLARE_SOA_COLUMN(XiMassCasc1, xiMassCasc1, float); -DECLARE_SOA_COLUMN(CosPACasc1, cosPACasc1, float); -DECLARE_SOA_COLUMN(DcaBachPVCasc1, dcaBachPVCasc1, float); -DECLARE_SOA_COLUMN(DcaV0BachCasc1, dcaV0BachCasc1, float); -DECLARE_SOA_COLUMN(NSigmaKBach1, nSigmaKBach1, float); - -DECLARE_SOA_COLUMN(PtCasc2, ptCasc2, float); -DECLARE_SOA_COLUMN(EtaCasc2, etaCasc2, float); -DECLARE_SOA_COLUMN(PhiCasc2, phiCasc2, float); -DECLARE_SOA_COLUMN(CascDecLength2, cascDecLength2, float); -DECLARE_SOA_COLUMN(OmegaMassCasc2, omegaMassCasc2, float); -DECLARE_SOA_COLUMN(XiMassCasc2, xiMassCasc2, float); -DECLARE_SOA_COLUMN(CosPACasc2, cosPACasc2, float); -DECLARE_SOA_COLUMN(DcaBachPVCasc2, dcaBachPVCasc2, float); -DECLARE_SOA_COLUMN(DcaV0BachCasc2, dcaV0BachCasc2, float); -DECLARE_SOA_COLUMN(NSigmaKBach2, nSigmaKBach2, float); - -DECLARE_SOA_COLUMN(DoubleOmegaMass, doubleOmegaMass, float); -} // namespace DoubleCascTables - -DECLARE_SOA_TABLE(DoubleCascTable, "AOD", "DOUBLECASCTABLE", - DoubleCascTables::PtCasc1, - DoubleCascTables::EtaCasc1, - DoubleCascTables::PhiCasc1, - DoubleCascTables::CascDecLength1, - DoubleCascTables::OmegaMassCasc1, - DoubleCascTables::XiMassCasc1, - DoubleCascTables::CosPACasc1, - DoubleCascTables::DcaBachPVCasc1, - DoubleCascTables::DcaV0BachCasc1, - DoubleCascTables::NSigmaKBach1, - DoubleCascTables::PtCasc2, - DoubleCascTables::EtaCasc2, - DoubleCascTables::PhiCasc2, - DoubleCascTables::CascDecLength2, - DoubleCascTables::OmegaMassCasc2, - DoubleCascTables::XiMassCasc2, - DoubleCascTables::CosPACasc2, - DoubleCascTables::DcaBachPVCasc2, - DoubleCascTables::DcaV0BachCasc2, - DoubleCascTables::NSigmaKBach2, - DoubleCascTables::DoubleOmegaMass); - -} // namespace o2::aod - -#endif // PWGLF_DATAMODEL_LFDOUBLECASCTABLES_H_ diff --git a/PWGLF/DataModel/LFDoubleOmegaTables.h b/PWGLF/DataModel/LFDoubleOmegaTables.h new file mode 100644 index 00000000000..dd14aa99004 --- /dev/null +++ b/PWGLF/DataModel/LFDoubleOmegaTables.h @@ -0,0 +1,112 @@ +// 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. + +#ifndef PWGLF_DATAMODEL_LFDOUBLEOMEGATABLES_H_ +#define PWGLF_DATAMODEL_LFDOUBLEOMEGATABLES_H_ + +#include + +#include + +namespace o2::aod +{ + +namespace DoubleOmegaTables +{ +DECLARE_SOA_COLUMN(PtCasc1, ptCasc1, float); // signed pt of the cascade +DECLARE_SOA_COLUMN(EtaCasc1, etaCasc1, float); +DECLARE_SOA_COLUMN(PhiCasc1, phiCasc1, float); +DECLARE_SOA_COLUMN(CascDecLength1, cascDecLength1, float); +DECLARE_SOA_COLUMN(OmegaMassCasc1, omegaMassCasc1, float); +DECLARE_SOA_COLUMN(XiMassCasc1, xiMassCasc1, float); +DECLARE_SOA_COLUMN(CosPACasc1, cosPACasc1, float); +DECLARE_SOA_COLUMN(DcaBachPVCasc1, dcaBachPVCasc1, float); +DECLARE_SOA_COLUMN(DcaV0BachCasc1, dcaV0BachCasc1, float); +DECLARE_SOA_COLUMN(NSigmaKBach1, nSigmaKBach1, float); + +DECLARE_SOA_COLUMN(PtLambda, ptLambda, float); +DECLARE_SOA_COLUMN(EtaLambda, etaLambda, float); +DECLARE_SOA_COLUMN(PhiLambda, phiLambda, float); +DECLARE_SOA_COLUMN(LambdaDecLength, lambdaDecLength, float); +DECLARE_SOA_COLUMN(LambdaMass, lambdaMass, float); +DECLARE_SOA_COLUMN(CosPALambda, cosPALambda, float); +DECLARE_SOA_COLUMN(DcaPosPVLambda, dcaPosPVLambda, float); +DECLARE_SOA_COLUMN(DcaNegPVLambda, dcaNegPVLambda, float); +DECLARE_SOA_COLUMN(DcaLambdaDaughters, dcaLambdaDaughters, float); +DECLARE_SOA_COLUMN(NSigmaPrLambda, nSigmaPrLambda, float); +DECLARE_SOA_COLUMN(NSigmaPiLambda, nSigmaPiLambda, float); + +DECLARE_SOA_COLUMN(PtKaon, ptKaon, float); +DECLARE_SOA_COLUMN(EtaKaon, etaKaon, float); +DECLARE_SOA_COLUMN(PhiKaon, phiKaon, float); +DECLARE_SOA_COLUMN(ChargeKaon, chargeKaon, int8_t); +DECLARE_SOA_COLUMN(DcaXYKaon, dcaXYKaon, float); +DECLARE_SOA_COLUMN(DcaZKaon, dcaZKaon, float); +DECLARE_SOA_COLUMN(NSigmaKKaon, nSigmaKKaon, float); + +DECLARE_SOA_COLUMN(DoubleOmegaMass, doubleOmegaMass, float); + +DECLARE_SOA_COLUMN(GenPt, genPt, float); +DECLARE_SOA_COLUMN(GenEta, genEta, float); +DECLARE_SOA_COLUMN(GenPhi, genPhi, float); +DECLARE_SOA_COLUMN(GenDecayLength, genDecayLength, float); +DECLARE_SOA_COLUMN(PdgDoubleOmega, pdgDoubleOmega, int); +DECLARE_SOA_COLUMN(IsReco, isReco, bool); +} // namespace DoubleOmegaTables + +#define DOUBLE_OMEGA_RECO_COLUMNS \ + DoubleOmegaTables::PtCasc1, \ + DoubleOmegaTables::EtaCasc1, \ + DoubleOmegaTables::PhiCasc1, \ + DoubleOmegaTables::CascDecLength1, \ + DoubleOmegaTables::OmegaMassCasc1, \ + DoubleOmegaTables::XiMassCasc1, \ + DoubleOmegaTables::CosPACasc1, \ + DoubleOmegaTables::DcaBachPVCasc1, \ + DoubleOmegaTables::DcaV0BachCasc1, \ + DoubleOmegaTables::NSigmaKBach1, \ + DoubleOmegaTables::PtLambda, \ + DoubleOmegaTables::EtaLambda, \ + DoubleOmegaTables::PhiLambda, \ + DoubleOmegaTables::LambdaDecLength, \ + DoubleOmegaTables::LambdaMass, \ + DoubleOmegaTables::CosPALambda, \ + DoubleOmegaTables::DcaPosPVLambda, \ + DoubleOmegaTables::DcaNegPVLambda, \ + DoubleOmegaTables::DcaLambdaDaughters, \ + DoubleOmegaTables::NSigmaPrLambda, \ + DoubleOmegaTables::NSigmaPiLambda, \ + DoubleOmegaTables::PtKaon, \ + DoubleOmegaTables::EtaKaon, \ + DoubleOmegaTables::PhiKaon, \ + DoubleOmegaTables::ChargeKaon, \ + DoubleOmegaTables::DcaXYKaon, \ + DoubleOmegaTables::DcaZKaon, \ + DoubleOmegaTables::NSigmaKKaon, \ + DoubleOmegaTables::DoubleOmegaMass + +DECLARE_SOA_TABLE(DoubleOmegaTable, "AOD", "DOUBLEOMEGATABLE", + DOUBLE_OMEGA_RECO_COLUMNS); + +DECLARE_SOA_TABLE(DoubleOmegaTableMC, "AOD", "DBLOMEGAMCTABLE", + DOUBLE_OMEGA_RECO_COLUMNS, + DoubleOmegaTables::GenPt, + DoubleOmegaTables::GenEta, + DoubleOmegaTables::GenPhi, + DoubleOmegaTables::GenDecayLength, + DoubleOmegaTables::PdgDoubleOmega, + DoubleOmegaTables::IsReco); + +#undef DOUBLE_OMEGA_RECO_COLUMNS + +} // namespace o2::aod + +#endif // PWGLF_DATAMODEL_LFDOUBLEOMEGATABLES_H_ diff --git a/PWGLF/TableProducer/Strangeness/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/CMakeLists.txt index 83b05e9a1bc..a6cd708d50f 100644 --- a/PWGLF/TableProducer/Strangeness/CMakeLists.txt +++ b/PWGLF/TableProducer/Strangeness/CMakeLists.txt @@ -56,8 +56,8 @@ o2physics_add_dpl_workflow(hstrangecorrelationfilter PUBLIC_LINK_LIBRARIES O2::DCAFitter O2Physics::AnalysisCore O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(double-casc-tree-creator - SOURCES doubleCascTreeCreator.cxx +o2physics_add_dpl_workflow(double-omega-tree-creator + SOURCES doubleOmegaTreeCreator.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) diff --git a/PWGLF/TableProducer/Strangeness/doubleCascTreeCreator.cxx b/PWGLF/TableProducer/Strangeness/doubleCascTreeCreator.cxx deleted file mode 100644 index 736cda47318..00000000000 --- a/PWGLF/TableProducer/Strangeness/doubleCascTreeCreator.cxx +++ /dev/null @@ -1,352 +0,0 @@ -// 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. - -#include "PWGLF/DataModel/LFDoubleCascTables.h" -#include "PWGLF/DataModel/LFStrangenessTables.h" - -#include "Common/CCDB/EventSelectionParams.h" -#include "Common/Core/Zorro.h" -#include "Common/Core/ZorroSummary.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/PIDResponseTPC.h" - -#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 Collisions = soa::Join::iterator; -using FullCascades = aod::CascDataExt; -using TracksFull = soa::Join; - -struct doubleCascCand { - float ptCasc1 = -999.f; // signed pt of the cascade - float etaCasc1 = -999.f; - float phiCasc1 = -999.f; - float cascDecLength1 = -999.f; - float omegaMassCasc1 = -999.f; - float xiMassCasc1 = -999.f; - float cosPACasc1 = -999.f; - float dcaBachPVCasc1 = -999.f; - float dcaV0BachCasc1 = -999.f; - float nSigmaKBach1 = -999.f; - - float ptCasc2 = -999.f; - float etaCasc2 = -999.f; - float phiCasc2 = -999.f; - float cascDecLength2 = -999.f; - float omegaMassCasc2 = -999.f; - float xiMassCasc2 = -999.f; - float cosPACasc2 = -999.f; - float dcaBachPVCasc2 = -999.f; - float dcaV0BachCasc2 = -999.f; - float nSigmaKBach2 = -999.f; - float doubleOmegaMass = -999.f; -}; - -struct doubleCascTreeCreator { - Produces doubleCascTable; - std::vector doubleCascCands; - Service ccdb; - o2::vertexing::DCAFitterN<2> fitter; - - int mRunNumber; - - Zorro zorro; - OutputObj zorroSummary{"zorroSummary"}; - - Configurable cfgSkimmedProcessing{"cfgSkimmedProcessing", false, "Skimmed dataset processing"}; - Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - - ConfigurableAxis centAxis{"centAxis", {106, 0, 106}, "binning for the centrality"}; - ConfigurableAxis zVtxAxis{"zVtxBins", {100, -20.f, 20.f}, "Binning for the vertex z in cm"}; - - // binning of (anti)lambda mass QA histograms - ConfigurableAxis massOmegaAxis{"massOmegaAxis", {400, o2::constants::physics::MassOmegaMinus - 0.05f, o2::constants::physics::MassOmegaMinus + 0.05}, "binning for the Omega invariant-mass"}; - ConfigurableAxis massXiAxis{"massXiAxis", {400, o2::constants::physics::MassXiMinus - 0.05f, o2::constants::physics::MassXiMinus + 0.05f}, "binning for the Xi invariant-mass"}; - - Configurable zVtxMax{"zVtxMax", 10.0f, "maximum z position of the primary vertex"}; - Configurable etaMax{"etaMax", 0.9f, "maximum eta"}; - ConfigurableAxis momAxis{"momAxisFine", {5.e2, 0.f, 5.f}, "momentum axis binning"}; - - Configurable cascPtMin{"cascPtMin", 1.f, "minimum (anti)casc pT (GeV/c)"}; - Configurable cascPtMax{"cascPtMax", 5.f, "maximum (anti)casc pT (GeV/c)"}; - - Configurable minNCrossedRows{"minNCrossedRows", 100, "Minimum number of crossed TPC rows"}; - Configurable minNITSClus{"minNITSClus", 0., "Minimum number of ITS clusters"}; - Configurable minNTPCClus{"minNTPCClus", 80, "Minimum number of TPC clusters"}; - Configurable maxNSharedTPCClus{"maxNSharedTPCClus", 5, "Maximum number of shared TPC clusters"}; - - Configurable minCascCosPA{"minCascCosPA", 0.99f, "Minimum cosine of the pointing angle of the cascade"}; - Configurable nSigmaTPCCut{"nSigmaTPCCut", 3.f, "Number of sigmas for the TPC PID"}; - Configurable dcaBachToPV{"dcaBachToPV", 0.05f, "DCA of the bachelor to the primary vertex"}; - Configurable dcaV0Bach{"dcaV0Bach", 1.f, "DCA between the V0 daughters"}; - Configurable mXiWindow{"mXiWindow", 0.02f, "mXiWindow"}; - Configurable mOmegaWindow{"mOmegaWindow", 0.01f, "mOmegaWindow"}; - - HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - - template - bool selectTrack(T const& track) - { - if (std::abs(track.eta()) > etaMax) { - return false; - } - if (track.itsNCls() < minNITSClus || - track.tpcNClsFound() < minNTPCClus || - track.tpcNClsCrossedRows() < minNCrossedRows || - track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable() || - track.tpcNClsShared() > maxNSharedTPCClus) { - return false; - } - return true; - } - - template - void initCCDB(Bc const& bc) - { - if (mRunNumber == bc.runNumber()) { - return; - } - auto timestamp = bc.timestamp(); - LOG(info) << "Retrieved GRP for timestamp " << timestamp; - mRunNumber = bc.runNumber(); - if (cfgSkimmedProcessing) { - zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), "fDoubleOmega,fOmegaXi"); - zorro.populateHistRegistry(histos, bc.runNumber()); - } - } - - void init(o2::framework::InitContext&) - { - mRunNumber = 0; - ccdb->setURL(ccdburl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - ccdb->setFatalWhenNull(false); - - zorroSummary.setObject(zorro.getZorroSummary()); - - // event QA - histos.add("QA/zVtx", ";#it{z}_{vtx} (cm);Entries", HistType::kTH1F, {zVtxAxis}); - histos.add("QA/massXi1", ";#it{p}_{T} (GeV/#it{c});#it{M}(#Lambda + #pi^{-}) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massXiAxis}); - histos.add("QA/massOmega1", ";#it{p}_{T} (GeV/#it{c});#it{M}(#Omega + #pi^{-}) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massOmegaAxis}); - histos.add("QA/massXi2", ";#it{p}_{T} (GeV/#it{c});#it{M}(#Lambda + #pi^{-}) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massXiAxis}); - histos.add("QA/massOmega2", ";#it{p}_{T} (GeV/#it{c});#it{M}(#Omega + #pi^{-}) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massOmegaAxis}); - } - - template - bool isSelectedCasc(C const& collision, T const&, FullCascades::iterator const& casc) - { - - auto bachelor = casc.bachelor_as(); - auto posDau = casc.posTrack_as(); - auto negDau = casc.negTrack_as(); - - if (!selectTrack(bachelor) || !selectTrack(posDau) || !selectTrack(negDau)) { - return false; - } - if (casc.sign() > 0) { - if (TMath::Abs(posDau.tpcNSigmaPi()) > nSigmaTPCCut || TMath::Abs(negDau.tpcNSigmaPr()) > nSigmaTPCCut) { - return false; - } - } else if (casc.sign() < 0) { - if (TMath::Abs(negDau.tpcNSigmaPi()) > nSigmaTPCCut || TMath::Abs(posDau.tpcNSigmaPr()) > nSigmaTPCCut) { - return false; - } - } - if (TMath::Abs(casc.dcabachtopv()) < dcaBachToPV) { - return false; - } - if (TMath::Abs(casc.dcacascdaughters()) > dcaV0Bach) { - return false; - } - if (casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()) < minCascCosPA) { - return false; - } - if (TMath::Abs(casc.eta()) > etaMax) { - return false; - } - // mass cuts - bool massInWindow = false; - if (casc.mOmega() > o2::constants::physics::MassOmegaMinus - mOmegaWindow && casc.mOmega() < o2::constants::physics::MassOmegaMinus + mOmegaWindow) { - massInWindow = true; - } - if (casc.mXi() > o2::constants::physics::MassXiMinus - mXiWindow && casc.mXi() < o2::constants::physics::MassXiMinus + mXiWindow) { - massInWindow = true; - } - if (!massInWindow) { - return false; - } - return true; - }; - - template - float doubleOmegaMass(T const&, FullCascades::iterator const& casc1, FullCascades::iterator const& casc2) - { - // the fake omega decay is the one with the smaller radius - auto& fakeOmega = casc1.cascradius() < casc2.cascradius() ? casc1 : casc2; - auto& realOmega = casc1.cascradius() < casc2.cascradius() ? casc2 : casc1; - auto kaon = fakeOmega.bachelor_as(); - float momKaon[3] = {kaon.px(), kaon.py(), kaon.pz()}; - float momLambda[3] = {fakeOmega.pxlambda(), fakeOmega.pylambda(), fakeOmega.pzlambda()}; - // now compute real Omega-lambda-kaon mass - float momTot[3] = {momKaon[0] + momLambda[0] + realOmega.px(), momKaon[1] + momLambda[1] + realOmega.py(), momKaon[2] + momLambda[2] + realOmega.pz()}; - float eK = std::sqrt(o2::constants::physics::MassKaonCharged * o2::constants::physics::MassKaonCharged + momKaon[0] * momKaon[0] + momKaon[1] * momKaon[1] + momKaon[2] * momKaon[2]); - float eL = std::sqrt(o2::constants::physics::MassLambda0 * o2::constants::physics::MassLambda0 + momLambda[0] * momLambda[0] + momLambda[1] * momLambda[1] + momLambda[2] * momLambda[2]); - float eO = std::sqrt(o2::constants::physics::MassOmegaMinus * o2::constants::physics::MassOmegaMinus + realOmega.px() * realOmega.px() + realOmega.py() * realOmega.py() + realOmega.pz() * realOmega.pz()); - float eTot = eK + eL + eO; - float mass = std::sqrt(eTot * eTot - momTot[0] * momTot[0] - momTot[1] * momTot[1] - momTot[2] * momTot[2]); - return mass; - } - - template - void fillDoubleCasc(C const& collision, T const& tracks, FullCascades const& cascades) - { - doubleCascCands.clear(); - - for (auto& casc1 : cascades) { - if (!isSelectedCasc(collision, tracks, casc1)) { - continue; - } - histos.fill(HIST("QA/massXi1"), casc1.pt(), casc1.mXi()); - histos.fill(HIST("QA/massOmega1"), casc1.pt(), casc1.mOmega()); - for (auto& casc2 : cascades) { - if (!isSelectedCasc(collision, tracks, casc2)) { - continue; - } - histos.fill(HIST("QA/massXi2"), casc2.pt(), casc2.mXi()); - histos.fill(HIST("QA/massOmega2"), casc2.pt(), casc2.mOmega()); - - // check that the cascades do not share any track - std::vector trackIdsCasc1 = {casc1.posTrackId(), casc1.negTrackId(), casc1.bachelorId()}; - std::vector trackIdsCasc2 = {casc2.posTrackId(), casc2.negTrackId(), casc2.bachelorId()}; - bool shareTrack = false; - for (auto id1 : trackIdsCasc1) { - for (auto id2 : trackIdsCasc2) { - if (id1 == id2) { - shareTrack = true; - break; - } - } - if (shareTrack) { - break; - } - } - if (shareTrack) { - continue; - } - - auto bach1 = casc1.bachelor_as(); - auto bach2 = casc2.bachelor_as(); - - doubleCascCand cand; - cand.ptCasc1 = casc1.pt(); - cand.etaCasc1 = casc1.eta(); - cand.phiCasc1 = casc1.phi(); - cand.cascDecLength1 = std::hypot(casc1.x() - collision.posX(), casc1.y() - collision.posY(), casc1.z() - collision.posZ()); - cand.omegaMassCasc1 = casc1.mOmega(); - cand.xiMassCasc1 = casc1.mXi(); - cand.cosPACasc1 = casc1.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); - cand.dcaBachPVCasc1 = casc1.dcabachtopv(); - cand.dcaV0BachCasc1 = casc1.dcacascdaughters(); - cand.nSigmaKBach1 = bach1.tpcNSigmaKa(); - - cand.ptCasc2 = casc2.pt(); - cand.etaCasc2 = casc2.eta(); - cand.phiCasc2 = casc2.phi(); - cand.cascDecLength2 = std::hypot(casc2.x() - collision.posX(), casc2.y() - collision.posY(), casc2.z() - collision.posZ()); - cand.omegaMassCasc2 = casc2.mOmega(); - cand.xiMassCasc2 = casc2.mXi(); - cand.cosPACasc2 = casc2.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); - cand.dcaBachPVCasc2 = casc2.dcabachtopv(); - cand.dcaV0BachCasc2 = casc2.dcacascdaughters(); - cand.nSigmaKBach2 = bach2.tpcNSigmaKa(); - - cand.doubleOmegaMass = doubleOmegaMass(tracks, casc1, casc2); - - doubleCascCands.push_back(cand); - } - } - }; - - void processData(Collisions const& collision, TracksFull const& tracks, FullCascades const& cascades, aod::BCsWithTimestamps const&) - { - auto bc = collision.bc_as(); - initCCDB(bc); - - if (!collision.sel8()) - return; - - if (std::abs(collision.posZ()) > zVtxMax) - return; - - if (!collision.selection_bit(aod::evsel::kNoITSROFrameBorder) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) - return; - - if (cfgSkimmedProcessing) { - zorro.isSelected(collision.bc_as().globalBC()); /// Just let Zorro do the accounting - } - histos.fill(HIST("QA/zVtx"), collision.posZ()); - fillDoubleCasc(collision, tracks, cascades); - - for (auto& cand : doubleCascCands) { - doubleCascTable( - cand.ptCasc1, - cand.etaCasc1, - cand.phiCasc1, - cand.cascDecLength1, - cand.omegaMassCasc1, - cand.xiMassCasc1, - cand.cosPACasc1, - cand.dcaBachPVCasc1, - cand.dcaV0BachCasc1, - cand.nSigmaKBach1, - cand.ptCasc2, - cand.etaCasc2, - cand.phiCasc2, - cand.cascDecLength2, - cand.omegaMassCasc2, - cand.xiMassCasc2, - cand.cosPACasc2, - cand.dcaBachPVCasc2, - cand.dcaV0BachCasc2, - cand.nSigmaKBach2, - cand.doubleOmegaMass); - } - } - PROCESS_SWITCH(doubleCascTreeCreator, processData, "process (Run 3)", false); -}; - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; -} diff --git a/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx b/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx new file mode 100644 index 00000000000..a377cd2bee1 --- /dev/null +++ b/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx @@ -0,0 +1,760 @@ +// 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. + +#include "PWGLF/DataModel/LFDoubleOmegaTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/Core/Zorro.h" +#include "Common/Core/ZorroSummary.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTPC.h" + +#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 CollisionsTable = soa::Join; +using Collisions = CollisionsTable::iterator; +using CollisionsMC = soa::Join; +using FullCascades = aod::CascDataExt; +using FullV0s = aod::V0Datas; +using TracksFull = soa::Join; + +struct DoubleOmegaCandidate { + float ptCasc1 = -999.f; + float etaCasc1 = -999.f; + float phiCasc1 = -999.f; + float cascDecLength1 = -999.f; + float omegaMassCasc1 = -999.f; + float xiMassCasc1 = -999.f; + float cosPACasc1 = -999.f; + float dcaBachPVCasc1 = -999.f; + float dcaV0BachCasc1 = -999.f; + float nSigmaKBach1 = -999.f; + + float ptLambda = -999.f; + float etaLambda = -999.f; + float phiLambda = -999.f; + float lambdaDecLength = -999.f; + float lambdaMass = -999.f; + float cosPALambda = -999.f; + float dcaPosPVLambda = -999.f; + float dcaNegPVLambda = -999.f; + float dcaLambdaDaughters = -999.f; + float nSigmaPrLambda = -999.f; + float nSigmaPiLambda = -999.f; + + float ptKaon = -999.f; + float etaKaon = -999.f; + float phiKaon = -999.f; + int8_t chargeKaon = 0; + float dcaXYKaon = -999.f; + float dcaZKaon = -999.f; + float nSigmaKKaon = -999.f; + + float doubleOmegaMass = -999.f; +}; + +struct DoubleOmegaMCInfo { + int64_t motherId = -1; + float pt = -999.f; + float eta = -999.f; + float phi = -999.f; + float decayLength = -999.f; + int pdgCode = 0; +}; + +struct LambdaCandidate { + float px = 0.f; + float py = 0.f; + float pz = 0.f; + float decayLength = -999.f; + float mass = -999.f; + float cosPA = -999.f; + float dcaPosToPV = -999.f; + float dcaNegToPV = -999.f; + float dcaDaughters = -999.f; + float nSigmaPr = -999.f; + float nSigmaPi = -999.f; +}; + +struct doubleOmegaTreeCreator { + Produces doubleOmegaTable; + Produces doubleOmegaTableMC; + Service ccdb; + + static constexpr int kDoubleOmegaPdg = 1060020020; + enum FindableDaughter : uint8_t { + kOmega, + kLambda, + kKaon, + kNFindableDaughters + }; + std::vector reconstructedDoubleOmegaIds; + + Preslice tracksPerCollision = aod::track::collisionId; + Preslice cascadesPerCollision = aod::cascade::collisionId; + Preslice v0sPerCollision = aod::v0data::collisionId; + + int mRunNumber = 0; + + Zorro zorro; + OutputObj zorroSummary{"zorroSummary"}; + + Configurable cfgSkimmedProcessing{"cfgSkimmedProcessing", false, "Skimmed dataset processing"}; + Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + + ConfigurableAxis zVtxAxis{"zVtxBins", {100, -20.f, 20.f}, "Binning for the vertex z in cm"}; + ConfigurableAxis massOmegaAxis{"massOmegaAxis", {400, o2::constants::physics::MassOmegaMinus - 0.05f, o2::constants::physics::MassOmegaMinus + 0.05f}, "binning for the Omega invariant-mass"}; + ConfigurableAxis massXiAxis{"massXiAxis", {400, o2::constants::physics::MassXiMinus - 0.05f, o2::constants::physics::MassXiMinus + 0.05f}, "binning for the Xi invariant-mass"}; + ConfigurableAxis massLambdaAxis{"massLambdaAxis", {400, o2::constants::physics::MassLambda0 - 0.05f, o2::constants::physics::MassLambda0 + 0.05f}, "binning for the Lambda invariant-mass"}; + ConfigurableAxis momAxis{"momAxisFine", {5.e2, 0.f, 5.f}, "momentum axis binning"}; + + Configurable zVtxMax{"zVtxMax", 10.0f, "maximum z position of the primary vertex"}; + Configurable etaMax{"etaMax", 0.9f, "maximum eta"}; + Configurable cascPtMin{"cascPtMin", 1.f, "minimum (anti)cascade pT (GeV/c)"}; + Configurable cascPtMax{"cascPtMax", 5.f, "maximum (anti)cascade pT (GeV/c)"}; + + Configurable minNCrossedRows{"minNCrossedRows", 100, "Minimum number of crossed TPC rows"}; + Configurable minNITSClus{"minNITSClus", 0., "Minimum number of ITS clusters"}; + Configurable minNTPCClus{"minNTPCClus", 80, "Minimum number of TPC clusters"}; + Configurable maxNSharedTPCClus{"maxNSharedTPCClus", 5, "Maximum number of shared TPC clusters"}; + + Configurable minCascCosPA{"minCascCosPA", 0.99f, "Minimum cosine of the pointing angle of the cascade"}; + Configurable minLambdaCosPA{"minLambdaCosPA", 0.995f, "Minimum cosine of the pointing angle of the Lambda"}; + Configurable nSigmaTPCCut{"nSigmaTPCCut", 3.f, "Number of sigmas for the TPC PID"}; + Configurable dcaBachToPV{"dcaBachToPV", 0.05f, "Minimum DCA of a cascade bachelor to the primary vertex"}; + Configurable dcaKaonToPV{"dcaKaonToPV", 0.05f, "Minimum transverse DCA of the kaon to the primary vertex"}; + Configurable dcaV0DauToPV{"dcaV0DauToPV", 0.05f, "Minimum DCA of Lambda daughters to the primary vertex"}; + Configurable dcaV0Bach{"dcaV0Bach", 1.f, "Maximum DCA between the V0 and cascade bachelor"}; + Configurable dcaLambdaDaughters{"dcaLambdaDaughters", 1.f, "Maximum DCA between Lambda daughters"}; + Configurable mXiWindow{"mXiWindow", 0.02f, "Xi mass window used by the cascade compatibility mode"}; + Configurable mOmegaWindow{"mOmegaWindow", 0.01f, "Omega mass window"}; + Configurable mLambdaWindow{"mLambdaWindow", 0.01f, "Lambda mass window"}; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + template + bool selectTrack(T const& track) + { + if (std::abs(track.eta()) > etaMax) { + return false; + } + if (track.itsNCls() < minNITSClus || + track.tpcNClsFound() < minNTPCClus || + track.tpcNClsCrossedRows() < minNCrossedRows || + track.tpcNClsCrossedRows() < 0.8f * track.tpcNClsFindable() || + track.tpcNClsShared() > maxNSharedTPCClus) { + return false; + } + return true; + } + + template + void initCCDB(Bc const& bc) + { + if (mRunNumber == bc.runNumber()) { + return; + } + LOG(info) << "Retrieved GRP for timestamp " << bc.timestamp(); + mRunNumber = bc.runNumber(); + if (cfgSkimmedProcessing) { + zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), "fDoubleOmega,fOmegaXi"); + zorro.populateHistRegistry(histos, bc.runNumber()); + } + } + + void init(InitContext const&) + { + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + zorroSummary.setObject(zorro.getZorroSummary()); + + histos.add("QA/zVtx", ";#it{z}_{vtx} (cm);Entries", HistType::kTH1F, {zVtxAxis}); + histos.add("QA/massXi", ";#it{p}_{T} (GeV/#it{c});#it{M}(#Lambda + #pi) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massXiAxis}); + histos.add("QA/massOmega", ";#it{p}_{T} (GeV/#it{c});#it{M}(#Lambda + K) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massOmegaAxis}); + histos.add("QA/massLambda", ";#it{p}_{T} (GeV/#it{c});#it{M}(p + #pi) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massLambdaAxis}); + histos.add("MC/findableDaughters", ";Findable daughter;Entries", HistType::kTH1F, {{kNFindableDaughters, -0.5f, static_cast(kNFindableDaughters) - 0.5f}}); + auto findableDaughters = histos.get(HIST("MC/findableDaughters")); + findableDaughters->GetXaxis()->SetBinLabel(kOmega + 1, "#Omega"); + findableDaughters->GetXaxis()->SetBinLabel(kLambda + 1, "#Lambda"); + findableDaughters->GetXaxis()->SetBinLabel(kKaon + 1, "K"); + } + + template + bool isSelectedOmega(C const& collision, T const&, Casc const& casc) + { + auto bachelor = casc.template bachelor_as(); + auto posDau = casc.template posTrack_as(); + auto negDau = casc.template negTrack_as(); + + if (!selectTrack(bachelor) || !selectTrack(posDau) || !selectTrack(negDau)) { + return false; + } + if (casc.sign() > 0) { + if (std::abs(posDau.tpcNSigmaPi()) > nSigmaTPCCut || std::abs(negDau.tpcNSigmaPr()) > nSigmaTPCCut) { + return false; + } + } else if (casc.sign() < 0) { + if (std::abs(negDau.tpcNSigmaPi()) > nSigmaTPCCut || std::abs(posDau.tpcNSigmaPr()) > nSigmaTPCCut) { + return false; + } + } else { + return false; + } + if (std::abs(bachelor.tpcNSigmaKa()) > nSigmaTPCCut || + std::abs(casc.dcabachtopv()) < dcaBachToPV || + std::abs(casc.dcacascdaughters()) > dcaV0Bach || + casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()) < minCascCosPA || + std::abs(casc.eta()) > etaMax || + std::abs(casc.mOmega() - o2::constants::physics::MassOmegaMinus) > mOmegaWindow) { + return false; + } + return true; + } + + template + bool isSelectedLambda(C const&, T const&, V0 const& v0, int8_t charge) + { + auto posDau = v0.template posTrack_as(); + auto negDau = v0.template negTrack_as(); + if (!selectTrack(posDau) || !selectTrack(negDau) || + std::abs(v0.eta()) > etaMax || + v0.v0cosPA() < minLambdaCosPA || + std::abs(v0.dcapostopv()) < dcaV0DauToPV || + std::abs(v0.dcanegtopv()) < dcaV0DauToPV || + std::abs(v0.dcaV0daughters()) > dcaLambdaDaughters) { + return false; + } + + const bool isMatter = charge < 0; + const float mass = isMatter ? v0.mLambda() : v0.mAntiLambda(); + const float protonNSigma = isMatter ? posDau.tpcNSigmaPr() : negDau.tpcNSigmaPr(); + const float pionNSigma = isMatter ? negDau.tpcNSigmaPi() : posDau.tpcNSigmaPi(); + return std::abs(mass - o2::constants::physics::MassLambda0) <= mLambdaWindow && + std::abs(protonNSigma) <= nSigmaTPCCut && + std::abs(pionNSigma) <= nSigmaTPCCut; + } + + template + bool isSelectedKaon(Track const& track, int8_t charge) + { + return track.sign() == charge && + selectTrack(track) && + std::abs(track.tpcNSigmaKa()) <= nSigmaTPCCut && + std::abs(track.dcaXY()) >= dcaKaonToPV; + } + + template + LambdaCandidate getLambdaCandidateFromV0(C const& collision, T const&, V0 const& v0, int8_t charge) + { + auto posDau = v0.template posTrack_as(); + auto negDau = v0.template negTrack_as(); + const bool isMatter = charge < 0; + return {v0.px(), + v0.py(), + v0.pz(), + std::hypot(v0.x() - collision.posX(), v0.y() - collision.posY(), v0.z() - collision.posZ()), + isMatter ? v0.mLambda() : v0.mAntiLambda(), + v0.v0cosPA(), + v0.dcapostopv(), + v0.dcanegtopv(), + v0.dcaV0daughters(), + isMatter ? posDau.tpcNSigmaPr() : negDau.tpcNSigmaPr(), + isMatter ? negDau.tpcNSigmaPi() : posDau.tpcNSigmaPi()}; + } + + template + LambdaCandidate getLambdaCandidateFromCascade(C const& collision, T const&, Casc const& casc) + { + auto posDau = casc.template posTrack_as(); + auto negDau = casc.template negTrack_as(); + const bool isMatter = casc.sign() < 0; + return {casc.pxlambda(), + casc.pylambda(), + casc.pzlambda(), + std::hypot(casc.xlambda() - collision.posX(), casc.ylambda() - collision.posY(), casc.zlambda() - collision.posZ()), + casc.mLambda(), + casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()), + casc.dcapostopv(), + casc.dcanegtopv(), + casc.dcaV0daughters(), + isMatter ? posDau.tpcNSigmaPr() : negDau.tpcNSigmaPr(), + isMatter ? negDau.tpcNSigmaPi() : posDau.tpcNSigmaPi()}; + } + + static float etaFromMomentum(float px, float py, float pz) + { + const float pt = std::hypot(px, py); + return pt > 0.f ? std::asinh(pz / pt) : 0.f; + } + + template + DoubleOmegaCandidate makeCandidate(C const& collision, Omega const& omega, LambdaCandidate const& lambda, Kaon const& kaon) + { + DoubleOmegaCandidate cand; + auto omegaBach = omega.template bachelor_as(); + + cand.ptCasc1 = omega.pt(); + cand.etaCasc1 = omega.eta(); + cand.phiCasc1 = omega.phi(); + cand.cascDecLength1 = std::hypot(omega.x() - collision.posX(), omega.y() - collision.posY(), omega.z() - collision.posZ()); + cand.omegaMassCasc1 = omega.mOmega(); + cand.xiMassCasc1 = omega.mXi(); + cand.cosPACasc1 = omega.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + cand.dcaBachPVCasc1 = omega.dcabachtopv(); + cand.dcaV0BachCasc1 = omega.dcacascdaughters(); + cand.nSigmaKBach1 = omegaBach.tpcNSigmaKa(); + + cand.ptLambda = std::hypot(lambda.px, lambda.py); + cand.etaLambda = etaFromMomentum(lambda.px, lambda.py, lambda.pz); + cand.phiLambda = std::atan2(lambda.py, lambda.px); + cand.lambdaDecLength = lambda.decayLength; + cand.lambdaMass = lambda.mass; + cand.cosPALambda = lambda.cosPA; + cand.dcaPosPVLambda = lambda.dcaPosToPV; + cand.dcaNegPVLambda = lambda.dcaNegToPV; + cand.dcaLambdaDaughters = lambda.dcaDaughters; + cand.nSigmaPrLambda = lambda.nSigmaPr; + cand.nSigmaPiLambda = lambda.nSigmaPi; + + cand.ptKaon = kaon.pt(); + cand.etaKaon = kaon.eta(); + cand.phiKaon = kaon.phi(); + cand.chargeKaon = kaon.sign(); + cand.dcaXYKaon = kaon.dcaXY(); + cand.dcaZKaon = kaon.dcaZ(); + cand.nSigmaKKaon = kaon.tpcNSigmaKa(); + + const std::array momentum{ + omega.px() + lambda.px + kaon.px(), + omega.py() + lambda.py + kaon.py(), + omega.pz() + lambda.pz + kaon.pz()}; + const float energyOmega = std::sqrt(o2::constants::physics::MassOmegaMinus * o2::constants::physics::MassOmegaMinus + omega.p() * omega.p()); + const float momentumLambda2 = lambda.px * lambda.px + lambda.py * lambda.py + lambda.pz * lambda.pz; + const float energyLambda = std::sqrt(o2::constants::physics::MassLambda0 * o2::constants::physics::MassLambda0 + momentumLambda2); + const float energyKaon = std::sqrt(o2::constants::physics::MassKaonCharged * o2::constants::physics::MassKaonCharged + kaon.p() * kaon.p()); + const float energy = energyOmega + energyLambda + energyKaon; + const float mass2 = energy * energy - momentum[0] * momentum[0] - momentum[1] * momentum[1] - momentum[2] * momentum[2]; + cand.doubleOmegaMass = std::sqrt(std::max(0.f, mass2)); + return cand; + } + + template + bool getMCInfo(DoubleOmegaMCInfo& mcInfo, Omega const& omega, V0 const& v0, Kaon const& kaon, + CascLabels const& cascLabels, V0Labels const& v0Labels, TrackLabels const& trackLabels) + { + if (omega.globalIndex() >= cascLabels.size() || + v0.globalIndex() >= v0Labels.size() || + kaon.globalIndex() >= trackLabels.size()) { + LOG(info) << "Skipping MC info retrieval for double Omega candidate with invalid indices: " + << "omega index = " << omega.globalIndex() << ", v0 index = " << v0.globalIndex() + << ", kaon index = " << kaon.globalIndex(); + return false; + } + + auto omegaLabel = cascLabels.rawIteratorAt(omega.globalIndex()); + auto v0Label = v0Labels.rawIteratorAt(v0.globalIndex()); + auto kaonLabel = trackLabels.rawIteratorAt(kaon.globalIndex()); + if (!omegaLabel.has_mcParticle() || !v0Label.has_mcParticle() || !kaonLabel.has_mcParticle()) { + return false; + } + + auto mcOmega = omegaLabel.template mcParticle_as(); + auto mcLambda = v0Label.template mcParticle_as(); + auto mcKaon = kaonLabel.template mcParticle_as(); + const bool isMatter = omega.sign() < 0; + if (mcOmega.pdgCode() != (isMatter ? 3334 : -3334) || + mcLambda.pdgCode() != (isMatter ? 3122 : -3122) || + mcKaon.pdgCode() != (isMatter ? -321 : 321)) { + return false; + } + for (const auto& omegaMother : mcOmega.template mothers_as()) { + if (omegaMother.pdgCode() != (isMatter ? kDoubleOmegaPdg : -kDoubleOmegaPdg)) { + continue; + } + for (const auto& lambdaMother : mcLambda.template mothers_as()) { + if (omegaMother.globalIndex() != lambdaMother.globalIndex()) { + continue; + } + for (const auto& kaonMother : mcKaon.template mothers_as()) { + if (omegaMother.globalIndex() != kaonMother.globalIndex()) { + continue; + } + mcInfo.motherId = omegaMother.globalIndex(); + mcInfo.pt = omegaMother.pt(); + mcInfo.eta = omegaMother.eta(); + mcInfo.phi = omegaMother.phi(); + mcInfo.decayLength = std::hypot(mcOmega.vx() - omegaMother.vx(), + mcOmega.vy() - omegaMother.vy(), + mcOmega.vz() - omegaMother.vz()); + mcInfo.pdgCode = omegaMother.pdgCode(); + return true; + } + } + } + return false; + } + + template + bool getGeneratedMCInfo(DoubleOmegaMCInfo& info, McParticle const& particle) + { + if (std::abs(particle.pdgCode()) != kDoubleOmegaPdg) { + return false; + } + + const bool isMatter = particle.pdgCode() > 0; + const int expectedOmegaPdg = isMatter ? 3334 : -3334; + const int expectedLambdaPdg = isMatter ? 3122 : -3122; + const int expectedKaonPdg = isMatter ? -321 : 321; + bool foundOmega = false; + bool foundLambda = false; + bool foundKaon = false; + std::array decayVertex{0.f, 0.f, 0.f}; + + for (const auto& daughter : particle.template daughters_as()) { + if (daughter.pdgCode() == expectedOmegaPdg) { + foundOmega = true; + decayVertex = {daughter.vx(), daughter.vy(), daughter.vz()}; + } else if (daughter.pdgCode() == expectedLambdaPdg) { + foundLambda = true; + } else if (daughter.pdgCode() == expectedKaonPdg) { + foundKaon = true; + } + } + if (!foundOmega || !foundLambda || !foundKaon) { + return false; + } + + info.motherId = particle.globalIndex(); + info.pt = particle.pt(); + info.eta = particle.eta(); + info.phi = particle.phi(); + info.decayLength = std::hypot(decayVertex[0] - particle.vx(), + decayVertex[1] - particle.vy(), + decayVertex[2] - particle.vz()); + info.pdgCode = particle.pdgCode(); + return true; + } + + template + void fillFindableDaughters(Candidates const& candidates, Labels const& labels, int expectedAbsPdg, FindableDaughter daughterType) + { + for (const auto& candidate : candidates) { + + if (expectedAbsPdg == 3334) { + LOG(info) << "Checking candidate with global index " << candidate.globalIndex() << " for findable daughter of type " << daughterType; + LOG(info) << "Total number of labels: " << labels.size(); + } + + if (candidate.globalIndex() >= labels.size()) { + continue; + } + auto label = labels.rawIteratorAt(candidate.globalIndex()); + if (!label.has_mcParticle()) { + continue; + } + auto mcParticle = label.template mcParticle_as(); + if (std::abs(mcParticle.pdgCode()) != expectedAbsPdg) { + continue; + } + if(expectedAbsPdg == 3334){ + LOG(info) << "Found a cascade with pdg code " << mcParticle.pdgCode() << " and global index " << mcParticle.globalIndex(); + for (const auto& mother : mcParticle.template mothers_as()) { + LOG(info) << "Mother pdg code: " << mother.pdgCode() << ", global index: " << mother.globalIndex(); + } + } + for (const auto& mother : mcParticle.template mothers_as()) { + if (std::abs(mother.pdgCode()) == kDoubleOmegaPdg) { + histos.fill(HIST("MC/findableDaughters"), daughterType); + break; + } + } + } + } + + template + void checkFindables(Cascades const& cascades, V0s const& v0s, Tracks const& tracks, + CascLabels const& cascLabels, V0Labels const& v0Labels, TrackLabels const& trackLabels) + { + fillFindableDaughters(cascades, cascLabels, 3334, kOmega); + fillFindableDaughters(v0s, v0Labels, 3122, kLambda); + fillFindableDaughters(tracks, trackLabels, 321, kKaon); + } + + void writeDataCandidate(DoubleOmegaCandidate const& cand) + { + doubleOmegaTable(cand.ptCasc1, + cand.etaCasc1, + cand.phiCasc1, + cand.cascDecLength1, + cand.omegaMassCasc1, + cand.xiMassCasc1, + cand.cosPACasc1, + cand.dcaBachPVCasc1, + cand.dcaV0BachCasc1, + cand.nSigmaKBach1, + cand.ptLambda, + cand.etaLambda, + cand.phiLambda, + cand.lambdaDecLength, + cand.lambdaMass, + cand.cosPALambda, + cand.dcaPosPVLambda, + cand.dcaNegPVLambda, + cand.dcaLambdaDaughters, + cand.nSigmaPrLambda, + cand.nSigmaPiLambda, + cand.ptKaon, + cand.etaKaon, + cand.phiKaon, + cand.chargeKaon, + cand.dcaXYKaon, + cand.dcaZKaon, + cand.nSigmaKKaon, + cand.doubleOmegaMass); + } + + void writeMCCandidate(DoubleOmegaCandidate const& cand, DoubleOmegaMCInfo const& mcInfo, bool isReco) + { + doubleOmegaTableMC(cand.ptCasc1, + cand.etaCasc1, + cand.phiCasc1, + cand.cascDecLength1, + cand.omegaMassCasc1, + cand.xiMassCasc1, + cand.cosPACasc1, + cand.dcaBachPVCasc1, + cand.dcaV0BachCasc1, + cand.nSigmaKBach1, + cand.ptLambda, + cand.etaLambda, + cand.phiLambda, + cand.lambdaDecLength, + cand.lambdaMass, + cand.cosPALambda, + cand.dcaPosPVLambda, + cand.dcaNegPVLambda, + cand.dcaLambdaDaughters, + cand.nSigmaPrLambda, + cand.nSigmaPiLambda, + cand.ptKaon, + cand.etaKaon, + cand.phiKaon, + cand.chargeKaon, + cand.dcaXYKaon, + cand.dcaZKaon, + cand.nSigmaKKaon, + cand.doubleOmegaMass, + mcInfo.pt, + mcInfo.eta, + mcInfo.phi, + mcInfo.decayLength, + mcInfo.pdgCode, + isReco); + } + + template + void fillFromV0s(C const& collision, T const& tracks, Cascades const& cascades, V0s const& v0s, McLabels const&... mcLabels) + { + for (const auto& omega : cascades) { + LOG(info) << "Processing omega with global index " << omega.globalIndex(); + if (!isSelectedOmega(collision, tracks, omega)) { + continue; + } + histos.fill(HIST("QA/massXi"), omega.pt(), omega.mXi()); + histos.fill(HIST("QA/massOmega"), omega.pt(), omega.mOmega()); + const std::array omegaTrackIds{omega.posTrackId(), omega.negTrackId(), omega.bachelorId()}; + for (const auto& v0 : v0s) { + const int8_t charge = omega.sign(); + if (!isSelectedLambda(collision, tracks, v0, charge)) { + continue; + } + if (std::find(omegaTrackIds.begin(), omegaTrackIds.end(), v0.posTrackId()) != omegaTrackIds.end() || + std::find(omegaTrackIds.begin(), omegaTrackIds.end(), v0.negTrackId()) != omegaTrackIds.end()) { + continue; + } + const auto lambda = getLambdaCandidateFromV0(collision, tracks, v0, charge); + histos.fill(HIST("QA/massLambda"), std::hypot(lambda.px, lambda.py), lambda.mass); + for (const auto& kaon : tracks) { + if (!isSelectedKaon(kaon, charge) || + kaon.globalIndex() == v0.posTrackId() || + kaon.globalIndex() == v0.negTrackId() || + std::find(omegaTrackIds.begin(), omegaTrackIds.end(), kaon.globalIndex()) != omegaTrackIds.end()) { + continue; + } + auto cand = makeCandidate(collision, omega, lambda, kaon); + if constexpr (isMC) { + static_assert(sizeof...(McLabels) == 3); + DoubleOmegaMCInfo mcInfo; + if (getMCInfo(mcInfo, omega, v0, kaon, mcLabels...)) { + reconstructedDoubleOmegaIds.push_back(mcInfo.motherId); + writeMCCandidate(cand, mcInfo, true); + } + } else { + writeDataCandidate(cand); + } + } + } + } + } + + template + void fillFromCascades(C const& collision, T const& tracks, Cascades const& cascades) + { + for (const auto& omega : cascades) { + if (!isSelectedOmega(collision, tracks, omega)) { + continue; + } + histos.fill(HIST("QA/massXi"), omega.pt(), omega.mXi()); + histos.fill(HIST("QA/massOmega"), omega.pt(), omega.mOmega()); + const std::array omegaTrackIds{omega.posTrackId(), omega.negTrackId(), omega.bachelorId()}; + + for (const auto& lambdaKaonSource : cascades) { + if (lambdaKaonSource.globalIndex() == omega.globalIndex() || + lambdaKaonSource.sign() != omega.sign() || + !isSelectedOmega(collision, tracks, lambdaKaonSource)) { + continue; + } + + const std::array sourceTrackIds{ + lambdaKaonSource.posTrackId(), lambdaKaonSource.negTrackId(), lambdaKaonSource.bachelorId()}; + bool sharesTrack = false; + for (const auto omegaTrackId : omegaTrackIds) { + if (std::find(sourceTrackIds.begin(), sourceTrackIds.end(), omegaTrackId) != sourceTrackIds.end()) { + sharesTrack = true; + break; + } + } + if (sharesTrack) { + continue; + } + + auto kaon = lambdaKaonSource.template bachelor_as(); + const auto lambda = getLambdaCandidateFromCascade(collision, tracks, lambdaKaonSource); + histos.fill(HIST("QA/massLambda"), std::hypot(lambda.px, lambda.py), lambda.mass); + writeDataCandidate(makeCandidate(collision, omega, lambda, kaon)); + } + } + } + + template + bool acceptCollision(C const& collision, aod::BCsWithTimestamps const&) + { + auto bc = collision.template bc_as(); + initCCDB(bc); + + if (!collision.sel8() || + std::abs(collision.posZ()) > zVtxMax || + !collision.selection_bit(aod::evsel::kNoITSROFrameBorder) || + !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { + return false; + } + if (cfgSkimmedProcessing) { + zorro.isSelected(bc.globalBC()); + } + histos.fill(HIST("QA/zVtx"), collision.posZ()); + return true; + } + + void processData(Collisions const& collision, + TracksFull const& tracks, + FullCascades const& cascades, + FullV0s const& v0s, + aod::BCsWithTimestamps const& bcs) + { + if (acceptCollision(collision, bcs)) { + fillFromV0s(collision, tracks, cascades, v0s); + } + } + PROCESS_SWITCH(doubleOmegaTreeCreator, processData, "Reconstruct Omega + Lambda(V0) + kaon(track)", true); + + void processDataFromCascades(Collisions const& collision, + TracksFull const& tracks, + FullCascades const& cascades, + aod::BCsWithTimestamps const& bcs) + { + if (acceptCollision(collision, bcs)) { + fillFromCascades(collision, tracks, cascades); + } + } + PROCESS_SWITCH(doubleOmegaTreeCreator, processDataFromCascades, "Reconstruct Omega + Lambda + kaon from two cascade rows", false); + + void processMC(CollisionsMC const& collisions, + TracksFull const& tracks, + FullCascades const& cascades, + FullV0s const& v0s, + aod::McTrackLabels const& trackLabels, + aod::McCascLabels const& cascLabels, + aod::McV0Labels const& v0Labels, + aod::McParticles const& mcParticles, + aod::BCsWithTimestamps const&) + { + reconstructedDoubleOmegaIds.clear(); + + for (const auto& collision : collisions) { + // if (!acceptCollision(collision, bcs)) { + // continue; + // } + auto tracksThisCollision = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); + auto cascadesThisCollision = cascades.sliceBy(cascadesPerCollision, collision.globalIndex()); + auto v0sThisCollision = v0s.sliceBy(v0sPerCollision, collision.globalIndex()); + cascadesThisCollision.bindExternalIndices(&tracks); + v0sThisCollision.bindExternalIndices(&tracks); + checkFindables(cascadesThisCollision, v0sThisCollision, tracksThisCollision, cascLabels, v0Labels, trackLabels); + fillFromV0s(collision, tracksThisCollision, cascadesThisCollision, v0sThisCollision, cascLabels, v0Labels, trackLabels); + } + + for (const auto& mcParticle : mcParticles) { + DoubleOmegaMCInfo mcInfo; + if (!getGeneratedMCInfo(mcInfo, mcParticle)) { + continue; + } + + if (std::find(reconstructedDoubleOmegaIds.begin(), reconstructedDoubleOmegaIds.end(), mcInfo.motherId) != reconstructedDoubleOmegaIds.end()) { + continue; + } + writeMCCandidate(DoubleOmegaCandidate{}, mcInfo, false); + } + } + PROCESS_SWITCH(doubleOmegaTreeCreator, processMC, "Reconstruct and truth-match Omega + Lambda + kaon in MC", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} From 99d656e089ac56946a10e18a6ba06497e723ad15 Mon Sep 17 00:00:00 2001 From: Francesco Mazzaschi Date: Fri, 21 Aug 2026 11:37:06 +0200 Subject: [PATCH 2/2] added findables + vertexing --- PWGLF/DataModel/LFDoubleOmegaTables.h | 99 +- .../Strangeness/doubleOmegaTreeCreator.cxx | 1149 +++++++++++------ 2 files changed, 816 insertions(+), 432 deletions(-) diff --git a/PWGLF/DataModel/LFDoubleOmegaTables.h b/PWGLF/DataModel/LFDoubleOmegaTables.h index dd14aa99004..1657f31110e 100644 --- a/PWGLF/DataModel/LFDoubleOmegaTables.h +++ b/PWGLF/DataModel/LFDoubleOmegaTables.h @@ -21,38 +21,24 @@ namespace o2::aod namespace DoubleOmegaTables { -DECLARE_SOA_COLUMN(PtCasc1, ptCasc1, float); // signed pt of the cascade -DECLARE_SOA_COLUMN(EtaCasc1, etaCasc1, float); -DECLARE_SOA_COLUMN(PhiCasc1, phiCasc1, float); -DECLARE_SOA_COLUMN(CascDecLength1, cascDecLength1, float); -DECLARE_SOA_COLUMN(OmegaMassCasc1, omegaMassCasc1, float); -DECLARE_SOA_COLUMN(XiMassCasc1, xiMassCasc1, float); -DECLARE_SOA_COLUMN(CosPACasc1, cosPACasc1, float); -DECLARE_SOA_COLUMN(DcaBachPVCasc1, dcaBachPVCasc1, float); -DECLARE_SOA_COLUMN(DcaV0BachCasc1, dcaV0BachCasc1, float); -DECLARE_SOA_COLUMN(NSigmaKBach1, nSigmaKBach1, float); - -DECLARE_SOA_COLUMN(PtLambda, ptLambda, float); -DECLARE_SOA_COLUMN(EtaLambda, etaLambda, float); -DECLARE_SOA_COLUMN(PhiLambda, phiLambda, float); -DECLARE_SOA_COLUMN(LambdaDecLength, lambdaDecLength, float); -DECLARE_SOA_COLUMN(LambdaMass, lambdaMass, float); -DECLARE_SOA_COLUMN(CosPALambda, cosPALambda, float); -DECLARE_SOA_COLUMN(DcaPosPVLambda, dcaPosPVLambda, float); -DECLARE_SOA_COLUMN(DcaNegPVLambda, dcaNegPVLambda, float); -DECLARE_SOA_COLUMN(DcaLambdaDaughters, dcaLambdaDaughters, float); -DECLARE_SOA_COLUMN(NSigmaPrLambda, nSigmaPrLambda, float); -DECLARE_SOA_COLUMN(NSigmaPiLambda, nSigmaPiLambda, float); - -DECLARE_SOA_COLUMN(PtKaon, ptKaon, float); -DECLARE_SOA_COLUMN(EtaKaon, etaKaon, float); -DECLARE_SOA_COLUMN(PhiKaon, phiKaon, float); -DECLARE_SOA_COLUMN(ChargeKaon, chargeKaon, int8_t); -DECLARE_SOA_COLUMN(DcaXYKaon, dcaXYKaon, float); -DECLARE_SOA_COLUMN(DcaZKaon, dcaZKaon, float); -DECLARE_SOA_COLUMN(NSigmaKKaon, nSigmaKKaon, float); - -DECLARE_SOA_COLUMN(DoubleOmegaMass, doubleOmegaMass, float); +DECLARE_SOA_COLUMN(PtDoubleOmega, ptDoubleOmega, float); +DECLARE_SOA_COLUMN(EtaDoubleOmega, etaDoubleOmega, float); +DECLARE_SOA_COLUMN(PhiDoubleOmega, phiDoubleOmega, float); +DECLARE_SOA_COLUMN(DecayVtxX, decayVtxX, float); +DECLARE_SOA_COLUMN(DecayVtxY, decayVtxY, float); +DECLARE_SOA_COLUMN(DecayVtxZ, decayVtxZ, float); +DECLARE_SOA_COLUMN(CosPAOmega, cosPAOmega, float); +DECLARE_SOA_COLUMN(CosPADirectLambda, cosPADirectLambda, float); +DECLARE_SOA_COLUMN(CosPADoubleOmega, cosPADoubleOmega, float); +DECLARE_SOA_COLUMN(DCAxyOmegaToPV, dcaXYOmegaToPV, float); +DECLARE_SOA_COLUMN(DCAzOmegaToPV, dcaZOmegaToPV, float); +DECLARE_SOA_COLUMN(DCAxyDirectLambdaToPV, dcaXYDirectLambdaToPV, float); +DECLARE_SOA_COLUMN(DCAzDirectLambdaToPV, dcaZDirectLambdaToPV, float); +DECLARE_SOA_COLUMN(DCAxyDirectKaonToPV, dcaXYDirectKaonToPV, float); +DECLARE_SOA_COLUMN(DCAzDirectKaonToPV, dcaZDirectKaonToPV, float); +DECLARE_SOA_COLUMN(MassDoubleOmega, massDoubleOmega, float); +DECLARE_SOA_COLUMN(MassOmega, massOmega, float); +DECLARE_SOA_COLUMN(MassXi, massXi, float); DECLARE_SOA_COLUMN(GenPt, genPt, float); DECLARE_SOA_COLUMN(GenEta, genEta, float); @@ -62,36 +48,25 @@ DECLARE_SOA_COLUMN(PdgDoubleOmega, pdgDoubleOmega, int); DECLARE_SOA_COLUMN(IsReco, isReco, bool); } // namespace DoubleOmegaTables -#define DOUBLE_OMEGA_RECO_COLUMNS \ - DoubleOmegaTables::PtCasc1, \ - DoubleOmegaTables::EtaCasc1, \ - DoubleOmegaTables::PhiCasc1, \ - DoubleOmegaTables::CascDecLength1, \ - DoubleOmegaTables::OmegaMassCasc1, \ - DoubleOmegaTables::XiMassCasc1, \ - DoubleOmegaTables::CosPACasc1, \ - DoubleOmegaTables::DcaBachPVCasc1, \ - DoubleOmegaTables::DcaV0BachCasc1, \ - DoubleOmegaTables::NSigmaKBach1, \ - DoubleOmegaTables::PtLambda, \ - DoubleOmegaTables::EtaLambda, \ - DoubleOmegaTables::PhiLambda, \ - DoubleOmegaTables::LambdaDecLength, \ - DoubleOmegaTables::LambdaMass, \ - DoubleOmegaTables::CosPALambda, \ - DoubleOmegaTables::DcaPosPVLambda, \ - DoubleOmegaTables::DcaNegPVLambda, \ - DoubleOmegaTables::DcaLambdaDaughters, \ - DoubleOmegaTables::NSigmaPrLambda, \ - DoubleOmegaTables::NSigmaPiLambda, \ - DoubleOmegaTables::PtKaon, \ - DoubleOmegaTables::EtaKaon, \ - DoubleOmegaTables::PhiKaon, \ - DoubleOmegaTables::ChargeKaon, \ - DoubleOmegaTables::DcaXYKaon, \ - DoubleOmegaTables::DcaZKaon, \ - DoubleOmegaTables::NSigmaKKaon, \ - DoubleOmegaTables::DoubleOmegaMass +#define DOUBLE_OMEGA_RECO_COLUMNS \ + DoubleOmegaTables::PtDoubleOmega, \ + DoubleOmegaTables::EtaDoubleOmega, \ + DoubleOmegaTables::PhiDoubleOmega, \ + DoubleOmegaTables::DecayVtxX, \ + DoubleOmegaTables::DecayVtxY, \ + DoubleOmegaTables::DecayVtxZ, \ + DoubleOmegaTables::CosPAOmega, \ + DoubleOmegaTables::CosPADirectLambda, \ + DoubleOmegaTables::CosPADoubleOmega, \ + DoubleOmegaTables::DCAxyOmegaToPV, \ + DoubleOmegaTables::DCAzOmegaToPV, \ + DoubleOmegaTables::DCAxyDirectLambdaToPV, \ + DoubleOmegaTables::DCAzDirectLambdaToPV, \ + DoubleOmegaTables::DCAxyDirectKaonToPV, \ + DoubleOmegaTables::DCAzDirectKaonToPV, \ + DoubleOmegaTables::MassDoubleOmega, \ + DoubleOmegaTables::MassOmega, \ + DoubleOmegaTables::MassXi DECLARE_SOA_TABLE(DoubleOmegaTable, "AOD", "DOUBLEOMEGATABLE", DOUBLE_OMEGA_RECO_COLUMNS); diff --git a/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx b/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx index a377cd2bee1..1b8a649d7cc 100644 --- a/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx +++ b/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx @@ -13,14 +13,20 @@ #include "PWGLF/DataModel/LFStrangenessTables.h" #include "Common/CCDB/EventSelectionParams.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/Zorro.h" #include "Common/Core/ZorroSummary.h" +#include "Common/Core/trackUtilities.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponseTPC.h" #include #include +#include +#include +#include +#include #include #include #include @@ -30,6 +36,7 @@ #include #include #include +#include #include #include @@ -38,6 +45,7 @@ #include #include #include +#include #include using namespace o2; @@ -48,42 +56,30 @@ using CollisionsTable = soa::Join; using FullCascades = aod::CascDataExt; -using FullV0s = aod::V0Datas; using TracksFull = soa::Join; +using TracksFullIU = soa::Join; +using TracksFullIUMC = soa::Join; +using TracksFullMC = soa::Join; struct DoubleOmegaCandidate { - float ptCasc1 = -999.f; - float etaCasc1 = -999.f; - float phiCasc1 = -999.f; - float cascDecLength1 = -999.f; - float omegaMassCasc1 = -999.f; - float xiMassCasc1 = -999.f; - float cosPACasc1 = -999.f; - float dcaBachPVCasc1 = -999.f; - float dcaV0BachCasc1 = -999.f; - float nSigmaKBach1 = -999.f; - - float ptLambda = -999.f; - float etaLambda = -999.f; - float phiLambda = -999.f; - float lambdaDecLength = -999.f; - float lambdaMass = -999.f; - float cosPALambda = -999.f; - float dcaPosPVLambda = -999.f; - float dcaNegPVLambda = -999.f; - float dcaLambdaDaughters = -999.f; - float nSigmaPrLambda = -999.f; - float nSigmaPiLambda = -999.f; - - float ptKaon = -999.f; - float etaKaon = -999.f; - float phiKaon = -999.f; - int8_t chargeKaon = 0; - float dcaXYKaon = -999.f; - float dcaZKaon = -999.f; - float nSigmaKKaon = -999.f; - - float doubleOmegaMass = -999.f; + float pt = -999.f; + float eta = -999.f; + float phi = -999.f; + float x = -999.f; + float y = -999.f; + float z = -999.f; + float cosPAOmega = -999.f; + float cosPADirectLambda = -999.f; + float cosPADoubleOmega = -999.f; + float dcaXYOmegaToPV = -999.f; + float dcaZOmegaToPV = -999.f; + float dcaXYDirectLambdaToPV = -999.f; + float dcaZDirectLambdaToPV = -999.f; + float dcaXYDirectKaonToPV = -999.f; + float dcaZDirectKaonToPV = -999.f; + float mass = -999.f; + float massOmega = -999.f; + float massXi = -999.f; }; struct DoubleOmegaMCInfo { @@ -99,14 +95,33 @@ struct LambdaCandidate { float px = 0.f; float py = 0.f; float pz = 0.f; - float decayLength = -999.f; float mass = -999.f; - float cosPA = -999.f; - float dcaPosToPV = -999.f; - float dcaNegToPV = -999.f; - float dcaDaughters = -999.f; - float nSigmaPr = -999.f; - float nSigmaPi = -999.f; +}; + +struct BuiltLambdaCandidate { + LambdaCandidate candidate; + o2::track::TrackParCov parentTrack; + std::array decayVertex{}; + int64_t v0Id = -1; + int64_t posTrackId = -1; + int64_t negTrackId = -1; +}; + +struct BuiltOmegaCandidate { + float px = 0.f; + float py = 0.f; + float pz = 0.f; + float x = 0.f; + float y = 0.f; + float z = 0.f; + float massOmega = -999.f; + float massXi = -999.f; + int8_t sign = 0; + int64_t cascadeId = -1; + int64_t posTrackId = -1; + int64_t negTrackId = -1; + int64_t bachelorId = -1; + o2::track::TrackParCov parentTrack; }; struct doubleOmegaTreeCreator { @@ -115,25 +130,33 @@ struct doubleOmegaTreeCreator { Service ccdb; static constexpr int kDoubleOmegaPdg = 1060020020; - enum FindableDaughter : uint8_t { - kOmega, - kLambda, - kKaon, - kNFindableDaughters + enum FindabilityStep : uint8_t { + kAllGenerated, + kOmegaLambdaKaon, + kTwoLambdasTwoKaons, + kFinalState, + kFindable, + kFindableSelectedTracks, + kFindableCascadeAndV0, + kNFindabilitySteps }; std::vector reconstructedDoubleOmegaIds; - Preslice tracksPerCollision = aod::track::collisionId; - Preslice cascadesPerCollision = aod::cascade::collisionId; - Preslice v0sPerCollision = aod::v0data::collisionId; + Preslice tracksIUPerCollision = aod::track::collisionId; + Preslice rawCascadesPerCollision = aod::cascade::collisionId; + Preslice rawV0sPerCollision = aod::v0::collisionId; int mRunNumber = 0; + float mBz = 0.f; + o2::vertexing::DCAFitterN<2> fitter; + o2::vertexing::DCAFitterN<3> fitter3Body; Zorro zorro; OutputObj zorroSummary{"zorroSummary"}; Configurable cfgSkimmedProcessing{"cfgSkimmedProcessing", false, "Skimmed dataset processing"}; Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable cfgMaterialCorrection{"cfgMaterialCorrection", static_cast(o2::base::Propagator::MatCorrType::USEMatCorrLUT), "Material correction for the raw V0/cascade fits"}; ConfigurableAxis zVtxAxis{"zVtxBins", {100, -20.f, 20.f}, "Binning for the vertex z in cm"}; ConfigurableAxis massOmegaAxis{"massOmegaAxis", {400, o2::constants::physics::MassOmegaMinus - 0.05f, o2::constants::physics::MassOmegaMinus + 0.05f}, "binning for the Omega invariant-mass"}; @@ -151,17 +174,21 @@ struct doubleOmegaTreeCreator { Configurable minNTPCClus{"minNTPCClus", 80, "Minimum number of TPC clusters"}; Configurable maxNSharedTPCClus{"maxNSharedTPCClus", 5, "Maximum number of shared TPC clusters"}; - Configurable minCascCosPA{"minCascCosPA", 0.99f, "Minimum cosine of the pointing angle of the cascade"}; - Configurable minLambdaCosPA{"minLambdaCosPA", 0.995f, "Minimum cosine of the pointing angle of the Lambda"}; Configurable nSigmaTPCCut{"nSigmaTPCCut", 3.f, "Number of sigmas for the TPC PID"}; Configurable dcaBachToPV{"dcaBachToPV", 0.05f, "Minimum DCA of a cascade bachelor to the primary vertex"}; - Configurable dcaKaonToPV{"dcaKaonToPV", 0.05f, "Minimum transverse DCA of the kaon to the primary vertex"}; + Configurable dcaKaonToPV{"dcaKaonToPV", 0.05f, "Minimum absolute transverse DCA of the direct kaon to the primary vertex"}; + Configurable dcaOmegaToPV{"dcaOmegaToPV", 0.f, "Minimum absolute transverse DCA of the Omega to the primary vertex"}; + Configurable dcaDirectLambdaToPV{"dcaDirectLambdaToPV", 0.f, "Minimum absolute transverse DCA of the direct Lambda to the primary vertex"}; Configurable dcaV0DauToPV{"dcaV0DauToPV", 0.05f, "Minimum DCA of Lambda daughters to the primary vertex"}; Configurable dcaV0Bach{"dcaV0Bach", 1.f, "Maximum DCA between the V0 and cascade bachelor"}; Configurable dcaLambdaDaughters{"dcaLambdaDaughters", 1.f, "Maximum DCA between Lambda daughters"}; Configurable mXiWindow{"mXiWindow", 0.02f, "Xi mass window used by the cascade compatibility mode"}; Configurable mOmegaWindow{"mOmegaWindow", 0.01f, "Omega mass window"}; Configurable mLambdaWindow{"mLambdaWindow", 0.01f, "Lambda mass window"}; + Configurable minCosPAOmega{"minCosPAOmega", -1.f, "Minimum Omega cosPA relative to the double-Omega decay vertex"}; + Configurable minCosPADirectLambda{"minCosPADirectLambda", -1.f, "Minimum direct-Lambda cosPA relative to the double-Omega decay vertex"}; + Configurable minCosPADoubleOmega{"minCosPADoubleOmega", -1.f, "Minimum double-Omega cosPA relative to the primary vertex"}; + Configurable minDoubleOmegaDecayRadius{"minDoubleOmegaDecayRadius", 1.f, "Minimum double-Omega transverse decay radius in cm"}; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -187,7 +214,26 @@ struct doubleOmegaTreeCreator { if (mRunNumber == bc.runNumber()) { return; } - LOG(info) << "Retrieved GRP for timestamp " << bc.timestamp(); + auto* grpmag = ccdb->getForTimeStamp("GLO/Config/GRPMagField", bc.timestamp()); + if (!grpmag) { + LOG(fatal) << "Could not retrieve GRPMagField for timestamp " << bc.timestamp(); + } + o2::base::Propagator::initFieldFromGRP(grpmag); + mBz = o2::base::Propagator::Instance()->getNominalBz(); + fitter.setBz(mBz); + fitter3Body.setBz(mBz); + + if (static_cast(cfgMaterialCorrection.value) == o2::base::Propagator::MatCorrType::USEMatCorrLUT) { + auto* lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->getForTimeStamp("GLO/Param/MatLUT", bc.timestamp())); + if (!lut) { + LOG(fatal) << "Could not retrieve material LUT for timestamp " << bc.timestamp(); + } + o2::base::Propagator::Instance()->setMatLUT(lut); + } + fitter.setMatCorrType(static_cast(cfgMaterialCorrection.value)); + fitter3Body.setMatCorrType(static_cast(cfgMaterialCorrection.value)); + + LOG(info) << "Retrieved GRP for timestamp " << bc.timestamp() << " with magnetic field " << mBz << " kG"; mRunNumber = bc.runNumber(); if (cfgSkimmedProcessing) { zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), "fDoubleOmega,fOmegaXi"); @@ -202,119 +248,262 @@ struct doubleOmegaTreeCreator { ccdb->setLocalObjectValidityChecking(); ccdb->setFatalWhenNull(false); + fitter.setPropagateToPCA(true); + fitter.setMaxR(200.); + fitter.setMaxDZIni(4.); + fitter.setMinParamChange(1.e-3); + fitter.setMinRelChi2Change(0.9); + fitter.setUseAbsDCA(true); + + fitter3Body.setPropagateToPCA(true); + fitter3Body.setMaxR(200.); + fitter3Body.setMaxDZIni(4.); + fitter3Body.setMinParamChange(1.e-3); + fitter3Body.setMinRelChi2Change(0.9); + fitter3Body.setUseAbsDCA(true); + zorroSummary.setObject(zorro.getZorroSummary()); histos.add("QA/zVtx", ";#it{z}_{vtx} (cm);Entries", HistType::kTH1F, {zVtxAxis}); histos.add("QA/massXi", ";#it{p}_{T} (GeV/#it{c});#it{M}(#Lambda + #pi) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massXiAxis}); histos.add("QA/massOmega", ";#it{p}_{T} (GeV/#it{c});#it{M}(#Lambda + K) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massOmegaAxis}); histos.add("QA/massLambda", ";#it{p}_{T} (GeV/#it{c});#it{M}(p + #pi) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massLambdaAxis}); - histos.add("MC/findableDaughters", ";Findable daughter;Entries", HistType::kTH1F, {{kNFindableDaughters, -0.5f, static_cast(kNFindableDaughters) - 0.5f}}); - auto findableDaughters = histos.get(HIST("MC/findableDaughters")); - findableDaughters->GetXaxis()->SetBinLabel(kOmega + 1, "#Omega"); - findableDaughters->GetXaxis()->SetBinLabel(kLambda + 1, "#Lambda"); - findableDaughters->GetXaxis()->SetBinLabel(kKaon + 1, "K"); + histos.add("MC/generatedAndFindable", ";Double-#Omega;Entries", HistType::kTH1F, {{kNFindabilitySteps, -0.5f, static_cast(kNFindabilitySteps) - 0.5f}}); + auto generatedAndFindable = histos.get(HIST("MC/generatedAndFindable")); + generatedAndFindable->GetXaxis()->SetBinLabel(kAllGenerated + 1, "All generated"); + generatedAndFindable->GetXaxis()->SetBinLabel(kOmegaLambdaKaon + 1, "#Omega + #Lambda + K"); + generatedAndFindable->GetXaxis()->SetBinLabel(kTwoLambdasTwoKaons + 1, "2#Lambda + 2K"); + generatedAndFindable->GetXaxis()->SetBinLabel(kFinalState + 1, "2p + 2#pi^{-} + 2K^{-} (and c.c.)"); + generatedAndFindable->GetXaxis()->SetBinLabel(kFindable + 1, "Findable"); + generatedAndFindable->GetXaxis()->SetBinLabel(kFindableSelectedTracks + 1, "Findable, selected tracks"); + generatedAndFindable->GetXaxis()->SetBinLabel(kFindableCascadeAndV0 + 1, "Findable cascade + V0"); + } + + template + bool isSelectedKaon(Track const& track, int8_t charge) + { + return track.sign() == charge && + selectTrack(track) && + std::abs(track.tpcNSigmaKa()) <= nSigmaTPCCut; } - template - bool isSelectedOmega(C const& collision, T const&, Casc const& casc) + static float invariantMass2Body(std::array const& momentum1, float mass1, + std::array const& momentum2, float mass2) { - auto bachelor = casc.template bachelor_as(); - auto posDau = casc.template posTrack_as(); - auto negDau = casc.template negTrack_as(); + const float momentum1Squared = momentum1[0] * momentum1[0] + momentum1[1] * momentum1[1] + momentum1[2] * momentum1[2]; + const float momentum2Squared = momentum2[0] * momentum2[0] + momentum2[1] * momentum2[1] + momentum2[2] * momentum2[2]; + const float energy = std::sqrt(momentum1Squared + mass1 * mass1) + std::sqrt(momentum2Squared + mass2 * mass2); + const std::array totalMomentum{ + momentum1[0] + momentum2[0], momentum1[1] + momentum2[1], momentum1[2] + momentum2[2]}; + const float massSquared = energy * energy - totalMomentum[0] * totalMomentum[0] - totalMomentum[1] * totalMomentum[1] - totalMomentum[2] * totalMomentum[2]; + return std::sqrt(std::max(0.f, massSquared)); + } - if (!selectTrack(bachelor) || !selectTrack(posDau) || !selectTrack(negDau)) { + template + bool buildLambda(T const&, V0 const& v0, int8_t charge, BuiltLambdaCandidate& builtLambda) + { + auto posTrack = v0.template posTrack_as(); + auto negTrack = v0.template negTrack_as(); + if (!selectTrack(posTrack) || !selectTrack(negTrack)) { return false; } - if (casc.sign() > 0) { - if (std::abs(posDau.tpcNSigmaPi()) > nSigmaTPCCut || std::abs(negDau.tpcNSigmaPr()) > nSigmaTPCCut) { - return false; - } - } else if (casc.sign() < 0) { - if (std::abs(negDau.tpcNSigmaPi()) > nSigmaTPCCut || std::abs(posDau.tpcNSigmaPr()) > nSigmaTPCCut) { - return false; - } - } else { + + const bool isMatter = charge < 0; + const float protonNSigma = isMatter ? posTrack.tpcNSigmaPr() : negTrack.tpcNSigmaPr(); + const float pionNSigma = isMatter ? negTrack.tpcNSigmaPi() : posTrack.tpcNSigmaPi(); + if (std::abs(protonNSigma) > nSigmaTPCCut || std::abs(pionNSigma) > nSigmaTPCCut || + std::abs(posTrack.dcaXY()) < dcaV0DauToPV || std::abs(negTrack.dcaXY()) < dcaV0DauToPV) { + return false; + } + + auto posTrackParCov = getTrackParCov(posTrack); + auto negTrackParCov = getTrackParCov(negTrack); + int nCandidates = 0; + try { + nCandidates = fitter.process(posTrackParCov, negTrackParCov); + } catch (...) { + LOG(error) << "Exception while fitting raw V0 " << v0.globalIndex(); return false; } - if (std::abs(bachelor.tpcNSigmaKa()) > nSigmaTPCCut || - std::abs(casc.dcabachtopv()) < dcaBachToPV || - std::abs(casc.dcacascdaughters()) > dcaV0Bach || - casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()) < minCascCosPA || - std::abs(casc.eta()) > etaMax || - std::abs(casc.mOmega() - o2::constants::physics::MassOmegaMinus) > mOmegaWindow) { + if (nCandidates == 0) { return false; } + + std::array posMomentum{}; + std::array negMomentum{}; + fitter.getTrack(0).getPxPyPzGlo(posMomentum); + fitter.getTrack(1).getPxPyPzGlo(negMomentum); + const std::array lambdaMomentum{ + posMomentum[0] + negMomentum[0], + posMomentum[1] + negMomentum[1], + posMomentum[2] + negMomentum[2]}; + const auto& fittedVertex = fitter.getPCACandidate(); + const std::array decayVertex{ + static_cast(fittedVertex[0]), static_cast(fittedVertex[1]), static_cast(fittedVertex[2])}; + const float mass = isMatter ? invariantMass2Body(posMomentum, o2::constants::physics::MassProton, + negMomentum, o2::constants::physics::MassPionCharged) + : invariantMass2Body(posMomentum, o2::constants::physics::MassPionCharged, + negMomentum, o2::constants::physics::MassProton); + const float dcaDaughters = std::sqrt(std::abs(fitter.getChi2AtPCACandidate())); + const float eta = etaFromMomentum(lambdaMomentum[0], lambdaMomentum[1], lambdaMomentum[2]); + if (std::abs(eta) > etaMax || + dcaDaughters > dcaLambdaDaughters || + std::abs(mass - o2::constants::physics::MassLambda0) > mLambdaWindow) { + return false; + } + + builtLambda.candidate = { + lambdaMomentum[0], + lambdaMomentum[1], + lambdaMomentum[2], + mass}; + builtLambda.parentTrack = fitter.createParentTrackParCov(0); + builtLambda.decayVertex = decayVertex; + builtLambda.v0Id = v0.globalIndex(); + builtLambda.posTrackId = v0.posTrackId(); + builtLambda.negTrackId = v0.negTrackId(); return true; } - template - bool isSelectedLambda(C const&, T const&, V0 const& v0, int8_t charge) + template + bool buildOmega(T const& tracks, Casc const& cascade, BuiltOmegaCandidate& builtOmega) { - auto posDau = v0.template posTrack_as(); - auto negDau = v0.template negTrack_as(); - if (!selectTrack(posDau) || !selectTrack(negDau) || - std::abs(v0.eta()) > etaMax || - v0.v0cosPA() < minLambdaCosPA || - std::abs(v0.dcapostopv()) < dcaV0DauToPV || - std::abs(v0.dcanegtopv()) < dcaV0DauToPV || - std::abs(v0.dcaV0daughters()) > dcaLambdaDaughters) { + auto bachelor = cascade.template bachelor_as(); + if (bachelor.sign() == 0 || + !selectTrack(bachelor) || + std::abs(bachelor.tpcNSigmaKa()) > nSigmaTPCCut || + std::abs(bachelor.dcaXY()) < dcaBachToPV) { return false; } - const bool isMatter = charge < 0; - const float mass = isMatter ? v0.mLambda() : v0.mAntiLambda(); - const float protonNSigma = isMatter ? posDau.tpcNSigmaPr() : negDau.tpcNSigmaPr(); - const float pionNSigma = isMatter ? negDau.tpcNSigmaPi() : posDau.tpcNSigmaPi(); - return std::abs(mass - o2::constants::physics::MassLambda0) <= mLambdaWindow && - std::abs(protonNSigma) <= nSigmaTPCCut && - std::abs(pionNSigma) <= nSigmaTPCCut; - } + auto v0 = cascade.template v0_as(); + BuiltLambdaCandidate lambda; + if (!buildLambda(tracks, v0, bachelor.sign(), lambda)) { + return false; + } - template - bool isSelectedKaon(Track const& track, int8_t charge) - { - return track.sign() == charge && - selectTrack(track) && - std::abs(track.tpcNSigmaKa()) <= nSigmaTPCCut && - std::abs(track.dcaXY()) >= dcaKaonToPV; - } + auto v0TrackParCov = lambda.parentTrack; + auto bachelorTrackParCov = getTrackParCov(bachelor); + int nCandidates = 0; + try { + nCandidates = fitter.process(v0TrackParCov, bachelorTrackParCov); + } catch (...) { + LOG(error) << "Exception while fitting raw cascade " << cascade.globalIndex(); + return false; + } + if (nCandidates == 0) { + return false; + } - template - LambdaCandidate getLambdaCandidateFromV0(C const& collision, T const&, V0 const& v0, int8_t charge) - { - auto posDau = v0.template posTrack_as(); - auto negDau = v0.template negTrack_as(); - const bool isMatter = charge < 0; - return {v0.px(), - v0.py(), - v0.pz(), - std::hypot(v0.x() - collision.posX(), v0.y() - collision.posY(), v0.z() - collision.posZ()), - isMatter ? v0.mLambda() : v0.mAntiLambda(), - v0.v0cosPA(), - v0.dcapostopv(), - v0.dcanegtopv(), - v0.dcaV0daughters(), - isMatter ? posDau.tpcNSigmaPr() : negDau.tpcNSigmaPr(), - isMatter ? negDau.tpcNSigmaPi() : posDau.tpcNSigmaPi()}; + std::array lambdaMomentum{}; + std::array bachelorMomentum{}; + fitter.getTrack(0).getPxPyPzGlo(lambdaMomentum); + fitter.getTrack(1).getPxPyPzGlo(bachelorMomentum); + const std::array omegaMomentum{ + lambdaMomentum[0] + bachelorMomentum[0], + lambdaMomentum[1] + bachelorMomentum[1], + lambdaMomentum[2] + bachelorMomentum[2]}; + const auto& fittedVertex = fitter.getPCACandidate(); + const std::array decayVertex{ + static_cast(fittedVertex[0]), static_cast(fittedVertex[1]), static_cast(fittedVertex[2])}; + const float massOmega = invariantMass2Body(lambdaMomentum, o2::constants::physics::MassLambda0, + bachelorMomentum, o2::constants::physics::MassKaonCharged); + const float massXi = invariantMass2Body(lambdaMomentum, o2::constants::physics::MassLambda0, + bachelorMomentum, o2::constants::physics::MassPionCharged); + const float dcaDaughters = std::sqrt(std::abs(fitter.getChi2AtPCACandidate())); + const float eta = etaFromMomentum(omegaMomentum[0], omegaMomentum[1], omegaMomentum[2]); + if (dcaDaughters > dcaV0Bach || + std::abs(eta) > etaMax || + std::abs(massOmega - o2::constants::physics::MassOmegaMinus) > mOmegaWindow) { + return false; + } + + builtOmega = { + omegaMomentum[0], + omegaMomentum[1], + omegaMomentum[2], + decayVertex[0], + decayVertex[1], + decayVertex[2], + massOmega, + massXi, + static_cast(bachelor.sign()), + cascade.globalIndex(), + v0.posTrackId(), + v0.negTrackId(), + cascade.bachelorId(), + fitter.createParentTrackParCov()}; + return true; } - template - LambdaCandidate getLambdaCandidateFromCascade(C const& collision, T const&, Casc const& casc) + template + bool buildOmegaFromBuilder(T const& tracks, Casc const& cascade, BuiltOmegaCandidate& builtOmega) { - auto posDau = casc.template posTrack_as(); - auto negDau = casc.template negTrack_as(); - const bool isMatter = casc.sign() < 0; - return {casc.pxlambda(), - casc.pylambda(), - casc.pzlambda(), - std::hypot(casc.xlambda() - collision.posX(), casc.ylambda() - collision.posY(), casc.zlambda() - collision.posZ()), - casc.mLambda(), - casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()), - casc.dcapostopv(), - casc.dcanegtopv(), - casc.dcaV0daughters(), - isMatter ? posDau.tpcNSigmaPr() : negDau.tpcNSigmaPr(), - isMatter ? negDau.tpcNSigmaPi() : posDau.tpcNSigmaPi()}; + auto bachelor = cascade.template bachelor_as(); + if (bachelor.sign() == 0 || + !selectTrack(bachelor) || + std::abs(bachelor.tpcNSigmaKa()) > nSigmaTPCCut || + std::abs(bachelor.dcaXY()) < dcaBachToPV) { + return false; + } + + BuiltLambdaCandidate lambda; + if (!buildLambda(tracks, cascade, bachelor.sign(), lambda)) { + return false; + } + + auto v0TrackParCov = lambda.parentTrack; + auto bachelorTrackParCov = getTrackParCov(bachelor); + int nCandidates = 0; + try { + nCandidates = fitter.process(v0TrackParCov, bachelorTrackParCov); + } catch (...) { + LOG(error) << "Exception while refitting cascade " << cascade.globalIndex(); + return false; + } + if (nCandidates == 0) { + return false; + } + + std::array lambdaMomentum{}; + std::array bachelorMomentum{}; + fitter.getTrack(0).getPxPyPzGlo(lambdaMomentum); + fitter.getTrack(1).getPxPyPzGlo(bachelorMomentum); + const std::array omegaMomentum{ + lambdaMomentum[0] + bachelorMomentum[0], + lambdaMomentum[1] + bachelorMomentum[1], + lambdaMomentum[2] + bachelorMomentum[2]}; + const auto& fittedVertex = fitter.getPCACandidate(); + const std::array decayVertex{ + static_cast(fittedVertex[0]), static_cast(fittedVertex[1]), static_cast(fittedVertex[2])}; + const float massOmega = invariantMass2Body(lambdaMomentum, o2::constants::physics::MassLambda0, + bachelorMomentum, o2::constants::physics::MassKaonCharged); + const float massXi = invariantMass2Body(lambdaMomentum, o2::constants::physics::MassLambda0, + bachelorMomentum, o2::constants::physics::MassPionCharged); + const float dcaDaughters = std::sqrt(std::abs(fitter.getChi2AtPCACandidate())); + const float eta = etaFromMomentum(omegaMomentum[0], omegaMomentum[1], omegaMomentum[2]); + if (dcaDaughters > dcaV0Bach || + std::abs(eta) > etaMax || + std::abs(massOmega - o2::constants::physics::MassOmegaMinus) > mOmegaWindow) { + return false; + } + + builtOmega = { + omegaMomentum[0], + omegaMomentum[1], + omegaMomentum[2], + decayVertex[0], + decayVertex[1], + decayVertex[2], + massOmega, + massXi, + static_cast(bachelor.sign()), + cascade.globalIndex(), + cascade.posTrackId(), + cascade.negTrackId(), + cascade.bachelorId(), + fitter.createParentTrackParCov()}; + return true; } static float etaFromMomentum(float px, float py, float pz) @@ -323,88 +512,174 @@ struct doubleOmegaTreeCreator { return pt > 0.f ? std::asinh(pz / pt) : 0.f; } - template - DoubleOmegaCandidate makeCandidate(C const& collision, Omega const& omega, LambdaCandidate const& lambda, Kaon const& kaon) + template + bool buildDoubleOmega(C const& collision, BuiltOmegaCandidate const& omega, BuiltLambdaCandidate const& directLambda, Kaon const& kaon, DoubleOmegaCandidate& cand) { - DoubleOmegaCandidate cand; - auto omegaBach = omega.template bachelor_as(); - - cand.ptCasc1 = omega.pt(); - cand.etaCasc1 = omega.eta(); - cand.phiCasc1 = omega.phi(); - cand.cascDecLength1 = std::hypot(omega.x() - collision.posX(), omega.y() - collision.posY(), omega.z() - collision.posZ()); - cand.omegaMassCasc1 = omega.mOmega(); - cand.xiMassCasc1 = omega.mXi(); - cand.cosPACasc1 = omega.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); - cand.dcaBachPVCasc1 = omega.dcabachtopv(); - cand.dcaV0BachCasc1 = omega.dcacascdaughters(); - cand.nSigmaKBach1 = omegaBach.tpcNSigmaKa(); - - cand.ptLambda = std::hypot(lambda.px, lambda.py); - cand.etaLambda = etaFromMomentum(lambda.px, lambda.py, lambda.pz); - cand.phiLambda = std::atan2(lambda.py, lambda.px); - cand.lambdaDecLength = lambda.decayLength; - cand.lambdaMass = lambda.mass; - cand.cosPALambda = lambda.cosPA; - cand.dcaPosPVLambda = lambda.dcaPosToPV; - cand.dcaNegPVLambda = lambda.dcaNegToPV; - cand.dcaLambdaDaughters = lambda.dcaDaughters; - cand.nSigmaPrLambda = lambda.nSigmaPr; - cand.nSigmaPiLambda = lambda.nSigmaPi; - - cand.ptKaon = kaon.pt(); - cand.etaKaon = kaon.eta(); - cand.phiKaon = kaon.phi(); - cand.chargeKaon = kaon.sign(); - cand.dcaXYKaon = kaon.dcaXY(); - cand.dcaZKaon = kaon.dcaZ(); - cand.nSigmaKKaon = kaon.tpcNSigmaKa(); - - const std::array momentum{ - omega.px() + lambda.px + kaon.px(), - omega.py() + lambda.py + kaon.py(), - omega.pz() + lambda.pz + kaon.pz()}; - const float energyOmega = std::sqrt(o2::constants::physics::MassOmegaMinus * o2::constants::physics::MassOmegaMinus + omega.p() * omega.p()); - const float momentumLambda2 = lambda.px * lambda.px + lambda.py * lambda.py + lambda.pz * lambda.pz; - const float energyLambda = std::sqrt(o2::constants::physics::MassLambda0 * o2::constants::physics::MassLambda0 + momentumLambda2); - const float energyKaon = std::sqrt(o2::constants::physics::MassKaonCharged * o2::constants::physics::MassKaonCharged + kaon.p() * kaon.p()); - const float energy = energyOmega + energyLambda + energyKaon; - const float mass2 = energy * energy - momentum[0] * momentum[0] - momentum[1] * momentum[1] - momentum[2] * momentum[2]; - cand.doubleOmegaMass = std::sqrt(std::max(0.f, mass2)); - return cand; + auto omegaTrackParCov = omega.parentTrack; + auto lambdaTrackParCov = directLambda.parentTrack; + auto kaonTrackParCov = getTrackParCov(kaon); + int nCandidates = 0; + try { + nCandidates = fitter3Body.process(omegaTrackParCov, lambdaTrackParCov, kaonTrackParCov); + } catch (...) { + LOG(error) << "Exception while fitting the double-Omega candidate"; + return false; + } + if (nCandidates == 0) { + return false; + } + + auto omegaTrackAtVertex = fitter3Body.getTrack(0); + auto lambdaTrackAtVertex = fitter3Body.getTrack(1); + auto kaonTrackAtVertex = fitter3Body.getTrack(2); + std::array omegaMomentum{}; + std::array lambdaMomentum{}; + std::array kaonMomentum{}; + omegaTrackAtVertex.getPxPyPzGlo(omegaMomentum); + lambdaTrackAtVertex.getPxPyPzGlo(lambdaMomentum); + kaonTrackAtVertex.getPxPyPzGlo(kaonMomentum); + + const std::array totalMomentum{ + omegaMomentum[0] + lambdaMomentum[0] + kaonMomentum[0], + omegaMomentum[1] + lambdaMomentum[1] + kaonMomentum[1], + omegaMomentum[2] + lambdaMomentum[2] + kaonMomentum[2]}; + const auto& fittedVertex = fitter3Body.getPCACandidate(); + const std::array decayVertex{ + static_cast(fittedVertex[0]), static_cast(fittedVertex[1]), static_cast(fittedVertex[2])}; + const std::array primaryVertex{collision.posX(), collision.posY(), collision.posZ()}; + const o2::math_utils::Point3D primaryVertexPoint{collision.posX(), collision.posY(), collision.posZ()}; + const std::array omegaDecayVertex{omega.x, omega.y, omega.z}; + + std::array dcaOmega{}; + std::array dcaDirectLambda{}; + std::array dcaDirectKaon{}; + const auto matCorr = static_cast(cfgMaterialCorrection.value); + if (!o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertexPoint, omegaTrackAtVertex, 2.f, matCorr, &dcaOmega) || + !o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertexPoint, lambdaTrackAtVertex, 2.f, matCorr, &dcaDirectLambda) || + !o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertexPoint, kaonTrackAtVertex, 2.f, matCorr, &dcaDirectKaon)) { + return false; + } + + const auto momentumSquared = [](std::array const& momentum) { + return momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2]; + }; + const float energy = std::sqrt(momentumSquared(omegaMomentum) + o2::constants::physics::MassOmegaMinus * o2::constants::physics::MassOmegaMinus) + + std::sqrt(momentumSquared(lambdaMomentum) + o2::constants::physics::MassLambda0 * o2::constants::physics::MassLambda0) + + std::sqrt(momentumSquared(kaonMomentum) + o2::constants::physics::MassKaonCharged * o2::constants::physics::MassKaonCharged); + const float massSquared = energy * energy - totalMomentum[0] * totalMomentum[0] - totalMomentum[1] * totalMomentum[1] - totalMomentum[2] * totalMomentum[2]; + + cand.pt = std::hypot(totalMomentum[0], totalMomentum[1]); + cand.eta = etaFromMomentum(totalMomentum[0], totalMomentum[1], totalMomentum[2]); + cand.phi = std::atan2(totalMomentum[1], totalMomentum[0]); + cand.x = decayVertex[0]; + cand.y = decayVertex[1]; + cand.z = decayVertex[2]; + cand.cosPAOmega = RecoDecay::cpa(decayVertex, omegaDecayVertex, omegaMomentum); + cand.cosPADirectLambda = RecoDecay::cpa(decayVertex, directLambda.decayVertex, lambdaMomentum); + cand.cosPADoubleOmega = RecoDecay::cpa(primaryVertex, decayVertex, totalMomentum); + cand.dcaXYOmegaToPV = dcaOmega[0]; + cand.dcaZOmegaToPV = dcaOmega[1]; + cand.dcaXYDirectLambdaToPV = dcaDirectLambda[0]; + cand.dcaZDirectLambdaToPV = dcaDirectLambda[1]; + cand.dcaXYDirectKaonToPV = dcaDirectKaon[0]; + cand.dcaZDirectKaonToPV = dcaDirectKaon[1]; + if (cand.cosPAOmega < minCosPAOmega || + cand.cosPADirectLambda < minCosPADirectLambda || + cand.cosPADoubleOmega < minCosPADoubleOmega || + std::hypot(decayVertex[0], decayVertex[1]) < minDoubleOmegaDecayRadius || + std::abs(cand.dcaXYOmegaToPV) < dcaOmegaToPV || + std::abs(cand.dcaXYDirectLambdaToPV) < dcaDirectLambdaToPV || + std::abs(cand.dcaXYDirectKaonToPV) < dcaKaonToPV) { + return false; + } + cand.mass = std::sqrt(std::max(0.f, massSquared)); + cand.massOmega = omega.massOmega; + cand.massXi = omega.massXi; + return true; } - template - bool getMCInfo(DoubleOmegaMCInfo& mcInfo, Omega const& omega, V0 const& v0, Kaon const& kaon, - CascLabels const& cascLabels, V0Labels const& v0Labels, TrackLabels const& trackLabels) + template + int64_t getLambdaMCLabel(Track const& posTrack, Track const& negTrack, int8_t sign) { - if (omega.globalIndex() >= cascLabels.size() || - v0.globalIndex() >= v0Labels.size() || - kaon.globalIndex() >= trackLabels.size()) { - LOG(info) << "Skipping MC info retrieval for double Omega candidate with invalid indices: " - << "omega index = " << omega.globalIndex() << ", v0 index = " << v0.globalIndex() - << ", kaon index = " << kaon.globalIndex(); - return false; + if (!posTrack.has_mcParticle() || !negTrack.has_mcParticle()) { + return -1; + } + + auto mcPosTrack = posTrack.template mcParticle_as(); + auto mcNegTrack = negTrack.template mcParticle_as(); + const int expectedPosPdg = sign < 0 ? 2212 : 211; + const int expectedNegPdg = sign < 0 ? -211 : -2212; + const int expectedLambdaPdg = sign < 0 ? 3122 : -3122; + if (mcPosTrack.pdgCode() != expectedPosPdg || mcNegTrack.pdgCode() != expectedNegPdg) { + return -1; + } + + for (const auto& posMother : mcPosTrack.template mothers_as()) { + if (posMother.pdgCode() != expectedLambdaPdg) { + continue; + } + for (const auto& negMother : mcNegTrack.template mothers_as()) { + if (posMother.globalIndex() == negMother.globalIndex()) { + return posMother.globalIndex(); + } + } + } + return -1; + } + + template + int64_t getOmegaMCLabel(Track const& posTrack, Track const& negTrack, Track const& bachelorTrack, + int8_t sign, McParticles const& mcParticles) + { + const int64_t lambdaLabel = getLambdaMCLabel(posTrack, negTrack, sign); + if (lambdaLabel < 0 || !bachelorTrack.has_mcParticle()) { + return -1; + } + + auto mcLambda = mcParticles.rawIteratorAt(lambdaLabel); + auto mcBachelor = bachelorTrack.template mcParticle_as(); + const int expectedBachelorPdg = sign < 0 ? -321 : 321; + const int expectedOmegaPdg = sign < 0 ? 3334 : -3334; + if (mcBachelor.pdgCode() != expectedBachelorPdg) { + return -1; } - auto omegaLabel = cascLabels.rawIteratorAt(omega.globalIndex()); - auto v0Label = v0Labels.rawIteratorAt(v0.globalIndex()); - auto kaonLabel = trackLabels.rawIteratorAt(kaon.globalIndex()); - if (!omegaLabel.has_mcParticle() || !v0Label.has_mcParticle() || !kaonLabel.has_mcParticle()) { + for (const auto& lambdaMother : mcLambda.template mothers_as()) { + if (lambdaMother.pdgCode() != expectedOmegaPdg) { + continue; + } + for (const auto& bachelorMother : mcBachelor.template mothers_as()) { + if (lambdaMother.globalIndex() == bachelorMother.globalIndex()) { + return lambdaMother.globalIndex(); + } + } + } + return -1; + } + + template + bool getMCInfo(DoubleOmegaMCInfo& mcInfo, + Track const& omegaPosTrack, Track const& omegaNegTrack, Track const& omegaBachelorTrack, + Track const& lambdaPosTrack, Track const& lambdaNegTrack, Track const& directKaonTrack, + int8_t sign, McParticles const& mcParticles) + { + const int64_t omegaLabel = getOmegaMCLabel(omegaPosTrack, omegaNegTrack, omegaBachelorTrack, sign, mcParticles); + const int64_t lambdaLabel = getLambdaMCLabel(lambdaPosTrack, lambdaNegTrack, sign); + if (omegaLabel < 0 || lambdaLabel < 0 || !directKaonTrack.has_mcParticle()) { return false; } - auto mcOmega = omegaLabel.template mcParticle_as(); - auto mcLambda = v0Label.template mcParticle_as(); - auto mcKaon = kaonLabel.template mcParticle_as(); - const bool isMatter = omega.sign() < 0; - if (mcOmega.pdgCode() != (isMatter ? 3334 : -3334) || - mcLambda.pdgCode() != (isMatter ? 3122 : -3122) || - mcKaon.pdgCode() != (isMatter ? -321 : 321)) { + auto mcOmega = mcParticles.rawIteratorAt(omegaLabel); + auto mcLambda = mcParticles.rawIteratorAt(lambdaLabel); + auto mcKaon = directKaonTrack.template mcParticle_as(); + const int expectedKaonPdg = sign < 0 ? -321 : 321; + const int expectedDoubleOmegaPdg = sign < 0 ? kDoubleOmegaPdg : -kDoubleOmegaPdg; + if (mcKaon.pdgCode() != expectedKaonPdg) { return false; } + for (const auto& omegaMother : mcOmega.template mothers_as()) { - if (omegaMother.pdgCode() != (isMatter ? kDoubleOmegaPdg : -kDoubleOmegaPdg)) { + if (omegaMother.pdgCode() != expectedDoubleOmegaPdg) { continue; } for (const auto& lambdaMother : mcLambda.template mothers_as()) { @@ -412,18 +687,17 @@ struct doubleOmegaTreeCreator { continue; } for (const auto& kaonMother : mcKaon.template mothers_as()) { - if (omegaMother.globalIndex() != kaonMother.globalIndex()) { - continue; + if (omegaMother.globalIndex() == kaonMother.globalIndex()) { + mcInfo.motherId = omegaMother.globalIndex(); + mcInfo.pt = omegaMother.pt(); + mcInfo.eta = omegaMother.eta(); + mcInfo.phi = omegaMother.phi(); + mcInfo.decayLength = std::hypot(mcOmega.vx() - omegaMother.vx(), + mcOmega.vy() - omegaMother.vy(), + mcOmega.vz() - omegaMother.vz()); + mcInfo.pdgCode = omegaMother.pdgCode(); + return true; } - mcInfo.motherId = omegaMother.globalIndex(); - mcInfo.pt = omegaMother.pt(); - mcInfo.eta = omegaMother.eta(); - mcInfo.phi = omegaMother.phi(); - mcInfo.decayLength = std::hypot(mcOmega.vx() - omegaMother.vx(), - mcOmega.vy() - omegaMother.vy(), - mcOmega.vz() - omegaMother.vz()); - mcInfo.pdgCode = omegaMother.pdgCode(); - return true; } } } @@ -471,115 +745,65 @@ struct doubleOmegaTreeCreator { return true; } - template - void fillFindableDaughters(Candidates const& candidates, Labels const& labels, int expectedAbsPdg, FindableDaughter daughterType) + template + bool getLambdaFinalStateIds(McParticle const& lambda, int sign, std::array& finalStateIds) { - for (const auto& candidate : candidates) { - - if (expectedAbsPdg == 3334) { - LOG(info) << "Checking candidate with global index " << candidate.globalIndex() << " for findable daughter of type " << daughterType; - LOG(info) << "Total number of labels: " << labels.size(); - } - - if (candidate.globalIndex() >= labels.size()) { - continue; - } - auto label = labels.rawIteratorAt(candidate.globalIndex()); - if (!label.has_mcParticle()) { - continue; - } - auto mcParticle = label.template mcParticle_as(); - if (std::abs(mcParticle.pdgCode()) != expectedAbsPdg) { - continue; - } - if(expectedAbsPdg == 3334){ - LOG(info) << "Found a cascade with pdg code " << mcParticle.pdgCode() << " and global index " << mcParticle.globalIndex(); - for (const auto& mother : mcParticle.template mothers_as()) { - LOG(info) << "Mother pdg code: " << mother.pdgCode() << ", global index: " << mother.globalIndex(); - } - } - for (const auto& mother : mcParticle.template mothers_as()) { - if (std::abs(mother.pdgCode()) == kDoubleOmegaPdg) { - histos.fill(HIST("MC/findableDaughters"), daughterType); - break; - } + bool foundProton = false; + bool foundPion = false; + for (const auto& daughter : lambda.template daughters_as()) { + if (daughter.pdgCode() == sign * 2212 && !foundProton) { + finalStateIds[0] = daughter.globalIndex(); + foundProton = true; + } else if (daughter.pdgCode() == -sign * 211 && !foundPion) { + finalStateIds[1] = daughter.globalIndex(); + foundPion = true; } } - } - - template - void checkFindables(Cascades const& cascades, V0s const& v0s, Tracks const& tracks, - CascLabels const& cascLabels, V0Labels const& v0Labels, TrackLabels const& trackLabels) - { - fillFindableDaughters(cascades, cascLabels, 3334, kOmega); - fillFindableDaughters(v0s, v0Labels, 3122, kLambda); - fillFindableDaughters(tracks, trackLabels, 321, kKaon); + return foundProton && foundPion; } void writeDataCandidate(DoubleOmegaCandidate const& cand) { - doubleOmegaTable(cand.ptCasc1, - cand.etaCasc1, - cand.phiCasc1, - cand.cascDecLength1, - cand.omegaMassCasc1, - cand.xiMassCasc1, - cand.cosPACasc1, - cand.dcaBachPVCasc1, - cand.dcaV0BachCasc1, - cand.nSigmaKBach1, - cand.ptLambda, - cand.etaLambda, - cand.phiLambda, - cand.lambdaDecLength, - cand.lambdaMass, - cand.cosPALambda, - cand.dcaPosPVLambda, - cand.dcaNegPVLambda, - cand.dcaLambdaDaughters, - cand.nSigmaPrLambda, - cand.nSigmaPiLambda, - cand.ptKaon, - cand.etaKaon, - cand.phiKaon, - cand.chargeKaon, - cand.dcaXYKaon, - cand.dcaZKaon, - cand.nSigmaKKaon, - cand.doubleOmegaMass); + doubleOmegaTable(cand.pt, + cand.eta, + cand.phi, + cand.x, + cand.y, + cand.z, + cand.cosPAOmega, + cand.cosPADirectLambda, + cand.cosPADoubleOmega, + cand.dcaXYOmegaToPV, + cand.dcaZOmegaToPV, + cand.dcaXYDirectLambdaToPV, + cand.dcaZDirectLambdaToPV, + cand.dcaXYDirectKaonToPV, + cand.dcaZDirectKaonToPV, + cand.mass, + cand.massOmega, + cand.massXi); } void writeMCCandidate(DoubleOmegaCandidate const& cand, DoubleOmegaMCInfo const& mcInfo, bool isReco) { - doubleOmegaTableMC(cand.ptCasc1, - cand.etaCasc1, - cand.phiCasc1, - cand.cascDecLength1, - cand.omegaMassCasc1, - cand.xiMassCasc1, - cand.cosPACasc1, - cand.dcaBachPVCasc1, - cand.dcaV0BachCasc1, - cand.nSigmaKBach1, - cand.ptLambda, - cand.etaLambda, - cand.phiLambda, - cand.lambdaDecLength, - cand.lambdaMass, - cand.cosPALambda, - cand.dcaPosPVLambda, - cand.dcaNegPVLambda, - cand.dcaLambdaDaughters, - cand.nSigmaPrLambda, - cand.nSigmaPiLambda, - cand.ptKaon, - cand.etaKaon, - cand.phiKaon, - cand.chargeKaon, - cand.dcaXYKaon, - cand.dcaZKaon, - cand.nSigmaKKaon, - cand.doubleOmegaMass, + doubleOmegaTableMC(cand.pt, + cand.eta, + cand.phi, + cand.x, + cand.y, + cand.z, + cand.cosPAOmega, + cand.cosPADirectLambda, + cand.cosPADoubleOmega, + cand.dcaXYOmegaToPV, + cand.dcaZOmegaToPV, + cand.dcaXYDirectLambdaToPV, + cand.dcaZDirectLambdaToPV, + cand.dcaXYDirectKaonToPV, + cand.dcaZDirectKaonToPV, + cand.mass, + cand.massOmega, + cand.massXi, mcInfo.pt, mcInfo.eta, mcInfo.phi, @@ -588,40 +812,57 @@ struct doubleOmegaTreeCreator { isReco); } - template - void fillFromV0s(C const& collision, T const& tracks, Cascades const& cascades, V0s const& v0s, McLabels const&... mcLabels) + template + void fillFromRawTables(C const& collision, T const& tracks, Cascades const& cascades, V0s const& v0s, + aod::McParticles const* mcParticles = nullptr) { - for (const auto& omega : cascades) { - LOG(info) << "Processing omega with global index " << omega.globalIndex(); - if (!isSelectedOmega(collision, tracks, omega)) { + for (const auto& rawCascade : cascades) { + BuiltOmegaCandidate omega; + if (!buildOmega(tracks, rawCascade, omega)) { continue; } - histos.fill(HIST("QA/massXi"), omega.pt(), omega.mXi()); - histos.fill(HIST("QA/massOmega"), omega.pt(), omega.mOmega()); - const std::array omegaTrackIds{omega.posTrackId(), omega.negTrackId(), omega.bachelorId()}; - for (const auto& v0 : v0s) { - const int8_t charge = omega.sign(); - if (!isSelectedLambda(collision, tracks, v0, charge)) { + const float omegaPt = std::hypot(omega.px, omega.py); + histos.fill(HIST("QA/massXi"), omegaPt, omega.massXi); + histos.fill(HIST("QA/massOmega"), omegaPt, omega.massOmega); + const std::array omegaTrackIds{omega.posTrackId, omega.negTrackId, omega.bachelorId}; + + for (const auto& rawV0 : v0s) { + BuiltLambdaCandidate lambda; + if (!buildLambda(tracks, rawV0, omega.sign, lambda)) { continue; } - if (std::find(omegaTrackIds.begin(), omegaTrackIds.end(), v0.posTrackId()) != omegaTrackIds.end() || - std::find(omegaTrackIds.begin(), omegaTrackIds.end(), v0.negTrackId()) != omegaTrackIds.end()) { + if (std::find(omegaTrackIds.begin(), omegaTrackIds.end(), lambda.posTrackId) != omegaTrackIds.end() || + std::find(omegaTrackIds.begin(), omegaTrackIds.end(), lambda.negTrackId) != omegaTrackIds.end()) { continue; } - const auto lambda = getLambdaCandidateFromV0(collision, tracks, v0, charge); - histos.fill(HIST("QA/massLambda"), std::hypot(lambda.px, lambda.py), lambda.mass); + histos.fill(HIST("QA/massLambda"), std::hypot(lambda.candidate.px, lambda.candidate.py), lambda.candidate.mass); + for (const auto& kaon : tracks) { - if (!isSelectedKaon(kaon, charge) || - kaon.globalIndex() == v0.posTrackId() || - kaon.globalIndex() == v0.negTrackId() || + if (!isSelectedKaon(kaon, omega.sign) || + kaon.globalIndex() == lambda.posTrackId || + kaon.globalIndex() == lambda.negTrackId || std::find(omegaTrackIds.begin(), omegaTrackIds.end(), kaon.globalIndex()) != omegaTrackIds.end()) { continue; } - auto cand = makeCandidate(collision, omega, lambda, kaon); + DoubleOmegaCandidate cand; + if (!buildDoubleOmega(collision, omega, lambda, kaon, cand)) { + continue; + } if constexpr (isMC) { - static_assert(sizeof...(McLabels) == 3); + if (mcParticles == nullptr) { + continue; + } + auto omegaV0 = rawCascade.template v0_as(); + auto omegaPosTrack = omegaV0.template posTrack_as(); + auto omegaNegTrack = omegaV0.template negTrack_as(); + auto omegaBachelorTrack = rawCascade.template bachelor_as(); + auto lambdaPosTrack = rawV0.template posTrack_as(); + auto lambdaNegTrack = rawV0.template negTrack_as(); DoubleOmegaMCInfo mcInfo; - if (getMCInfo(mcInfo, omega, v0, kaon, mcLabels...)) { + if (getMCInfo(mcInfo, + omegaPosTrack, omegaNegTrack, omegaBachelorTrack, + lambdaPosTrack, lambdaNegTrack, kaon, + omega.sign, *mcParticles)) { reconstructedDoubleOmegaIds.push_back(mcInfo.motherId); writeMCCandidate(cand, mcInfo, true); } @@ -636,18 +877,19 @@ struct doubleOmegaTreeCreator { template void fillFromCascades(C const& collision, T const& tracks, Cascades const& cascades) { - for (const auto& omega : cascades) { - if (!isSelectedOmega(collision, tracks, omega)) { + for (const auto& omegaRow : cascades) { + BuiltOmegaCandidate omega; + if (!buildOmegaFromBuilder(tracks, omegaRow, omega)) { continue; } - histos.fill(HIST("QA/massXi"), omega.pt(), omega.mXi()); - histos.fill(HIST("QA/massOmega"), omega.pt(), omega.mOmega()); - const std::array omegaTrackIds{omega.posTrackId(), omega.negTrackId(), omega.bachelorId()}; + const float omegaPt = std::hypot(omega.px, omega.py); + histos.fill(HIST("QA/massXi"), omegaPt, omega.massXi); + histos.fill(HIST("QA/massOmega"), omegaPt, omega.massOmega); + const std::array omegaTrackIds{omega.posTrackId, omega.negTrackId, omega.bachelorId}; for (const auto& lambdaKaonSource : cascades) { - if (lambdaKaonSource.globalIndex() == omega.globalIndex() || - lambdaKaonSource.sign() != omega.sign() || - !isSelectedOmega(collision, tracks, lambdaKaonSource)) { + if (lambdaKaonSource.globalIndex() == omegaRow.globalIndex() || + lambdaKaonSource.sign() != omega.sign) { continue; } @@ -665,9 +907,18 @@ struct doubleOmegaTreeCreator { } auto kaon = lambdaKaonSource.template bachelor_as(); - const auto lambda = getLambdaCandidateFromCascade(collision, tracks, lambdaKaonSource); - histos.fill(HIST("QA/massLambda"), std::hypot(lambda.px, lambda.py), lambda.mass); - writeDataCandidate(makeCandidate(collision, omega, lambda, kaon)); + if (!isSelectedKaon(kaon, omega.sign)) { + continue; + } + BuiltLambdaCandidate lambda; + if (!buildLambda(tracks, lambdaKaonSource, omega.sign, lambda)) { + continue; + } + histos.fill(HIST("QA/massLambda"), std::hypot(lambda.candidate.px, lambda.candidate.py), lambda.candidate.mass); + DoubleOmegaCandidate cand; + if (buildDoubleOmega(collision, omega, lambda, kaon, cand)) { + writeDataCandidate(cand); + } } } } @@ -692,16 +943,16 @@ struct doubleOmegaTreeCreator { } void processData(Collisions const& collision, - TracksFull const& tracks, - FullCascades const& cascades, - FullV0s const& v0s, + TracksFullIU const& tracks, + aod::V0s const& v0s, + aod::Cascades const& cascades, aod::BCsWithTimestamps const& bcs) { if (acceptCollision(collision, bcs)) { - fillFromV0s(collision, tracks, cascades, v0s); + fillFromRawTables(collision, tracks, cascades, v0s); } } - PROCESS_SWITCH(doubleOmegaTreeCreator, processData, "Reconstruct Omega + Lambda(V0) + kaon(track)", true); + PROCESS_SWITCH(doubleOmegaTreeCreator, processData, "Reconstruct Omega + Lambda + kaon from raw V0/cascade indices", true); void processDataFromCascades(Collisions const& collision, TracksFull const& tracks, @@ -715,30 +966,29 @@ struct doubleOmegaTreeCreator { PROCESS_SWITCH(doubleOmegaTreeCreator, processDataFromCascades, "Reconstruct Omega + Lambda + kaon from two cascade rows", false); void processMC(CollisionsMC const& collisions, - TracksFull const& tracks, - FullCascades const& cascades, - FullV0s const& v0s, - aod::McTrackLabels const& trackLabels, - aod::McCascLabels const& cascLabels, - aod::McV0Labels const& v0Labels, + TracksFullIUMC const& tracks, + aod::V0s const& v0s, + aod::Cascades const& cascades, aod::McParticles const& mcParticles, - aod::BCsWithTimestamps const&) + aod::BCsWithTimestamps const& bcs) { reconstructedDoubleOmegaIds.clear(); for (const auto& collision : collisions) { - // if (!acceptCollision(collision, bcs)) { - // continue; - // } - auto tracksThisCollision = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); - auto cascadesThisCollision = cascades.sliceBy(cascadesPerCollision, collision.globalIndex()); - auto v0sThisCollision = v0s.sliceBy(v0sPerCollision, collision.globalIndex()); - cascadesThisCollision.bindExternalIndices(&tracks); + auto bc = collision.template bc_as(); + initCCDB(bc); + if (!acceptCollision(collision, bcs)) { + continue; + } + auto tracksThisCollision = tracks.sliceBy(tracksIUPerCollision, collision.globalIndex()); + auto cascadesThisCollision = cascades.sliceBy(rawCascadesPerCollision, collision.globalIndex()); + auto v0sThisCollision = v0s.sliceBy(rawV0sPerCollision, collision.globalIndex()); v0sThisCollision.bindExternalIndices(&tracks); - checkFindables(cascadesThisCollision, v0sThisCollision, tracksThisCollision, cascLabels, v0Labels, trackLabels); - fillFromV0s(collision, tracksThisCollision, cascadesThisCollision, v0sThisCollision, cascLabels, v0Labels, trackLabels); + cascadesThisCollision.bindExternalIndices(&tracks); + cascadesThisCollision.bindExternalIndices(&v0s); + fillFromRawTables(collision, tracksThisCollision, cascadesThisCollision, v0sThisCollision, &mcParticles); } - + LOG(info) << "Found " << reconstructedDoubleOmegaIds.size() << " reconstructed double-Omega candidates in MC"; for (const auto& mcParticle : mcParticles) { DoubleOmegaMCInfo mcInfo; if (!getGeneratedMCInfo(mcInfo, mcParticle)) { @@ -752,6 +1002,165 @@ struct doubleOmegaTreeCreator { } } PROCESS_SWITCH(doubleOmegaTreeCreator, processMC, "Reconstruct and truth-match Omega + Lambda + kaon in MC", false); + + void processFindableTracks(TracksFullMC const& tracks, + aod::Cascades const& cascades, + aod::V0s const& v0s, + aod::McParticles const& mcParticles) + { + for (const auto& mcParticle : mcParticles) { + if (std::abs(mcParticle.pdgCode()) != kDoubleOmegaPdg) { + continue; + } + histos.fill(HIST("MC/generatedAndFindable"), kAllGenerated); + + const int sign = mcParticle.pdgCode() > 0 ? 1 : -1; + int64_t omegaId = -1; + int64_t directLambdaId = -1; + int64_t directKaonId = -1; + + // First stage: double-Omega -> Omega + Lambda + kaon. + for (const auto& daughter : mcParticle.daughters_as()) { + if (daughter.pdgCode() == sign * 3334 && omegaId < 0) { + omegaId = daughter.globalIndex(); + } else if (daughter.pdgCode() == sign * 3122 && directLambdaId < 0) { + directLambdaId = daughter.globalIndex(); + } else if (daughter.pdgCode() == -sign * 321 && directKaonId < 0) { + directKaonId = daughter.globalIndex(); + } + } + if (omegaId < 0 || directLambdaId < 0 || directKaonId < 0) { + continue; + } + histos.fill(HIST("MC/generatedAndFindable"), kOmegaLambdaKaon); + + // Second stage: Omega -> Lambda + kaon, giving two Lambdas and two kaons. + const auto omega = mcParticles.rawIteratorAt(omegaId); + int64_t omegaLambdaId = -1; + int64_t omegaKaonId = -1; + for (const auto& daughter : omega.daughters_as()) { + if (daughter.pdgCode() == sign * 3122 && omegaLambdaId < 0) { + omegaLambdaId = daughter.globalIndex(); + } else if (daughter.pdgCode() == -sign * 321 && omegaKaonId < 0) { + omegaKaonId = daughter.globalIndex(); + } + } + if (omegaLambdaId < 0 || omegaKaonId < 0) { + continue; + } + histos.fill(HIST("MC/generatedAndFindable"), kTwoLambdasTwoKaons); + + // Third stage: both Lambdas -> proton + pion. + const auto omegaLambda = mcParticles.rawIteratorAt(omegaLambdaId); + const auto directLambda = mcParticles.rawIteratorAt(directLambdaId); + std::array omegaLambdaFinalStateIds{}; + std::array directLambdaFinalStateIds{}; + if (!getLambdaFinalStateIds(omegaLambda, sign, omegaLambdaFinalStateIds)) { + continue; + } + if (!getLambdaFinalStateIds(directLambda, sign, directLambdaFinalStateIds)) { + continue; + } + histos.fill(HIST("MC/generatedAndFindable"), kFinalState); + + const std::array finalStateIds{ + omegaKaonId, + omegaLambdaFinalStateIds[0], + omegaLambdaFinalStateIds[1], + directLambdaFinalStateIds[0], + directLambdaFinalStateIds[1], + directKaonId}; + std::array reconstructedDaughters{}; + std::array selectedReconstructedDaughters{}; + std::array daughterTrackIds{}; + daughterTrackIds.fill(-1); + for (const auto& track : tracks) { + if (!track.has_mcParticle()) { + continue; + } + const auto daughterId = std::find(finalStateIds.begin(), finalStateIds.end(), track.mcParticleId()); + if (daughterId != finalStateIds.end()) { + const auto daughterIndex = std::distance(finalStateIds.begin(), daughterId); + reconstructedDaughters[daughterIndex] = true; + if (daughterTrackIds[daughterIndex] < 0) { + daughterTrackIds[daughterIndex] = track.globalIndex(); + } + if (selectTrack(track)) { + selectedReconstructedDaughters[daughterIndex] = true; + } + } + } + bool allDaughtersReconstructed = true; + bool allDaughtersSelected = true; + for (size_t iDaughter = 0; iDaughter < reconstructedDaughters.size(); ++iDaughter) { + allDaughtersReconstructed &= reconstructedDaughters[iDaughter]; + allDaughtersSelected &= selectedReconstructedDaughters[iDaughter]; + } + if (allDaughtersReconstructed) { + histos.fill(HIST("MC/generatedAndFindable"), kFindable); + bool omegaInCascade = false; + for (const auto& cascade : cascades) { + const auto cascadeV0 = cascade.v0(); + const std::array cascadeTrackIds{ + cascadeV0.posTrackId(), cascadeV0.negTrackId(), cascade.bachelorId()}; + bool allOmegaTracksInCascade = true; + for (size_t iDaughter = 0; iDaughter < 3; ++iDaughter) { + if (std::find(cascadeTrackIds.begin(), cascadeTrackIds.end(), daughterTrackIds[iDaughter]) == cascadeTrackIds.end()) { + allOmegaTracksInCascade = false; + break; + } + } + if (allOmegaTracksInCascade) { + omegaInCascade = true; + break; + } + } + + bool directLambdaInV0 = false; + for (const auto& v0 : v0s) { + const std::array v0TrackIds{v0.posTrackId(), v0.negTrackId()}; + bool allDirectLambdaTracksInV0 = true; + for (size_t iDaughter = 3; iDaughter < 5; ++iDaughter) { + if (std::find(v0TrackIds.begin(), v0TrackIds.end(), daughterTrackIds[iDaughter]) == v0TrackIds.end()) { + allDirectLambdaTracksInV0 = false; + break; + } + } + if (allDirectLambdaTracksInV0) { + directLambdaInV0 = true; + break; + } + } + + if (omegaInCascade && directLambdaInV0) { + histos.fill(HIST("MC/generatedAndFindable"), kFindableCascadeAndV0); + } + } + if (allDaughtersSelected) { + histos.fill(HIST("MC/generatedAndFindable"), kFindableSelectedTracks); + LOG(debug) << "----------------------------------------"; + for (size_t iDaughter = 0; iDaughter < daughterTrackIds.size(); ++iDaughter) { + LOG(debug) << "+++++++"; + const auto daughterTrack = tracks.rawIteratorAt(daughterTrackIds[iDaughter]); + const auto daughterMCParticle = daughterTrack.mcParticle_as(); + LOG(debug) << "Dau" << iDaughter + 1 + << ": hasITS = " << daughterTrack.hasITS() + << ", hasTPC = " << daughterTrack.hasTPC() + << ", hasTOF = " << daughterTrack.hasTOF() + << ", hasTRD = " << daughterTrack.hasTRD() + << ", isITSAfterburner = " << daughterTrack.isITSAfterburner(); + LOG(debug) << "Dau" << iDaughter + 1 + << ": eta = " << daughterTrack.eta() + << ", pt = " << daughterTrack.pt() + << ", CollisionID = " << daughterTrack.collisionId() + << ", production point = (" << daughterMCParticle.vx() + << ", " << daughterMCParticle.vy() + << ", " << daughterMCParticle.vz() << ") cm"; + } + } + } + } + PROCESS_SWITCH(doubleOmegaTreeCreator, processFindableTracks, "Count generated and track-findable double-Omega decays", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)