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..1657f31110e --- /dev/null +++ b/PWGLF/DataModel/LFDoubleOmegaTables.h @@ -0,0 +1,87 @@ +// 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(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); +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::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); + +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..1b8a649d7cc --- /dev/null +++ b/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx @@ -0,0 +1,1169 @@ +// 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/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 +#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 TracksFull = soa::Join; +using TracksFullIU = soa::Join; +using TracksFullIUMC = soa::Join; +using TracksFullMC = soa::Join; + +struct DoubleOmegaCandidate { + 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 { + 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 mass = -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 { + Produces doubleOmegaTable; + Produces doubleOmegaTableMC; + Service ccdb; + + static constexpr int kDoubleOmegaPdg = 1060020020; + enum FindabilityStep : uint8_t { + kAllGenerated, + kOmegaLambdaKaon, + kTwoLambdasTwoKaons, + kFinalState, + kFindable, + kFindableSelectedTracks, + kFindableCascadeAndV0, + kNFindabilitySteps + }; + std::vector reconstructedDoubleOmegaIds; + + 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"}; + 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 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 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}; + + 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; + } + 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"); + zorro.populateHistRegistry(histos, bc.runNumber()); + } + } + + void init(InitContext const&) + { + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + 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/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; + } + + static float invariantMass2Body(std::array const& momentum1, float mass1, + std::array const& momentum2, float mass2) + { + 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)); + } + + 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; + } + + 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 (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 buildOmega(T const& tracks, Casc const& cascade, BuiltOmegaCandidate& builtOmega) + { + auto bachelor = cascade.template bachelor_as(); + if (bachelor.sign() == 0 || + !selectTrack(bachelor) || + std::abs(bachelor.tpcNSigmaKa()) > nSigmaTPCCut || + std::abs(bachelor.dcaXY()) < dcaBachToPV) { + return false; + } + + auto v0 = cascade.template v0_as(); + BuiltLambdaCandidate lambda; + if (!buildLambda(tracks, v0, 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 fitting raw 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(), + v0.posTrackId(), + v0.negTrackId(), + cascade.bachelorId(), + fitter.createParentTrackParCov()}; + return true; + } + + template + bool buildOmegaFromBuilder(T const& tracks, Casc const& cascade, BuiltOmegaCandidate& builtOmega) + { + 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) + { + const float pt = std::hypot(px, py); + return pt > 0.f ? std::asinh(pz / pt) : 0.f; + } + + template + bool buildDoubleOmega(C const& collision, BuiltOmegaCandidate const& omega, BuiltLambdaCandidate const& directLambda, Kaon const& kaon, DoubleOmegaCandidate& 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 + int64_t getLambdaMCLabel(Track const& posTrack, Track const& negTrack, int8_t sign) + { + 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; + } + + 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 = 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() != expectedDoubleOmegaPdg) { + 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()) { + 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 + bool getLambdaFinalStateIds(McParticle const& lambda, int sign, std::array& finalStateIds) + { + 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; + } + } + return foundProton && foundPion; + } + + void writeDataCandidate(DoubleOmegaCandidate const& cand) + { + 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.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, + mcInfo.decayLength, + mcInfo.pdgCode, + isReco); + } + + template + void fillFromRawTables(C const& collision, T const& tracks, Cascades const& cascades, V0s const& v0s, + aod::McParticles const* mcParticles = nullptr) + { + for (const auto& rawCascade : cascades) { + BuiltOmegaCandidate omega; + if (!buildOmega(tracks, rawCascade, omega)) { + continue; + } + 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(), lambda.posTrackId) != omegaTrackIds.end() || + std::find(omegaTrackIds.begin(), omegaTrackIds.end(), lambda.negTrackId) != omegaTrackIds.end()) { + continue; + } + histos.fill(HIST("QA/massLambda"), std::hypot(lambda.candidate.px, lambda.candidate.py), lambda.candidate.mass); + + for (const auto& kaon : tracks) { + if (!isSelectedKaon(kaon, omega.sign) || + kaon.globalIndex() == lambda.posTrackId || + kaon.globalIndex() == lambda.negTrackId || + std::find(omegaTrackIds.begin(), omegaTrackIds.end(), kaon.globalIndex()) != omegaTrackIds.end()) { + continue; + } + DoubleOmegaCandidate cand; + if (!buildDoubleOmega(collision, omega, lambda, kaon, cand)) { + continue; + } + if constexpr (isMC) { + 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, + omegaPosTrack, omegaNegTrack, omegaBachelorTrack, + lambdaPosTrack, lambdaNegTrack, kaon, + omega.sign, *mcParticles)) { + 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& omegaRow : cascades) { + BuiltOmegaCandidate omega; + if (!buildOmegaFromBuilder(tracks, omegaRow, omega)) { + continue; + } + 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() == omegaRow.globalIndex() || + lambdaKaonSource.sign() != omega.sign) { + 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(); + 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); + } + } + } + } + + 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, + TracksFullIU const& tracks, + aod::V0s const& v0s, + aod::Cascades const& cascades, + aod::BCsWithTimestamps const& bcs) + { + if (acceptCollision(collision, bcs)) { + fillFromRawTables(collision, tracks, cascades, v0s); + } + } + PROCESS_SWITCH(doubleOmegaTreeCreator, processData, "Reconstruct Omega + Lambda + kaon from raw V0/cascade indices", 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, + TracksFullIUMC const& tracks, + aod::V0s const& v0s, + aod::Cascades const& cascades, + aod::McParticles const& mcParticles, + aod::BCsWithTimestamps const& bcs) + { + reconstructedDoubleOmegaIds.clear(); + + for (const auto& collision : collisions) { + 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); + 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)) { + 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); + + 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) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}