diff --git a/PWGHF/Core/SelectorCuts.h b/PWGHF/Core/SelectorCuts.h index 0351681ce23..1bb1495a6d6 100644 --- a/PWGHF/Core/SelectorCuts.h +++ b/PWGHF/Core/SelectorCuts.h @@ -755,7 +755,7 @@ static const std::vector labelsCutVar = {"pT Ka from Omegac"}; namespace hf_cuts_xic_to_xi_pi { -static constexpr int NBinsPt = 11; +static constexpr int NBinsPt = 12; static constexpr int NCutVars = 28; // default values for the pT bin edges (can be used to configure histogram axis) // offset by 1 from the bin numbers in cuts array @@ -771,22 +771,24 @@ constexpr double BinsPt[NBinsPt + 1] = { 10.0, 12.0, 16.0, - 24.0}; + 24.0, + 50.0}; const auto vecBinsPt = std::vector{BinsPt, BinsPt + NBinsPt + 1}; // default values for the cuts -constexpr double Cuts[NBinsPt][NCutVars] = {{0.2, 0.99, 0.97, 0.99, 0.99, 0.1, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 0 < pt < 1 */ - {0.5, 0.99, 0.97, 0.99, 0.99, 0.1, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 1 < pt < 2 */ - {0.5, 0.99, 0.97, 0.99, 0.99, 0.1, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 2 < pt < 3 */ - {0.5, 0.99, 0.97, 0.99, 0.99, 0.1, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 3 < pt < 4 */ - {0.5, 0.99, 0.97, 0.99, 0.99, 0.1, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 4 < pt < 5 */ - {0.5, 0.99, 0.97, 0.99, 0.99, 0.1, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 5 < pt < 6 */ - {0.5, 0.99, 0.97, 0.99, 0.99, 0.1, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 6 < pt < 8 */ - {0.5, 0.99, 0.97, 0.99, 0.99, 0.1, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 8 < pt < 10 */ - {0.5, 0.99, 0.97, 0.99, 0.99, 0.1, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 10 < pt < 12 */ - {0.5, 0.99, 0.97, 0.99, 0.99, 0.1, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 12 < pt < 16 */ - {0.5, 0.99, 0.97, 0.99, 0.99, 0.1, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}}; /* 16 < pt < 24 */ +constexpr double Cuts[NBinsPt][NCutVars] = {{0.2, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 0 < pt < 1 */ + {0.5, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 1 < pt < 2 */ + {0.5, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 2 < pt < 3 */ + {0.5, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 3 < pt < 4 */ + {0.5, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 4 < pt < 5 */ + {0.5, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 5 < pt < 6 */ + {0.5, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 6 < pt < 8 */ + {0.5, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 8 < pt < 10 */ + {0.5, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 10 < pt < 12 */ + {0.5, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 12 < pt < 16 */ + {0.5, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}, /* 16 < pt < 24 */ + {0.5, 0.97, 0.97, 0.99, 0.99, 2.0, 0.2, 1.0, 0.04, 0.06, 0.06, 0.05, 0.3, 70, 60, 100, 120, 250, 250, 0.4, 100, 300, 0., 0., 1.5, 0., 0., 0.4}}; /* 24 < pt < 50 -> Implemented just for test*/ // row labels static const std::vector labelsPt = { @@ -800,14 +802,16 @@ static const std::vector labelsPt = { "pT bin 7", "pT bin 8", "pT bin 9", - "pT bin 10"}; + "pT bin 10", + "pT bin 11"}; // column labels static const std::vector labelsCutVar = {"ptPiFromCharmBaryon", "cosPACasc", "cosPAV0", "cosPaCascToXic", "cosPaV0ToCasc", - "dcaCharmBaryonDau", "dcaCascDau", "dcaV0Dau", "dcaXYToPvCascDau", "dcaXYToPvV0Dau0", "dcaXYToPvV0Dau1", "kfDcaXYPiFromXic", "kfDcaXYCascToPv", - "chi2GeoXic", "chi2GeoCasc", "chi2GeoV0", - "chi2TopoXicToPv", "chi2TopoPiFromXicToPv", "chi2TopoCascToPv", "chi2TopoV0ToPv", "chi2TopoV0ToCasc", "chi2TopoCascToXic", - "cascldl", "v0ldl", "decayLenXYXic", "decayLenXYCasc", "decayLenXYLambda", "cTauXic"}; + "dcaCharmBaryonDau", "dcaCascDau", "dcaV0Dau", "dcaXYToPvCascDau", "dcaXYToPvV0Dau0", + "dcaXYToPvV0Dau1", "kfDcaXYPiFromXic", "kfDcaXYCascToPv", "chi2GeoXic", "chi2GeoCasc", + "chi2GeoV0", "chi2TopoXicToPv", "chi2TopoPiFromXicToPv", "chi2TopoCascToPv", "chi2TopoV0ToPv", + "chi2TopoV0ToCasc", "chi2TopoCascToXic", "cascldl", "v0ldl", "decayLenXYXic", + "decayLenXYCasc", "decayLenXYLambda", "cTauXic"}; } // namespace hf_cuts_xic_to_xi_pi namespace hf_cuts_xic_to_p_k_pi diff --git a/PWGHF/D2H/Tasks/taskXic0ToXiPi.cxx b/PWGHF/D2H/Tasks/taskXic0ToXiPi.cxx index c9c09f94859..510a5045916 100644 --- a/PWGHF/D2H/Tasks/taskXic0ToXiPi.cxx +++ b/PWGHF/D2H/Tasks/taskXic0ToXiPi.cxx @@ -165,6 +165,15 @@ struct HfTaskXic0ToXiPi { registry.add("hMassVsPtVsYVsCentVsPtPion", "Thn for Xic0 candidates with Cent&pTpi", HistType::kTHnSparseD, axesWithCent); registry.get(HIST("hBdtScoreVsMassVsPtVsYVsCentVsPtPion"))->Sumw2(); registry.get(HIST("hMassVsPtVsYVsCentVsPtPion"))->Sumw2(); + } else { + const AxisSpec thnAxisPromptScore{thnConfigAxisPromptScore, "BDT score prompt."}; + const AxisSpec thnAxisPtPion{thnConfigAxisPtPion, "Pt of Pion from Xic0."}; + std::vector const axesWithBdtWithoutCent = {thnAxisPromptScore, thnAxisMass, thnAxisPt, thnAxisY, thnAxisPtPion, thnConfigAxisNumPvContr}; + std::vector const axesWithoutCent = {thnAxisMass, thnAxisPt, thnAxisY, thnAxisPtPion, thnConfigAxisNumPvContr}; + registry.add("hBdtScoreVsMassVsPtVsYVsPtPion", "Thn for Xic0 candidates with BDT&Cent&pTpi", HistType::kTHnSparseD, axesWithBdtWithoutCent); + registry.add("hMassVsPtVsYVsPtPion", "Thn for Xic0 candidates with Cent&pTpi", HistType::kTHnSparseD, axesWithoutCent); + registry.get(HIST("hBdtScoreVsMassVsPtVsYVsPtPion"))->Sumw2(); + registry.get(HIST("hMassVsPtVsYVsPtPion"))->Sumw2(); } } @@ -192,22 +201,41 @@ struct HfTaskXic0ToXiPi { double const ptXic = RecoDecay::pt(candidate.pxCharmBaryon(), candidate.pyCharmBaryon()); double const ptPiFromXic = RecoDecay::pt(candidate.pxBachFromCharmBaryon(), candidate.pyBachFromCharmBaryon()); if constexpr (ApplyMl) { - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsCentVsPtPion"), - candidate.mlProbToXiPi()[0], - candidate.invMassCharmBaryon(), - ptXic, - yCharmBaryon, - centrality, - ptPiFromXic, - numPvContributors); + if constexpr (UseCentrality) { + registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsCentVsPtPion"), + candidate.mlProbToXiPi()[0], + candidate.invMassCharmBaryon(), + ptXic, + yCharmBaryon, + centrality, + ptPiFromXic, + numPvContributors); + } else { + registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsPtPion"), + candidate.mlProbToXiPi()[0], + candidate.invMassCharmBaryon(), + ptXic, + yCharmBaryon, + ptPiFromXic, + numPvContributors); + } } else { - registry.fill(HIST("hMassVsPtVsYVsCentVsPtPion"), - candidate.invMassCharmBaryon(), - ptXic, - yCharmBaryon, - centrality, - ptPiFromXic, - numPvContributors); + if constexpr (UseCentrality) { + registry.fill(HIST("hMassVsPtVsYVsCentVsPtPion"), + candidate.invMassCharmBaryon(), + ptXic, + yCharmBaryon, + centrality, + ptPiFromXic, + numPvContributors); + } else { + registry.fill(HIST("hMassVsPtVsYVsPtPion"), + candidate.invMassCharmBaryon(), + ptXic, + yCharmBaryon, + ptPiFromXic, + numPvContributors); + } } } @@ -348,7 +376,7 @@ struct HfTaskXic0ToXiPi { } } } - PROCESS_SWITCH(HfTaskXic0ToXiPi, processDataWithKFParticle, "process HfTaskXic0ToXiPi with KFParticle", true); + PROCESS_SWITCH(HfTaskXic0ToXiPi, processDataWithKFParticle, "process HfTaskXic0ToXiPi with KFParticle", false); void processDataWithDCAFitterMl(Xic0CandsMl const& candidates, CollisionsWithEvSels const& collisions) diff --git a/PWGHF/TableProducer/CMakeLists.txt b/PWGHF/TableProducer/CMakeLists.txt index 26b3400ca28..6a6dde9c3b6 100644 --- a/PWGHF/TableProducer/CMakeLists.txt +++ b/PWGHF/TableProducer/CMakeLists.txt @@ -90,6 +90,11 @@ o2physics_add_dpl_workflow(candidate-creator-xic0-omegac0 PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle O2Physics::AnalysisCCDB O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(candidate-creator-xic0-omegac0-qa + SOURCES candidateCreatorXic0Omegac0Qa.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(candidate-creator-xic-to-xi-pi-pi SOURCES candidateCreatorXicToXiPiPi.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle O2Physics::AnalysisCCDB O2Physics::EventFilteringUtils @@ -187,6 +192,11 @@ o2physics_add_dpl_workflow(candidate-selector-to-xi-pi PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(candidate-selector-to-xi-pi-qa + SOURCES candidateSelectorToXiPiQa.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(candidate-selector-xic-to-p-k-pi SOURCES candidateSelectorXicToPKPi.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore @@ -279,6 +289,11 @@ o2physics_add_dpl_workflow(tree-creator-to-xi-pi PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(tree-creator-to-xi-pi-qa + SOURCES treeCreatorToXiPiQa.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(tree-creator-xic0-to-xi-pi-kf SOURCES treeCreatorXic0ToXiPiKf.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGHF/TableProducer/candidateCreatorXic0Omegac0Qa.cxx b/PWGHF/TableProducer/candidateCreatorXic0Omegac0Qa.cxx new file mode 100644 index 00000000000..78eb469c92c --- /dev/null +++ b/PWGHF/TableProducer/candidateCreatorXic0Omegac0Qa.cxx @@ -0,0 +1,2128 @@ +// 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. + +/// \file candidateCreatorXic0Omegac0Qa.cxx +/// \brief Reconstruction of Xic0 and Xicp candiates with hadronic decay chain +/// +/// \author Jinhyun Park , Pusan National University +/// \author Krista Smith , Pusan National University + +#ifndef HomogeneousField +#define HomogeneousField // o2-linter: disable=name/macro (required by KFParticle) +#endif + +#include "PWGHF/Core/DecayChannelsLegacy.h" +#include "PWGHF/DataModel/AliasTables.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/DataModel/TrackIndexSkimmingTables.h" +#include "PWGHF/Utils/utilsBfieldCCDB.h" +#include "PWGHF/Utils/utilsEvSelHf.h" +#include "PWGHF/Utils/utilsTrkCandHf.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/mcCentrality.h" +#include "PWGLF/Utils/strangenessBuilderHelper.h" // -> Added to test removal of strangeness builder workflow + +#include "Common/Core/ZorroSummary.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/CollisionAssociationTables.h" +#include "Common/DataModel/EventSelection.h" +#include "Tools/KFparticle/KFUtilities.h" + +#include "CommonConstants/PhysicsConstants.h" +#include "DCAFitter/DCAFitterN.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/DCA.h" + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::hf_centrality; +using namespace o2::hf_occupancy; +using namespace o2::constants::physics; +using namespace o2::hf_evsel; + +struct HfCandidateCreatorXic0Omegac0Qa { + + // Cursor to fill tables + struct : ProducesGroup { + // Candidates created with DCAFitter + Produces rowCandToXiPi; + Produces rowCandToOmegaPi; + Produces rowCandToOmegaKa; + // Candidate created with KFParticle + Produces rowCandToXiPiKf; + Produces rowCandToOmegaPiKf; + Produces rowCandToOmegaKaKf; + } cursors; + + // Configurables + struct : ConfigurableGroup { + // Options for internal V0 building + // ...Initial values taken from PWGLF/Utiles/strangenessBuilderModule.h + // --------------------------------------------------------------------- + Configurable minCrossedRowsFromLF{"minCrossedRowsFromLF", 50, "minimun TPC crossed rows for daughter tracks. Used for internal V0 Building"}; + Configurable dcanegtopvFromLF{"dcanegtopvFromLF", .1, "DCV Neg to PV"}; + Configurable dcapostopvFromLF{"dcapostopvFromLF", .1, "DCV Pos To PV"}; + Configurable v0cospaFromLF{"v0cospaFromLF", 0.95, "V0 CosPA"}; + Configurable dcav0dauFromLF{"dcav0dauFromLF", 1.0, "DCA V0 Daughters"}; + Configurable v0radiusFromLF{"v0radiusFromLF", 0.9, "v0radius"}; + Configurable maxDaughterEtaFromLF{"maxDaughterEtaFromLF", 5.0, "Maximun daughter eta (in abs value)"}; + + // Options for internal cascade building + // ...Initial values taken from PWGLF/Utiles/strangenessBuilderModule.h + // -------------------------------------------------------------------- + Configurable dcabachtopvFromLF{"dcabachtopvFromLF", .1, "DCV Bach to PV"}; + Configurable cascradiusFromLF{"cascradiusFromLF", .1, "DCV Bach to PV"}; + Configurable casccospaFromLF{"casccospaFromLF", 0.95, "Cascade CosPA"}; + Configurable dcacascdauFromLF{"dcacascdauFromLF", 1.0, "DCA cascade daughters"}; + Configurable lambdaMassWindowFromLF{"lambdaMassWindowFromLF", 0.10, "Distance from Lambda mass(does not apply to KF path)"}; + + // Options for internal cascade building - KF Building specifics + // ...Initial values taken from PWGLF/Utiles/strangenessBuilderModule.h + // -------------------------------------------------------------------- + Configurable kfTuneForOmegaFromLF{"kfTuneForOmegaFromLF", false, "if enabled, take main cascade properties from omega fit instread of Xi fit(=default)"}; + Configurable kfConstructMethodFromLF{"kfConstructMethodFromLF", 2, "2 : Daughter particle masses stay fixed in construction process"}; + Configurable kfUseV0MassConstraintFromLF{"kfUseV0MassConstraintFromLF", false, "KF : Use Lambda mass constraint"}; + Configurable kfUseCascadeMassConstraintFromLF{"kfUseCascadeMassConstraintFromLF", false, "KF : Use Cascade mass constraint - WARNING : Not adequate for inv mass analysis of Xi"}; + Configurable kfDoDCAFitterPreMinimV0FromLF{"kfDoDCAFitterPreMinimV0FromLF", true, "KF : do DCAFitter pre-optimization before KF fit to include material correction for V0"}; + Configurable kfDoDCAFitterPreMinimCascFromLF{"kfDoDCAFitterPreMinimCascFromLF", true, "KF : do DCAFitter pre-optimization before KF fit to include material correction for Xi"}; + } LFConfigs; + + // For cascade building using LF strangeness builder + o2::pwglf::strangenessBuilderHelper straHelper; + + // Configurables + struct : ConfigurableGroup { + // Switch for filling histograms + // ----------------------------- + Configurable fillHistograms{"fillHistograms", true, "fill validation plots"}; + // Magnetic field setting from CCDB + // -------------------------------- + Configurable isRun2{"isRun2", false, "enable Run2 or Run3 GRP objects for magnetic field"}; + Configurable ccdbUrl{"ccdbUrl", "https://alice-ccdb.cern.ch", "url of the ccdb object"}; + Configurable ccdbPathLut{"ccdbPathLut", "GLO/Param/MatLUT", "Path for LUT parameterization"}; + Configurable ccdbPathGrp{"ccdbPathGrp", "GLO/GRP/GRP", "path of the group file (Run2)"}; + Configurable ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object (Run3)"}; + + // Cascade pre selection + // -------------------- + Configurable doCascadePreselection{"doCascadePreselection", true, "Use invariant mass and dcaXY cuts to preselect cascade candidates"}; + Configurable massToleranceCascade{"massToleranceCascade", 0.01, "Invariant mass tolerance for cascades"}; + Configurable dcaXYToPVCascadeMax{"dcaXYToPVCascadeMax", 3, "Max cascade DCA to PV in XY plane"}; + Configurable dcaV0DaughtersMax{"dcaV0DaughtersMax", 1.0, "Max DCA of V0 daughter"}; + Configurable dcaCascDaughtersMax{"dcaCascDaughtersMax", 1.0, "Max DCA of cascade daughter"}; + + // Options for DCAFitter + // --------------------- + Configurable propagateToPCA{"propagateToPCA", true, "Create tracks version propagated to PCA"}; + Configurable maxR{"maxR", 200., "Reject PCA's above this radius"}; + Configurable maxDZIni{"maxDZIni", 4., "Reject (if>0) PCA candidate if tracks DZ exceeds this threshold"}; + Configurable maxDXYIni{"maxDXYIni", 4., "Reject (if>0) PCA candidate if tracks DXY exceeds this threshold"}; + Configurable minParamChange{"minParamChange", 1.e-3, "Stop iteration if largest change of any X is smaller than this"}; + Configurable minRelChi2Change{"minRelChi2Change", 0.9, "Stop iteration if Chi2/Chi2old > this"}; + Configurable maxChi2{"maxChi2", 100, "Discard vertices with Chi2/Nprongs > this(or sum {DCAi^2}/Nprongs for abs. distance minimization)"}; + Configurable useAbsDCA{"useAbsDCA", true, "Minimise abs. distance rather than chi2"}; + Configurable useWeightedFinalPCA{"useWeightedFinalPCA", true, "Recalculate vertex position using track covariance, effective only if useAbsDCA is true"}; + + // Options for KFParticle + // ---------------------- + // V0 cuts + Configurable lambdaMassWindow{"lambdaMassWindow", 0.0075, "Distance from LambdaMass"}; + // For KF Particle operation + Configurable kfConstructMethod{"kfConstructMethod", 2, "KF Construct method"}; + Configurable kfUseV0MassConstraint{"kfUseV0MassConstraint", false, "KF: Use lambda mass constraint"}; + Configurable kfUseCascadeMassConstraint{"kfUseCascadeMassConstraint", false, "KF: Use lambda mass constraint"}; + + // Options for QA histogram binning + // ----------------------------- + + // For Cascade + Configurable nBinMassCasc{"nBinMassCasc", 1000, "nBinCascMass"}; + Configurable minMassCasc{"minMassCasc", 1.0, "xiMassMin"}; + Configurable maxMassCasc{"maxMassCasc", 2.0, "xiMassMax"}; + Configurable nBinPtCasc{"nBinPtCasc", 100, "nBinPtXi"}; + Configurable minPtCasc{"minPtCasc", 0.0, "minimun value of cascade"}; + Configurable maxPtCasc{"maxPtCasc", 20.0, "maximum value of cascade"}; + + // For Prong0 + Configurable nBinMassProng0{"nBinMassProng0", 100, "nBinMassPi"}; + Configurable minMassProng0{"minMassProng0", 0.0, "ptPiMin"}; + Configurable maxMassProng0{"maxMassProng0", 20.0, "ptPiMax"}; + Configurable nBinPtProng0{"nBinPtProng0", 100, "nBinPtProng0"}; + Configurable minPtProng0{"minPtProng0", 0.0, "ptPiMin"}; + Configurable maxPtProng0{"maxPtProng0", 20.0, "ptPiMax"}; + + // For Charm Baryon + Configurable nBinMassCharmBaryon{"nBinMassCharmBaryon", 3000, "nBinXic0Mass"}; + Configurable minMassCharmBaryon{"minMassCharmBaryon", 1.0, "xic0MassMin"}; + Configurable maxMassCharmBaryon{"maxMassCharmBaryon", 4.0, "xic0MassMax"}; + Configurable nBinPtCharmBaryon{"nBinPtCharmBaryon", 100, "nBinXic0Pt"}; + Configurable minPtCharmBaryon{"minPtCharmBaryon", 0.0, "xic0PtMin"}; + Configurable maxPtCharmBaryon{"maxPtCharmBaryon", 20.0, "xic0PtMax"}; + + // Etc + Configurable nBinCpa2Prong{"nBinCpa2Prong", 240, "nBinCpa2Prong"}; + Configurable cpa2ProngMin{"cpa2ProngMin", -1.2, "cpa2ProngMin"}; + Configurable cpa2ProngMax{"cpa2ProngMax", 1.2, "cpa2ProngMax"}; + Configurable nBinImpParXYXi{"nBinImpParXYXi", 30, "nBinImpParXYXi"}; + Configurable impParXYXiMin{"impParXYXiMin", -1.5, "impParXYXiMin"}; + Configurable impParXYXiMax{"impParXYXiMax", 1.5, "impParXYXiMax"}; + Configurable nBinImpParXYPi{"nBinImpParXYPi", 30, "nBinImpParXYPi"}; + Configurable impParXYPiMin{"impParXYPiMin", -1.5, "impParXYPiMin"}; + Configurable impParXYPiMax{"impParXYPiMax", 1.5, "impParXYPiMax"}; + } configs; + + // For magnetic field + Service ccdb; + o2::base::MatLayerCylSet* lut; + o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; + + // DCAFitter + o2::vertexing::DCAFitterN<2> df; + + int runNumber{0}; + double magneticField{0.}; + + enum CharmBaryonCandCounter { All = 0, + HfFlagPass, + CascReconstructed, + VertexFit }; + + // Table aliases + using SelectedCollisions = soa::Join; + using TracksWCovIU = soa::Join; + using TracksWCovDcaExtraPidPrPiKa = soa::Join; + using TracksWCovExtraPidIU = soa::Join; + + HistogramRegistry registry{"hists"}; + OutputObj zorroSummary{"zorroSummary"}; + HfEventSelection hfEvSel; + + // PDG Id of daughter tracks & V0s & cascades & charm baryons - Used in KFParticle + int pdgIdOfV0DauPos, pdgIdOfV0DauNeg, pdgIdOfBach, pdgIdOfCharmBach; + int pdgIdOfAntiV0DauPos, pdgIdOfAntiV0DauNeg, pdgIdOfAntiBach, pdgIdOfAntiCharmBach; + int pdgIdOfV0, pdgIdOfCascade, pdgIdOfCharmBaryon; + + // Track PID - Used in DCAFitter + int trackPidOfCascade; + + // Mass of daughter tracks & V0s & cascades & charm baryons; + int massOfV0DauPos, massOfV0DauNeg, massOfBach, massOfCharmBach; + int massOfV0, massOfCascade, massOfCharmBaryon; + + // Pointer of histograms for QA + std::shared_ptr hInvMassCharmBaryonToXiPi, hInvMassCharmBaryonToOmegaPi, hInvMassCharmBaryonToOmegaKa; + std::shared_ptr hCandidateCounterToXiPi, hCandidateCounterToOmegaPi, hCandidateCounterToOmegaKa; + + void init(InitContext const&) + { + std::vector processesToXiPiDca{doprocessToXiPiWithDCAFitterNoCent, /*doprocessToXiPiWithDCAFitterNoCentWithTrackedCasc,*/ doprocessToXiPiWithDCAFitterCentFT0C, doprocessToXiPiWithDCAFitterCentFT0M}; + std::vector processesToOmegaPiDca{doprocessToOmegaPiWithDCAFitterNoCent, doprocessToOmegaPiWithDCAFitterCentFT0C, doprocessToOmegaPiWithDCAFitterCentFT0M}; + std::vector processesToOmegaKaDca{doprocessToOmegaKaWithDCAFitterNoCent, doprocessToOmegaKaWithDCAFitterCentFT0C, doprocessToOmegaKaWithDCAFitterCentFT0M}; + std::vector processesToXiPiKf{doprocessToXiPiWithKFParticleNoCent, doprocessToXiPiWithKFParticleCentFT0C, doprocessToXiPiWithKFParticleCentFT0M}; + std::vector processesToOmegaPiKf{doprocessToOmegaPiWithKFParticleNoCent, doprocessToOmegaPiWithKFParticleCentFT0C, doprocessToOmegaPiWithKFParticleCentFT0M}; + std::vector processesToOmegaKaKf{doprocessToOmegaKaWithKFParticleNoCent, doprocessToOmegaKaWithKFParticleCentFT0C, doprocessToOmegaKaWithKFParticleCentFT0M}; + std::vector processesCollMonitoring{doprocessCollisionsNoCent, doprocessCollisionsCentFT0C, doprocessCollisionsCentFT0M}; + + int xipiEnabledDca = std::accumulate(processesToXiPiDca.begin(), processesToXiPiDca.end(), 0); + int xipiEnabledKf = std::accumulate(processesToXiPiKf.begin(), processesToXiPiKf.end(), 0); + int omegapiEnabledDca = std::accumulate(processesToOmegaPiDca.begin(), processesToOmegaPiDca.end(), 0); + int omegapiEnabledKf = std::accumulate(processesToOmegaPiKf.begin(), processesToOmegaPiKf.end(), 0); + int omegakaEnabledDca = std::accumulate(processesToOmegaKaDca.begin(), processesToOmegaKaDca.end(), 0); + int omegakaEnabledKf = std::accumulate(processesToOmegaKaKf.begin(), processesToOmegaKaKf.end(), 0); + + // Exit if workflow is not configured correctly - More than one process function enabled for candidate to XiPi + if ((xipiEnabledDca > 0) && (xipiEnabledKf > 0)) { + LOGP(fatal, "More than one process function enabled for candidte decaying to xi pi"); + } + + // Exit if workflow is not configured correctly - More than one process function enabled for candidate to OmegaPi + if ((omegapiEnabledDca > 0) && (omegapiEnabledKf > 0)) { + LOGP(fatal, "More than one process function enabled for candidte decaying to omega pi"); + } + + // Exit if workflow is not configured correctly - More than one process function enabled for candidate to OmegaKa + if ((omegakaEnabledDca > 0) && (omegakaEnabledKf > 0)) { + LOGP(fatal, "More than one process function enabled for candidte decaying to omega ka"); + } + + // Exit if workflow is not configured correctly - More than one process enabled for collision monitoring + if (std::accumulate(processesCollMonitoring.begin(), processesCollMonitoring.end(), 0) > 1) { + LOGP(fatal, "More than one process fucntion for CollMonitoring was enabled. Please choose only one process function"); + } + + // Assign pdg & mass hypothesis for each decay channel + if (xipiEnabledDca || xipiEnabledKf) { + pdgIdOfV0DauPos = kProton; + pdgIdOfV0DauNeg = kPiMinus; + pdgIdOfBach = kPiMinus; + pdgIdOfCharmBach = kPiPlus; + + pdgIdOfAntiV0DauPos = kPiPlus; + pdgIdOfAntiV0DauNeg = kProton; + pdgIdOfAntiBach = kPiPlus; + pdgIdOfAntiCharmBach = kPiMinus; + + pdgIdOfV0 = kLambda0; + pdgIdOfCascade = kXiMinus; + pdgIdOfCharmBaryon = kXiC0; + + trackPidOfCascade = o2::track::PID::XiMinus; + massOfCharmBach = o2::constants::physics::MassPiPlus; + massOfV0 = o2::constants::physics::MassLambda; + massOfCascade = o2::constants::physics::MassXiMinus; + } else if (omegapiEnabledDca || omegapiEnabledKf) { + pdgIdOfV0DauPos = kProton; + pdgIdOfV0DauNeg = kPiMinus; + pdgIdOfBach = kKMinus; + pdgIdOfCharmBach = kPiPlus; + + pdgIdOfAntiV0DauPos = kPiPlus; + pdgIdOfAntiV0DauNeg = kProton; + pdgIdOfAntiBach = kKPlus; + pdgIdOfAntiCharmBach = kPiMinus; + + pdgIdOfV0 = kLambda0; + pdgIdOfCascade = kOmegaMinus; + pdgIdOfCharmBaryon = kOmegaC0; + + trackPidOfCascade = o2::track::PID::OmegaMinus; + massOfCharmBach = o2::constants::physics::MassPiPlus; + massOfV0 = o2::constants::physics::MassLambda; + massOfCascade = o2::constants::physics::MassOmegaMinus; + } else if (omegakaEnabledDca || omegakaEnabledKf) { + pdgIdOfV0DauPos = kProton; + pdgIdOfV0DauNeg = kPiMinus; + pdgIdOfBach = kKMinus; + pdgIdOfCharmBach = kKPlus; + + pdgIdOfAntiV0DauPos = kPiPlus; + pdgIdOfAntiV0DauNeg = kProton; + pdgIdOfAntiBach = kKPlus; + pdgIdOfAntiCharmBach = kKMinus; + + pdgIdOfV0 = kLambda0; + pdgIdOfCascade = kOmegaMinus; + pdgIdOfCharmBaryon = kOmegaC0; + + trackPidOfCascade = o2::track::PID::OmegaMinus; + massOfCharmBach = o2::constants::physics::MassKPlus; + massOfV0 = o2::constants::physics::MassLambda; + massOfCascade = o2::constants::physics::MassOmegaMinus; + } + LOGF(info, "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"); + LOGF(info, "PDG ID of V0 positive daughter: %d", pdgIdOfV0DauPos); + LOGF(info, "PDG ID of V0 negative daughter: %d", pdgIdOfV0DauNeg); + LOGF(info, "PDG ID of Bachelor: %d", pdgIdOfBach); + LOGF(info, "PDG ID of Charm Bachelor: %d", pdgIdOfCharmBach); + LOGF(info, "----------"); + LOGF(info, "PDG ID of anti V0 positive daughter: %d", pdgIdOfAntiV0DauPos); + LOGF(info, "PDG ID of anti V0 negative daughter: %d", pdgIdOfAntiV0DauNeg); + LOGF(info, "PDG ID of anti Bachelor: %d", pdgIdOfAntiBach); + LOGF(info, "PDG ID of anti Charm Bachelor: %d", pdgIdOfAntiCharmBach); + LOGF(info, "----------"); + LOGF(info, "PDG ID of V0: %d", pdgIdOfV0); + LOGF(info, "PDG ID of Cascade: %d", pdgIdOfCascade); + LOGF(info, "PDG ID of Charm Baryon: %d", pdgIdOfCharmBaryon); + LOGF(info, "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"); + + // Add histogram to indicate which sv method was used + registry.add("hVertexerType", "Use KF or DCAFitterN;Vertexer type;entries", {kTH1F, {{2, 0.0, 2.0}}}); + registry.get(HIST("hVertexerType"))->GetXaxis()->SetBinLabel(1 + aod::hf_cand::VertexerType::DCAFitter, "DCAFitter"); + registry.get(HIST("hVertexerType"))->GetXaxis()->SetBinLabel(1 + aod::hf_cand::VertexerType::KfParticle, "KFParticle"); + + // initialize ccdb + // --------------- + ccdb->setURL(configs.ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(configs.ccdbPathLut)); + straHelper.lut = lut; + runNumber = 0; + + // Initilization for strangeness builder helper + // -------------------------------------------- + + // Settings for internal V0 building + straHelper.v0selections.minCrossedRows = LFConfigs.minCrossedRowsFromLF; + straHelper.v0selections.dcanegtopv = LFConfigs.dcanegtopvFromLF; + straHelper.v0selections.dcapostopv = LFConfigs.dcapostopvFromLF; + straHelper.v0selections.v0cospa = LFConfigs.v0cospaFromLF; + straHelper.v0selections.dcav0dau = LFConfigs.dcav0dauFromLF; + straHelper.v0selections.v0radius = LFConfigs.v0radiusFromLF; + straHelper.v0selections.maxDaughterEta = LFConfigs.maxDaughterEtaFromLF; + + // Settings for internal Cascade building + straHelper.cascadeselections.minCrossedRows = LFConfigs.minCrossedRowsFromLF; + straHelper.cascadeselections.dcabachtopv = LFConfigs.dcabachtopvFromLF; + straHelper.cascadeselections.cascradius = LFConfigs.cascradiusFromLF; + straHelper.cascadeselections.casccospa = LFConfigs.casccospaFromLF; + straHelper.cascadeselections.dcacascdau = LFConfigs.dcacascdauFromLF; + straHelper.cascadeselections.lambdaMassWindow = LFConfigs.lambdaMassWindowFromLF; + straHelper.cascadeselections.maxDaughterEta = LFConfigs.maxDaughterEtaFromLF; + + // Fitter setting + straHelper.fitter.setPropagateToPCA(configs.propagateToPCA); + straHelper.fitter.setMaxR(configs.maxR); + straHelper.fitter.setMaxDZIni(configs.maxDZIni); + straHelper.fitter.setMinParamChange(configs.minParamChange); + straHelper.fitter.setUseAbsDCA(configs.useAbsDCA); + straHelper.fitter.setWeightedFinalPCA(configs.useWeightedFinalPCA); + + // Extra initialization for DCAFitter + // ---------------------------------- + if (xipiEnabledDca == 1 || omegapiEnabledDca == 1 || omegakaEnabledDca == 1) { + registry.get(HIST("hVertexerType"))->Fill(aod::hf_cand::VertexerType::DCAFitter); + df.setPropagateToPCA(configs.propagateToPCA); + df.setMaxR(configs.maxR); + df.setMaxDZIni(configs.maxDZIni); + df.setMaxDXYIni(configs.maxDXYIni); + df.setMinParamChange(configs.minParamChange); + df.setMinRelChi2Change(configs.minRelChi2Change); + df.setMaxChi2(configs.maxChi2); + df.setUseAbsDCA(configs.useAbsDCA); + df.setWeightedFinalPCA(configs.useWeightedFinalPCA); + } + + // Extra initialization for KFParticle + // ----------------------------------- + if (xipiEnabledKf == 1 || omegapiEnabledKf == 1 || omegakaEnabledKf == 1) { + registry.get(HIST("hVertexerType"))->Fill(aod::hf_cand::VertexerType::KfParticle); + } + + // initailize HF event selection helper + // ------------------------------------ + hfEvSel.init(registry, &zorroSummary); + + // Histograms for QA + // ----------------- + registry.add("ReconstructedDecayChannel", "DecayChannel", {kTH1F, {{3, 0.0, 3.0}}}); + registry.get(HIST("ReconstructedDecayChannel"))->GetXaxis()->SetBinLabel(1 + hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi, "To #Xi #pi"); + registry.get(HIST("ReconstructedDecayChannel"))->GetXaxis()->SetBinLabel(1 + hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaPi, "To #Omega #pi"); + registry.get(HIST("ReconstructedDecayChannel"))->GetXaxis()->SetBinLabel(1 + hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaK, "To #Omega K"); + + if (xipiEnabledDca != 0 || xipiEnabledKf != 0) { + hInvMassCharmBaryonToXiPi = registry.add("hInvMassCharmBaryonToXiPi", "Charm baryon invariant mass - #Xi #pi decay;inv. mass (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{configs.nBinMassCharmBaryon, configs.minMassCharmBaryon, configs.maxMassCharmBaryon}}}); + hCandidateCounterToXiPi = registry.add("hCandidateCounterToXiPi", "Candidate counter wrt derived data - #Xi #pi decay;status;entries", {HistType::kTH1D, {{4, -0.5, 3.5}}}); + registry.get(HIST("hCandidateCounterToXiPi"))->GetXaxis()->SetBinLabel(1 + All, "Total"); + registry.get(HIST("hCandidateCounterToXiPi"))->GetXaxis()->SetBinLabel(1 + HfFlagPass, "HfFlagPass"); + registry.get(HIST("hCandidateCounterToXiPi"))->GetXaxis()->SetBinLabel(1 + CascReconstructed, "CascReconstructed"); + registry.get(HIST("hCandidateCounterToXiPi"))->GetXaxis()->SetBinLabel(1 + VertexFit, "VertexFit"); + registry.get(HIST("ReconstructedDecayChannel"))->Fill(hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi); + } + + if (omegapiEnabledDca != 0 || omegapiEnabledKf != 0) { + hInvMassCharmBaryonToOmegaPi = registry.add("hInvMassCharmBaryonToOmegaPi", "Charm baryon invariant mass - #Omega #pi decay;inv. mass (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{configs.nBinMassCharmBaryon, configs.minMassCharmBaryon, configs.maxMassCharmBaryon}}}); + hCandidateCounterToOmegaPi = registry.add("hCandidateCounterToOmegaPi", "Candidate counter wrt derived data - #Omega #pi decay;status;entries", {HistType::kTH1D, {{4, -0.5, 3.5}}}); + registry.get(HIST("hCandidateCounterToOmegaPi"))->GetXaxis()->SetBinLabel(1 + All, "Total"); + registry.get(HIST("hCandidateCounterToOmegaPi"))->GetXaxis()->SetBinLabel(1 + HfFlagPass, "HfFlagPass"); + registry.get(HIST("hCandidateCounterToOmegaPi"))->GetXaxis()->SetBinLabel(1 + CascReconstructed, "CascReconstructed"); + registry.get(HIST("hCandidateCounterToOmegaPi"))->GetXaxis()->SetBinLabel(1 + VertexFit, "VertexFit"); + registry.get(HIST("ReconstructedDecayChannel"))->Fill(hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaPi); + } + + if (omegakaEnabledDca != 0 || omegakaEnabledKf != 0) { + hInvMassCharmBaryonToOmegaKa = registry.add("hInvMassCharmBaryonToOmegaKa", "Charm baryon invariant mass - #Omega K decay;inv. mass (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{configs.nBinMassCharmBaryon, configs.minMassCharmBaryon, configs.maxMassCharmBaryon}}}); + hCandidateCounterToOmegaKa = registry.add("hCandidateCounterToOmegaKa", "Candidate counter wrt derived data - #Omega K decay;status;entries", {HistType::kTH1D, {{4, -0.5, 3.5}}}); + registry.get(HIST("hCandidateCounterToOmegaKa"))->GetXaxis()->SetBinLabel(1 + All, "Total"); + registry.get(HIST("hCandidateCounterToOmegaKa"))->GetXaxis()->SetBinLabel(1 + HfFlagPass, "HfFlagPass"); + registry.get(HIST("hCandidateCounterToOmegaKa"))->GetXaxis()->SetBinLabel(1 + CascReconstructed, "CascReconstructed"); + registry.get(HIST("hCandidateCounterToOmegaKa"))->GetXaxis()->SetBinLabel(1 + VertexFit, "VertexFit"); + registry.get(HIST("ReconstructedDecayChannel"))->Fill(hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaK); + } + + registry.add("hCascMass", "Inv mass of reconstructed cascade;Inv mass;Entries", {HistType::kTH1F, {{configs.nBinMassCasc, configs.minMassCasc, configs.maxMassCasc}}}); + registry.add("hCascPt", "Pt of reconstructed cascade;pT;Entries", {HistType::kTH1F, {{configs.nBinPtCasc, configs.minPtCasc, configs.maxPtCasc}}}); + + } // end of initialization + + //////////////////////////////////////////////////////////// + // // + // Candidate reconstruction with DCAFitter // + // // + //////////////////////////////////////////////////////////// + + // template function for running charm baryon reconstruction with DCAFitter + /// \brief centEstimator is for different centrality estimators + /// \brief decayChannel is for different decay channels. 0 for XiczeroOmegaczeroToXiPi, 1 for OmegaczeroToOmegaPi, 2 for OmegaczeroToOmeagaK + /// \brief Colls is for collision tables joined with different centraltiy estimators + /// \brief Hist is for QA histograms + template + void runCreatorWithDCAFitter(Colls const&, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const&, // -> Internal cascade building + aod::V0s const&, // -> Internal v0 building + aod::TracksWCovDca const& tracks, + TracksWCovIU const& lfTracks, + aod::BCsWithTimestamps const&, + Hist& hInvMassCharmBaryon, + Hist& hCandCounter) + { + // arrays which holds values for different decay channel reconstruction + + // Loop over candidate + for (auto const& cand : candidates) { + + // Fill cascandidates before selection + if (configs.fillHistograms) { + hCandCounter->Fill(All); + } + + // Apply hfflag selection for different candidate reconstruction + if (!TESTBIT(cand.hfflag(), decayChannel)) { + continue; + } else { + if (configs.fillHistograms) { + hCandCounter->Fill(HfFlagPass); + } + } + + // Event selection + auto collision = cand.collision_as(); + float centrality{-1.f}; + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); + if (rejectionMask != 0) { // None of the event selection satisfied -> Reject this candidate + continue; + } + + //------------------------------Set Magnetic field------------------------------ + auto bc = collision.template bc_as(); + if (runNumber != bc.runNumber()) { + LOG(info) << ">>>>>>>>>> Current run Number : " << runNumber; + initCCDB(bc, runNumber, ccdb, configs.isRun2 ? configs.ccdbPathGrp : configs.ccdbPathGrpMag, lut, configs.isRun2); + magneticField = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << ">>>>>>>>>> Magnetic field: " << magneticField; + } + straHelper.fitter.setBz(magneticField); // -> Magnetic field setting for internal cascade building + df.setBz(magneticField); // -> Magnetic field setting for charm baryon building + + //------------------Intenal Cascade building------------------ + auto cascAodElement = cand.cascade_as(); + auto v0AodElement = cascAodElement.v0_as(); + auto posTrack = lfTracks.rawIteratorAt(v0AodElement.posTrackId()); + auto negTrack = lfTracks.rawIteratorAt(v0AodElement.negTrackId()); + auto bachTrack = lfTracks.rawIteratorAt(cascAodElement.bachelorId()); + + // Make cascade starting from V0 + // If success, fill Cascade and V0 information for reconstruction + if (!straHelper.buildCascadeCandidate(collision.globalIndex(), + collision.posX(), collision.posY(), collision.posZ(), + posTrack, + negTrack, + bachTrack, + false, // calculateBachelorBaryonVariable + false, // useCascadeMomentumAtPV + true)) { + // LOG(info) << "!This cascade cannot be rebuilt(cascade ID : " << cand.cascadeId() << "/ collision ID : " << collision.globalIndex(); + continue; + } else { + float storeMass = (decayChannel == 0) ? straHelper.cascade.massXi : straHelper.cascade.massOmega; + float storePt = RecoDecay::pt(straHelper.cascade.cascadeMomentum); + registry.fill(HIST("hCascMass"), storeMass); + registry.fill(HIST("hCascPt"), storePt); + if (configs.fillHistograms) { + hCandCounter->Fill(CascReconstructed); + } + } + + //------------------------------Info of V0 and cascade------------------------------ + // V0 quantities from LF strangeness builder + std::array vertexV0 = {straHelper.cascade.v0Position[0], straHelper.cascade.v0Position[1], straHelper.cascade.v0Position[2]}; + std::array pVecV0 = {straHelper.cascade.v0Momentum[0], straHelper.cascade.v0Momentum[1], straHelper.cascade.v0Momentum[2]}; + std::array pVecV0DauPos = {straHelper.cascade.positiveMomentum[0], straHelper.cascade.positiveMomentum[1], straHelper.cascade.positiveMomentum[2]}; + std::array pVecV0DauNeg = {straHelper.cascade.negativeMomentum[0], straHelper.cascade.negativeMomentum[1], straHelper.cascade.negativeMomentum[2]}; + + // pseudo rapidity - V0 daughters + // float pseudorapV0Dau0 = RecoDecay::eta(pVecV0DauPos); + // float pseudorapV0Dau1 = RecoDecay::eta(pVecV0DauNeg); + + // Cascade quantities from LF strangeness builder + int chargeCasc = straHelper.cascade.charge; + std::array vertexCasc(straHelper.cascade.cascadePosition); + std::array const pVecCasc(straHelper.cascade.cascadeMomentum); + std::array covCasc = {0.}; + constexpr int NumCovElement = 6; + constexpr int MomInd[NumCovElement] = {9, 13, 14, 18, 19, 20}; + for (int i = 0; i < NumCovElement; i++) { + covCasc[MomInd[i]] = straHelper.cascade.covariance[MomInd[i]]; + covCasc[i] = straHelper.cascade.covariance[i]; + } + + // pseudo rapidity - cascade bachelor + // float pseudorapCascBachelor = RecoDecay::eta(straHelper.cascade.bachelorMomentum); + + //------------------------------Create cascade track------------------------------ + + o2::track::TrackParCov trackCasc; + if (chargeCasc < 0) { // Xi- or Omega- + trackCasc = o2::track::TrackParCov(vertexCasc, pVecCasc, covCasc, -1, true); + } else if (chargeCasc > 0) { // Xi+ or Omega+ + trackCasc = o2::track::TrackParCov(vertexCasc, pVecCasc, covCasc, 1, true); + } else { + continue; + } + + trackCasc.setAbsCharge(1); + trackCasc.setPID(trackPidOfCascade); + + //------------------------------Fit SV & Create Charm Baryon track------------------------------ + + // Perform secondary vertex fitting + auto trackCharmBachelor = tracks.rawIteratorAt(cand.prong0Id()); + auto trackParCovCharmBachelor = getTrackParCov(trackCharmBachelor); + try { + if (df.process(trackCasc, trackParCovCharmBachelor) == 0) { + continue; + } + } catch (std::runtime_error& e) { + LOG(error) << "Execption caught in charm DCA Fitter process call: " << e.what(); + continue; + } + + if (configs.fillHistograms) { + hCandCounter->Fill(VertexFit); + } + + //------------------------------Calculate physical properties----------------------------- + + // get track momenta + std::array pVecCascAsD, pVecCharmBachAsD; + df.getTrack(0).getPxPyPzGlo(pVecCascAsD); + df.getTrack(1).getPxPyPzGlo(pVecCharmBachAsD); + + std::array pVecCharmBaryon = {pVecCascAsD[0] + pVecCharmBachAsD[0], pVecCascAsD[1] + pVecCharmBachAsD[1], pVecCascAsD[2] + pVecCharmBachAsD[2]}; + std::array pVecBach = {straHelper.cascade.bachelorMomentum[0], straHelper.cascade.bachelorMomentum[1], straHelper.cascade.bachelorMomentum[2]}; + + // get PV Properties + auto primaryVertex = getPrimaryVertex(collision); + auto covMatrixPV = primaryVertex.getCov(); + std::array pvCoord = {collision.posX(), collision.posY(), collision.posZ()}; + + // get SV Properties + auto const& secondaryVertex = df.getPCACandidate(); + auto covMatrixSV = df.calcPCACovMatrixFlat(); + + // DCAxy and DCAz. Computed with propagatToDCABxByBz method + auto trackParCovV0DauPos = getTrackParCov(posTrack); + auto trackParCovV0DauNeg = getTrackParCov(negTrack); + auto trackParCovBach = getTrackParCov(bachTrack); + + o2::dataformats::DCA impactParameterV0DauPos, impactParameterV0DauNeg, impactParameterBach; + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovV0DauPos, 2.f, matCorr, &impactParameterV0DauPos); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovV0DauNeg, 2.f, matCorr, &impactParameterV0DauNeg); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovBach, 2.f, matCorr, &impactParameterBach); + + float dcaxyV0DauPos = impactParameterV0DauPos.getY(); + float dcaxyV0DauNeg = impactParameterV0DauNeg.getY(); + float dcaxyBach = impactParameterBach.getY(); + float dcazV0DauPos = impactParameterV0DauPos.getZ(); + float dcazV0DauNeg = impactParameterV0DauNeg.getZ(); + float dcazBach = impactParameterBach.getZ(); + + // get impact parameter. Compute with propagateToDCABxByBz method + o2::dataformats::DCA impactParameterCasc, impactParameterCharmBach; + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackCasc, 2.f, matCorr, &impactParameterCasc); // trackCasc is TrackParCov object + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovCharmBachelor, 2.f, matCorr, &impactParameterCharmBach); + // float impactParCharmBachFromCharmBayonXY = impactParameterCharmBach.getY(); + // float impactParCharmBachFromCharmBayonZ = impactParameterCharmBach.getZ(); + + // get v0 invariant mass - from LF Table + float mLambda = straHelper.v0.massLambda; // from LF Table + + // get Casc mass - from LF Table + float mCasc = 0.; + if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi) { + mCasc = straHelper.cascade.massXi; + } else if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaPi) { + mCasc = straHelper.cascade.massOmega; + } else if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaK) { + mCasc = straHelper.cascade.massOmega; + } + + // get Charm baryon invariant mass + auto arrMomenta = std::array{pVecCascAsD, pVecCharmBachAsD}; + float massCharmBaryonCand = RecoDecay::m(arrMomenta, std::array{massOfCascade, massOfCharmBach}); + if (configs.fillHistograms) { + hInvMassCharmBaryon->Fill(massCharmBaryonCand); + } + + // calculate cosine of pointing angle + std::array vtxCoordCharmBaryon = df.getPCACandidatePos(); + float cpaV0 = RecoDecay::cpa(pvCoord, vertexV0, pVecV0); + float cpaCasc = RecoDecay::cpa(pvCoord, vertexCasc, pVecCasc); + float cpaCharmBaryon = RecoDecay::cpa(pvCoord, vtxCoordCharmBaryon, pVecCharmBaryon); + float cpaxyV0 = RecoDecay::cpaXY(pvCoord, vertexV0, pVecV0); + float cpaxyCasc = RecoDecay::cpaXY(pvCoord, vertexCasc, pVecCasc); + float cpaxyCharmBaryon = RecoDecay::cpaXY(pvCoord, vtxCoordCharmBaryon, pVecCharmBaryon); + + // calculate decay length + float decLenV0 = RecoDecay::distance(vertexCasc, vertexV0); + float decLenCasc = RecoDecay::distance(vtxCoordCharmBaryon, vertexCasc); + float decLenCharmBaryon = RecoDecay::distance(pvCoord, vtxCoordCharmBaryon); + float phi, theta; + getPointDirection(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta); + auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixSV, phi, theta)); + auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixSV, phi, 0.)); + + // calcuate ctau + float ctV0 = RecoDecay::ct(pVecV0, decLenV0, MassLambda0); + float ctCasc = RecoDecay::ct(pVecCasc, decLenCasc, massOfCascade); + float ctOmegac0 = RecoDecay::ct(pVecCharmBaryon, decLenCharmBaryon, MassOmegaC0); + float ctXic0 = RecoDecay::ct(pVecCharmBaryon, decLenCharmBaryon, MassXiC0); + + // get eta + float etaV0DauPos = posTrack.eta(); + float etaV0DauNeg = negTrack.eta(); + float etaBach = bachTrack.eta(); + float etaCharmBach = trackCharmBachelor.eta(); + float etaV0 = RecoDecay::eta(pVecV0); + float etaCasc = RecoDecay::eta(pVecCasc); + float etaCharmBaryon = RecoDecay::eta(pVecCharmBaryon); + + // DCA between daughters + float dcaV0Dau = straHelper.cascade.v0DaughterDCA; + float dcaCascDau = straHelper.cascade.cascadeDaughterDCA; + float dcaCharmBaryonDau = std::sqrt(df.getChi2AtPCACandidate()); + + //------------------------------Fill QA histograms----------------------------- + + //------------------------------Fill the table----------------------------- + if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi) { + cursors.rowCandToXiPi(collision.globalIndex(), + pvCoord[0], pvCoord[1], pvCoord[2], + secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], + vertexCasc[0], vertexCasc[1], vertexCasc[2], + vertexV0[0], vertexV0[1], vertexV0[2], + bachTrack.sign(), + covMatrixSV[0], covMatrixSV[1], covMatrixSV[2], covMatrixSV[3], covMatrixSV[4], covMatrixSV[5], + pVecCharmBaryon[0], pVecCharmBaryon[1], pVecCharmBaryon[2], + pVecCascAsD[0], pVecCascAsD[1], pVecCascAsD[2], + pVecCharmBachAsD[0], pVecCharmBachAsD[1], pVecCharmBachAsD[2], + pVecV0[0], pVecV0[1], pVecV0[2], + pVecBach[0], pVecBach[1], pVecBach[2], + pVecV0DauPos[0], pVecV0DauPos[1], pVecV0DauPos[2], + pVecV0DauNeg[0], pVecV0DauNeg[1], pVecV0DauNeg[2], + impactParameterCasc.getY(), impactParameterCharmBach.getY(), + impactParameterCasc.getZ(), impactParameterCharmBach.getZ(), + std::sqrt(impactParameterCasc.getSigmaY2()), std::sqrt(impactParameterCharmBach.getSigmaY2()), + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), + cand.cascadeId(), trackCharmBachelor.globalIndex(), cand.prong0Id(), + mLambda, mCasc, massCharmBaryonCand, + cpaV0, cpaCharmBaryon, cpaCasc, + cpaxyV0, cpaxyCharmBaryon, cpaxyCasc, + ctOmegac0, ctCasc, ctV0, ctXic0, + etaV0DauPos, etaV0DauNeg, etaBach, etaCharmBach, etaCharmBaryon, etaCasc, etaV0, + dcaxyV0DauPos, dcaxyV0DauNeg, dcaxyBach, + dcazV0DauPos, dcazV0DauNeg, dcazBach, + dcaCascDau, dcaV0Dau, dcaCharmBaryonDau, + decLenCharmBaryon, decLenCasc, decLenV0, errorDecayLength, errorDecayLengthXY); + } else if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaPi) { + cursors.rowCandToOmegaPi(collision.globalIndex(), + pvCoord[0], pvCoord[1], pvCoord[2], + secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], + vertexCasc[0], vertexCasc[1], vertexCasc[2], + vertexV0[0], vertexV0[1], vertexV0[2], + bachTrack.sign(), + covMatrixSV[0], covMatrixSV[1], covMatrixSV[2], covMatrixSV[3], covMatrixSV[4], covMatrixSV[5], + pVecCharmBaryon[0], pVecCharmBaryon[1], pVecCharmBaryon[2], + pVecCascAsD[0], pVecCascAsD[1], pVecCascAsD[2], + pVecCharmBachAsD[0], pVecCharmBachAsD[1], pVecCharmBachAsD[2], + pVecV0[0], pVecV0[1], pVecV0[2], + pVecBach[0], pVecBach[1], pVecBach[2], + pVecV0DauPos[0], pVecV0DauPos[1], pVecV0DauPos[2], + pVecV0DauNeg[0], pVecV0DauNeg[1], pVecV0DauNeg[2], + impactParameterCasc.getY(), impactParameterCharmBach.getY(), + impactParameterCasc.getZ(), impactParameterCharmBach.getZ(), + std::sqrt(impactParameterCasc.getSigmaY2()), std::sqrt(impactParameterCharmBach.getSigmaY2()), + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), + cand.cascadeId(), trackCharmBachelor.globalIndex(), cand.prong0Id(), + mLambda, mCasc, massCharmBaryonCand, + cpaV0, cpaCharmBaryon, cpaCasc, + cpaxyV0, cpaxyCharmBaryon, cpaxyCasc, + ctOmegac0, ctCasc, ctV0, + etaV0DauPos, etaV0DauNeg, etaBach, etaCharmBach, etaCharmBaryon, etaCasc, etaV0, + dcaxyV0DauPos, dcaxyV0DauNeg, dcaxyBach, + dcazV0DauPos, dcazV0DauNeg, dcazBach, + dcaCascDau, dcaV0Dau, dcaCharmBaryonDau, + decLenCharmBaryon, decLenCasc, decLenV0, errorDecayLength, errorDecayLengthXY, cand.hfflag()); + } else { + cursors.rowCandToOmegaKa(collision.globalIndex(), + pvCoord[0], pvCoord[1], pvCoord[2], + secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], + vertexCasc[0], vertexCasc[1], vertexCasc[2], + vertexV0[0], vertexV0[1], vertexV0[2], + bachTrack.sign(), + covMatrixSV[0], covMatrixSV[1], covMatrixSV[2], covMatrixSV[3], covMatrixSV[4], covMatrixSV[5], + pVecCharmBaryon[0], pVecCharmBaryon[1], pVecCharmBaryon[2], + pVecCascAsD[0], pVecCascAsD[1], pVecCascAsD[2], + pVecCharmBachAsD[0], pVecCharmBachAsD[1], pVecCharmBachAsD[2], + pVecV0[0], pVecV0[1], pVecV0[2], + pVecBach[0], pVecBach[1], pVecBach[2], + pVecV0DauPos[0], pVecV0DauPos[1], pVecV0DauPos[2], + pVecV0DauNeg[0], pVecV0DauNeg[1], pVecV0DauNeg[2], + impactParameterCasc.getY(), impactParameterCharmBach.getY(), + impactParameterCasc.getZ(), impactParameterCharmBach.getZ(), + std::sqrt(impactParameterCasc.getSigmaY2()), std::sqrt(impactParameterCharmBach.getSigmaY2()), + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), + cand.cascadeId(), trackCharmBachelor.globalIndex(), cand.prong0Id(), + mLambda, mCasc, massCharmBaryonCand, + cpaV0, cpaCharmBaryon, cpaCasc, + cpaxyV0, cpaxyCharmBaryon, cpaxyCasc, + ctOmegac0, ctCasc, ctV0, + etaV0DauPos, etaV0DauNeg, etaBach, etaCharmBach, etaCharmBaryon, etaCasc, etaV0, + dcaxyV0DauPos, dcaxyV0DauNeg, dcaxyBach, + dcazV0DauPos, dcazV0DauNeg, dcazBach, + dcaCascDau, dcaV0Dau, dcaCharmBaryonDau, + decLenCharmBaryon, decLenCasc, decLenV0, errorDecayLength, errorDecayLengthXY); + } + } // candidate loop + } // end of run function + + // template function for running Charm Baryon reconstruction via KFParticle method + /// \brief centEstimator is for different centrality estimators + /// \brief decayChannel is for different decay channels. 0 for XiczeroOmegaczeroToXiPi, 1 for OmegaczeroToOmegaPi, 2 for OmegaczeroToOmegaK + /// \brief Colls is for collision tables joined with different centrality estimators + /// \brief Hist is for QA histograms + template + void runCreatorWithKfParticle(Colls const&, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const&, // -> Implemented for internal cascade building + aod::V0s const&, // -> Implemented for internal cascade building + TracksWCovDcaExtraPidPrPiKa const& tracks, + TracksWCovExtraPidIU const& lfTracks, + aod::BCsWithTimestamps const&, + Hist& hInvMassCharmBaryon, + Hist& hCandCounter) + { + // Loop over candidates + for (auto const& cand : candidates) { + + // Fill cascandidates before selection + if (configs.fillHistograms) { + hCandCounter->Fill(All); + } + + // Apply hfflag selection + if (!TESTBIT(cand.hfflag(), decayChannel)) { + continue; + } else { + if (configs.fillHistograms) { + hCandCounter->Fill(HfFlagPass); + } + } + + // Event selection + auto collision = cand.collision_as(); + float centrality{-1.f}; + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); + if (rejectionMask != 0) { // None of the event seletion satisfied -> Reject this candidate + continue; + } + + //------------------------------Set Magnetic field------------------------------ + auto bc = collision.template bc_as(); + if (runNumber != bc.runNumber()) { + LOG(info) << ">>>>>>>>>> Current run Number : " << runNumber; + initCCDB(bc, runNumber, ccdb, configs.isRun2 ? configs.ccdbPathGrp : configs.ccdbPathGrpMag, lut, configs.isRun2); + magneticField = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << ">>>>>>>>>> Magnetic field: " << magneticField; + } + // magnetic field setting for KFParticle + straHelper.fitter.setBz(magneticField); // -> Manetic field setting for internal cascade building + KFParticle::SetField(magneticField); // -> Magnetic field setting for CharmBaryon building + + //------------------Intenal Cascade building------------------ + auto cascAodElement = cand.cascade_as(); + auto v0AodElement = cascAodElement.v0_as(); + auto posTrack = lfTracks.rawIteratorAt(v0AodElement.posTrackId()); + auto negTrack = lfTracks.rawIteratorAt(v0AodElement.negTrackId()); + auto bachTrack = lfTracks.rawIteratorAt(cascAodElement.bachelorId()); + + // Make cascade starting from V0 + // If success, fill Cascade and V0 information for reconstruction + if (!straHelper.buildCascadeCandidateWithKF(collision.globalIndex(), + collision.posX(), collision.posY(), collision.posZ(), + posTrack, + negTrack, + bachTrack, + false, // calculateBachelorBaryonVariables + LFConfigs.kfConstructMethodFromLF, + LFConfigs.kfTuneForOmegaFromLF, + LFConfigs.kfUseV0MassConstraintFromLF, + LFConfigs.kfUseCascadeMassConstraintFromLF, + LFConfigs.kfDoDCAFitterPreMinimV0FromLF, + LFConfigs.kfDoDCAFitterPreMinimCascFromLF)) { + LOG(info) << "This cascade cannot be rebuilt"; + continue; + } else { + float storeMass = (decayChannel == 0) ? straHelper.cascade.massXi : straHelper.cascade.massOmega; + float storePt = RecoDecay::pt(straHelper.cascade.cascadeMomentum); + registry.fill(HIST("hCascMass"), storeMass); + registry.fill(HIST("hCascPt"), storePt); + if (configs.fillHistograms) { + hCandCounter->Fill(CascReconstructed); + } + } + + //------------------------------Cascade pre-selection------------------------------ + // ! only for Xic0, Omegac0 -> Omega K + if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaK) { + if (configs.doCascadePreselection) { + if (std::abs(straHelper.cascade.cascadeDCAxy) > configs.dcaXYToPVCascadeMax) { + continue; + } + // pre selection on dcaV0Daughters + if (std::abs(straHelper.cascade.v0DaughterDCA) > configs.dcaV0DaughtersMax) { + continue; + } + // pre selection on dcaCascDaughters + if (std::abs(straHelper.cascade.cascadeDaughterDCA) > configs.dcaCascDaughtersMax) { + continue; + } + // pre selection on invariantmass + if (std::abs(straHelper.cascade.massOmega - massOfCascade) > configs.massToleranceCascade) { + continue; + } + } + } + + //----------Create charm bayron as KF Partible object starting from V0---------- + + // Create KFParticle object of V0 Daughter & Bachelor + const KFPTrack kfTrack0 = createKFPTrackFromTrack(posTrack); + const KFPTrack kfTrack1 = createKFPTrackFromTrack(negTrack); + const KFPTrack kfTrackBach = createKFPTrackFromTrack(bachTrack); + + bool isAnti = (bachTrack.signed1Pt() > 0 ? true : false); + + KFParticle kfPos(kfTrack0, (isAnti ? pdgIdOfAntiV0DauPos : pdgIdOfV0DauPos)); + KFParticle kfNeg(kfTrack1, (isAnti ? pdgIdOfAntiV0DauNeg : pdgIdOfV0DauNeg)); + KFParticle kfBach(kfTrackBach, (isAnti ? pdgIdOfAntiBach : pdgIdOfBach)); + KFParticle kfBachRej(kfTrackBach, (isAnti ? pdgIdOfAntiBach : pdgIdOfBach)); // Rej -> Used for Omegac0->OmegaPi only + + // ~~~~~~~Construct V0 with KF~~~~~~~ + const KFParticle* v0Daughters[2] = {&kfPos, &kfNeg}; + KFParticle kfV0; + kfV0.SetConstructMethod(configs.kfConstructMethod); + try { + kfV0.Construct(v0Daughters, 2); + } catch (std::runtime_error& e) { + LOG(debug) << "Failed to construct v0" << e.what(); + continue; + } + + // Require lambda pre-selection before mass constraint + float massLam, sigMassLam; + kfV0.GetMass(massLam, sigMassLam); + if (std::abs(massLam - MassLambda0) > configs.lambdaMassWindow) { + continue; + } + if (sigMassLam <= 0) { + continue; + } + if ((kfV0.GetNDF() <= 0) || (kfV0.GetChi2() <= 0)) { + continue; + } + + // Set mass constraint to lambda + KFParticle kfV0MassConstrained = kfV0; + kfV0MassConstrained.SetNonlinearMassConstraint(massOfV0); + if (configs.kfUseV0MassConstraint) { + kfV0 = kfV0MassConstrained; + } + kfV0.TransportToDecayVertex(); + + //~~~~~~~Construct cascade with KF~~~~~~~ + const KFParticle* cascDaughters[2] = {&kfBach, &kfV0}; + const KFParticle* cascDaughtersRej[2] = {&kfBachRej, &kfV0}; + KFParticle kfCasc, kfCascRej; + + kfCasc.SetConstructMethod(configs.kfConstructMethod); + try { + kfCasc.Construct(cascDaughters, 2); + } catch (std::runtime_error& e) { + LOG(debug) << "Failed to construct Cascade: " << e.what(); + } + + float massCasc, sigMassCasc, massCascRej, sigMassCascRej; + kfCasc.GetMass(massCasc, sigMassCasc); + + if (sigMassCasc <= 0) { + continue; + } + if (std::abs(massCasc - massCasc) > configs.massToleranceCascade) { + continue; + } + if (kfCasc.GetNDF() <= 0 || kfCasc.GetChi2() <= 0) { + continue; + } + + // perform cascade building on casc_rej - only for Omega + if constexpr (decayChannel != hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi) { + kfCascRej.SetConstructMethod(configs.kfConstructMethod); + try { + kfCascRej.Construct(cascDaughtersRej, 2); + } catch (std::runtime_error& e) { + LOG(debug) << "Failed to construct Cascade_rej: " << e.what(); + } + + kfCascRej.GetMass(massCascRej, sigMassCascRej); + } + + // Set mass constraint to cascade + KFParticle kfCascMassConstrained = kfCasc; + kfCascMassConstrained.SetNonlinearMassConstraint(massOfCascade); + if (configs.kfUseCascadeMassConstraint) { + kfCasc = kfCascMassConstrained; + } + kfCasc.TransportToDecayVertex(); + + //~~~~~~~Construct Charm Baryon with KF~~~~~~~ + auto trackCharmBachelor = tracks.rawIteratorAt(cand.prong0Id()); + const KFPTrack kfTrackCharmBach = createKFPTrackFromTrack(trackCharmBachelor); + const KFParticle kfCharmBach(kfTrackCharmBach, (isAnti ? pdgIdOfAntiCharmBach : pdgIdOfCharmBach)); + const KFParticle* charmBaryonDaughters[2] = {&kfCharmBach, &kfCasc}; + + KFParticle kfCharmBaryon; + kfCharmBaryon.SetConstructMethod(configs.kfConstructMethod); + try { + kfCharmBaryon.Construct(charmBaryonDaughters, 2); + } catch (std::runtime_error& e) { + LOG(debug) << "Failed to construct Charm baryon: " << e.what(); + } + + float massCharmBaryon, sigMassCharmBaryon; + kfCharmBaryon.GetMass(massCharmBaryon, sigMassCharmBaryon); + if (sigMassCharmBaryon <= 0) { + continue; + } + if (kfCharmBaryon.GetNDF() <= 0 || kfCharmBaryon.GetChi2() <= 0) { + continue; + } + kfCharmBaryon.TransportToDecayVertex(); + if (configs.fillHistograms) { + hCandCounter->Fill(VertexFit); + hInvMassCharmBaryon->Fill(massCharmBaryon); + } + + // Set production vertex + // PV + const KFPVertex kfVertex = createKFPVertexFromCollision(collision); + const KFParticle kfPv(kfVertex); + const KFParticle kfPosOrigin = kfPos; + const KFParticle kfNegOrigin = kfNeg; + + // To V0 + kfPos.SetProductionVertex(kfV0); + kfNeg.SetProductionVertex(kfV0); + + // To Casc + KFParticle kfBachToCasc = kfBach; + KFParticle kfV0ToCasc = kfV0; + kfBach.SetProductionVertex(kfCasc); + kfV0ToCasc.SetProductionVertex(kfCasc); + + // To Charm baryon + KFParticle kfCascToCharmBaryon = kfCasc; + KFParticle kfCharmBachToCharmBaryon = kfCharmBach; + kfCascToCharmBaryon.SetProductionVertex(kfCharmBaryon); + kfCharmBachToCharmBaryon.SetProductionVertex(kfCharmBaryon); + + // To Pv + KFParticle kfV0ToPv = kfV0; + KFParticle kfCascToPv = kfCasc; + KFParticle kfCharmBachToPv = kfCharmBach; + KFParticle kfCharmBaryonToPv = kfCharmBaryon; + kfV0ToPv.SetProductionVertex(kfPv); + kfCascToPv.SetProductionVertex(kfPv); + kfCharmBachToPv.SetProductionVertex(kfPv); + kfCharmBaryonToPv.SetProductionVertex(kfPv); + + //----------Reconstruct information after vertex fit---------- + std::array vertexV0 = {kfV0.GetX(), kfV0.GetY(), kfV0.GetZ()}; + std::array vertexCasc = {kfCasc.GetX(), kfCasc.GetY(), kfCasc.GetZ()}; + + std::array pVecV0DauPos = {kfPos.GetPx(), kfPos.GetPy(), kfPos.GetPz()}; + std::array pVecV0DauNeg = {kfNeg.GetPx(), kfNeg.GetPy(), kfNeg.GetPz()}; + std::array pVecV0 = {kfV0.GetPx(), kfV0.GetPy(), kfV0.GetPz()}; + std::array pVecBach = {kfBachToCasc.GetPx(), kfBachToCasc.GetPy(), kfBachToCasc.GetPz()}; + std::array pVecCharmBachelorAsD = {kfCharmBachToCharmBaryon.GetPx(), kfCharmBachToCharmBaryon.GetPy(), kfCharmBachToCharmBaryon.GetPz()}; + std::array pVecCharmBaryon = {kfCharmBaryon.GetPx(), kfCharmBaryon.GetPy(), kfCharmBaryon.GetPy()}; + + auto* covVtxCharmBaryon = kfCharmBaryon.CovarianceMatrix(); + float covMatrixPv[6]; + kfVertex.GetCovarianceMatrix(covMatrixPv); + + std::array pvCoord = {collision.posX(), collision.posY(), collision.posZ()}; + + auto trackParCovV0DauPos = getTrackParCovFromKFP(kfPos, kfPos.GetPDG(), 1); + auto trackParCovV0DauNeg = getTrackParCovFromKFP(kfNeg, kfNeg.GetPDG(), -1); + auto trackParCovBach = getTrackParCovFromKFP(kfBachToCasc, kfBachToCasc.GetPDG(), (isAnti ? 1 : -1)); + auto trackParCovCharmBach = getTrackParCovFromKFP(kfCharmBachToCharmBaryon, kfCharmBachToCharmBaryon.GetPDG(), (isAnti ? -1 : 1)); + auto trackParCovCasc = getTrackParCovFromKFP(kfCascToCharmBaryon, kfCascToCharmBaryon.GetPDG(), (isAnti ? 1 : -1)); + trackParCovV0DauPos.setAbsCharge(1); + trackParCovV0DauNeg.setAbsCharge(1); + trackParCovBach.setAbsCharge(1); + trackParCovCharmBach.setAbsCharge(1); + trackParCovCasc.setAbsCharge(1); + + //----------Calculate physical quantities and fill candidate table---------- + + // impact parameters + std::array impactParameterV0DauPos; + std::array impactParameterV0DauNeg; + std::array impactParameterBach; + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCovV0DauPos, 2.f, matCorr, &impactParameterV0DauPos); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCovV0DauNeg, 2.f, matCorr, &impactParameterV0DauNeg); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCovBach, 2.f, matCorr, &impactParameterBach); + float dcaxyV0DauPos = impactParameterV0DauPos[0]; + float dcaxyV0DauNeg = impactParameterV0DauNeg[0]; + float dcaxyBach = impactParameterBach[0]; + float dcazV0DauPos = impactParameterV0DauPos[1]; + float dcazV0DauNeg = impactParameterV0DauNeg[1]; + float dcazBach = impactParameterBach[1]; + + o2::dataformats::DCA impactParameterCasc, impactParameterCharmBachelor; + auto primaryVertex = getPrimaryVertex(collision); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovCasc, 2.f, matCorr, &impactParameterCasc); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovCharmBach, 2.f, matCorr, &impactParameterCharmBachelor); + + // get Chi2Topo/NDF + float chi2NdfTopoV0ToPv = kfV0ToPv.GetChi2() / kfV0ToPv.GetNDF(); + float chi2NdfTopoCascToPv = kfCascToPv.GetChi2() / kfCascToPv.GetNDF(); + float chi2NdfTopoCharmBachToPv = kfCharmBachToPv.GetChi2() / kfCharmBachToPv.GetNDF(); + float chi2NdfTopoCharmBaryonToPv = kfCharmBaryonToPv.GetChi2() / kfCharmBaryonToPv.GetNDF(); + float chi2NdfTopoBachToCasc = kfBachToCasc.GetChi2() / kfBachToCasc.GetNDF(); + float chi2NdfTopoV0ToCasc = kfV0ToCasc.GetChi2() / kfV0ToCasc.GetNDF(); + float chi2NdfTopoCharmBachToCharmBaryon = kfCharmBachToCharmBaryon.GetChi2() / kfCharmBachToCharmBaryon.GetChi2(); + float chi2NdfTopoCascToCharmBaryon = kfCascToCharmBaryon.GetChi2() / kfCascToCharmBaryon.GetChi2(); + + // get ldl + float ldlV0 = ldlFromKF(kfV0, kfPv); + float ldlCasc = ldlFromKF(kfCasc, kfPv); + float ldlCharmBaryon = ldlFromKF(kfCharmBaryon, kfPv); + + // get DCAs + float kfDcaV0Daughters = kfNeg.GetDistanceFromParticle(kfPos); + float kfDcaCascDaughters = kfBachToCasc.GetDistanceFromParticle(kfV0ToCasc); + float kfDcaCharmBaryonDaughters = kfCharmBachToCharmBaryon.GetDistanceFromParticle(kfCascToCharmBaryon); + float kfDcaXYCharmBachelorToPv = kfCharmBachToCharmBaryon.GetDistanceFromVertexXY(kfPv); + float kfDcaXYCascToPv = kfCascToCharmBaryon.GetDistanceFromVertexXY(kfPv); + + // get decay length - In XY + float decayLXYV0, errDecayLXYV0; + kfV0ToCasc.GetDecayLengthXY(decayLXYV0, errDecayLXYV0); + + float decayLXYCasc, errDecayLXYCasc; + kfCascToCharmBaryon.GetDecayLengthXY(decayLXYCasc, errDecayLXYCasc); + + float decayLXYCharmBaryon, errDecayLXYCharmBaryon; + kfCharmBaryonToPv.GetDecayLengthXY(decayLXYCharmBaryon, errDecayLXYCharmBaryon); + + // get decay length - In XYZ + float decayLV0 = RecoDecay::distance(std::array{kfCasc.GetX(), kfCasc.GetY(), kfCasc.GetZ()}, std::array{kfV0.GetX(), kfV0.GetY(), kfV0.GetZ()}); + float decayLCasc = RecoDecay::distance(std::array{kfCharmBaryon.GetX(), kfCharmBaryon.GetY(), kfCharmBaryon.GetZ()}, std::array{kfCasc.GetX(), kfCasc.GetY(), kfCasc.GetZ()}); + float decayLCharmBaryon = RecoDecay::distance(std::array{collision.posX(), collision.posY(), collision.posZ()}, std::array{kfCharmBaryon.GetX(), kfCharmBaryon.GetY(), kfCharmBaryon.GetZ()}); + + double phiCharmBaryon, thetaCharmBaryon; + getPointDirection(std::array{kfV0.GetX(), kfV0.GetY(), kfV0.GetZ()}, std::array{kfCharmBaryon.GetX(), kfCharmBaryon.GetY(), kfCharmBaryon.GetZ()}, phiCharmBaryon, thetaCharmBaryon); + float errDecayLCharmBaryon = std::sqrt(getRotatedCovMatrixXX(covMatrixPv, phiCharmBaryon, thetaCharmBaryon) + getRotatedCovMatrixXX(covVtxCharmBaryon, phiCharmBaryon, thetaCharmBaryon)); + + // get cosine of pointing angle + float cosPaV0ToPv = cpaFromKF(kfV0, kfPv); + float cosPaCascToPv = cpaFromKF(kfCasc, kfPv); + float cosPaCharmBaryonToPv = cpaFromKF(kfCharmBaryon, kfPv); + + float cosPaXYV0ToPv = RecoDecay::cpaXY(pvCoord, vertexV0, pVecV0); + float cosPaXYCascToPv = cpaXYFromKF(kfCasc, kfPv); + float cosPaXYCharmBaryonToPv = cpaXYFromKF(kfCharmBaryon, kfPv); + + float cosPaV0ToCasc = cpaFromKF(kfV0, kfCasc); + float cosPaCascToCharmBaryon = cpaFromKF(kfCasc, kfCharmBaryon); + float cosPaXYV0ToCasc = RecoDecay::cpaXY(vertexCasc, vertexV0, pVecV0); + float cosPaXYCascToCharmBaryon = cpaXYFromKF(kfCasc, kfCharmBaryon); + + float deviationCharmBachToPv = kfCalculateChi2ToPrimaryVertex(kfCharmBaryon, kfPv); // -> For Omegac0 + + // KF pT, eta + float ptCasc = kfCascToCharmBaryon.GetPt(); + float ptCharmBachelor = kfCharmBachToCharmBaryon.GetPt(); + float ptCharmBaryon = kfCharmBaryon.GetPt(); + float yCharmBaryon = kfCharmBaryon.GetRapidity(); + + // get KF cosThetaStar + float cosThetaStarCharmBachelor; + if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi) { + cosThetaStarCharmBachelor = cosThetaStarFromKF(0, pdgIdOfCharmBaryon, pdgIdOfCharmBach, pdgIdOfCascade, kfCharmBachToCharmBaryon, kfCascToCharmBaryon); + } + float cosThetaStarKaFromOmegac0, cosThetaStarKaFromXic0; // -> Only requeted to be calculated and filled for Xic0/Omegac0 -> Omega K + if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaK) { + cosThetaStarKaFromOmegac0 = cosThetaStarFromKF(0, 4332, 321, 3334, kfCharmBachToCharmBaryon, kfCascToCharmBaryon); + cosThetaStarKaFromXic0 = cosThetaStarFromKF(0, 4132, 321, 3334, kfCharmBachToCharmBaryon, kfCascToCharmBaryon); + } + + // KF ct + float ctV0 = kfV0ToCasc.GetLifeTime(); + float ctCasc = kfCascToCharmBaryon.GetLifeTime(); + float ctCharmBaryon = kfCharmBaryonToPv.GetLifeTime(); + + //------------------------------Calculate physical quantities and fill candidate table------------------------------ + + //------------------------------Fill the table------------------------------ + if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi) { + cursors.rowCandToXiPiKf(collision.globalIndex(), // Global index of collision + pvCoord[0], pvCoord[1], pvCoord[2], // coordination of PV + vertexCasc[0], vertexCasc[1], vertexCasc[2], // Decay position of kfCasc + vertexV0[0], vertexV0[1], vertexV0[2], + bachTrack.sign(), + covVtxCharmBaryon[0], covVtxCharmBaryon[1], covVtxCharmBaryon[2], covVtxCharmBaryon[3], covVtxCharmBaryon[4], covVtxCharmBaryon[5], + pVecCharmBaryon[0], pVecCharmBaryon[1], pVecCharmBaryon[2], // x, y, z momentum of charm baryon + kfCascToCharmBaryon.GetPx(), kfCascToCharmBaryon.GetPy(), kfCascToCharmBaryon.GetPz(), + pVecCharmBachelorAsD[0], pVecCharmBachelorAsD[1], pVecCharmBachelorAsD[2], + pVecV0[0], pVecV0[1], pVecV0[2], + pVecBach[0], pVecBach[1], pVecBach[2], + pVecV0DauPos[0], pVecV0DauPos[1], pVecV0DauPos[2], + pVecV0DauNeg[0], pVecV0DauNeg[1], pVecV0DauNeg[2], + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), + cand.cascadeId(), trackCharmBachelor.globalIndex(), cand.prong0Id(), + massLam, massCasc, massCharmBaryon, + cosPaV0ToPv, cosPaCascToPv, + ctCasc, ctV0, ctCharmBaryon, + kfPos.GetEta(), kfNeg.GetEta(), kfBach.GetEta(), kfCharmBachToCharmBaryon.GetEta(), + kfCharmBaryon.GetEta(), kfCasc.GetEta(), kfV0.GetEta(), + dcaxyV0DauPos, dcaxyV0DauNeg, dcaxyBach, + kfDcaCascDaughters, kfDcaV0Daughters, kfDcaCharmBaryonDaughters, + kfDcaXYCharmBachelorToPv, kfDcaXYCascToPv, + kfV0.GetChi2(), kfCasc.GetChi2(), kfCharmBaryon.GetChi2(), kfV0MassConstrained.GetChi2(), kfCascMassConstrained.GetChi2(), + ldlV0, ldlCasc, // ldlCharmBaryon, + chi2NdfTopoV0ToPv, chi2NdfTopoCascToPv, chi2NdfTopoCharmBachToPv, chi2NdfTopoCharmBaryonToPv, + chi2NdfTopoV0ToCasc, chi2NdfTopoCascToCharmBaryon, + decayLXYV0, decayLXYCasc, decayLXYCharmBaryon, + cosPaV0ToCasc, cosPaCascToCharmBaryon, // cosPaXYV0ToCasc, cosPaXYCascToCharmBaryon, + yCharmBaryon, // ptCharmBachelor, ptCharmBaryon, + cosThetaStarCharmBachelor, + kfV0.GetNDF(), kfCasc.GetNDF(), kfCharmBaryon.GetNDF(), kfV0MassConstrained.GetNDF(), kfCascMassConstrained.GetNDF(), + kfV0.GetChi2() / kfV0.GetNDF(), kfCasc.GetChi2() / kfCasc.GetNDF(), kfCharmBaryon.GetChi2() / kfCharmBaryon.GetNDF(), kfV0MassConstrained.GetChi2() / kfV0MassConstrained.GetNDF(), kfCascMassConstrained.GetChi2() / kfCascMassConstrained.GetNDF()); + } else if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaPi) { + cursors.rowCandToOmegaPi(collision.globalIndex(), + pvCoord[0], pvCoord[1], pvCoord[2], + 0.f, 0.f, 0.f, // -> vertexCharmBaryonFromFitter. For KF, this is 0 + vertexCasc[0], vertexCasc[1], vertexCasc[2], + vertexV0[0], vertexV0[1], vertexV0[2], + trackCharmBachelor.sign(), + covVtxCharmBaryon[0], covVtxCharmBaryon[1], covVtxCharmBaryon[2], covVtxCharmBaryon[3], covVtxCharmBaryon[4], covVtxCharmBaryon[5], + kfCharmBaryon.GetPx(), kfCharmBaryon.GetPy(), kfCharmBaryon.GetPz(), + kfCascToCharmBaryon.GetPx(), kfCascToCharmBaryon.GetPy(), kfCascToCharmBaryon.GetPz(), + kfCharmBachToCharmBaryon.GetPx(), kfCharmBachToCharmBaryon.GetPy(), kfCharmBachToCharmBaryon.GetPz(), + pVecV0[0], pVecV0[1], pVecV0[2], + pVecBach[0], pVecBach[1], pVecBach[2], + pVecV0DauPos[0], pVecV0DauPos[1], pVecV0DauPos[2], + pVecV0DauNeg[0], pVecV0DauNeg[1], pVecV0DauNeg[2], + impactParameterCasc.getY(), impactParameterCharmBachelor.getY(), + impactParameterCasc.getZ(), impactParameterCharmBachelor.getZ(), + std::sqrt(impactParameterCasc.getSigmaY2()), std::sqrt(impactParameterCharmBachelor.getSigmaY2()), + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), + cand.cascadeId(), trackCharmBachelor.globalIndex(), cand.prong0Id(), + massLam, massCasc, massCharmBaryon, + cosPaV0ToPv, cosPaCharmBaryonToPv, cosPaCascToPv, cosPaXYV0ToPv, cosPaXYCharmBaryonToPv, cosPaXYCascToPv, + ctCharmBaryon, ctCasc, ctV0, + kfPos.GetEta(), kfNeg.GetEta(), kfBach.GetEta(), kfCharmBachToCharmBaryon.GetEta(), + kfCharmBaryon.GetEta(), kfCasc.GetEta(), kfV0.GetEta(), + dcaxyV0DauPos, dcaxyV0DauNeg, dcaxyBach, + dcazV0DauPos, dcazV0DauNeg, dcazBach, + kfDcaCascDaughters, straHelper.cascade.v0DaughterDCA, kfDcaCharmBaryonDaughters, + decayLCharmBaryon, decayLCasc, decayLV0, errDecayLCharmBaryon, errDecayLXYCharmBaryon, cand.hfflag()); + + cursors.rowCandToOmegaPiKf(kfDcaXYCharmBachelorToPv, kfDcaXYCascToPv, + kfV0.GetChi2(), kfCasc.GetChi2(), kfCharmBaryon.GetChi2(), kfV0MassConstrained.GetChi2(), kfCascMassConstrained.GetChi2(), + ldlV0, ldlCasc, ldlCharmBaryon, + chi2NdfTopoV0ToPv, chi2NdfTopoCascToPv, chi2NdfTopoCharmBachToPv, chi2NdfTopoCharmBaryonToPv, deviationCharmBachToPv, + chi2NdfTopoV0ToCasc, chi2NdfTopoCascToCharmBaryon, + decayLXYV0, decayLXYCasc, decayLXYCharmBaryon, + cosPaV0ToCasc, cosPaCascToCharmBaryon, cosPaXYV0ToCasc, cosPaXYCascToCharmBaryon, + kfCharmBaryon.GetRapidity(), ptCharmBachelor, ptCharmBaryon, + cosThetaStarKaFromOmegac0, + kfV0.GetNDF(), kfCasc.GetNDF(), kfCharmBaryon.GetNDF(), kfV0MassConstrained.GetNDF(), kfCascMassConstrained.GetNDF(), + kfV0.GetChi2() / kfV0.GetNDF(), kfCasc.Chi2() / kfCasc.GetNDF(), kfCharmBaryon.GetChi2() / kfCharmBaryon.GetNDF(), + kfV0MassConstrained.GetChi2() / kfV0.GetNDF(), kfCascMassConstrained.Chi2() / kfCasc.GetNDF(), + massCascRej); + + } else { + cursors.rowCandToOmegaKaKf(collision.globalIndex(), + collision.posX(), collision.posY(), collision.posZ(), + kfPv.GetX(), kfPv.GetY(), kfPv.GetZ(), + straHelper.cascade.v0Position[0], straHelper.cascade.v0Position[1], straHelper.cascade.v0Position[2], + straHelper.cascade.v0Momentum[0], straHelper.cascade.v0Momentum[1], straHelper.cascade.v0Momentum[2], + straHelper.cascade.cascadePosition[0], straHelper.cascade.cascadePosition[1], straHelper.cascade.cascadePosition[2], + straHelper.cascade.cascadeMomentum[0], straHelper.cascade.cascadeMomentum[1], straHelper.cascade.cascadeMomentum[2], + kfV0.GetX(), kfV0.GetY(), kfV0.GetZ(), + kfV0.GetPx(), kfV0.GetPy(), kfV0.GetPz(), + kfCasc.GetX(), kfCasc.GetY(), kfCasc.GetZ(), + kfCasc.GetPx(), kfCasc.GetPy(), kfCasc.GetPz(), + kfCharmBaryon.GetX(), kfCharmBaryon.GetY(), kfCharmBaryon.GetZ(), + kfCharmBaryon.GetPx(), kfCharmBaryon.GetPy(), kfCharmBaryon.GetPz(), + straHelper.cascade.charge, + kfPos.GetEta(), kfNeg.GetEta(), kfBach.GetEta(), kfCharmBach.GetEta(), kfV0.GetEta(), kfCasc.GetEta(), kfCharmBaryon.GetEta(), kfCharmBaryon.GetRapidity(), + impactParameterCharmBachelor.getY(), std::sqrt(impactParameterCharmBachelor.getSigmaY2()), impactParameterCasc.getY(), std::sqrt(impactParameterCasc.getSigmaY2()), + kfDcaV0Daughters, kfDcaCascDaughters, kfDcaCharmBaryonDaughters, + cosPaV0ToPv, cosPaCascToPv, cosPaCharmBaryonToPv, cosPaXYV0ToPv, cosPaXYCascToPv, cosPaXYCharmBaryonToPv, cosPaV0ToCasc, cosPaCascToCharmBaryon, cosPaXYV0ToCasc, cosPaXYCascToCharmBaryon, + kfV0.GetChi2() / kfV0.GetNDF(), kfCasc.GetChi2() / kfCasc.GetNDF(), kfCharmBaryon.GetChi2() / kfCharmBaryon.GetNDF(), + kfV0MassConstrained.GetChi2() / kfV0MassConstrained.GetNDF(), kfCascMassConstrained.GetChi2() / kfCascMassConstrained.GetNDF(), + chi2NdfTopoV0ToCasc, chi2NdfTopoBachToCasc, chi2NdfTopoCharmBachToCharmBaryon, chi2NdfTopoCascToCharmBaryon, // Topological constraints to mother + chi2NdfTopoV0ToPv, chi2NdfTopoCascToPv, chi2NdfTopoCharmBachToPv, chi2NdfTopoCharmBaryonToPv, // Topological constraints to PV + ldlV0, ldlCasc, ldlCharmBaryon, + decayLXYV0, decayLXYCasc, decayLXYCharmBaryon, + massLam, sigMassLam, massCasc, sigMassCasc, massCascRej, sigMassCascRej, massCharmBaryon, sigMassCharmBaryon, + ptCharmBaryon, ptCharmBachelor, ptCasc, + cosThetaStarKaFromOmegac0, cosThetaStarKaFromXic0, ctV0, ctCasc, ctCharmBaryon, + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), cand.cascadeId(), cand.prong0Id(), trackCharmBachelor.globalIndex()); + } + } // end candidate loop + } // end of runCreator + + ///////////////////////////////////////////////////// + /// /// + /// Process functions with DCAFitter /// + /// /// + ///////////////////////////////////////////////////// + + /*~~~~~~~~~~~~~~*/ + /*~~~To Xi Pi~~~*/ + /*~~~~~~~~~~~~~~*/ + void processToXiPiWithDCAFitterNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + aod::TracksWCovDca const& tracks, + TracksWCovIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithDCAFitterNoCent, "Charm candidte reconstruction with Xi Pi via DcaFitter method, no centrality", true); + +#if 0 + void processToXiPiWithDCAFitterNoCentWithTrackedCasc(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::TrackedCascades const& cascades, + aod::V0s const& v0s, + TrackedCascLinked const& trackedCascLinked, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, trackedCascFull, trackedCascLinked, tracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithDCAFitterNoCentWithTrackedCasc, "Charm candidte reconstruction with Xi Pi via DcaFitter method with tracked cascade, no centrality", false); +#endif + + void processToXiPiWithDCAFitterCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + aod::TracksWCovDca const& tracks, + TracksWCovIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithDCAFitterCentFT0C, "Charm candidate reconstruction with Xi Pi via DcaFitter method, centrality selection on FT0C", false); + + void processToXiPiWithDCAFitterCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + aod::TracksWCovDca const& tracks, + TracksWCovIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithDCAFitterCentFT0M, "Charm candidate reconstruction with Xi Pi via DcaFitter method, centrality selection on FT0M", false); + + /*~~~~~~~~~~~~~~~~~*/ + /*~~~To Omega Pi~~~*/ + /*~~~~~~~~~~~~~~~~~*/ + void processToOmegaPiWithDCAFitterNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + aod::TracksWCovDca const& tracks, + TracksWCovIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithDCAFitterNoCent, "Charm candidte reconstruction with Omega Pi via DcaFitter method, no centrality", false); + + void processToOmegaPiWithDCAFitterCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + aod::TracksWCovDca const& tracks, + TracksWCovIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithDCAFitterCentFT0C, "Charm candidate reconstruction with Omega Pi via DcaFitter method, centrality selection on FT0C", false); + + void processToOmegaPiWithDCAFitterCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + aod::TracksWCovDca const& tracks, + TracksWCovIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hInvMassCharmBaryonToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithDCAFitterCentFT0M, "Charm candidate reconstruction with Omega Pi via DcaFitter method, centrality selection on FT0M", false); + + /*~~~~~~~~~~~~~~~~~*/ + /*~~~To Omega Ka~~~*/ + /*~~~~~~~~~~~~~~~~~*/ + void processToOmegaKaWithDCAFitterNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + aod::TracksWCovDca const& tracks, + TracksWCovIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaKa); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithDCAFitterNoCent, "Charm candidte reconstruction with Omega Ka via DcaFitter method, no centrality", false); + + void processToOmegaKaWithDCAFitterCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + aod::TracksWCovDca const& tracks, + TracksWCovIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaKa); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithDCAFitterCentFT0C, "Charm candidate reconstruction with Omega Ka via DcaFitter method, centrality selection on FT0C", false); + + void processToOmegaKaWithDCAFitterCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + aod::TracksWCovDca const& tracks, + TracksWCovIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithDCAFitterCentFT0M, "Charm candidate reconstruction with Omega Ka via DcaFitter method, centrality selection on FT0M", false); + + ///////////////////////////////////////////////////// + /// /// + /// Process functions with KFParticle /// + /// /// + ///////////////////////////////////////////////////// + + /*~~~~~~~~~~~~~~*/ + /*~~~To Xi Pi~~~*/ + /*~~~~~~~~~~~~~~*/ + void processToXiPiWithKFParticleNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + TracksWCovExtraPidIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithKFParticleNoCent, "Charm Baryon decaying to Xi Pi reconstruction via KFParticle method, no centrality", false); + + void processToXiPiWithKFParticleCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + TracksWCovExtraPidIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithKFParticleCentFT0C, "Charm Baryon decaying to Xi Pi reconstruction via KFParticle method, centrality on FT0C", false); + + void processToXiPiWithKFParticleCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + TracksWCovExtraPidIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithKFParticleCentFT0M, "Charm Baryon decaying to Xi Pireconstruction via KFParticle method, centrality on FT0M", false); + + /*~~~~~~~~~~~~~~~~~*/ + /*~~~To Omega Pi~~~*/ + /*~~~~~~~~~~~~~~~~~*/ + void processToOmegaPiWithKFParticleNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + TracksWCovExtraPidIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithKFParticleNoCent, "Charm Baryon decaying to Omega Pi reconstruction via KFParticle method, no centrality", false); + + void processToOmegaPiWithKFParticleCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + TracksWCovExtraPidIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithKFParticleCentFT0C, "Charm Baryon decaying to Omega Pi reconstruction via KFParticle method, centrality on FT0C", false); + + void processToOmegaPiWithKFParticleCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + TracksWCovExtraPidIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithKFParticleCentFT0M, "Charm Baryong decaying to Omega Pi reconstruction via KFParticle method, centrality on FT0M", false); + + /*~~~~~~~~~~~~~~~~~*/ + /*~~~To Omega Ka~~~*/ + /*~~~~~~~~~~~~~~~~~*/ + void processToOmegaKaWithKFParticleNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + TracksWCovExtraPidIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaKa); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithKFParticleNoCent, "Charm Baryon decaying to Omega Ka reconstruction via KFParticle method, no centrality", false); + + void processToOmegaKaWithKFParticleCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + TracksWCovExtraPidIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaKa); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithKFParticleCentFT0C, "Charm Baryon decaying to Omega Ka reconstruction via KFParticle method, centrality on FT0C", false); + + void processToOmegaKaWithKFParticleCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + TracksWCovExtraPidIU const& lfTracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, lfTracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaKa); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithKFParticleCentFT0M, "Charm Baryong decaying to Omega Ka reconstruction via KFParticle method, centrality on FT0M", false); + + /////////////////////////////////////////////////////////////// + /// /// + /// Process functions for Collision monitoring /// + /// /// + /////////////////////////////////////////////////////////////// + + void processCollisionsNoCent(soa::Join const& collisions, + aod::BCsWithTimestamps const&) + { + for (const auto& collision : collisions) { + + // bitmask with event selection info + float centrality{-1.f}; + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); + + // monitor the satisfied event selection + hfEvSel.fillHistograms(collision, rejectionMask, centrality, occupancy); + } + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processCollisionsNoCent, "Collision monitoring - No Centrality", true); + + void processCollisionsCentFT0C(soa::Join const& collisions, + aod::BCsWithTimestamps const&) + { + for (const auto& collision : collisions) { + + // bitmask with event selection info + float centrality{-1.f}; + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); + + // monitor the satisfied event selection + hfEvSel.fillHistograms(collision, rejectionMask, centrality, occupancy); + } + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processCollisionsCentFT0C, "Collision monitoring - Centrality selection with FT0C", false); + + void processCollisionsCentFT0M(soa::Join const& collisions, + aod::BCsWithTimestamps const&) + { + for (const auto& collision : collisions) { + + // bitmask with event selection info + float centrality{-1.f}; + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); + + // monitor the satisfied event selection + hfEvSel.fillHistograms(collision, rejectionMask, centrality, occupancy); + } + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processCollisionsCentFT0M, "Collision monitoring - Centrality selection with FT0M", false); +}; + +struct HfCandidateCreatorXic0Omegac0QaMc { + + // Cursor to fill tables + struct : ProducesGroup { + + Produces rowMcMatchRecXicToXiPi; + Produces rowMcMatchGenXicToXiPi; + Produces rowMcMatchRecOmegacToXiPi; + Produces rowMcMatchGenOmegacToXiPi; + Produces rowMcMatchRecToOmegaPi; + Produces rowMcMatchGenToOmegaPi; + Produces rowMcMatchRecToOmegaKa; + Produces rowMcMatchGenToOmegaKa; + + } cursors; + + // Configurables + struct : ConfigurableGroup { + + Configurable rejectBackground{"rejectBackground", true, "Reject particles from background events"}; // -> Used for only Xic0 + Configurable acceptTrackInteractionWithMaterial{"acceptTrackInteractionWithMaterial", false, "Accept candidates with final daughters interacting with materials"}; + Configurable fillMcHistograms{"fillMcHistograms", true, "Fill validation plots"}; + Configurable fillResidualTable{"fillResidualTable", true, "Fill table contaning residuals and pulls of PV and SV"}; + // Configurable matchDecayedPions{"matchedDecayedPions", true, "Match also candidates with daughter pion tracks that decay with kinked toploogy"}; + + } configs; + + enum McMatchFlag : uint8_t { + None = 0, + CharmBaryonUnmatched, + CascUnmatched, + V0Unmatched, + NumberMcMatchFlag + }; + + std::array pdgOfCharmBaryon{+kXiC0, +kOmegaC0, +kOmegaC0, +kOmegaC0}; + std::array pdgOfCascade{+kXiMinus, +kXiMinus, +kOmegaMinus, +kOmegaMinus}; + std::array pdgOfCharmBachelor{+kPiPlus, +kPiPlus, +kPiPlus, +kKPlus}; + std::array pdgOfBachelor{+kPiMinus, +kPiMinus, +kKMinus, +kKMinus}; + + // Table aliases + using TracksWMcIU = soa::Join; + using McCollisionsNoCents = soa::Join; + using McCollisionsFT0Cs = soa::Join; + using McCollisionsFT0Ms = soa::Join; + using McCollisionsCentFT0Ms = soa::Join; // -> Used for subscription for process functions of centrality with FT0Ms + using BCsInfo = soa::Join; + + Preslice mcParticlesPerMcCollision = aod::mcparticle::mcCollisionId; + PresliceUnsorted colPerMcCollision = aod::mccollisionlabel::mcCollisionId; // -> Why use unsorted?? + PresliceUnsorted colPerMcCollisionFT0C = aod::mccollisionlabel::mcCollisionId; + PresliceUnsorted colPerMcCollisionFT0M = aod::mccollisionlabel::mcCollisionId; + + HistogramRegistry registry{"registry"}; + HfEventSelectionMc hfEvSelMc; + + void init(InitContext& initContext) + { + const auto& workflows = initContext.services().get(); + for (const DeviceSpec& device : workflows.devices) { + if (device.name.compare("hf-candidate-creator-xic0-omegac0-qa") == 0) { + hfEvSelMc.init(device, registry); + break; + } + } + + // Add histograms for QA + if (configs.fillMcHistograms) { + registry.add("hGenCharmBaryonPtRapidityTight", "Generated charm baryon #it{p}_{T};#if{p}_{T} (GeV/#it{c});entries", {HistType::kTH1D, {{20, 0.0, 20.0}}}); + registry.add("hGenCharmBaryonPtRapidityLoose", "Generated charm baryon #it{p}_{T};#if{p}_{T} (GeV/#it{c});entries", {HistType::kTH1D, {{20, 0.0, 20.0}}}); + } + } + + // Helper function to fill in table for MC Rec + ///@brief decay Channel is the chosen reconstruction channel + ///@brief flag + ///@brief debug + ///@brief origin + ///@brief collisionMatched + ///@brief pt + ///@brief pdgcode + template + void fillRecoMcTableByDecayChannel(int8_t flag, int8_t debug, int8_t origin, bool collisionMatched, float pt, int pdgCode) + { + if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::XiczeroToXiPi) { + cursors.rowMcMatchRecXicToXiPi(flag, debug, origin, collisionMatched, pt, pdgCode); + } else if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::OmegaczeroToXiPi) { + cursors.rowMcMatchRecOmegacToXiPi(flag, debug, origin, collisionMatched, pt, pdgCode); + } else if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::OmegaczeroToOmegaPi) { + cursors.rowMcMatchRecToOmegaPi(flag, debug, origin, collisionMatched, pt, pdgCode); + } else { + cursors.rowMcMatchRecToOmegaKa(flag, debug, origin, collisionMatched, pt, pdgCode); + } + } + + // Helper function to fill in table for MC Gen + ///@brief decay Channel is the chosen reconstruction channel + ///@brief flag + ///@brief debug + ///@brief origin + ///@brief collisionMatched + ///@brief pt + ///@brief pdgcode + template + void fillGenMcTableByDecayChannel(int8_t flag, int8_t debugGenCharmBaryon, int8_t debugGenCascade, int8_t debugGenLambda, float ptCharmBaryon, float yCharmBaryon, int8_t origin, int idxMother) + { + if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::XiczeroToXiPi) { + cursors.rowMcMatchGenXicToXiPi(flag, debugGenCharmBaryon, debugGenCascade, debugGenLambda, ptCharmBaryon, yCharmBaryon, origin, idxMother); + } else if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::OmegaczeroToXiPi) { + cursors.rowMcMatchGenOmegacToXiPi(flag, debugGenCharmBaryon, debugGenCascade, debugGenLambda, ptCharmBaryon, yCharmBaryon, origin, idxMother); + } else if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::OmegaczeroToOmegaPi) { + cursors.rowMcMatchGenToOmegaPi(flag, debugGenCharmBaryon, debugGenCascade, debugGenLambda, ptCharmBaryon, yCharmBaryon, origin, idxMother); + } else { + cursors.rowMcMatchGenToOmegaKa(flag, debugGenCharmBaryon, debugGenCascade, debugGenLambda, ptCharmBaryon, yCharmBaryon, origin, idxMother); + } + } + + template + void runXic0Omegac0Mc(TRecoCand const& candidates, + TracksWMcIU const&, + aod::McParticles const& mcParticles, + Colls const& collsWithMcLabels, + McCollisions const& mcCollisions, + BCsInfo const&) + { + int indexRec{-1}; + int indexRecCharmBaryon{-1}; + int8_t sign{-9}; + int8_t signCasc{-9}; + int8_t signV0{-9}; + int8_t flag{0}; + int8_t origin{0}; + int8_t debug{0}; + int8_t debugGenCharmBaryon{0}; + int8_t debugGenCasc{0}; + int8_t debugGenV0{0}; + bool collisionMatched = false; + float ptCharmBaryonGen = -999.; + float yCharmBaryonGen = -999.; + + //////////////////////////////////// + // Match reconstructed candidates // + //////////////////////////////////// + + for (const auto& candidate : candidates) { + + flag = 0; + origin = RecoDecay::OriginType::None; + debug = McMatchFlag::None; + collisionMatched = false; + std::vector idxBhadMothers{}; + + auto arrayDaughters = std::array{candidate.template bachelorFromCharmBaryon_as(), + candidate.template bachelor_as(), + candidate.template posTrack_as(), + candidate.template negTrack_as()}; + auto arrayDaughtersCasc = std::array{candidate.template bachelor_as(), + candidate.template posTrack_as(), + candidate.template negTrack_as()}; + auto arrayDaughtersV0 = std::array{candidate.template posTrack_as(), + candidate.template negTrack_as()}; + + // Reject particles from background events + if (configs.rejectBackground) { + bool fromBkg{false}; + for (auto const& daughter : arrayDaughters) { + if (daughter.has_mcParticle()) { + auto mcParticle = daughter.mcParticle(); + if (mcParticle.fromBackgroundEvent()) { + fromBkg = true; + break; + } + } + } + if (fromBkg) { + // fill the tables + fillRecoMcTableByDecayChannel(flag, debug, origin, collisionMatched, -1.f, 0); + continue; + } + } + + // CharmBaryon -> Charm bachelor + Cascade + indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, pdgOfCharmBaryon[decayChannel], std::array{pdgOfCharmBachelor[decayChannel], pdgOfBachelor[decayChannel], +kProton, +kPiMinus}, true, &sign, 3); + indexRecCharmBaryon = indexRec; + if (indexRec == -1) { // Xic0 not reconstructed + debug = McMatchFlag::CharmBaryonUnmatched; + } + if (indexRec > -1) { + // Cascade -> Bachelor + V0 + indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughtersCasc, pdgOfCascade[decayChannel], std::array{pdgOfBachelor[decayChannel], +kProton, +kPiMinus}, true, &signCasc, 2); + if (indexRec == -1) { // Xi- not reconstructed + debug = McMatchFlag::CascUnmatched; + } + if (indexRec > -1) { + // V0 -> Pos + Neg + indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughtersV0, +kLambda0, std::array{+kProton, +kPiMinus}, true, &signV0, 1); + if (indexRec == -1) { // V0 not reconstructed + debug = McMatchFlag::V0Unmatched; + } + if (indexRec > -1) { + flag = sign * (1 << decayChannel); + collisionMatched = candidate.template collision_as().mcCollisionId() == mcParticles.iteratorAt(indexRecCharmBaryon).mcCollisionId(); + } + } + } + + // Check if Xic0 is from b-hadron decay(prompt vs non-prompt) + if (flag != 0) { + auto particle = mcParticles.rawIteratorAt(indexRecCharmBaryon); + origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, false, &idxBhadMothers); + } + if (origin == RecoDecay::OriginType::NonPrompt) { + auto bHadMother = mcParticles.rawIteratorAt(idxBhadMothers[0]); + fillRecoMcTableByDecayChannel(flag, debug, origin, collisionMatched, bHadMother.pt(), bHadMother.pdgCode()); + } else { + fillRecoMcTableByDecayChannel(flag, debug, origin, collisionMatched, -1.f, 0); + } + if (debug == McMatchFlag::CascUnmatched || debug == McMatchFlag::V0Unmatched) { + LOGF(info, "WARNING: Charm baryon decays in the expected final state but the condition on the intermediate states are not fullfilled"); + } + } // candidate loop + + /////////////////////////////// + // Match generated particles // + /////////////////////////////// + + for (const auto& mcCollision : mcCollisions) { + + const auto& mcParticlesPerMcColl = mcParticles.sliceBy(mcParticlesPerMcCollision, mcCollision.globalIndex()); + + float centrality{-1.f}; + uint16_t rejectionMask{0}; + int nSplitColl{0}; + + if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::None) { + const auto collSlice = collsWithMcLabels.sliceBy(colPerMcCollision, mcCollision.globalIndex()); + rejectionMask = hfEvSelMc.getHfMcCollisionRejectionMask(mcCollision, collSlice, centrality); + nSplitColl = collSlice.size(); + } else if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FT0C) { + const auto collSlice = collsWithMcLabels.sliceBy(colPerMcCollisionFT0C, mcCollision.globalIndex()); + rejectionMask = hfEvSelMc.getHfMcCollisionRejectionMask(mcCollision, collSlice, centrality); + nSplitColl = collSlice.size(); + } else if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FT0M) { + const auto collSlice = collsWithMcLabels.sliceBy(colPerMcCollisionFT0M, mcCollision.globalIndex()); + rejectionMask = hfEvSelMc.getHfMcCollisionRejectionMask(mcCollision, collSlice, centrality); + nSplitColl = collSlice.size(); + } + + hfEvSelMc.fillHistograms(mcCollision, rejectionMask, nSplitColl); + + if (rejectionMask != 0) { // At least one event selection not satisfied --> Reject all particles from this event + for (unsigned int i = 0; i < mcParticlesPerMcColl.size(); ++i) { + fillGenMcTableByDecayChannel(0, 0, 0, 0, -999., -999., RecoDecay::OriginType::None, -1); + } + continue; + } + + // Match generated particles + for (auto const& particle : mcParticlesPerMcColl) { + ptCharmBaryonGen = -999.; + yCharmBaryonGen = -999.; + flag = 0; + sign = 0; + debugGenCharmBaryon = 0; + debugGenCasc = 0; + debugGenV0 = 0; + origin = RecoDecay::OriginType::None; + std::vector idxBhadMothers{}; + float kYCutTight = 0.5; + float kYCutLoose = 0.8; + + // Reject particles from background events + if (particle.fromBackgroundEvent() && configs.rejectBackground) { + fillGenMcTableByDecayChannel(flag, debugGenCharmBaryon, debugGenCasc, debugGenV0, ptCharmBaryonGen, yCharmBaryonGen, origin, -1); + continue; + } + + // Charm Baryon -> Cascade + Charm bachelor + if (RecoDecay::isMatchedMCGen(mcParticles, particle, pdgOfCharmBaryon[decayChannel], std::array{pdgOfCascade[decayChannel], pdgOfCharmBachelor[decayChannel]}, true, &sign)) { + debugGenCharmBaryon = 1; + ptCharmBaryonGen = particle.pt(); + yCharmBaryonGen = particle.y(); + debug = 1; // -> Matched Xic0 + + for (auto const& daughterCharm : particle.template daughters_as()) { + if (std::abs(daughterCharm.pdgCode()) != pdgOfCascade[decayChannel]) { + continue; + } + // Xi -> Lambda + pi + if (RecoDecay::isMatchedMCGen(mcParticles, daughterCharm, pdgOfCascade[decayChannel], std::array{+kLambda0, pdgOfBachelor[decayChannel]}, true)) { + debugGenCasc = 1; // -> Matched Xi- + for (auto const& daughterCascade : daughterCharm.template daughters_as()) { + if (std::abs(daughterCascade.pdgCode()) != +kLambda0) { + continue; + } + + // Lambda -> p + pi + if (RecoDecay::isMatchedMCGen(mcParticles, daughterCascade, +kLambda0, std::array{+kProton, +kPiMinus}, true)) { + debugGenV0 = 1; // -> Matched Lambda0 + flag = sign * (1 << decayChannel); + } + } // V0 daughter loop + } // cascade daughter loop + } + } // charm daughter loop + + if (flag != 0) { + origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, false, &idxBhadMothers); + if (std::abs(yCharmBaryonGen) < kYCutTight) { + registry.fill(HIST("hGenCharmBaryonPtRapidityTight"), ptCharmBaryonGen); + } + if (std::abs(yCharmBaryonGen) < kYCutLoose) { + registry.fill(HIST("hGenCharmBaryonPtRapidityLoose"), ptCharmBaryonGen); + } + } + + // Check if charm is prompt or non-prompt + if (origin == RecoDecay::OriginType::NonPrompt) { + fillGenMcTableByDecayChannel(flag, debugGenCharmBaryon, debugGenCasc, debugGenV0, ptCharmBaryonGen, yCharmBaryonGen, origin, idxBhadMothers[0]); + } else { + fillGenMcTableByDecayChannel(flag, debugGenCharmBaryon, debugGenCasc, debugGenV0, ptCharmBaryonGen, yCharmBaryonGen, origin, -1); + } + + } // particle loop + + } // end of collision loop + + } // template run function + + void processMcEmpty(aod::Collisions const&) + { + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcEmpty, "Empty process function to prevent workflow from getting stuck", true); + + ///////////////////////////////////////////////////// + /// /// + /// Process functions with DCAFitter /// + /// /// + ///////////////////////////////////////////////////// + + //~~~~~~~~~~~~~~~~// + //~~~~To Xi Pi~~~~// + //~~~~~~~~~~~~~~~~// + void processMcXicToXiPiWithDCAFitterNoCent(aod::HfCandToXiPi const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsNoCents const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcXicToXiPiWithDCAFitterNoCent, "Perform MC matching of DCAFitter reconstructed Xic0 to Xi Pi. No cents", false); + + void processMcXicToXiPiWithDCAFitterCentFT0C(aod::HfCandToXiPi const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsFT0Cs const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcXicToXiPiWithDCAFitterCentFT0C, "Perform MC matching of DCAFitter reconstructed Xic0 to Xi Pi. Cents with FT0C", false); + + void processMcXicToXiPiWithDCAFitterCentFT0M(aod::HfCandToXiPi const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + McCollisionsCentFT0Ms const& mcCollisions, + McCollisionsFT0Ms const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcXicToXiPiWithDCAFitterCentFT0M, "Perform MC matching of DCAFitter reconstructed Xic0 to Xi Pi. Cents with FT0M", false); + + void processMcOmegacToXiPiWithDCAFitterNoCent(aod::HfCandToXiPi const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsNoCents const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToXiPiWithDCAFitterNoCent, "Perform MC matching of DCAFitter reconstructed Omegac0 to Xi Pi. No cents", false); + + void processMcOmegacToXiPiWithDCAFitterCentFT0C(aod::HfCandToXiPi const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsFT0Cs const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToXiPiWithDCAFitterCentFT0C, "Perform MC matching of DCAFitter reconstructed Omeagc0 to Xi Pi. Cents with FT0C", false); + + void processMcOmegacToXiPiWithDCAFitterCentFT0M(aod::HfCandToXiPi const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + McCollisionsCentFT0Ms const& mcCollisions, + McCollisionsFT0Ms const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToXiPiWithDCAFitterCentFT0M, "Perform MC matching of DCAFitter reconstructed Omegac0 to Xi Pi. Cents with FT0M", false); + + //~~~~~~~~~~~~~~~~~~~// + //~~~~To Omega Pi~~~~// + //~~~~~~~~~~~~~~~~~~~// + void processMcOmegacToOmegaPiWithDCAFitterNoCent(aod::HfCandToOmegaPi const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsNoCents const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaPiWithDCAFitterNoCent, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Pi. No cents", false); + + void processMcOmegacToOmegaPiWithDCAFitterCentFT0C(aod::HfCandToOmegaPi const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsFT0Cs const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaPiWithDCAFitterCentFT0C, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Pi. Cents with FT0C", false); + + void processMcOmegacToOmegaPiWithDCAFitterCentFT0M(aod::HfCandToOmegaPi const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + McCollisionsCentFT0Ms const& mcCollisions, + McCollisionsFT0Ms const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaPiWithDCAFitterCentFT0M, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Pi. Cents with FT0M", false); + + //~~~~~~~~~~~~~~~~~~// + //~~~~To Omega Ka~~~~// + //~~~~~~~~~~~~~~~~~~// + void processMcOmegacToOmegaKaWithDCAFitterNoCent(aod::HfCandToOmegaK const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsNoCents const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaKaWithDCAFitterNoCent, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Ka. No cents", false); + + void processMcOmegacToOmegaKaWithDCAFitterCentFT0C(aod::HfCandToOmegaK const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsFT0Cs const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaKaWithDCAFitterCentFT0C, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Ka. Cents with FT0C", false); + + void processMcOmegacToOmegaKaWithDCAFitterCentFT0M(aod::HfCandToOmegaK const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + McCollisionsCentFT0Ms const& mcCollisions, + McCollisionsFT0Ms const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaKaWithDCAFitterCentFT0M, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Ka. Cents with FT0M", false); + + ///////////////////////////////////////////////////// + /// /// + /// Process functions with KFParticle /// + /// /// + ///////////////////////////////////////////////////// + void processMcXicToXiPiWithKFParticleNoCent(aod::HfCandToXiPiKf const& candidates, + TracksWMcIU const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsNoCents const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcXicToXiPiWithKFParticleNoCent, "Perform MC matching of KFParticle reconstructed Xic0 to Xi Pi. No cents", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/TableProducer/candidateSelectorToXiPiQa.cxx b/PWGHF/TableProducer/candidateSelectorToXiPiQa.cxx new file mode 100644 index 00000000000..1fadcabb6ae --- /dev/null +++ b/PWGHF/TableProducer/candidateSelectorToXiPiQa.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. + +/// \file candidateSelectorToXiPiQa.cxx +/// \brief Selection of Xic0 and Xicp candidates +/// +/// \author Jinhyun Park , Pusan National University +/// \author Krista Smith , Pusan National University + +#include "PWGHF/Core/HfMlResponseXic0ToXiPi.h" +#include "PWGHF/Core/HfMlResponseXic0ToXiPiKf.h" +#include "PWGHF/Core/SelectorCuts.h" +#include "PWGHF/DataModel/AliasTables.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/Utils/utilsAnalysis.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TrackSelectorPID.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::analysis; + +enum PidInfoStored { + PiFromLam = 0, + PrFromLam, + PiFromCasc, + PiFromCharm +}; + +enum { + doDcaFitter = 0, + doKfParticle +}; + +/// Struct for applying Omegac0/Xic0 selection cuts +struct HfCandidateSelectorToXiPiQa { + + // DCAFitter + Produces hfSelToXiPi; + // KFParticle + Produces hfSelToXiPiKf; + // ML selection - Filled with both DCAFitter and KFParticle + Produces hfMlToXiPi; + + // cuts from SelectorCuts.h - pT dependent cuts + Configurable> binsPt{"binsPt", std::vector{hf_cuts_xic_to_xi_pi::vecBinsPt}, "pT bin limits"}; + Configurable> cuts{"cuts", {hf_cuts_xic_to_xi_pi::Cuts[0], hf_cuts_xic_to_xi_pi::NBinsPt, hf_cuts_xic_to_xi_pi::NCutVars, hf_cuts_xic_to_xi_pi::labelsPt, hf_cuts_xic_to_xi_pi::labelsCutVar}, "Xic0 candidate selection per Pt Bin"}; + // ML inference + Configurable applyMl{"applyMl", false, "Flag to apply ML selections"}; + Configurable> binsPtMl{"binsPtMl", std::vector{hf_cuts_ml::vecBinsPt}, "pT bin limits for ML application"}; + Configurable> cutDirMl{"cutDirMl", std::vector{hf_cuts_ml::vecCutDir}, "Whether to reject score values greater or smaller than the threshold"}; + Configurable> cutsMl{"cutsMl", {hf_cuts_ml::Cuts[0], hf_cuts_ml::NBinsPt, hf_cuts_ml::NCutScores, hf_cuts_ml::labelsPt, hf_cuts_ml::labelsCutScore}, "ML selections per pT bin"}; + Configurable nClassesMl{"nClassesMl", static_cast(hf_cuts_ml::NCutScores), "Number of classes in ML model"}; + Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"feature1", "feature2"}, "Names of ML model input features"}; + // CCDB configuration + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable> modelPathsCCDB{"modelPathsCCDB", std::vector{"EventFiltering/PWGHF/BDTXic0ToXipiKf"}, "Paths of models on CCDB"}; + Configurable> onnxFileNames{"onnxFileNames", std::vector{"ModelHandler_onnx_Xic0ToXipiKf.onnx"}, "ONNX file names for each pT bin (if not from CCDB full path)"}; + Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB"}; + Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; + // LF analysis selections + Configurable etaTrackCharmBachMax{"etaTrackCharmBachMax", 0.8, "Max absolute value of eta for charm baryon bachelor"}; + Configurable etaTrackLFDauMax{"etaTrackLFDauMax", 1.0, "Max absolute value of eta for V0 and cascade daughters"}; + Configurable ptPiFromCascMin{"ptPiFromCascMin", 0.15, "Min pT pi <--Casc"}; + Configurable ptPiFromCharmBaryonMin{"ptPiFromCharmBaryonMin", 0.2, "Min pT pi <--Casc"}; + Configurable radiusCascMin{"radiusCascMin", 0.5, "Min Cascade radius"}; + Configurable radiusV0Min{"radiusV0Min", 1.1, "Min V0 radius"}; + Configurable impactParXYPiFromCharmBaryonMin{"impactParXYPiFromCharmBaryonMin", 0., "Min dcaxy pi from charm baryon track to pV"}; + Configurable impactParXYPiFromCharmBaryonMax{"impactParXYPiFromCharmBaryonMax", 10., "Max dcaxy pi from charm baryon track to pV"}; + Configurable impactParXYCascMin{"impactParXYCascMin", 0., "Min dcaxy casc track to pV"}; + Configurable impactParXYCascMax{"impactParXYCascMax", 10., "Max dcaxy casc track to pV"}; + Configurable impactParZPiFromCharmBaryonMin{"impactParZPiFromCharmBaryonMin", 0., "Min dcaz pi from charm baryon track to pV"}; + Configurable impactParZPiFromCharmBaryonMax{"impactParZPiFromCharmBaryonMax", 10., "Max dcaz pi from charm baryon track to pV"}; + Configurable impactParZCascMin{"impactParZCascMin", 0., "Min dcaz casc track to pV"}; + Configurable impactParZCascMax{"impactParZCascMax", 10., "Max dcaz casc track to pV"}; + Configurable applyTrkSelLf{"applyTrkSelLf", true, "Apply track selection for LF daughters"}; + // Mass window + Configurable v0MassWindow{"v0MassWindow", 0.01, "v0 mass window"}; + Configurable cascMassWindow{"cascMassWindow", 0.01, "cascade mass window"}; + Configurable invMassCharmBaryonMin{"invMassCharmBaryonMin", 2.0, "Lower limit of invariant mass spectrum charm baryon"}; + Configurable invMassCharmBaryonMax{"invMassCharmBaryonMax", 3.1, "Lower limit of invariant mass spectrum charm baryon"}; + // PID options + Configurable usePidTpcOnly{"usePidTpcOnly", false, "Perform PID using only TPC"}; + Configurable usePidTpcTofCombined{"usePidTpcTofCombined", true, "Perform PID using TPC & TOF"}; + // PID - TPC selections + Configurable ptPiPidTpcMin{"ptPiPidTpcMin", -1, "Lower bound of track pT for TPC PID for pion selection"}; + Configurable ptPiPidTpcMax{"ptPiPidTpcMax", 9999.9, "Upper bound of track pT for TPC PID for pion selection"}; + Configurable nSigmaTpcPiMax{"nSigmaTpcPiMax", 3., "Nsigma cut on TPC only for pion selection"}; + Configurable nSigmaTpcCombinedPiMax{"nSigmaTpcCombinedPiMax", 0., "Nsigma cut on TPC combined with TOF for pion selection"}; + Configurable ptPrPidTpcMin{"ptPrPidTpcMin", -1, "Lower bound of track pT for TPC PID for proton selection"}; + Configurable ptPrPidTpcMax{"ptPrPidTpcMax", 9999.9, "Upper bound of track pT for TPC PID for proton selection"}; + Configurable nSigmaTpcPrMax{"nSigmaTpcPrMax", 3., "Nsigma cut on TPC only for proton selection"}; + Configurable nSigmaTpcCombinedPrMax{"nSigmaTpcCombinedPrMax", 0., "Nsigma cut on TPC combined with TOF for proton selection"}; + // PID - TOF selections + Configurable ptPiPidTofMin{"ptPiPidTofMin", -1, "Lower bound of track pT for TOF PID for pion selection"}; + Configurable ptPiPidTofMax{"ptPiPidTofMax", 9999.9, "Upper bound of track pT for TOF PID for pion selection"}; + Configurable nSigmaTofPiMax{"nSigmaTofPiMax", 3., "Nsigma cut on TOF only for pion selection"}; + Configurable nSigmaTofCombinedPiMax{"nSigmaTofCombinedPiMax", 0., "Nsigma cut on TOF combined with TPC for pion selection"}; + Configurable ptPrPidTofMin{"ptPrPidTofMin", -1, "Lower bound of track pT for TOF PID for proton selection"}; + Configurable ptPrPidTofMax{"ptPrPidTofMax", 9999.9, "Upper bound of track pT for TOF PID for proton selection"}; + Configurable nSigmaTofPrMax{"nSigmaTofPrMax", 5., "Nsigma cut on TOF only for proton selection"}; + Configurable nSigmaTofCombinedPrMax{"nSigmaTofCombinedPrMax", 0., "Nsigma cut on TOF combined with TPC for proton selection"}; + // detector track quality selections + Configurable nClustersTpcMin{"nClustersTpcMin", 70, "Minimum number of TPC clusters requirement"}; + Configurable nTpcCrossedRowsMin{"nTpcCrossedRowsMin", 70, "Minimum number of TPC crossed rows requirement"}; + Configurable tpcCrossedRowsOverFindableClustersRatioMin{"tpcCrossedRowsOverFindableClustersRatioMin", 0.8, "Minimum ratio TPC crossed rows over findable clusters requirement"}; + Configurable tpcChi2PerClusterMax{"tpcChi2PerClusterMax", 4, "Maximum value of chi2 fit over TPC clusters"}; + Configurable nClustersItsMin{"nClustersItsMin", 3, "Minimum number of ITS clusters requirement for pi <- charm baryon"}; + Configurable nClustersItsInnBarrMin{"nClustersItsInnBarrMin", 1, "Minimum number of ITS clusters in inner barrel requirement for pi <- charm baryon"}; + Configurable itsChi2PerClusterMax{"itsChi2PerClusterMax", 36, "Maximum value of chi2 fit over ITS clusters for pi <- charm baryon"}; + + o2::analysis::HfMlResponseXic0ToXiPi hfMlResponseDca; + o2::analysis::HfMlResponseXic0ToXiPiKf hfMlResponseKf; + std::vector outputMlXic0ToXiPi = {}; + o2::ccdb::CcdbApi ccdbApi; + + TrackSelectorPi selectorPion; + TrackSelectorPr selectorProton; + + using TracksSel = soa::Join; + // using TracksSelLf = soa::Join; + + HistogramRegistry registry{"registry"}; // for QA of selections + + void init(InitContext const&) + { + selectorPion.setRangePtTpc(ptPiPidTpcMin, ptPiPidTpcMax); + selectorPion.setRangeNSigmaTpc(-nSigmaTpcPiMax, nSigmaTpcPiMax); + selectorPion.setRangeNSigmaTpcCondTof(-nSigmaTpcCombinedPiMax, nSigmaTpcCombinedPiMax); + selectorPion.setRangePtTof(ptPiPidTofMin, ptPiPidTofMax); + selectorPion.setRangeNSigmaTof(-nSigmaTofPiMax, nSigmaTofPiMax); + selectorPion.setRangeNSigmaTofCondTpc(-nSigmaTofCombinedPiMax, nSigmaTofCombinedPiMax); + + selectorProton.setRangePtTpc(ptPrPidTpcMin, ptPrPidTpcMax); + selectorProton.setRangeNSigmaTpc(-nSigmaTpcPrMax, nSigmaTpcPrMax); + selectorProton.setRangeNSigmaTpcCondTof(-nSigmaTpcCombinedPrMax, nSigmaTpcCombinedPrMax); + selectorProton.setRangePtTof(ptPrPidTofMin, ptPrPidTofMax); + selectorProton.setRangeNSigmaTof(-nSigmaTofPrMax, nSigmaTofPrMax); + selectorProton.setRangeNSigmaTofCondTpc(-nSigmaTofCombinedPrMax, nSigmaTofCombinedPrMax); + + const AxisSpec axisSel{2, -0.5, 1.5, "status"}; + const AxisSpec axisSelOnLfDca{14, -0.5, 13.5, "status"}; + const AxisSpec axisSelOnLfKf{23, -0.5, 22.5, "status"}; + const AxisSpec axisSelOnHfDca{6, -0.5, 5.5, "status"}; + const AxisSpec axisSelOnHfKf{11, -0.5, 10.5, "status"}; + + registry.add("hSelSignDec", "hSelSignDec;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelStatusCluster", "hSelStatusCluster:# of events Passed;;", {HistType::kTH1F, {{6, -0.5, 5.5}}}); + registry.get(HIST("hSelStatusCluster"))->GetXaxis()->SetBinLabel(1, "All"); + registry.get(HIST("hSelStatusCluster"))->GetXaxis()->SetBinLabel(2, "TpcCluster PiFromV0"); + registry.get(HIST("hSelStatusCluster"))->GetXaxis()->SetBinLabel(3, "TpcCluster PrFromV0"); + registry.get(HIST("hSelStatusCluster"))->GetXaxis()->SetBinLabel(4, "TpcCluster PiFromCasc"); + registry.get(HIST("hSelStatusCluster"))->GetXaxis()->SetBinLabel(5, "TpcCluster PiFromCharm"); + registry.get(HIST("hSelStatusCluster"))->GetXaxis()->SetBinLabel(6, "ItsCluster PiFromCharm"); + + registry.add("hSelStatusPID", "hSelStatusPID;# of events Passed;;", {HistType::kTH1F, {{4, -0.5, 3.5}}}); + registry.get(HIST("hSelStatusPID"))->GetXaxis()->SetBinLabel(1, "All"); + registry.get(HIST("hSelStatusPID"))->GetXaxis()->SetBinLabel(2, "Lambda"); + registry.get(HIST("hSelStatusPID"))->GetXaxis()->SetBinLabel(3, "Cascade"); + registry.get(HIST("hSelStatusPID"))->GetXaxis()->SetBinLabel(4, "CharmBaryon"); + + // For QA of LF & HF selection + if (doprocessSelectionDCAFitter) { + registry.add("hSelStatusLf", "hSelStatusLf;# of candidate passed;", {HistType::kTH1F, {axisSelOnLfDca}}); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(1, "All"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(2, "etaV0Dau"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(3, "radiusV0"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(4, "radiusCasc"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(5, "cosPAV0"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(6, "cosPACasc"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(7, "dcaV0Dau"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(8, "dcaCascDau"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(9, "dcaXYToPvV0Dau0"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(10, "dcaXYToPvV0Dau1"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(11, "dcaXYToPvCascDau"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(12, "ptPiFromCasc"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(13, "impactParCascXY"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(14, "impactParCascZ"); + + registry.add("hSelStatusHf", "hSelStatusHf;# of candidate passed;", {HistType::kTH1F, {axisSelOnHfDca}}); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(1, "All"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(2, "etaTrackCharmBach"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(3, "dcaCharmBaryonDau"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(4, "ptPiFromCharmBaryon"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(5, "impactParBachFromCharmXY"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(6, "impactParBachFromCharmZ"); + } + + if (doprocessSelectionKFParticle) { + registry.add("hSelStatusLf", "hSelStatusLf;# of candidate passed;", {HistType::kTH1F, {axisSelOnLfKf}}); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(1, "All"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(2, "etaV0Dau"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(3, "radiusV0"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(4, "radiusCasc"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(5, "cosPAV0"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(6, "cosPACasc"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(7, "dcaV0Dau"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(8, "dcaCascDau"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(9, "dcaXYToPvV0Dau0"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(10, "dcaXYToPvV0Dau1"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(11, "dcaXYToPvCascDau"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(12, "ptPiFromCasc"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(13, "cosPAV0ToCasc"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(14, "kfDcaXYCascToPv"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(15, "chi2GeoV0"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(16, "chi2GeoCasc"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(17, "chi2TopoV0ToPv"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(18, "chi2TopoCascToPv"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(19, "chi2TopoV0ToCasc"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(20, "v0ldl"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(21, "cascldl"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(22, "decayLenXYLambda"); + registry.get(HIST("hSelStatusLf"))->GetXaxis()->SetBinLabel(23, "decayLenXYCasc"); + + registry.add("hSelStatusHf", "hSelStatusHf;# of candidate passed;", {HistType::kTH1F, {axisSelOnHfKf}}); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(1, "All"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(2, "etaTrackCharmBach"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(3, "dcaCharmBaryonDau"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(4, "ptPiFromCharmBaryon"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(5, "cosPaCascToXic"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(6, "kfDcaXYPiFromXic"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(7, "chi2GeoXic"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(8, "chi2TopoPiFromXicToPv"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(9, "chi2TopoCascToXic"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(10, "decayLenXYXic"); + registry.get(HIST("hSelStatusHf"))->GetXaxis()->SetBinLabel(11, "cTauXic"); + } + + // invarinat mass histograms + registry.add("hInvMassCharmBaryon", "Charm baryon invariant mass; int mass; entries", {HistType::kTH1F, {{1500, 1.5, 4.5}}}); + registry.add("hInvMassCharmBaryonBkg", "Charm baryon invariant mass, rejected; int mass; entries", {HistType::kTH1F, {{1500, 1.5, 4.5}}}); + + // HfMlResponse initialization + if (applyMl) { + if (doprocessSelectionDCAFitter) { + hfMlResponseDca.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); + if (loadModelsFromCCDB) { + ccdbApi.init(ccdbUrl); + hfMlResponseDca.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB, timestampCCDB); + } else { + hfMlResponseDca.setModelPathsLocal(onnxFileNames); + } + hfMlResponseDca.cacheInputFeaturesIndices(namesInputFeatures); + hfMlResponseDca.init(); + } else { + hfMlResponseKf.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); + if (loadModelsFromCCDB) { + ccdbApi.init(ccdbUrl); + hfMlResponseKf.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB, timestampCCDB); + } else { + hfMlResponseKf.setModelPathsLocal(onnxFileNames); + } + hfMlResponseKf.cacheInputFeaturesIndices(namesInputFeatures); + hfMlResponseKf.init(); + } + } + } + + // LF cuts - Cuts on LF tracks reco + // Selection on LF related informations + // returns true if all cuts are passed + template + bool SelectOnLF(const T& candidate, const int& inputPtBin) + { + + registry.fill(HIST("hSelStatusLf"), 0.0); + + // Eta selection of V0, Cascade daughters + double etaV0PosDau = candidate.etaV0PosDau(); + double etaV0NegDau = candidate.etaV0NegDau(); + double etaPiFromCasc = candidate.etaBachFromCasc(); + + if (std::abs(etaV0PosDau) > etaTrackLFDauMax || std::abs(etaV0NegDau) > etaTrackLFDauMax || std::abs(etaPiFromCasc) > etaTrackLFDauMax) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 1.0); + + // Minimum radius cut + double radiusV0 = RecoDecay::sqrtSumOfSquares(candidate.xDecayVtxV0(), candidate.yDecayVtxV0()); + double radiusCasc = RecoDecay::sqrtSumOfSquares(candidate.xDecayVtxCascade(), candidate.yDecayVtxCascade()); + + if (radiusV0 < radiusV0Min) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 2.0); + if (radiusCasc < radiusCascMin) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 3.0); + + // Cosine of pointing angle + if (candidate.cosPAV0() < cuts->get(inputPtBin, "cosPAV0")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 4.0); + if (candidate.cosPACasc() < cuts->get(inputPtBin, "cosPACasc")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 5.0); + + // Distance of Closest Approach(DCA) + if (candidate.dcaV0Dau() > cuts->get(inputPtBin, "dcaV0Dau")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 6.0); + + if (candidate.dcaCascDau() > cuts->get(inputPtBin, "dcaCascDau")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 7.0); + + if (std::abs(candidate.dcaXYToPvV0Dau0()) < cuts->get(inputPtBin, "dcaXYToPvV0Dau0")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 8.0); + + if (std::abs(candidate.dcaXYToPvV0Dau1()) < cuts->get(inputPtBin, "dcaXYToPvV0Dau1")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 9.0); + + if (std::abs(candidate.dcaXYToPvCascDau()) < cuts->get(inputPtBin, "dcaXYToPvCascDau")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 10.0); + + // pT: Bachelor + double ptPiFromCasc = RecoDecay::sqrtSumOfSquares(candidate.pxBachFromCasc(), candidate.pyBachFromCasc()); + if (ptPiFromCasc < ptPiFromCascMin) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 11.0); + + // Extra cuts for KFParticle + if constexpr (svReco == doKfParticle) { + // Cosine of Pointing angle + if (candidate.cosPaV0ToCasc() < cuts->get(inputPtBin, "cosPaV0ToCasc")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 12.0); + + // DCA + if (candidate.kfDcaXYCascToPv() > cuts->get(inputPtBin, "kfDcaXYCascToPv")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 13.0); + + // Chi2 + if (candidate.chi2GeoV0() < 0 || candidate.chi2GeoV0() > cuts->get(inputPtBin, "chi2GeoV0")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 14.0); + if (candidate.chi2GeoCasc() < 0 || candidate.chi2GeoCasc() > cuts->get(inputPtBin, "chi2GeoCasc")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 15.0); + if (candidate.chi2TopoV0ToPv() > 0 && candidate.chi2TopoV0ToPv() < cuts->get(inputPtBin, "chi2TopoV0ToPv")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 16.0); + if (candidate.chi2TopoCascToPv() < 0 || candidate.chi2TopoCascToPv() > cuts->get(inputPtBin, "chi2TopoCascToPv")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 17.0); + if (candidate.chi2TopoV0ToCasc() < 0 || candidate.chi2TopoV0ToCasc() > cuts->get(inputPtBin, "chi2TopoV0ToCasc")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 18.0); + + // ldl + if (candidate.v0ldl() < cuts->get(inputPtBin, "v0ldl")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 19.0); + if (candidate.cascldl() < cuts->get(inputPtBin, "cascldl")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 20.0); + + // Decay length + if (std::abs(candidate.decayLenXYLambda()) < cuts->get(inputPtBin, "decayLenXYLambda")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 21.0); + if (std::abs(candidate.decayLenXYCasc()) < cuts->get(inputPtBin, "decayLenXYCasc")) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 22.0); + + } else { + // Impact parameter(DCA?) + if (std::abs(candidate.impactParCascXY()) < impactParXYCascMin || std::abs(candidate.impactParCascXY()) > impactParXYCascMax) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 12.0); + if (std::abs(candidate.impactParCascZ()) < impactParZCascMin || std::abs(candidate.impactParCascZ()) > impactParZCascMax) { + return false; + } + registry.fill(HIST("hSelStatusLf"), 13.0); + } + + // If passes all cuts, return true + return true; + } + + // HF cuts - Cuts on Charm baryon reco + // Apply cuts with charm baryon & charm bachelor related informations + // returns true if all cuts are passed + template + bool SelectOnHF(const T& candidate, const int& inputPtBin) + { + + registry.fill(HIST("hSelStatusHf"), 0.0); + + // eta selection on charm bayron bachelor + if (candidate.etaBachFromCharmBaryon() > etaTrackCharmBachMax) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 1.0); + // Distance of Closest Approach(DCA) + if (candidate.dcaCharmBaryonDau() > cuts->get(inputPtBin, "dcaCharmBaryonDau")) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 2.0); + + // pT: Charm Bachelor + double ptPiFromCharmBaryon = RecoDecay::pt(candidate.pxBachFromCharmBaryon(), candidate.pyBachFromCharmBaryon()); + if (ptPiFromCharmBaryon < cuts->get("ptPiFromCharmBaryon")) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 3.0); + + // specific selections with KFParticle output + if constexpr (svReco == doKfParticle) { + // Cosine of pointing angle + if (candidate.cosPaCascToXic() < cuts->get(inputPtBin, "cosPaCascToXic")) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 4.0); + + // DCA + if (std::abs(candidate.kfDcaXYPiFromXic()) > cuts->get(inputPtBin, "kfDcaXYPiFromXic")) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 5.0); + + // Chi2 + if (candidate.chi2GeoXic() < 0 || candidate.chi2GeoXic() > cuts->get(inputPtBin, "chi2GeoXic")) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 6.0); + if (candidate.chi2TopoPiFromXicToPv() < 0 || candidate.chi2TopoPiFromXicToPv() > cuts->get(inputPtBin, "chi2TopoXicToPv")) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 7.0); + if (candidate.chi2TopoCascToXic() < 0 || candidate.chi2TopoCascToXic() > cuts->get(inputPtBin, "chi2TopoCascToXic")) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 8.0); + + // Decay Length + if (std::abs(candidate.decayLenXYXic()) > cuts->get(inputPtBin, "decayLenXYXic")) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 9.0); + + // ctau + if (std::abs(candidate.cTauXic()) > cuts->get(inputPtBin, "cTauXic")) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 10.0); + } else { + // Impact parameter(DCA?) + if (std::abs(candidate.impactParBachFromCharmBaryonXY()) < impactParXYPiFromCharmBaryonMin || std::abs(candidate.impactParBachFromCharmBaryonXY()) > impactParXYPiFromCharmBaryonMax) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 4.0); + if (std::abs(candidate.impactParBachFromCharmBaryonZ()) < impactParZPiFromCharmBaryonMin || std::abs(candidate.impactParBachFromCharmBaryonZ()) > impactParZPiFromCharmBaryonMax) { + return false; + } + registry.fill(HIST("hSelStatusHf"), 5.0); + } + + // If passes all cuts, return true + return true; + } + + template + void runSelection(TCandTable const& candidates, + TracksSel const&) + { + // looping over charm baryon candidates + for (const auto& candidate : candidates) { + + bool resultSelections = true; // True if the candidate passes all the selections, False otherwise + outputMlXic0ToXiPi.clear(); + + auto trackV0PosDau = candidate.template posTrack_as(); + auto trackV0NegDau = candidate.template negTrack_as(); + auto trackPiFromCasc = candidate.template bachelor_as(); + auto trackPiFromCharm = candidate.template bachelorFromCharmBaryon_as(); + + auto trackPiFromLam = trackV0NegDau; + auto trackPrFromLam = trackV0PosDau; + + int8_t signDecay = candidate.signDecay(); // sign of pi <- cascade + + if (signDecay > 0) { + trackPiFromLam = trackV0PosDau; + trackPrFromLam = trackV0NegDau; + registry.fill(HIST("hSelSignDec"), 1); // anti-particle decay + } else { + registry.fill(HIST("hSelSignDec"), 0); // particle decay + } + + // pT selection + auto ptCandXic0 = RecoDecay::pt(candidate.pxCharmBaryon(), candidate.pyCharmBaryon()); + int pTBin = findBin(binsPt, ptCandXic0); + if (pTBin == -1) { + resultSelections = false; + } + + // Topological selection + const bool selectionResOnLF = SelectOnLF(candidate, pTBin); + const bool selectionResOnHF = SelectOnHF(candidate, pTBin); + if (!selectionResOnLF || !selectionResOnHF) { + resultSelections = false; + } + + // TPC clusters selections + if (resultSelections) { + registry.fill(HIST("hSelStatusCluster"), 0.0); + } + if (applyTrkSelLf) { + if (!isSelectedTrackTpcQuality(trackPiFromLam, nClustersTpcMin, nTpcCrossedRowsMin, tpcCrossedRowsOverFindableClustersRatioMin, tpcChi2PerClusterMax)) { + resultSelections = false; + } else { + if (resultSelections) { + registry.fill(HIST("hSelStatusCluster"), 1.0); + } + } + + if (!isSelectedTrackTpcQuality(trackPrFromLam, nClustersTpcMin, nTpcCrossedRowsMin, tpcCrossedRowsOverFindableClustersRatioMin, tpcChi2PerClusterMax)) { + resultSelections = false; + } else { + if (resultSelections) { + registry.fill(HIST("hSelStatusCluster"), 2.0); + } + } + + if (!isSelectedTrackTpcQuality(trackPiFromCasc, nClustersTpcMin, nTpcCrossedRowsMin, tpcCrossedRowsOverFindableClustersRatioMin, tpcChi2PerClusterMax)) { + resultSelections = false; + } else { + if (resultSelections) { + registry.fill(HIST("hSelStatusCluster"), 3.0); + } + } + } + + if (!isSelectedTrackTpcQuality(trackPiFromCharm, nClustersTpcMin, nTpcCrossedRowsMin, tpcCrossedRowsOverFindableClustersRatioMin, tpcChi2PerClusterMax)) { + resultSelections = false; + } else { + if (resultSelections) { + registry.fill(HIST("hSelStatusCluster"), 4.0); + } + } + + // ITS clusters selection + if (!isSelectedTrackItsQuality(trackPiFromCharm, nClustersItsMin, itsChi2PerClusterMax) || trackPiFromCharm.itsNClsInnerBarrel() < nClustersItsInnBarrMin) { + resultSelections = false; + } else { + if (resultSelections) { + registry.fill(HIST("hSelStatusCluster"), 5.0); + } + } + + // Track level PID selection + if (resultSelections) { + registry.fill(HIST("hSelStatusPID"), 0.0); + } + int statusPidPrFromLam = -999; + int statusPidPiFromLam = -999; + int statusPidPiFromCasc = -999; + int statusPidPiFromCharmBaryon = -999; + + int infoTpcStored = 0; + int infoTofStored = 0; + + if (usePidTpcOnly == usePidTpcTofCombined) { + LOGF(fatal, "Check the PID configurables, usePidTpcOnly and usePidTpcTofCombined can't have the same value"); + } + + if (trackPiFromLam.hasTPC()) { + SETBIT(infoTpcStored, PiFromLam); + } + if (trackPrFromLam.hasTPC()) { + SETBIT(infoTpcStored, PrFromLam); + } + if (trackPiFromCasc.hasTPC()) { + SETBIT(infoTpcStored, PiFromCasc); + } + if (trackPiFromCharm.hasTPC()) { + SETBIT(infoTpcStored, PiFromCharm); + } + if (trackPiFromLam.hasTOF()) { + SETBIT(infoTofStored, PiFromLam); + } + if (trackPrFromLam.hasTOF()) { + SETBIT(infoTofStored, PrFromLam); + } + if (trackPiFromCasc.hasTOF()) { + SETBIT(infoTofStored, PiFromCasc); + } + if (trackPiFromCharm.hasTOF()) { + SETBIT(infoTofStored, PiFromCharm); + } + + if (usePidTpcOnly) { + statusPidPrFromLam = selectorProton.statusTpc(trackPrFromLam); + statusPidPiFromLam = selectorPion.statusTpc(trackPiFromLam); + statusPidPiFromCasc = selectorPion.statusTpc(trackPiFromCasc); + statusPidPiFromCharmBaryon = selectorPion.statusTpc(trackPiFromCharm); + } else if (usePidTpcTofCombined) { + statusPidPrFromLam = selectorProton.statusTpcOrTof(trackPrFromLam); + statusPidPiFromLam = selectorPion.statusTpcOrTof(trackPiFromLam); + statusPidPiFromCasc = selectorPion.statusTpcOrTof(trackPiFromCasc); + statusPidPiFromCharmBaryon = selectorPion.statusTpcOrTof(trackPiFromCharm); + } + + bool statusPidLambda = (statusPidPrFromLam == TrackSelectorPID::Accepted) && (statusPidPiFromLam == TrackSelectorPID::Accepted); + if (statusPidLambda && resultSelections) { + registry.fill(HIST("hSelStatusPID"), 1.0); + } + bool statusPidCascade = (statusPidLambda && statusPidPiFromCasc == TrackSelectorPID::Accepted); + if (statusPidCascade && resultSelections) { + registry.fill(HIST("hSelStatusPID"), 2.0); + } + bool statusPidCharmBaryon = (statusPidCascade && statusPidPiFromCharmBaryon == TrackSelectorPID::Accepted); + if (statusPidCharmBaryon) { + if (resultSelections) { + registry.fill(HIST("hSelStatusPID"), 3.0); + } + } else { + resultSelections = false; + } + + // invariant mass cuts + bool statusInvMassLambda = true; + bool statusInvMassCascade = true; + bool statusInvMassCharmBaryon = true; + + double invMassLambda = candidate.invMassLambda(); + double invMassCascade = candidate.invMassCascade(); + double invMassCharmBaryon = candidate.invMassCharmBaryon(); + + if ((invMassLambda - o2::constants::physics::MassLambda0) > v0MassWindow) { + statusInvMassLambda = false; + } + if ((invMassCascade - o2::constants::physics::MassXiMinus) > cascMassWindow) { + statusInvMassCascade = false; + } + if ((invMassCharmBaryon < invMassCharmBaryonMin) || (invMassCharmBaryon > invMassCharmBaryonMax)) { + statusInvMassCharmBaryon = false; + } + + // ML BDT selection + if (applyMl) { + bool isSelectedMlXic0 = false; + std::vector inputFeaturesXic0 = {}; + if constexpr (svReco == doDcaFitter) { + inputFeaturesXic0 = hfMlResponseDca.getInputFeatures(candidate, trackPiFromLam, trackPiFromCasc, trackPiFromCharm); + isSelectedMlXic0 = hfMlResponseDca.isSelectedMl(inputFeaturesXic0, ptCandXic0, outputMlXic0ToXiPi); + } else { + inputFeaturesXic0 = hfMlResponseKf.getInputFeatures(candidate, trackPiFromLam, trackPiFromCasc, trackPiFromCharm); + isSelectedMlXic0 = hfMlResponseKf.isSelectedMl(inputFeaturesXic0, ptCandXic0, outputMlXic0ToXiPi); + } + if (!isSelectedMlXic0) { + continue; + } + hfMlToXiPi(outputMlXic0ToXiPi); + } + + // Fill in selection result + if constexpr (svReco == doDcaFitter) { + hfSelToXiPi(statusPidLambda, statusPidCascade, statusPidCharmBaryon, statusInvMassLambda, statusInvMassCascade, statusInvMassCharmBaryon, resultSelections, infoTpcStored, infoTofStored, + trackPiFromCharm.tpcNSigmaPi(), trackPiFromCasc.tpcNSigmaPi(), trackPiFromLam.tpcNSigmaPi(), trackPrFromLam.tpcNSigmaPr(), + trackPiFromCharm.tofNSigmaPi(), trackPiFromCasc.tofNSigmaPi(), trackPiFromLam.tofNSigmaPi(), trackPrFromLam.tofNSigmaPr()); + } else { + hfSelToXiPiKf(resultSelections, + trackPiFromCharm.tpcNSigmaPi(), trackPiFromCasc.tpcNSigmaPi(), trackPiFromLam.tpcNSigmaPi(), trackPrFromLam.tpcNSigmaPr(), + trackPiFromCharm.tofNSigmaPi(), trackPiFromCasc.tofNSigmaPi(), trackPiFromLam.tofNSigmaPi(), trackPrFromLam.tofNSigmaPr()); + } + + // Fill in invariant mass histogram + if (statusPidLambda && statusPidCascade && statusPidCharmBaryon && statusInvMassLambda && statusInvMassCascade && statusInvMassCharmBaryon && resultSelections) { + registry.fill(HIST("hInvMassCharmBaryon"), invMassCharmBaryon); + } else { + registry.fill(HIST("hInvMassCharmBaryonBkg"), invMassCharmBaryon); + } + + } // end of candidate loop + } // end run fuction + + /////////////////////////////////// + /// Process with DCAFitter // + /////////////////////////////////// + void processSelectionDCAFitter(aod::HfCandToXiPi const& candidates, TracksSel const& tracks) + { + runSelection(candidates, tracks); + } + PROCESS_SWITCH(HfCandidateSelectorToXiPiQa, processSelectionDCAFitter, "Xic0 candidate selection with DCAFitter output", true); + + //////////////////////////////////// + /// Process with KFParticle // + //////////////////////////////////// + void processSelectionKFParticle(aod::HfCandToXiPiKf const& candidates, TracksSel const& tracks) + { + runSelection(candidates, tracks); + } + PROCESS_SWITCH(HfCandidateSelectorToXiPiQa, processSelectionKFParticle, "Xic0 candidate selection with KFParticle output", false); + +}; // end struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/TableProducer/treeCreatorToXiPiQa.cxx b/PWGHF/TableProducer/treeCreatorToXiPiQa.cxx new file mode 100644 index 00000000000..72f2b9b7c71 --- /dev/null +++ b/PWGHF/TableProducer/treeCreatorToXiPiQa.cxx @@ -0,0 +1,977 @@ +// 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. + +/// \file treeCreatorToXiPiQa.cxx +/// \brief Writer of the omegac0 or xic0 to Xi Pi candidates in the form of flat tables to be stored in TTrees. +/// In this file are defined and filled the output tables +/// +/// \author Jinhyun Park , Pusan National University +/// \author Krista Smith , Pusan National University + +#include "PWGHF/Core/CentralityEstimation.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace o2; +using namespace o2::framework; + +// SV Reco method +enum { + DCAFITTER = 0, + KFPARTICLE +}; + +// Table size +enum { + FULL = 0, + LITE +}; + +namespace o2::aod +{ +namespace full +{ +// collision info +DECLARE_SOA_COLUMN(IsEventSel8, isEventSel8, bool); +DECLARE_SOA_COLUMN(IsEventSelZ, isEventSelZ, bool); +DECLARE_SOA_COLUMN(Centrality, centrality, float); +// from creator +DECLARE_SOA_COLUMN(XPv, xPv, float); +DECLARE_SOA_COLUMN(YPv, yPv, float); +DECLARE_SOA_COLUMN(ZPv, zPv, float); +DECLARE_SOA_COLUMN(XDecayVtxCharmBaryon, xDecayVtxCharmBaryon, float); +DECLARE_SOA_COLUMN(YDecayVtxCharmBaryon, yDecayVtxCharmBaryon, float); +DECLARE_SOA_COLUMN(ZDecayVtxCharmBaryon, zDecayVtxCharmBaryon, float); +DECLARE_SOA_COLUMN(XDecayVtxCascade, xDecayVtxCascade, float); +DECLARE_SOA_COLUMN(YDecayVtxCascade, yDecayVtxCascade, float); +DECLARE_SOA_COLUMN(ZDecayVtxCascade, zDecayVtxCascade, float); +DECLARE_SOA_COLUMN(XDecayVtxV0, xDecayVtxV0, float); +DECLARE_SOA_COLUMN(YDecayVtxV0, yDecayVtxV0, float); +DECLARE_SOA_COLUMN(ZDecayVtxV0, zDecayVtxV0, float); +DECLARE_SOA_COLUMN(SignDecay, signDecay, int8_t); // sign of pi <- xi +DECLARE_SOA_COLUMN(CovVtxCharmBaryonXX, covVtxCharmBaryonXX, float); +DECLARE_SOA_COLUMN(CovVtxCharmBaryonYY, covVtxCharmBaryonYY, float); +DECLARE_SOA_COLUMN(CovVtxCharmBaryonZZ, covVtxCharmBaryonZZ, float); +DECLARE_SOA_COLUMN(PxCharmBaryon, pxCharmBaryon, float); +DECLARE_SOA_COLUMN(PyCharmBaryon, pyCharmBaryon, float); +DECLARE_SOA_COLUMN(PzCharmBaryon, pzCharmBaryon, float); +DECLARE_SOA_COLUMN(PxCasc, pxCasc, float); +DECLARE_SOA_COLUMN(PyCasc, pyCasc, float); +DECLARE_SOA_COLUMN(PzCasc, pzCasc, float); +DECLARE_SOA_COLUMN(PxPiFromCharmBaryon, pxPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(PyPiFromCharmBaryon, pyPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(PzPiFromCharmBaryon, pzPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(PxLambda, pxLambda, float); +DECLARE_SOA_COLUMN(PyLambda, pyLambda, float); +DECLARE_SOA_COLUMN(PzLambda, pzLambda, float); +DECLARE_SOA_COLUMN(PxPiFromCasc, pxPiFromCasc, float); +DECLARE_SOA_COLUMN(PyPiFromCasc, pyPiFromCasc, float); +DECLARE_SOA_COLUMN(PzPiFromCasc, pzPiFromCasc, float); +DECLARE_SOA_COLUMN(PxPosV0Dau, pxPosV0Dau, float); +DECLARE_SOA_COLUMN(PyPosV0Dau, pyPosV0Dau, float); +DECLARE_SOA_COLUMN(PzPosV0Dau, pzPosV0Dau, float); +DECLARE_SOA_COLUMN(PxNegV0Dau, pxNegV0Dau, float); +DECLARE_SOA_COLUMN(PyNegV0Dau, pyNegV0Dau, float); +DECLARE_SOA_COLUMN(PzNegV0Dau, pzNegV0Dau, float); +DECLARE_SOA_COLUMN(ImpactParCascXY, impactParCascXY, float); +DECLARE_SOA_COLUMN(ImpactParPiFromCharmBaryonXY, impactParPiFromCharmBaryonXY, float); +DECLARE_SOA_COLUMN(ImpactParCascZ, impactParCascZ, float); +DECLARE_SOA_COLUMN(ImpactParPiFromCharmBaryonZ, impactParPiFromCharmBaryonZ, float); +DECLARE_SOA_COLUMN(ErrImpactParCascXY, errImpactParCascXY, float); +DECLARE_SOA_COLUMN(ErrImpactParPiFromCharmBaryonXY, errImpactParPiFromCharmBaryonXY, float); +DECLARE_SOA_COLUMN(InvMassLambda, invMassLambda, float); +DECLARE_SOA_COLUMN(InvMassCascade, invMassCascade, float); +DECLARE_SOA_COLUMN(InvMassCharmBaryon, invMassCharmBaryon, float); +DECLARE_SOA_COLUMN(CosPAV0, cosPAV0, float); +DECLARE_SOA_COLUMN(CosPACharmBaryon, cosPACharmBaryon, float); +DECLARE_SOA_COLUMN(CosPACasc, cosPACasc, float); +DECLARE_SOA_COLUMN(CosPAXYV0, cosPAXYV0, float); +DECLARE_SOA_COLUMN(CosPAXYCharmBaryon, cosPAXYCharmBaryon, float); +DECLARE_SOA_COLUMN(CosPAXYCasc, cosPAXYCasc, float); +DECLARE_SOA_COLUMN(CTauOmegac, cTauOmegac, float); +DECLARE_SOA_COLUMN(CTauCascade, cTauCascade, float); +DECLARE_SOA_COLUMN(CTauV0, cTauV0, float); +DECLARE_SOA_COLUMN(CTauXic, cTauXic, float); +DECLARE_SOA_COLUMN(EtaV0PosDau, etaV0PosDau, float); +DECLARE_SOA_COLUMN(EtaV0NegDau, etaV0NegDau, float); +DECLARE_SOA_COLUMN(EtaPiFromCasc, etaPiFromCasc, float); +DECLARE_SOA_COLUMN(EtaPiFromCharmBaryon, etaPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(EtaCharmBaryon, etaCharmBaryon, float); +DECLARE_SOA_COLUMN(EtaCascade, etaCascade, float); +DECLARE_SOA_COLUMN(EtaV0, etaV0, float); +DECLARE_SOA_COLUMN(DcaXYToPvV0Dau0, dcaXYToPvV0Dau0, float); +DECLARE_SOA_COLUMN(DcaXYToPvV0Dau1, dcaXYToPvV0Dau1, float); +DECLARE_SOA_COLUMN(DcaXYToPvCascDau, dcaXYToPvCascDau, float); +DECLARE_SOA_COLUMN(DcaZToPvV0Dau0, dcaZToPvV0Dau0, float); +DECLARE_SOA_COLUMN(DcaZToPvV0Dau1, dcaZToPvV0Dau1, float); +DECLARE_SOA_COLUMN(DcaZToPvCascDau, dcaZToPvCascDau, float); +DECLARE_SOA_COLUMN(DcaCascDau, dcaCascDau, float); +DECLARE_SOA_COLUMN(DcaV0Dau, dcaV0Dau, float); +DECLARE_SOA_COLUMN(DcaCharmBaryonDau, dcaCharmBaryonDau, float); +DECLARE_SOA_COLUMN(DecLenCharmBaryon, decLenCharmBaryon, float); +DECLARE_SOA_COLUMN(DecLenCascade, decLenCascade, float); +DECLARE_SOA_COLUMN(DecLenV0, decLenV0, float); +DECLARE_SOA_COLUMN(ErrorDecayLengthCharmBaryon, errorDecayLengthCharmBaryon, float); +DECLARE_SOA_COLUMN(ErrorDecayLengthXYCharmBaryon, errorDecayLengthXYCharmBaryon, float); +DECLARE_SOA_COLUMN(NormImpParCascade, normImpParCascade, double); +DECLARE_SOA_COLUMN(NormImpParPiFromCharmBar, normImpParPiFromCharmBar, double); +DECLARE_SOA_COLUMN(NormDecayLenCharmBar, normDecayLenCharmBar, double); +DECLARE_SOA_COLUMN(IsPionGlbTrkWoDca, isPionGlbTrkWoDca, bool); +DECLARE_SOA_COLUMN(PionItsNCls, pionItsNCls, uint8_t); +DECLARE_SOA_COLUMN(NTpcRowsPion, nTpcRowsPion, int16_t); +DECLARE_SOA_COLUMN(NTpcRowsPiFromCasc, nTpcRowsPiFromCasc, int16_t); +DECLARE_SOA_COLUMN(NTpcRowsPosV0Dau, nTpcRowsPosV0Dau, int16_t); +DECLARE_SOA_COLUMN(NTpcRowsNegV0Dau, nTpcRowsNegV0Dau, int16_t); +// from creator KF +DECLARE_SOA_COLUMN(KfDcaXYPiFromXic, kfDcaXYPiFromXic, float); +DECLARE_SOA_COLUMN(KfDcaXYCascToPv, kfDcaXYCascToPv, float); +DECLARE_SOA_COLUMN(Chi2GeoV0, chi2GeoV0, float); +DECLARE_SOA_COLUMN(Chi2GeoCasc, chi2GeoCasc, float); +DECLARE_SOA_COLUMN(Chi2GeoXic, chi2GeoXic, float); +DECLARE_SOA_COLUMN(Chi2MassV0, chi2MassV0, float); +DECLARE_SOA_COLUMN(Chi2MassCasc, chi2MassCasc, float); +DECLARE_SOA_COLUMN(V0ldl, v0ldl, float); +DECLARE_SOA_COLUMN(Cascldl, cascldl, float); +DECLARE_SOA_COLUMN(Xicldl, xicldl, float); +DECLARE_SOA_COLUMN(Chi2TopoV0ToPv, chi2TopoV0ToPv, float); +DECLARE_SOA_COLUMN(Chi2TopoCascToPv, chi2TopoCascToPv, float); +DECLARE_SOA_COLUMN(Chi2TopoPiFromXicToPv, chi2TopoPiFromXicToPv, float); +DECLARE_SOA_COLUMN(Chi2TopoXicToPv, chi2TopoXicToPv, float); +DECLARE_SOA_COLUMN(Chi2TopoV0ToCasc, chi2TopoV0ToCasc, float); +DECLARE_SOA_COLUMN(Chi2TopoCascToXic, chi2TopoCascToXic, float); +DECLARE_SOA_COLUMN(DecayLenXYLambda, decayLenXYLambda, float); +DECLARE_SOA_COLUMN(DecayLenXYCasc, decayLenXYCasc, float); +DECLARE_SOA_COLUMN(DecayLenXYXic, decayLenXYXic, float); +DECLARE_SOA_COLUMN(CosPaV0ToCasc, cosPaV0ToCasc, float); +DECLARE_SOA_COLUMN(CosPaV0ToPv, cosPaV0ToPv, float); +DECLARE_SOA_COLUMN(CosPaCascToXic, cosPaCascToXic, float); +DECLARE_SOA_COLUMN(CosPaCascToPv, cosPaCascToPv, float); +DECLARE_SOA_COLUMN(CosPaXicToPv, cosPaXicToPv, float); +DECLARE_SOA_COLUMN(KfRapXic, kfRapXic, float); +DECLARE_SOA_COLUMN(KfptPiFromXic, kfptPiFromXic, float); +DECLARE_SOA_COLUMN(KfptXic, kfptXic, float); +DECLARE_SOA_COLUMN(CosThetaStarPiFromXic, cosThetaStarPiFromXic, float); +DECLARE_SOA_COLUMN(CtXic, ctXic, float); +DECLARE_SOA_COLUMN(EtaXic, etaXic, float); +DECLARE_SOA_COLUMN(V0Ndf, v0Ndf, float); +DECLARE_SOA_COLUMN(CascNdf, cascNdf, float); +DECLARE_SOA_COLUMN(XicNdf, xicNdf, float); +DECLARE_SOA_COLUMN(MassV0Ndf, massV0Ndf, float); +DECLARE_SOA_COLUMN(MassCascNdf, massCascNdf, float); +DECLARE_SOA_COLUMN(V0Chi2OverNdf, v0Chi2OverNdf, float); +DECLARE_SOA_COLUMN(CascChi2OverNdf, cascChi2OverNdf, float); +DECLARE_SOA_COLUMN(XicChi2OverNdf, xicChi2OverNdf, float); +DECLARE_SOA_COLUMN(MassV0Chi2OverNdf, massV0Chi2OverNdf, float); +DECLARE_SOA_COLUMN(MassCascChi2OverNdf, massCascChi2OverNdf, float); +// from creator - MC +DECLARE_SOA_COLUMN(FlagMcMatchRec, flagMcMatchRec, int8_t); // reconstruction level +DECLARE_SOA_COLUMN(DebugMcRec, debugMcRec, int8_t); // debug flag for mis-association reconstruction level +DECLARE_SOA_COLUMN(OriginRec, originRec, int8_t); +DECLARE_SOA_COLUMN(CollisionMatched, collisionMatched, bool); +// from selector +DECLARE_SOA_COLUMN(StatusPidLambda, statusPidLambda, bool); +DECLARE_SOA_COLUMN(StatusPidCascade, statusPidCascade, bool); +DECLARE_SOA_COLUMN(StatusPidCharmBaryon, statusPidCharmBaryon, bool); +DECLARE_SOA_COLUMN(StatusInvMassLambda, statusInvMassLambda, bool); +DECLARE_SOA_COLUMN(StatusInvMassCascade, statusInvMassCascade, bool); +DECLARE_SOA_COLUMN(StatusInvMassCharmBaryon, statusInvMassCharmBaryon, bool); +DECLARE_SOA_COLUMN(ResultSelections, resultSelections, bool); +DECLARE_SOA_COLUMN(PidTpcInfoStored, pidTpcInfoStored, int); +DECLARE_SOA_COLUMN(PidTofInfoStored, pidTofInfoStored, int); +DECLARE_SOA_COLUMN(TpcNSigmaPiFromCharmBaryon, tpcNSigmaPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(TpcNSigmaPiFromCasc, tpcNSigmaPiFromCasc, float); +DECLARE_SOA_COLUMN(TpcNSigmaPiFromLambda, tpcNSigmaPiFromLambda, float); +DECLARE_SOA_COLUMN(TpcNSigmaPrFromLambda, tpcNSigmaPrFromLambda, float); +DECLARE_SOA_COLUMN(TofNSigmaPiFromCharmBaryon, tofNSigmaPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(TofNSigmaPiFromCasc, tofNSigmaPiFromCasc, float); +DECLARE_SOA_COLUMN(TofNSigmaPiFromLambda, tofNSigmaPiFromLambda, float); +DECLARE_SOA_COLUMN(TofNSigmaPrFromLambda, tofNSigmaPrFromLambda, float); +} // namespace full + +DECLARE_SOA_TABLE(HfToXiPiEvs, "AOD", "HFTOXIPIEV", + full::IsEventSel8, full::IsEventSelZ); + +DECLARE_SOA_TABLE(HfToXiPiFulls, "AOD", "HFTOXIPIFULL", + full::XPv, full::YPv, full::ZPv, full::Centrality, collision::NumContrib, collision::Chi2, + full::XDecayVtxCharmBaryon, full::YDecayVtxCharmBaryon, full::ZDecayVtxCharmBaryon, + full::XDecayVtxCascade, full::YDecayVtxCascade, full::ZDecayVtxCascade, + full::XDecayVtxV0, full::YDecayVtxV0, full::ZDecayVtxV0, + full::SignDecay, + full::CovVtxCharmBaryonXX, full::CovVtxCharmBaryonYY, full::CovVtxCharmBaryonZZ, + full::PxCharmBaryon, full::PyCharmBaryon, full::PzCharmBaryon, + full::PxCasc, full::PyCasc, full::PzCasc, + full::PxPiFromCharmBaryon, full::PyPiFromCharmBaryon, full::PzPiFromCharmBaryon, + full::PxLambda, full::PyLambda, full::PzLambda, + full::PxPiFromCasc, full::PyPiFromCasc, full::PzPiFromCasc, + full::PxPosV0Dau, full::PyPosV0Dau, full::PzPosV0Dau, + full::PxNegV0Dau, full::PyNegV0Dau, full::PzNegV0Dau, + full::ImpactParCascXY, full::ImpactParPiFromCharmBaryonXY, + full::ImpactParCascZ, full::ImpactParPiFromCharmBaryonZ, + full::ErrImpactParCascXY, full::ErrImpactParPiFromCharmBaryonXY, + full::InvMassLambda, full::InvMassCascade, full::InvMassCharmBaryon, + full::CosPAV0, full::CosPACharmBaryon, full::CosPACasc, full::CosPAXYV0, full::CosPAXYCharmBaryon, full::CosPAXYCasc, + full::CTauOmegac, full::CTauCascade, full::CTauV0, full::CTauXic, + full::EtaV0PosDau, full::EtaV0NegDau, full::EtaPiFromCasc, full::EtaPiFromCharmBaryon, + full::EtaCharmBaryon, full::EtaCascade, full::EtaV0, + full::DcaXYToPvV0Dau0, full::DcaXYToPvV0Dau1, full::DcaXYToPvCascDau, + full::DcaZToPvV0Dau0, full::DcaZToPvV0Dau1, full::DcaZToPvCascDau, + full::DcaCascDau, full::DcaV0Dau, full::DcaCharmBaryonDau, + full::DecLenCharmBaryon, full::DecLenCascade, full::DecLenV0, full::ErrorDecayLengthCharmBaryon, full::ErrorDecayLengthXYCharmBaryon, + full::NormImpParCascade, full::NormImpParPiFromCharmBar, full::NormDecayLenCharmBar, full::IsPionGlbTrkWoDca, full::PionItsNCls, + full::NTpcRowsPion, full::NTpcRowsPiFromCasc, full::NTpcRowsPosV0Dau, full::NTpcRowsNegV0Dau, + full::StatusPidLambda, full::StatusPidCascade, full::StatusPidCharmBaryon, + full::StatusInvMassLambda, full::StatusInvMassCascade, full::StatusInvMassCharmBaryon, full::ResultSelections, + full::PidTpcInfoStored, full::PidTofInfoStored, + full::TpcNSigmaPiFromCharmBaryon, full::TpcNSigmaPiFromCasc, full::TpcNSigmaPiFromLambda, full::TpcNSigmaPrFromLambda, + full::TofNSigmaPiFromCharmBaryon, full::TofNSigmaPiFromCasc, full::TofNSigmaPiFromLambda, full::TofNSigmaPrFromLambda, + full::FlagMcMatchRec, full::DebugMcRec, full::OriginRec, full::CollisionMatched); + +DECLARE_SOA_TABLE(HfToXiPiLites, "AOD", "HFTOXIPILITE", + full::XPv, full::YPv, full::ZPv, full::Centrality, collision::NumContrib, collision::Chi2, + full::XDecayVtxCharmBaryon, full::YDecayVtxCharmBaryon, full::ZDecayVtxCharmBaryon, + full::XDecayVtxCascade, full::YDecayVtxCascade, full::ZDecayVtxCascade, + full::XDecayVtxV0, full::YDecayVtxV0, full::ZDecayVtxV0, + full::SignDecay, + full::PxCharmBaryon, full::PyCharmBaryon, full::PzCharmBaryon, + full::PxPiFromCharmBaryon, full::PyPiFromCharmBaryon, full::PzPiFromCharmBaryon, + full::PxPiFromCasc, full::PyPiFromCasc, full::PzPiFromCasc, + full::PxPosV0Dau, full::PyPosV0Dau, full::PzPosV0Dau, + full::PxNegV0Dau, full::PyNegV0Dau, full::PzNegV0Dau, + full::ImpactParCascXY, full::ImpactParPiFromCharmBaryonXY, + full::ErrImpactParCascXY, full::ErrImpactParPiFromCharmBaryonXY, + full::InvMassLambda, full::InvMassCascade, full::InvMassCharmBaryon, + full::EtaV0PosDau, full::EtaV0NegDau, full::EtaPiFromCasc, full::EtaPiFromCharmBaryon, + full::DcaXYToPvV0Dau0, full::DcaXYToPvV0Dau1, full::DcaXYToPvCascDau, + full::DcaCascDau, full::DcaV0Dau, full::DcaCharmBaryonDau, + full::ErrorDecayLengthCharmBaryon, full::NormImpParCascade, full::NormImpParPiFromCharmBar, + full::IsPionGlbTrkWoDca, full::PionItsNCls, + full::NTpcRowsPion, full::NTpcRowsPiFromCasc, full::NTpcRowsPosV0Dau, full::NTpcRowsNegV0Dau, + full::StatusPidLambda, full::StatusPidCascade, full::StatusPidCharmBaryon, + full::StatusInvMassLambda, full::StatusInvMassCascade, full::StatusInvMassCharmBaryon, full::ResultSelections, + full::PidTpcInfoStored, full::PidTofInfoStored, + full::TpcNSigmaPiFromCharmBaryon, full::TpcNSigmaPiFromCasc, full::TpcNSigmaPiFromLambda, full::TpcNSigmaPrFromLambda, + full::TofNSigmaPiFromCharmBaryon, full::TofNSigmaPiFromCasc, full::TofNSigmaPiFromLambda, full::TofNSigmaPrFromLambda, + full::FlagMcMatchRec, full::OriginRec, full::CollisionMatched); + +DECLARE_SOA_TABLE(HfKfXicFulls, "AOD", "HFKFXICFULL", + full::Centrality, + // full::StatusPidLambda, full::StatusPidCascade, full::StatusPidCharmBaryon, + // full::StatusInvMassLambda, full::StatusInvMassCascade, full::StatusInvMassCharmBaryon, + full::ResultSelections, + full::TpcNSigmaPiFromCharmBaryon, full::TofNSigmaPiFromCharmBaryon, full::TpcNSigmaPiFromCasc, full::TofNSigmaPiFromCasc, + full::TpcNSigmaPiFromLambda, full::TofNSigmaPiFromLambda, full::TpcNSigmaPrFromLambda, full::TofNSigmaPrFromLambda, + full::KfDcaXYPiFromXic, full::DcaCascDau, full::DcaCharmBaryonDau, full::KfDcaXYCascToPv, + full::DcaXYToPvV0Dau0, full::DcaXYToPvV0Dau1, full::DcaXYToPvCascDau, + full::Chi2GeoV0, full::Chi2GeoCasc, full::Chi2GeoXic, + full::Chi2MassV0, full::Chi2MassCasc, + full::V0ldl, full::Cascldl, // full::Xicldl, + full::Chi2TopoV0ToPv, full::Chi2TopoCascToPv, full::Chi2TopoPiFromXicToPv, full::Chi2TopoXicToPv, + full::Chi2TopoV0ToCasc, full::Chi2TopoCascToXic, + full::DecayLenXYLambda, full::DecayLenXYCasc, full::DecayLenXYXic, + full::CosPaV0ToCasc, full::CosPaV0ToPv, full::CosPaCascToXic, full::CosPaCascToPv, // full::CosPaXicToPv, + full::InvMassLambda, full::InvMassCascade, full::InvMassCharmBaryon, + full::KfRapXic, // full::KfptPiFromXic, full::KfptXic, + full::CosThetaStarPiFromXic, full::CtXic, full::EtaXic, + full::V0Ndf, full::CascNdf, full::XicNdf, + full::MassV0Ndf, full::MassCascNdf, + full::V0Chi2OverNdf, full::CascChi2OverNdf, full::XicChi2OverNdf, + full::MassV0Chi2OverNdf, full::MassCascChi2OverNdf, + full::FlagMcMatchRec, full::DebugMcRec, full::OriginRec, full::CollisionMatched); + +} // namespace o2::aod + +/// Writes the full information in an output TTree +struct HfTreeCreatorToXiPiQa { + + Produces rowCandidateFull; + Produces rowCandidateLite; + Produces rowKfCandidate; + Produces rowEv; + + Configurable zPvCut{"zPvCut", 10., "Cut on absolute value of primary vertex z coordinate"}; + + using MyTrackTable = soa::Join; + using MyEventTable = soa::Join; + using MyEventTableWithFT0C = soa::Join; + using MyEventTableWithFT0M = soa::Join; + using MyEventTableWithNTracksPV = soa::Join; + + void init(InitContext const&) + { + if ((doprocessMcLiteXic0 && doprocessMcLiteOmegac0) || (doprocessMcFullXic0 && doprocessMcFullOmegac0)) { + LOGF(fatal, "Both Xic0 and Omegac0 MC processes enabled, please choose ONLY one!"); + } + } + + ////////////////////////////////////////////////////// + // // + // Fill functions to fill in the tables // + // // + ////////////////////////////////////////////////////// + + template + void fillEvent(const T& collision, float cutZPv) + { + rowEv(collision.sel8(), std::abs(collision.posZ()) < cutZPv); + } + + template + void fillCandidate(const T& candidate, int8_t flagMc, int8_t debugMc, int8_t originMc, bool collisionMatched) + { + // save all candidate information + float centrality = -999.f; + if constexpr (useCentrality) { + const auto& collision = candidate.template collision_as(); + centrality = o2::hf_centrality::getCentralityColl(collision); + } + + if constexpr (svReco == DCAFITTER) { + if constexpr (tableSize == LITE) { + rowCandidateLite(candidate.xPv(), + candidate.yPv(), + candidate.zPv(), + centrality, + candidate.template collision_as().numContrib(), + candidate.template collision_as().chi2(), + candidate.xDecayVtxCharmBaryon(), + candidate.yDecayVtxCharmBaryon(), + candidate.zDecayVtxCharmBaryon(), + candidate.xDecayVtxCascade(), + candidate.yDecayVtxCascade(), + candidate.zDecayVtxCascade(), + candidate.xDecayVtxV0(), + candidate.yDecayVtxV0(), + candidate.zDecayVtxV0(), + candidate.signDecay(), + candidate.pxCharmBaryon(), + candidate.pyCharmBaryon(), + candidate.pzCharmBaryon(), + candidate.pxBachFromCharmBaryon(), + candidate.pyBachFromCharmBaryon(), + candidate.pzBachFromCharmBaryon(), + candidate.pxBachFromCasc(), + candidate.pyBachFromCasc(), + candidate.pzBachFromCasc(), + candidate.pxPosV0Dau(), + candidate.pyPosV0Dau(), + candidate.pzPosV0Dau(), + candidate.pxNegV0Dau(), + candidate.pyNegV0Dau(), + candidate.pzNegV0Dau(), + candidate.impactParCascXY(), + candidate.impactParBachFromCharmBaryonXY(), + candidate.errImpactParCascXY(), + candidate.errImpactParBachFromCharmBaryonXY(), + candidate.invMassLambda(), + candidate.invMassCascade(), + candidate.invMassCharmBaryon(), + candidate.etaV0PosDau(), + candidate.etaV0NegDau(), + candidate.etaBachFromCasc(), + candidate.etaBachFromCharmBaryon(), + candidate.dcaXYToPvV0Dau0(), + candidate.dcaXYToPvV0Dau1(), + candidate.dcaXYToPvCascDau(), + candidate.dcaCascDau(), + candidate.dcaV0Dau(), + candidate.dcaCharmBaryonDau(), + candidate.errorDecayLengthCharmBaryon(), + candidate.impactParCascXY() / candidate.errImpactParCascXY(), + candidate.impactParBachFromCharmBaryonXY() / candidate.errImpactParBachFromCharmBaryonXY(), + candidate.template bachelorFromCharmBaryon_as().isGlobalTrackWoDCA(), + candidate.template bachelorFromCharmBaryon_as().itsNCls(), + candidate.template bachelorFromCharmBaryon_as().tpcNClsCrossedRows(), + candidate.template bachelor_as().tpcNClsCrossedRows(), + candidate.template posTrack_as().tpcNClsCrossedRows(), + candidate.template negTrack_as().tpcNClsCrossedRows(), + candidate.statusPidLambda(), + candidate.statusPidCascade(), + candidate.statusPidCharmBaryon(), + candidate.statusInvMassLambda(), + candidate.statusInvMassCascade(), + candidate.statusInvMassCharmBaryon(), + candidate.resultSelections(), + candidate.pidTpcInfoStored(), + candidate.pidTofInfoStored(), + candidate.tpcNSigmaPiFromCharmBaryon(), + candidate.tpcNSigmaPiFromCasc(), + candidate.tpcNSigmaPiFromLambda(), + candidate.tpcNSigmaPrFromLambda(), + candidate.tofNSigmaPiFromCharmBaryon(), + candidate.tofNSigmaPiFromCasc(), + candidate.tofNSigmaPiFromLambda(), + candidate.tofNSigmaPrFromLambda(), + flagMc, + originMc, + collisionMatched); + } else { + rowCandidateFull(candidate.xPv(), + candidate.yPv(), + candidate.zPv(), + centrality, + candidate.template collision_as().numContrib(), + candidate.template collision_as().chi2(), + candidate.xDecayVtxCharmBaryon(), + candidate.yDecayVtxCharmBaryon(), + candidate.zDecayVtxCharmBaryon(), + candidate.xDecayVtxCascade(), + candidate.yDecayVtxCascade(), + candidate.zDecayVtxCascade(), + candidate.xDecayVtxV0(), + candidate.yDecayVtxV0(), + candidate.zDecayVtxV0(), + candidate.signDecay(), + candidate.covVtxCharmBaryon0(), + candidate.covVtxCharmBaryon3(), + candidate.covVtxCharmBaryon5(), + candidate.pxCharmBaryon(), + candidate.pyCharmBaryon(), + candidate.pzCharmBaryon(), + candidate.pxCasc(), + candidate.pyCasc(), + candidate.pzCasc(), + candidate.pxBachFromCharmBaryon(), + candidate.pyBachFromCharmBaryon(), + candidate.pzBachFromCharmBaryon(), + candidate.pxLambda(), + candidate.pyLambda(), + candidate.pzLambda(), + candidate.pxBachFromCasc(), + candidate.pyBachFromCasc(), + candidate.pzBachFromCasc(), + candidate.pxPosV0Dau(), + candidate.pyPosV0Dau(), + candidate.pzPosV0Dau(), + candidate.pxNegV0Dau(), + candidate.pyNegV0Dau(), + candidate.pzNegV0Dau(), + candidate.impactParCascXY(), + candidate.impactParBachFromCharmBaryonXY(), + candidate.impactParCascZ(), + candidate.impactParBachFromCharmBaryonZ(), + candidate.errImpactParCascXY(), + candidate.errImpactParBachFromCharmBaryonXY(), + candidate.invMassLambda(), + candidate.invMassCascade(), + candidate.invMassCharmBaryon(), + candidate.cosPAV0(), + candidate.cosPACharmBaryon(), + candidate.cosPACasc(), + candidate.cosPAXYV0(), + candidate.cosPAXYCharmBaryon(), + candidate.cosPAXYCasc(), + candidate.cTauOmegac(), + candidate.cTauCascade(), + candidate.cTauV0(), + candidate.cTauXic(), + candidate.etaV0PosDau(), + candidate.etaV0NegDau(), + candidate.etaBachFromCasc(), + candidate.etaBachFromCharmBaryon(), + candidate.etaCharmBaryon(), + candidate.etaCascade(), + candidate.etaV0(), + candidate.dcaXYToPvV0Dau0(), + candidate.dcaXYToPvV0Dau1(), + candidate.dcaXYToPvCascDau(), + candidate.dcaZToPvV0Dau0(), + candidate.dcaZToPvV0Dau1(), + candidate.dcaZToPvCascDau(), + candidate.dcaCascDau(), + candidate.dcaV0Dau(), + candidate.dcaCharmBaryonDau(), + candidate.decLenCharmBaryon(), + candidate.decLenCascade(), + candidate.decLenV0(), + candidate.errorDecayLengthCharmBaryon(), + candidate.errorDecayLengthXYCharmBaryon(), + candidate.impactParCascXY() / candidate.errImpactParCascXY(), + candidate.impactParBachFromCharmBaryonXY() / candidate.errImpactParBachFromCharmBaryonXY(), + candidate.decLenCharmBaryon() / candidate.errorDecayLengthCharmBaryon(), + candidate.template bachelorFromCharmBaryon_as().isGlobalTrackWoDCA(), + candidate.template bachelorFromCharmBaryon_as().itsNCls(), + candidate.template bachelorFromCharmBaryon_as().tpcNClsCrossedRows(), + candidate.template bachelor_as().tpcNClsCrossedRows(), + candidate.template posTrack_as().tpcNClsCrossedRows(), + candidate.template negTrack_as().tpcNClsCrossedRows(), + candidate.statusPidLambda(), + candidate.statusPidCascade(), + candidate.statusPidCharmBaryon(), + candidate.statusInvMassLambda(), + candidate.statusInvMassCascade(), + candidate.statusInvMassCharmBaryon(), + candidate.resultSelections(), + candidate.pidTpcInfoStored(), + candidate.pidTofInfoStored(), + candidate.tpcNSigmaPiFromCharmBaryon(), + candidate.tpcNSigmaPiFromCasc(), + candidate.tpcNSigmaPiFromLambda(), + candidate.tpcNSigmaPrFromLambda(), + candidate.tofNSigmaPiFromCharmBaryon(), + candidate.tofNSigmaPiFromCasc(), + candidate.tofNSigmaPiFromLambda(), + candidate.tofNSigmaPrFromLambda(), + flagMc, + debugMc, + originMc, + collisionMatched); + } + } else { + if constexpr (tableSize == LITE) { + // currently, no lite sized table for KFParticle + } else { + rowKfCandidate(centrality, + candidate.resultSelections(), + candidate.tpcNSigmaPiFromCharmBaryon(), + candidate.tofNSigmaPiFromCharmBaryon(), + candidate.tpcNSigmaPiFromCasc(), + candidate.tofNSigmaPiFromCasc(), + candidate.tpcNSigmaPiFromLambda(), + candidate.tofNSigmaPiFromLambda(), + candidate.tpcNSigmaPrFromLambda(), + candidate.tofNSigmaPrFromLambda(), + candidate.kfDcaXYPiFromXic(), + candidate.dcaCascDau(), + candidate.dcaCharmBaryonDau(), + candidate.kfDcaXYCascToPv(), + candidate.dcaXYToPvV0Dau0(), + candidate.dcaXYToPvV0Dau1(), + candidate.dcaXYToPvCascDau(), + candidate.chi2GeoV0(), + candidate.chi2GeoCasc(), + candidate.chi2GeoXic(), + candidate.chi2MassV0(), + candidate.chi2MassCasc(), + candidate.v0ldl(), + candidate.cascldl(), + candidate.chi2TopoV0ToPv(), + candidate.chi2TopoCascToPv(), + candidate.chi2TopoPiFromXicToPv(), + candidate.chi2TopoXicToPv(), + candidate.chi2TopoV0ToCasc(), + candidate.chi2TopoCascToXic(), + candidate.decayLenXYLambda(), + candidate.decayLenXYCasc(), + candidate.decayLenXYXic(), + candidate.cosPaV0ToCasc(), + candidate.cosPAV0(), + candidate.cosPaCascToXic(), + candidate.cosPACasc(), + candidate.invMassLambda(), + candidate.invMassCascade(), + candidate.invMassCharmBaryon(), + candidate.kfRapXic(), + candidate.cosThetaStarPiFromXic(), + candidate.cTauXic(), + candidate.etaCharmBaryon(), + candidate.v0Ndf(), + candidate.cascNdf(), + candidate.xicNdf(), + candidate.massV0Ndf(), + candidate.massCascNdf(), + candidate.v0Chi2OverNdf(), + candidate.cascChi2OverNdf(), + candidate.xicChi2OverNdf(), + candidate.massV0Chi2OverNdf(), + candidate.massCascChi2OverNdf(), + flagMc, + debugMc, + originMc, + collisionMatched); + } + } + } + + //////////////////////////////////// + // // + // Process functions // + // // + //////////////////////////////////// + + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + //*~~~~~~~Data with DCAFitter~~~~~~~~*// + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + + void processDataFull(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateFull.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + void processDataLite(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + void processDataLiteWithFT0M(MyEventTableWithFT0M const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + void processDataLiteWithFT0C(MyEventTableWithFT0C const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + void processDataLiteWithNTracksPV(MyEventTableWithNTracksPV const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processDataFull, "Process data with full information w/o centrality", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processDataLite, "Process data and produce lite table version", true); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processDataLiteWithFT0M, "Process data and produce lite table version with FT0M", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processDataLiteWithFT0C, "Process data and produce lite table version with FT0C", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processDataLiteWithNTracksPV, "Process data and produce lite table version with NTracksPV", false); + + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + //*~~~~~~~Data with KFParticle~~~~~~~~*// + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + void processKfData(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + void processKfDataWithFT0M(MyEventTableWithFT0M const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + void processKfDataWithFT0C(MyEventTableWithFT0C const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + void processKfDataWithNTracksPV(MyEventTableWithNTracksPV const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfData, "Process KF data, no cent", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfDataWithFT0M, "Process KF data, with FT0M", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfDataWithFT0C, "Process KF data, with FT0C", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfDataWithNTracksPV, "Process KF data, with NTracksPV", false); + + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + //*~~~~~~~MC with DCAFitter~~~~~~~~*// + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + + void processMcFullXic0(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateFull.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcFullOmegac0(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateFull.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcLiteXic0(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), -7, candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcLiteXic0WithFT0C(MyEventTableWithFT0C const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), -7, candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcLiteXic0WithFT0M(MyEventTableWithFT0M const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), -7, candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcLiteXic0WithNTracksPV(MyEventTableWithNTracksPV const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), -7, candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcLiteOmegac0(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), -7, candidate.originMcRec(), candidate.collisionMatched()); + } + } + + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcFullXic0, "Process MC with full information for xic0 w/o centrality", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcFullOmegac0, "Process MC with full information for omegac0", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcLiteXic0, "Process MC and produce lite table version for xic0", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcLiteXic0WithFT0C, "Process MC and produce lite table version for Xic0 with FT0C", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcLiteXic0WithFT0M, "Process MC and produce lite table version for Xic0 with FT0M", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcLiteXic0WithNTracksPV, "Process MC and produce lite table version for Xic0 with NTracksPV", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcLiteOmegac0, "Process MC and produce lite table version for omegac0", false); + + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + //*~~~~~~~MC with KFParticle~~~~~~~~*// + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + void processKfMcXic0(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processKfMcXic0WithFT0C(MyEventTableWithFT0C const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processKfMcXic0WithFT0M(MyEventTableWithFT0M const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processKfMcXic0WithNTracksPV(MyEventTableWithNTracksPV const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate table + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfMcXic0, "Process MC with information for xic0", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfMcXic0WithFT0C, "Process MC with information for xic0 at FT0C", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfMcXic0WithFT0M, "Process MC with information for xic0 at FT0M", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfMcXic0WithNTracksPV, "Process MC with information for xic0 at Ntrack", false); +}; // end of struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}