diff --git a/PWGHF/Core/SelectorCuts.h b/PWGHF/Core/SelectorCuts.h index b655a243890..92f62d72cca 100644 --- a/PWGHF/Core/SelectorCuts.h +++ b/PWGHF/Core/SelectorCuts.h @@ -700,6 +700,55 @@ static const std::vector labelsPt = { static const std::vector labelsCutVar = {"m", "pT p", "pT K", "pT Pi", "chi2PCA", "decay length", "cos pointing angle", "decLengthXY", "normDecLXY", "ct", "impParXY"}; } // namespace hf_cuts_xic_to_p_k_pi +namespace hf_cuts_xic_to_xi_pi_pi +{ +static const int nBinsPt = 10; +static const int nCutVars = 12; +// default values for the pT bin edges (can be used to configure histogram axis) +// offset by 1 from the bin numbers in cuts array +constexpr double binsPt[nBinsPt + 1] = { + 0., + 1., + 2., + 3., + 4., + 5., + 6., + 8., + 12., + 24., + 36.}; +auto vecBinsPt = std::vector{binsPt, binsPt + nBinsPt + 1}; + +// default values for the cuts m ptXi ptPi0 ptPi1 chi2PCA dL dLXY cosp cospXY impParXY Xi Pi0 Pi1 +constexpr double cuts[nBinsPt][nCutVars] = {{0.4, 0.4, 0.4, 0.4, 1e-5, 0.5, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1}, /* 0 < pT < 1 */ + {0.4, 0.4, 0.4, 0.4, 1e-5, 0.5, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1}, /* 1 < pT < 2 */ + {0.4, 0.4, 0.4, 0.4, 1e-5, 0.5, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1}, /* 2 < pT < 3 */ + {0.4, 0.4, 0.4, 0.4, 1e-5, 0.5, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1}, /* 3 < pT < 4 */ + {0.4, 0.4, 0.4, 0.4, 1e-5, 0.5, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1}, /* 4 < pT < 5 */ + {0.4, 0.4, 0.4, 0.4, 1e-5, 0.5, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1}, /* 5 < pT < 6 */ + {0.4, 0.4, 0.4, 0.4, 1e-5, 0.5, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1}, /* 6 < pT < 8 */ + {0.4, 0.4, 0.4, 0.4, 1e-5, 0.5, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1}, /* 8 < pT < 10 */ + {0.4, 0.4, 0.4, 0.4, 1e-5, 0.5, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1}, /* 12 < pT < 24 */ + {0.4, 0.4, 0.4, 0.4, 1e-5, 0.5, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1}}; /* 24 < pT < 36 */ + +// row labels +static const std::vector labelsPt = { + "pT bin 0", + "pT bin 1", + "pT bin 2", + "pT bin 3", + "pT bin 4", + "pT bin 5", + "pT bin 6", + "pT bin 7", + "pT bin 8", + "pT bin 9"}; + +// column labels +static const std::vector labelsCutVar = {"m", "pT Xi", "pT Pi0", "pT Pi1", "chi2PCA", "max decay length", "max decay length XY", "cos pointing angle", "cos pointing angle XY", "max impParXY Xi", "max impParXY Pi0", "max impParXY Pi1"}; +} // namespace hf_cuts_xic_to_xi_pi_pi + namespace hf_cuts_xicc_to_p_k_pi_pi { static const int nBinsPt = 10; diff --git a/PWGHF/D2H/Tasks/CMakeLists.txt b/PWGHF/D2H/Tasks/CMakeLists.txt index 4a533812364..48c93bbebca 100644 --- a/PWGHF/D2H/Tasks/CMakeLists.txt +++ b/PWGHF/D2H/Tasks/CMakeLists.txt @@ -89,6 +89,11 @@ o2physics_add_dpl_workflow(task-xic PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(task-xic-to-xi-pi-pi + SOURCES taskXicToXiPiPi.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(task-xicc SOURCES taskXicc.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGHF/D2H/Tasks/taskXicToXiPiPi.cxx b/PWGHF/D2H/Tasks/taskXicToXiPiPi.cxx new file mode 100644 index 00000000000..3a5acb6dd6e --- /dev/null +++ b/PWGHF/D2H/Tasks/taskXicToXiPiPi.cxx @@ -0,0 +1,353 @@ +// 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 taskXicToXiPiPi.cxx +/// \brief Ξc± → (Ξ∓ → (Λ → p π∓) π∓) π± π± analysis task +/// \note adapted from taskBs.cxx +/// +/// \author Phil Stahlhut + +#include "CommonConstants/PhysicsConstants.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include "Framework/runDataProcessing.h" + +#include "PWGHF/Core/SelectorCuts.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" + +using namespace o2; +using namespace o2::aod; +using namespace o2::analysis; +using namespace o2::framework; +using namespace o2::framework::expressions; + +/// Xic analysis task +struct HfTaskXicToXiPiPi { + Configurable selectionFlagXic{"selectionFlagXic", 1, "Selection Flag for Xic"}; + Configurable yCandGenMax{"yCandGenMax", 0.5, "max. gen particle rapidity"}; + Configurable yCandRecoMax{"yCandRecoMax", 0.8, "max. cand. rapidity"}; + Configurable etaTrackMax{"etaTrackMax", 0.8, "max. track pseudo-rapidity"}; + Configurable ptTrackMin{"ptTrackMin", 0.1, "min. track transverse momentum"}; + Configurable> binsPt{"binsPt", std::vector{hf_cuts_xic_to_xi_pi_pi::vecBinsPt}, "pT bin limits"}; + // MC checks + Configurable checkDecayTypeMc{"checkDecayTypeMc", false, "Flag to enable DecayType histogram"}; + + Service pdg; + + Filter filterSelectCandidates = (aod::hf_sel_candidate_xic::isSelXicToXiPiPi >= selectionFlagXic); + + HistogramRegistry registry{ + "registry", + {{"hPtCand", "#Xi^{#plus}_{c} candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{400, 0., 40.}}}}, + {"hPtProng0", "#Xi^{#plus}_{c} candidates;prong 0 (#Xi^{#minus}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}}, + {"hPtProng1", "#Xi^{#plus}_{c} candidates;prong 1 (#pi^{#plus}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 20.}}}}, + {"hPtProng2", "#Xi^{#plus}_{c} candidates;prong 2 (#pi^{#plus}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 20.}}}}}}; + + void init(InitContext const&) + { + static const AxisSpec axisMassXic = {300, 1.8, 3.0, "inv. mass (GeV/#it{c}^{2})"}; + static const AxisSpec axisMassXiRes = {300, 1.0, 2.0, "inv. mass (GeV/#it{c}^{2})"}; + static const AxisSpec axisPt = {(std::vector)binsPt, "#it{p}_{T} (GeV/#it{c})"}; + + registry.add("hEta", "#Xi^{#plus}_{c} candidates;#Xi^{#plus}_{c} candidate #it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, axisPt}}); + registry.add("hRapidity", "#Xi^{#plus}_{c} candidates;#Xi^{#plus}_{c} candidate #it{y};entries", {HistType::kTH2F, {{100, -2., 2.}, axisPt}}); + registry.add("hCPA", "#Xi^{#plus}_{c} candidates;#Xi^{#plus}_{c} candidate cosine of pointing angle;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hCPAxy", "#Xi^{#plus}_{c} candidates;#Xi^{#plus}_{c} candidate cosine of pointing angle xy;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hMass", "#Xi^{#plus}_{c} candidates;inv. mass #Xi^{#mp} #pi^{#pm} #pi^{#pm} (GeV/#it{c}^{2});#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {axisMassXic, axisPt}}); + registry.add("hDecLength", "#Xi^{#plus}_{c} candidates;decay length (cm);entries", {HistType::kTH2F, {{200, 0., 0.4}, axisPt}}); + registry.add("hErrDecLength", "#Xi^{#plus}_{c} candidates;#Xi^{#plus}_{c} candidate decay length error (cm);entries", {HistType::kTH2F, {{100, 0., 1.}, axisPt}}); + registry.add("hDecLengthXY", "#Xi^{#plus}_{c} candidates;decay length xy (cm);entries", {HistType::kTH2F, {{200, 0., 0.4}, axisPt}}); + registry.add("hErrDecLengthXY", "#Xi^{#plus}_{c} candidates;#Xi^{#plus}_{c} candidate decay length xy error (cm);entries", {HistType::kTH2F, {{100, 0., 1.}, axisPt}}); + registry.add("hd0Prong0", "#Xi^{#plus}_{c} candidates;prong 0 (#Xi^{#mp}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{100, -0.05, 0.05}, axisPt}}); + registry.add("hd0Prong1", "#Xi^{#plus}_{c} candidates;prong 1 (#pi^{#pm}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{100, -0.05, 0.05}, axisPt}}); + registry.add("hd0Prong2", "#Xi^{#plus}_{c} candidates;prong 2 (#pi^{#pm}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{100, -0.05, 0.05}, axisPt}}); + registry.add("hImpParErr", "#Xi^{#plus}_{c} candidates;prongs impact parameter error (cm);entries", {HistType::kTH2F, {{200, -0.1, 0.1}, axisPt}}); + registry.add("hChi2PCA", "#Xi^{#plus}_{c} candidates (matched);sum of distances of the secondary vertex to its prongs;entries", {HistType::kTH2F, {{240, -0.01, 0.5}, axisPt}}); + registry.add("hCPAXi", "#Xi^{#plus}_{c} candidates;#Xi^{#minus} candidate cosine of pointing angle;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hCPAxyXi", "#Xi^{#plus}_{c} candidates;#Xi^{#minus} candidate cosine of pointing angle xy;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hCPALambda", "#Xi^{#plus}_{c} candidates;#Lambda candidate cosine of pointing angle;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hCPAxyLambda", "#Xi^{#plus}_{c} candidates;#Lambda candidate cosine of pointing angle xy;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hMassXiPi1", "#Xi^{#plus}_{c} candidates;inv. mass #Xi^{#mp} #pi^{#pm} (prong 1) (GeV/#it{c}^{2});#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {axisMassXiRes, axisPt}}); + registry.add("hMassXiPi2", "#Xi^{#plus}_{c} candidates;inv. mass #Xi^{#mp} #pi^{#pm} (prong 2) (GeV/#it{c}^{2});#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {axisMassXiRes, axisPt}}); + + if (doprocessMc) { + // MC reconstructed + registry.add("hPtGenSig", "#Xi^{#plus}_{c} candidates (gen+rec);candidate #it{p}_{T}^{gen.} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}); + registry.add("hPtRecSig", "#Xi^{#plus}_{c} candidates (matched);candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}); + registry.add("hPtRecBg", "#Xi^{#plus}_{c} candidates (unmatched);candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}); + registry.add("hPtProng0RecSig", "#Xi^{#plus}_{c} candidates (matched);prong 0 (#Xi^{#mp}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}); + registry.add("hPtProng0RecBg", "#Xi^{#plus}_{c} candidates (unmatched);prong 0 (#Xi^{#mp}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}); + registry.add("hPtProng1RecSig", "#Xi^{#plus}_{c} candidates (matched);prong 1 (#pi^{#pm}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 20.}}}); + registry.add("hPtProng1RecBg", "#Xi^{#plus}_{c} candidates (unmatched);prong 1 (#pi^{#pm}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 20.}}}); + registry.add("hPtProng2RecSig", "#Xi^{#plus}_{c} candidates (matched);prong 2 (#pi^{#pm}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 20.}}}); + registry.add("hPtProng2RecBg", "#Xi^{#plus}_{c} candidates (unmatched);prong 2 (#pi^{#pm}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 20.}}}); + registry.add("hEtaRecSig", "#Xi^{#plus}_{c} candidates (matched);#Xi^{#plus}_{c} candidate #it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, axisPt}}); + registry.add("hEtaRecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Xi^{#plus}_{c} candidate #it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, axisPt}}); + registry.add("hRapidityRecSig", "#Xi^{#plus}_{c} candidates (matched);#Xi^{#plus}_{c} candidate #it{y};entries", {HistType::kTH2F, {{100, -2., 2.}, axisPt}}); + registry.add("hRapidityRecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Xi^{#plus}_{c} candidate #it{y};entries", {HistType::kTH2F, {{100, -2., 2.}, axisPt}}); + registry.add("hCPARecSig", "#Xi^{#plus}_{c} candidates (matched);#Xi^{#plus}_{c} candidate cosine of pointing angle;entries", {HistType::kTH2F, {{220, -1.1, 1.1}, axisPt}}); + registry.add("hCPARecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Xi^{#plus}_{c} candidate cosine of pointing angle;entries", {HistType::kTH2F, {{220, -1.1, 1.1}, axisPt}}); + registry.add("hCPAxyRecSig", "#Xi^{#plus}_{c} candidates (matched);#Xi^{#plus}_{c} candidate CPAxy;entries", {HistType::kTH2F, {{220, -1.1, 1.1}, axisPt}}); + registry.add("hCPAxyRecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Xi^{#plus}_{c} candidate CPAxy;entries", {HistType::kTH2F, {{220, -1.1, 1.1}, axisPt}}); + registry.add("hMassRecSig", "#Xi^{#plus}_{c} candidates (matched);inv. mass #Xi^{#mp} #pi^{#pm} #pi^{#pm} (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{300, 1.8, 3.0}, axisPt}}); + registry.add("hMassRecBg", "#Xi^{#plus}_{c} candidates (unmatched);inv. mass #Xi^{#mp} #pi^{#pm} #pi^{#pm} (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{300, 1.8, 3.0}, axisPt}}); + registry.add("hDecLengthRecSig", "#Xi^{#plus}_{c} candidates (matched);#Xi^{#plus}_{c} candidate decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, axisPt}}); + registry.add("hDecLengthRecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Xi^{#plus}_{c} candidate decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, axisPt}}); + registry.add("hErrDecLengthRecSig", "#Xi^{#plus}_{c} candidates (matched);#Xi^{#plus}_{c} candidate decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, axisPt}}); + registry.add("hErrDecLengthRecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Xi^{#plus}_{c} candidate decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, axisPt}}); + registry.add("hDecLengthXYRecSig", "#Xi^{#plus}_{c} candidates (matched);#Xi^{#plus}_{c} candidate decay length xy (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, axisPt}}); + registry.add("hDecLengthXYRecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Xi^{#plus}_{c} candidate decay length xy(cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, axisPt}}); + registry.add("hErrDecLengthXYRecSig", "#Xi^{#plus}_{c} candidates (matched);#Xi^{#plus}_{c} candidate decay length xy (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, axisPt}}); + registry.add("hErrDecLengthXYRecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Xi^{#plus}_{c} candidate decay length xy(cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, axisPt}}); + registry.add("hd0Prong0RecSig", "#Xi^{#plus}_{c} candidates (matched);prong 0 (#Xi^{#mp}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{200, -0.05, 0.05}, axisPt}}); + registry.add("hd0Prong0RecBg", "#Xi^{#plus}_{c} candidates (unmatched);prong 0 (#Xi^{#mp}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{200, -0.05, 0.05}, axisPt}}); + registry.add("hd0Prong1RecSig", "#Xi^{#plus}_{c} candidates (matched);prong 1 (#pi^{#pm}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{200, -0.05, 0.05}, axisPt}}); + registry.add("hd0Prong1RecBg", "#Xi^{#plus}_{c} candidates (unmatched);prong 1 (#pi^{#pm}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{200, -0.05, 0.05}, axisPt}}); + registry.add("hd0Prong2RecSig", "#Xi^{#plus}_{c} candidates (matched);prong 2 (#pi^{#pm}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{200, -0.05, 0.05}, axisPt}}); + registry.add("hd0Prong2RecBg", "#Xi^{#plus}_{c} candidates (unmatched);prong 2 (#pi^{#pm}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{200, -0.05, 0.05}, axisPt}}); + registry.add("hImpParErrRecSig", "#Xi^{#plus}_{c} candidates (matched);prongs impact parameter error (cm);entries", {HistType::kTH2F, {{100, -0.1, 0.1}, axisPt}}); + registry.add("hImpParErrRecBg", "#Xi^{#plus}_{c} candidates (unmatched);prongs impact parameter error (cm);entries", {HistType::kTH2F, {{100, -0.1, 0.1}, axisPt}}); + registry.add("hChi2PCARecSig", "#Xi^{#plus}_{c} candidates (matched);sum of distances of the secondary vertex to its prongs;entries", {HistType::kTH2F, {{240, -0.01, 0.1}, axisPt}}); + registry.add("hChi2PCARecBg", "#Xi^{#plus}_{c} candidates (unmatched);sum of distances of the secondary vertex to its prongs;entries", {HistType::kTH2F, {{240, -0.01, 0.1}, axisPt}}); + registry.add("hCPAXiRecSig", "#Xi^{#plus}_{c} candidates (matched);#Xi^{#minus} cosine of pointing angle;entries", {HistType::kTH2F, {{220, -1.1, 1.1}, axisPt}}); + registry.add("hCPAXiRecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Xi^{#minus} cosine of pointing angle;entries", {HistType::kTH2F, {{220, -1.1, 1.1}, axisPt}}); + registry.add("hCPAxyXiRecSig", "#Xi^{#plus}_{c} candidates (matched);#Xi^{#minus} candidate cosine of pointing angle xy;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hCPAxyXiRecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Xi^{#minus} candidate cosine of pointing angle xy;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hCPALambdaRecSig", "#Xi^{#plus}_{c} candidates (matched);#Lambda candidate cosine of pointing angle;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hCPALambdaRecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Lambda candidate cosine of pointing angle;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hCPAxyLambdaRecSig", "#Xi^{#plus}_{c} candidates (matched);#Lambda candidate cosine of pointing angle xy;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hCPAxyLambdaRecBg", "#Xi^{#plus}_{c} candidates (unmatched);#Lambda candidate cosine of pointing angle xy;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, axisPt}}); + registry.add("hMassXiPi1RecSig", "#Xi^{#plus}_{c} candidates (matched);inv. mass #Xi^{#mp} #pi^{#pm} (prong 1) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{300, 1.0, 2.0}, axisPt}}); + registry.add("hMassXiPi1RecBg", "#Xi^{#plus}_{c} candidates (unmatched);inv. mass #Xi^{#mp} #pi^{#pm} (prong 1) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{300, 1.0, 2.0}, axisPt}}); + registry.add("hMassXiPi2RecSig", "#Xi^{#plus}_{c} candidates (matched);inv. mass #Xi^{#mp} #pi^{#pm} (prong 2) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{300, 1.0, 2.0}, axisPt}}); + registry.add("hMassXiPi2RecBg", "#Xi^{#plus}_{c} candidates (unmatched);inv. mass #Xi^{#mp} #pi^{#pm} (prong 2) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{300, 1.0, 2.0}, axisPt}}); + // MC generated + registry.add("hPtProng0Gen", "MC particles (generated);prong 0 (#Xi^{#mp}) #it{p}_{T}^{gen} (GeV/#it{c});entries", {HistType::kTH2F, {{300, 0., 30.}, axisPt}}); + registry.add("hPtProng1Gen", "MC particles (generated);prong 1 (#pi^{#pm}) #it{p}_{T}^{gen} (GeV/#it{c});entries", {HistType::kTH2F, {{200, 0., 20.}, axisPt}}); + registry.add("hPtProng2Gen", "MC particles (generated);prong 2 (#pi^{#pm}) #it{p}_{T}^{gen} (GeV/#it{c});entries", {HistType::kTH2F, {{200, 0., 20.}, axisPt}}); + registry.add("hEtaProng0Gen", "MC particles (generated);prong 0 (#Xi^{#mp}) #it{#eta}^{gen};entries", {HistType::kTH2F, {{100, -2, 2}, axisPt}}); + registry.add("hEtaProng1Gen", "MC particles (generated);prong 1 (#pi^{#pm}) #it{#eta}^{gen};entries", {HistType::kTH2F, {{100, -2, 2}, axisPt}}); + registry.add("hEtaProng2Gen", "MC particles (generated);prong 2 (#pi^{#pm}) #it{#eta}^{gen};entries", {HistType::kTH2F, {{100, -2, 2}, axisPt}}); + registry.add("hYProng0Gen", "MC particles (generated);prong 0 (#Xi^{#mp}) #it{y}^{gen};entries", {HistType::kTH2F, {{100, -2, 2}, axisPt}}); + registry.add("hYProng1Gen", "MC particles (generated);prong 1 (#pi^{#pm}) #it{y}^{gen};entries", {HistType::kTH2F, {{100, -2, 2}, axisPt}}); + registry.add("hYProng2Gen", "MC particles (generated);prong 2 (#pi^{#pm}) #it{y}^{gen};entries", {HistType::kTH2F, {{100, -2, 2}, axisPt}}); + registry.add("hPtGen", "MC particles (generated);candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}); + registry.add("hEtaGen", "MC particles (generated);#Xi^{#plus}_{c} candidate #it{#eta}^{gen};entries", {HistType::kTH2F, {{100, -2., 2.}, axisPt}}); + registry.add("hYGen", "MC particles (generated);#Xi^{#plus}_{c} candidate #it{y}^{gen};entries", {HistType::kTH2F, {{100, -2., 2.}, axisPt}}); + registry.add("hPtGenWithProngsInAcceptance", "MC particles (generated-daughters in acceptance);candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}); + registry.add("hEtaGenWithProngsInAcceptance", "MC particles (generated-daughters in acceptance);#Xi^{#plus}_{c} candidate #it{#eta}^{gen};entries", {HistType::kTH2F, {{100, -2., 2.}, axisPt}}); + registry.add("hYGenWithProngsInAcceptance", "MC particles (generated-daughters in acceptance);#Xi^{#plus}_{c} candidate #it{y}^{gen};entries", {HistType::kTH2F, {{100, -2., 2.}, axisPt}}); + } + + if (checkDecayTypeMc) { + constexpr uint8_t kNBinsDecayTypeMc = hf_cand_xic_to_xi_pi_pi::DecayType::NDecayType + 1; + TString labels[kNBinsDecayTypeMc]; + labels[hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiPiPi] = "#Xi^{+}_{c} #rightarrow #Xi^{#minus} #pi^{#plus}) #pi^{#plus}"; + labels[hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiResPiToXiPiPi] = "#Xi^{+}_{c} #rightarrow #Xi(1530)^{0} #pi^{#plus} #rightarrow #Xi^{#minus} #pi^{#plus}) #pi^{#plus}"; + labels[hf_cand_xic_to_xi_pi_pi::DecayType::NDecayType] = "Other decays"; + static const AxisSpec axisDecayType = {kNBinsDecayTypeMc, 0.5, kNBinsDecayTypeMc + 0.5, ""}; + registry.add("hDecayTypeMc", "DecayType", {HistType::kTH3F, {axisDecayType, axisMassXic, axisPt}}); + for (uint8_t iBin = 0; iBin < kNBinsDecayTypeMc; ++iBin) { + registry.get(HIST("hDecayTypeMc"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin]); + } + } + } + + /// Selection of Xic daughter in geometrical acceptance + /// \param etaProng is the pseudorapidity of Xic prong + /// \param ptProng is the pT of Xic prong + /// \return true if prong is in geometrical acceptance + template + bool isProngInAcceptance(const T& etaProng, const T& ptProng) + { + return std::abs(etaProng) <= etaTrackMax && ptProng >= ptTrackMin; + } + + void process(soa::Filtered> const& candidates) + { + for (const auto& candidate : candidates) { + auto yCandXic = candidate.y(o2::constants::physics::MassXiCPlus); + if (yCandRecoMax >= 0. && std::abs(yCandXic) > yCandRecoMax) { + continue; + } + + auto ptCandXic = candidate.pt(); + + registry.fill(HIST("hPtCand"), ptCandXic); + registry.fill(HIST("hPtProng0"), candidate.ptProng0()); + registry.fill(HIST("hPtProng1"), candidate.ptProng1()); + registry.fill(HIST("hPtProng2"), candidate.ptProng2()); + registry.fill(HIST("hEta"), candidate.eta(), ptCandXic); + registry.fill(HIST("hRapidity"), yCandXic, ptCandXic); + registry.fill(HIST("hCPA"), candidate.cpa(), ptCandXic); + registry.fill(HIST("hCPAxy"), candidate.cpaXY(), ptCandXic); + registry.fill(HIST("hMass"), candidate.invMassXic(), ptCandXic); + registry.fill(HIST("hDecLength"), candidate.decayLength(), ptCandXic); + registry.fill(HIST("hErrDecLength"), candidate.errorDecayLength(), ptCandXic); + registry.fill(HIST("hDecLengthXY"), candidate.decayLengthXY(), ptCandXic); + registry.fill(HIST("hErrDecLengthXY"), candidate.errorDecayLengthXY(), ptCandXic); + registry.fill(HIST("hd0Prong0"), candidate.impactParameter0(), ptCandXic); + registry.fill(HIST("hd0Prong1"), candidate.impactParameter1(), ptCandXic); + registry.fill(HIST("hd0Prong2"), candidate.impactParameter2(), ptCandXic); + registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter0(), ptCandXic); + registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter1(), ptCandXic); + registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter2(), ptCandXic); + registry.fill(HIST("hChi2PCA"), candidate.chi2PCA(), ptCandXic); + registry.fill(HIST("hCPAXi"), candidate.cosPaXi(), ptCandXic); + registry.fill(HIST("hCPAxyXi"), candidate.cosPaXYXi(), ptCandXic); + registry.fill(HIST("hCPALambda"), candidate.cosPaLambda(), ptCandXic); + registry.fill(HIST("hCPAxyLambda"), candidate.cosPaLambda(), ptCandXic); + registry.fill(HIST("hMassXiPi1"), candidate.invMassXiPi0(), ptCandXic); + registry.fill(HIST("hMassXiPi2"), candidate.invMassXiPi1(), ptCandXic); + } // candidate loop + } // process + + /// MC analysis and fill histograms + void processMc(soa::Filtered> const& candidates, + soa::Join const& mcParticles, + aod::TracksWMc const&) + { + std::vector arrDaughIndex; + + // MC rec + for (const auto& candidate : candidates) { + auto yCandXic = candidate.y(o2::constants::physics::MassXiCPlus); + if (yCandRecoMax >= 0. && std::abs(yCandXic) > yCandRecoMax) { + continue; + } + + auto ptCandXic = candidate.pt(); + int flagMcMatchRecXic = std::abs(candidate.flagMcMatchRec()); + + if (TESTBIT(flagMcMatchRecXic, hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiPiPi) || TESTBIT(flagMcMatchRecXic, hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiResPiToXiPiPi)) { + auto indexMother = RecoDecay::getMother(mcParticles, candidate.pi0_as().mcParticle_as>(), o2::constants::physics::Pdg::kXiCPlus, true); + auto particleMother = mcParticles.rawIteratorAt(indexMother); + + registry.fill(HIST("hPtGenSig"), particleMother.pt()); + registry.fill(HIST("hPtRecSig"), ptCandXic); + registry.fill(HIST("hPtProng0RecSig"), candidate.ptProng0()); + registry.fill(HIST("hPtProng1RecSig"), candidate.ptProng1()); + registry.fill(HIST("hPtProng2RecSig"), candidate.ptProng2()); + registry.fill(HIST("hEtaRecSig"), candidate.eta(), ptCandXic); + registry.fill(HIST("hRapidityRecSig"), yCandXic, ptCandXic); + registry.fill(HIST("hCPARecSig"), candidate.cpa(), ptCandXic); + registry.fill(HIST("hCPAxyRecSig"), candidate.cpaXY(), ptCandXic); + registry.fill(HIST("hMassRecSig"), candidate.invMassXic(), ptCandXic); + registry.fill(HIST("hDecLengthRecSig"), candidate.decayLength(), ptCandXic); + registry.fill(HIST("hErrDecLengthRecSig"), candidate.errorDecayLength(), ptCandXic); + registry.fill(HIST("hDecLengthXYRecSig"), candidate.decayLengthXY(), ptCandXic); + registry.fill(HIST("hErrDecLengthXYRecSig"), candidate.errorDecayLengthXY(), ptCandXic); + registry.fill(HIST("hd0Prong0RecSig"), candidate.impactParameter0(), ptCandXic); + registry.fill(HIST("hd0Prong1RecSig"), candidate.impactParameter1(), ptCandXic); + registry.fill(HIST("hd0Prong2RecSig"), candidate.impactParameter2(), ptCandXic); + registry.fill(HIST("hImpParErrRecSig"), candidate.errorImpactParameter0(), ptCandXic); + registry.fill(HIST("hImpParErrRecSig"), candidate.errorImpactParameter1(), ptCandXic); + registry.fill(HIST("hImpParErrRecSig"), candidate.errorImpactParameter2(), ptCandXic); + registry.fill(HIST("hChi2PCARecSig"), candidate.chi2PCA(), ptCandXic); + registry.fill(HIST("hCPAXiRecSig"), candidate.cosPaXi(), ptCandXic); + registry.fill(HIST("hCPAxyXiRecSig"), candidate.cosPaXYXi(), ptCandXic); + registry.fill(HIST("hCPALambdaRecSig"), candidate.cosPaLambda(), ptCandXic); + registry.fill(HIST("hCPAxyLambdaRecSig"), candidate.cosPaLambda(), ptCandXic); + } else { + registry.fill(HIST("hPtRecBg"), ptCandXic); + registry.fill(HIST("hPtProng0RecBg"), candidate.ptProng0()); + registry.fill(HIST("hPtProng1RecBg"), candidate.ptProng1()); + registry.fill(HIST("hPtProng2RecBg"), candidate.ptProng2()); + registry.fill(HIST("hEtaRecBg"), candidate.eta(), ptCandXic); + registry.fill(HIST("hRapidityRecBg"), yCandXic, ptCandXic); + registry.fill(HIST("hCPARecBg"), candidate.cpa(), ptCandXic); + registry.fill(HIST("hCPAxyRecBg"), candidate.cpaXY(), ptCandXic); + registry.fill(HIST("hMassRecBg"), candidate.invMassXic(), ptCandXic); + registry.fill(HIST("hDecLengthRecBg"), candidate.decayLength(), ptCandXic); + registry.fill(HIST("hErrDecLengthRecBg"), candidate.errorDecayLength(), ptCandXic); + registry.fill(HIST("hDecLengthXYRecBg"), candidate.decayLengthXY(), ptCandXic); + registry.fill(HIST("hErrDecLengthXYRecBg"), candidate.errorDecayLengthXY(), ptCandXic); + registry.fill(HIST("hd0Prong0RecBg"), candidate.impactParameter0(), ptCandXic); + registry.fill(HIST("hd0Prong1RecBg"), candidate.impactParameter1(), ptCandXic); + registry.fill(HIST("hd0Prong2RecBg"), candidate.impactParameter2(), ptCandXic); + registry.fill(HIST("hImpParErrRecBg"), candidate.errorImpactParameter0(), ptCandXic); + registry.fill(HIST("hImpParErrRecBg"), candidate.errorImpactParameter1(), ptCandXic); + registry.fill(HIST("hImpParErrRecBg"), candidate.errorImpactParameter2(), ptCandXic); + registry.fill(HIST("hChi2PCARecBg"), candidate.chi2PCA(), ptCandXic); + registry.fill(HIST("hCPAXiRecBg"), candidate.cosPaXi(), ptCandXic); + registry.fill(HIST("hCPAxyXiRecBg"), candidate.cosPaXYXi(), ptCandXic); + registry.fill(HIST("hCPALambdaRecBg"), candidate.cosPaLambda(), ptCandXic); + registry.fill(HIST("hCPAxyLambdaRecBg"), candidate.cosPaLambda(), ptCandXic); + } + + if (checkDecayTypeMc) { + if (TESTBIT(flagMcMatchRecXic, hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiPiPi)) { + registry.fill(HIST("hDecayTypeMc"), 1 + hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiPiPi, candidate.invMassXic(), ptCandXic); + } else if (TESTBIT(flagMcMatchRecXic, hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiResPiToXiPiPi)) { + registry.fill(HIST("hDecayTypeMc"), 1 + hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiResPiToXiPiPi, candidate.invMassXic(), ptCandXic); + } else { + registry.fill(HIST("hDecayTypeMc"), 1 + hf_cand_xic_to_xi_pi_pi::DecayType::NDecayType, candidate.invMassXic(), ptCandXic); + } + } + } // rec + + // MC gen. level + for (const auto& particle : mcParticles) { + if (TESTBIT(std::abs(particle.flagMcMatchGen()), hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiPiPi) || TESTBIT(std::abs(particle.flagMcMatchGen()), hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiResPiToXiPiPi)) { + arrDaughIndex.clear(); + + auto ptParticle = particle.pt(); + auto yParticle = RecoDecay::y(particle.pVector(), o2::constants::physics::MassXiCPlus); + if (yCandGenMax >= 0. && std::abs(yParticle) > yCandGenMax) { + continue; + } + + // get kinematic variables of Ξ π π + std::array ptProngs; + std::array yProngs; + std::array etaProngs; + int counter = 0; + RecoDecay::getDaughters(particle, &arrDaughIndex, std::array{+kXiMinus, +kPiPlus, +kPiPlus}, 2); + for (auto iProng = 0u; iProng < arrDaughIndex.size(); ++iProng) { + auto daughI = mcParticles.rawIteratorAt(arrDaughIndex[iProng]); + ptProngs[counter] = daughI.pt(); + etaProngs[counter] = daughI.eta(); + yProngs[counter] = RecoDecay::y(daughI.pVector(), pdg->Mass(daughI.pdgCode())); + counter++; + } + + registry.fill(HIST("hPtProng0Gen"), ptProngs[0], ptParticle); + registry.fill(HIST("hPtProng1Gen"), ptProngs[1], ptParticle); + registry.fill(HIST("hPtProng2Gen"), ptProngs[2], ptParticle); + registry.fill(HIST("hEtaProng0Gen"), etaProngs[0], ptParticle); + registry.fill(HIST("hEtaProng1Gen"), etaProngs[1], ptParticle); + registry.fill(HIST("hEtaProng2Gen"), etaProngs[2], ptParticle); + registry.fill(HIST("hYProng0Gen"), yProngs[0], ptParticle); + registry.fill(HIST("hYProng1Gen"), yProngs[1], ptParticle); + registry.fill(HIST("hYProng2Gen"), yProngs[2], ptParticle); + registry.fill(HIST("hPtGen"), ptParticle); + registry.fill(HIST("hYGen"), yParticle, ptParticle); + registry.fill(HIST("hEtaGen"), particle.eta(), ptParticle); + + // reject Xic daughters that are not in geometrical acceptance + if (!isProngInAcceptance(etaProngs[0], ptProngs[0]) || !isProngInAcceptance(etaProngs[1], ptProngs[1]) || !isProngInAcceptance(etaProngs[2], ptProngs[2])) { + continue; + } + registry.fill(HIST("hPtGenWithProngsInAcceptance"), ptParticle); + registry.fill(HIST("hEtaGenWithProngsInAcceptance"), particle.eta(), ptParticle); + registry.fill(HIST("hYGenWithProngsInAcceptance"), yParticle, ptParticle); + } + } // gen + } // process + PROCESS_SWITCH(HfTaskXicToXiPiPi, processMc, "Process MC", false); +}; // struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index 6c1b4d50356..e1b28dfaab4 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -1320,6 +1320,120 @@ DECLARE_SOA_TABLE(HfToOmegaKMCGen, "AOD", "HFTOOMEKMCGEN", //! hf_cand_xic0_omegac0::FlagMcMatchGen, hf_cand_xic0_omegac0::DebugGenCharmBar, hf_cand_xic0_omegac0::DebugGenCasc, hf_cand_xic0_omegac0::DebugGenLambda, hf_cand_xic0_omegac0::PtCharmBaryonGen, hf_cand_xic0_omegac0::RapidityCharmBaryonGen, hf_cand_xic0_omegac0::OriginGen, hf_cand::IdxBhadMotherPart, o2::soa::Marker<4>); +// specific Xic to Xi Pi Pi candidate properties +namespace hf_cand_xic_to_xi_pi_pi +{ +DECLARE_SOA_INDEX_COLUMN_FULL(Pi0, pi0, int, Tracks, "_pi0"); +DECLARE_SOA_INDEX_COLUMN_FULL(Pi1, pi1, int, Tracks, "_pi1"); +DECLARE_SOA_COLUMN(XPvErr, xPvErr, float); +DECLARE_SOA_COLUMN(YPvErr, yPvErr, float); +DECLARE_SOA_COLUMN(ZPvErr, zPvErr, float); +DECLARE_SOA_COLUMN(XSvErr, xSvErr, float); +DECLARE_SOA_COLUMN(YSvErr, ySvErr, float); +DECLARE_SOA_COLUMN(ZSvErr, zSvErr, float); +DECLARE_SOA_COLUMN(XDecayVtxXi, xDecayVtxXi, float); +DECLARE_SOA_COLUMN(YDecayVtxXi, yDecayVtxXi, float); +DECLARE_SOA_COLUMN(ZDecayVtxXi, zDecayVtxXi, float); +DECLARE_SOA_COLUMN(XDecayVtxLambda, xDecayVtxLambda, float); +DECLARE_SOA_COLUMN(YDecayVtxLambda, yDecayVtxLambda, float); +DECLARE_SOA_COLUMN(ZDecayVtxLambda, zDecayVtxLambda, float); +DECLARE_SOA_COLUMN(CosPaXi, cosPaXi, float); +DECLARE_SOA_COLUMN(CosPaXYXi, cosPaXYXi, float); +DECLARE_SOA_COLUMN(CosPaLambda, cosPaLambda, float); +DECLARE_SOA_COLUMN(CosPaXYLambda, cosPaXYLambda, float); +DECLARE_SOA_COLUMN(InvMassXic, invMassXic, float); +DECLARE_SOA_COLUMN(Sign, sign, float); +DECLARE_SOA_COLUMN(InvMassXiPi0, invMassXiPi0, float); +DECLARE_SOA_COLUMN(InvMassXiPi1, invMassXiPi1, float); +// KF specific columns +DECLARE_SOA_COLUMN(DcaPi0Pi1, dcaPi0Pi1, float); +DECLARE_SOA_COLUMN(DcaPi0Xi, dcaPi0Xi, float); +DECLARE_SOA_COLUMN(DcaPi1Xi, dcaPi1Xi, float); +// MC matching result: +DECLARE_SOA_COLUMN(FlagMcMatchRec, flagMcMatchRec, int8_t); // reconstruction level +DECLARE_SOA_COLUMN(FlagMcMatchGen, flagMcMatchGen, int8_t); // generator level +DECLARE_SOA_COLUMN(DebugMcRec, debugMcRec, int8_t); // debug flag for mis-association reconstruction level +DECLARE_SOA_COLUMN(DebugMcGen, debugMcGen, int8_t); +DECLARE_SOA_COLUMN(OriginRec, originRec, int8_t); +DECLARE_SOA_COLUMN(OriginGen, originGen, int8_t); + +// mapping of decay types +enum DecayType { XicToXiPiPi = 0, // Ξc± → Ξ∓ π± π± + XicToXiResPiToXiPiPi, // Ξc± → Ξ(1530) π± → Ξ∓ π± π± + NDecayType }; +} // end of namespace hf_cand_xic_to_xi_pi_pi + +// declare dedicated Xic to Xi Pi Pi candidate table +DECLARE_SOA_TABLE(HfCandXicBase, "AOD", "HFCANDXICBASE", + hf_cand::CollisionId, + collision::PosX, collision::PosY, collision::PosZ, + hf_cand_xic_to_xi_pi_pi::XPvErr, hf_cand_xic_to_xi_pi_pi::YPvErr, hf_cand_xic_to_xi_pi_pi::ZPvErr, + // 3-prong specific columns + cascdata::CascadeId, hf_cand_xic_to_xi_pi_pi::Pi0Id, hf_cand_xic_to_xi_pi_pi::Pi1Id, + cascdata::BachelorId, cascdata::PosTrackId, cascdata::NegTrackId, + hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, + hf_cand_xic_to_xi_pi_pi::XSvErr, hf_cand_xic_to_xi_pi_pi::YSvErr, hf_cand_xic_to_xi_pi_pi::ZSvErr, + hf_cand::ErrorDecayLength, hf_cand::ErrorDecayLengthXY, + hf_cand::Chi2PCA, hf_cand_xic_to_xi_pi_pi::InvMassXic, hf_cand_xic_to_xi_pi_pi::Sign, + hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, + hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, + hf_cand::PxProng2, hf_cand::PyProng2, hf_cand::PzProng2, + hf_cand::ImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ImpactParameter2, + hf_cand::ErrorImpactParameter0, hf_cand::ErrorImpactParameter1, hf_cand::ErrorImpactParameter2, + // cascade specific columns + hf_cand_xic_to_xi_pi_pi::XDecayVtxXi, hf_cand_xic_to_xi_pi_pi::YDecayVtxXi, hf_cand_xic_to_xi_pi_pi::ZDecayVtxXi, + hf_cand_xic_to_xi_pi_pi::XDecayVtxLambda, hf_cand_xic_to_xi_pi_pi::YDecayVtxLambda, hf_cand_xic_to_xi_pi_pi::ZDecayVtxLambda, + hf_cand_xic_to_xi_pi_pi::CosPaXi, hf_cand_xic_to_xi_pi_pi::CosPaXYXi, hf_cand_xic_to_xi_pi_pi::CosPaLambda, hf_cand_xic_to_xi_pi_pi::CosPaXYLambda, + hf_cand_xic_to_xi_pi_pi::InvMassXiPi0, hf_cand_xic_to_xi_pi_pi::InvMassXiPi1, + /* dynamic columns */ + hf_cand::RSecondaryVertex, + hf_cand::DecayLength, + hf_cand::DecayLengthXY, + hf_cand::DecayLengthNormalised, + hf_cand::DecayLengthXYNormalised, + hf_cand::ImpactParameterNormalised0, + hf_cand::ImpactParameterNormalised1, + hf_cand::ImpactParameterNormalised2, + /* dynamic columns that use daughter momentum components */ + hf_cand::PtProng0, + hf_cand::PtProng1, + hf_cand::PtProng2, + /* dynamic columns that use candidate momentum components */ + hf_cand::Pt, + hf_cand::P, + hf_cand::PVector, + hf_cand::CPA, + hf_cand::CPAXY, + hf_cand::Ct, + hf_cand::ImpactParameterXY, + hf_cand_3prong::MaxNormalisedDeltaIP, + hf_cand::Eta, + hf_cand::Phi, + hf_cand::Y, + hf_cand::E); + +// extended table with expression columns that can be used as arguments of dynamic columns +DECLARE_SOA_EXTENDED_TABLE_USER(HfCandXicExt, HfCandXicBase, "HFCANDXICEXT", + hf_cand_3prong::Px, hf_cand_3prong::Py, hf_cand_3prong::Pz); + +using HfCandXic = HfCandXicExt; + +DECLARE_SOA_TABLE(HfCandXicKF, "AOD", "HFCANDXICKF", + cascdata::KFCascadeChi2, cascdata::KFV0Chi2, + hf_cand_xic_to_xi_pi_pi::DcaPi0Pi1, hf_cand_xic_to_xi_pi_pi::DcaPi0Xi, hf_cand_xic_to_xi_pi_pi::DcaPi1Xi); + +// table with results of reconstruction level MC matching +DECLARE_SOA_TABLE(HfCandXicMcRec, "AOD", "HFCANDXICMCREC", //! + hf_cand_xic_to_xi_pi_pi::FlagMcMatchRec, + hf_cand_xic_to_xi_pi_pi::DebugMcRec, + hf_cand_xic_to_xi_pi_pi::OriginRec); + +// table with results of generator level MC matching +DECLARE_SOA_TABLE(HfCandXicMcGen, "AOD", "HFCANDXICMCGEN", //! + hf_cand_xic_to_xi_pi_pi::FlagMcMatchGen, + hf_cand_xic_to_xi_pi_pi::DebugMcGen, + hf_cand_xic_to_xi_pi_pi::OriginGen); + // specific chic candidate properties namespace hf_cand_chic { diff --git a/PWGHF/DataModel/CandidateSelectionTables.h b/PWGHF/DataModel/CandidateSelectionTables.h index fa9bc07a981..e5b38f39493 100644 --- a/PWGHF/DataModel/CandidateSelectionTables.h +++ b/PWGHF/DataModel/CandidateSelectionTables.h @@ -301,16 +301,22 @@ DECLARE_SOA_TABLE(HfSelChicToJpsiGamma, "AOD", "HFSELCHIC", //! namespace hf_sel_candidate_xic { +// XicPlus to P K Pi DECLARE_SOA_COLUMN(IsSelXicToPKPi, isSelXicToPKPi, int); //! DECLARE_SOA_COLUMN(IsSelXicToPiKP, isSelXicToPiKP, int); //! DECLARE_SOA_COLUMN(MlProbXicToPKPi, mlProbXicToPKPi, std::vector); //! DECLARE_SOA_COLUMN(MlProbXicToPiKP, mlProbXicToPiKP, std::vector); //! +// XicPlus to Xi Pi Pi +DECLARE_SOA_COLUMN(IsSelXicToXiPiPi, isSelXicToXiPiPi, int); //! } // namespace hf_sel_candidate_xic DECLARE_SOA_TABLE(HfSelXicToPKPi, "AOD", "HFSELXIC", //! hf_sel_candidate_xic::IsSelXicToPKPi, hf_sel_candidate_xic::IsSelXicToPiKP); DECLARE_SOA_TABLE(HfMlXicToPKPi, "AOD", "HFMLXIC", //! hf_sel_candidate_xic::MlProbXicToPKPi, hf_sel_candidate_xic::MlProbXicToPiKP); +// XicPlus to Xi Pi Pi +DECLARE_SOA_TABLE(HfSelXicToXiPiPi, "AOD", "HFSELXICTOXI2PI", //! + hf_sel_candidate_xic::IsSelXicToXiPiPi); namespace hf_sel_candidate_xicc { diff --git a/PWGHF/TableProducer/CMakeLists.txt b/PWGHF/TableProducer/CMakeLists.txt index cbb9580bbf4..205e5c8c01b 100644 --- a/PWGHF/TableProducer/CMakeLists.txt +++ b/PWGHF/TableProducer/CMakeLists.txt @@ -80,6 +80,11 @@ o2physics_add_dpl_workflow(candidate-creator-xic0-omegac0 PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter 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 + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(candidate-creator-xicc SOURCES candidateCreatorXicc.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter @@ -162,6 +167,11 @@ o2physics_add_dpl_workflow(candidate-selector-xic-to-p-k-pi PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(candidate-selector-xic-to-xi-pi-pi + SOURCES candidateSelectorXicToXiPiPi.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(candidate-selector-xicc-to-p-k-pi-pi SOURCES candidateSelectorXiccToPKPiPi.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore @@ -239,6 +249,11 @@ o2physics_add_dpl_workflow(tree-creator-xic-to-p-k-pi PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(tree-creator-xic-to-xi-pi-pi + SOURCES treeCreatorXicToXiPiPi.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(tree-creator-xicc-to-p-k-pi-pi SOURCES treeCreatorXiccToPKPiPi.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGHF/TableProducer/candidateCreatorXicToXiPiPi.cxx b/PWGHF/TableProducer/candidateCreatorXicToXiPiPi.cxx new file mode 100644 index 00000000000..36f9c083625 --- /dev/null +++ b/PWGHF/TableProducer/candidateCreatorXicToXiPiPi.cxx @@ -0,0 +1,656 @@ +// 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 candidateCreatorXicToXiPiPi.cxx +/// \brief Reconstruction of Ξc± → (Ξ∓ → (Λ → p π∓) π∓) π± π± candidates +/// +/// \author Phil Lennart Stahlhut , Heidelberg University +/// \author Jinjoo Seo , Heidelberg University + +#ifndef HomogeneousField +#define HomogeneousField +#endif + +#include +#include +#include +#include +#include + +#include + +#include "CommonConstants/PhysicsConstants.h" +#include "DCAFitter/DCAFitterN.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/DCA.h" + +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/CollisionAssociationTables.h" +#include "Tools/KFparticle/KFUtilities.h" + +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/Utils/utilsBfieldCCDB.h" + +using namespace o2; +using namespace o2::analysis; +using namespace o2::aod::hf_cand_xic_to_xi_pi_pi; +using namespace o2::constants::physics; +using namespace o2::framework; + +/// Reconstruction of heavy-flavour 3-prong decay candidates +struct HfCandidateCreatorXicToXiPiPi { + Produces rowCandidateBase; + Produces rowCandidateKF; + + Configurable fillHistograms{"fillHistograms", true, "do validation plots"}; + // magnetic field setting from CCDB + Configurable isRun2{"isRun2", false, "enable Run 2 or Run 3 GRP objects for magnetic field"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable ccdbPathLut{"ccdbPathLut", "GLO/Param/MatLUT", "Path for LUT parametrization"}; + Configurable ccdbPathGrp{"ccdbPathGrp", "GLO/GRP/GRP", "Path of the grp file (Run 2)"}; + Configurable ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object (Run 3)"}; + // cascade preselections + Configurable doCascadePreselection{"doCascadePreselection", true, "Use invariant mass and dcaXY cuts to preselect cascade candidates"}; + Configurable massToleranceCascade{"massToleranceCascade", 0.01, "Invariant mass tolerance for cascade"}; + Configurable dcaXYToPVCascadeMax{"dcaXYToPVCascadeMax", 3, "Max cascade DCA to PV in xy plane"}; + // DCA fitter + 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 threshold"}; + Configurable minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"}; + Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; + Configurable useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"}; + Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; + // KFParticle + Configurable constrainXicPlusToPv{"constrainXicPlusToPv", false, "Constrain XicPlus to PV"}; + Configurable kfConstructMethod{"kfConstructMethod", 0, "Construct method of XicPlus: 0 no mass constraint, 2 mass constraint"}; + Configurable rejDiffCollTrack{"rejDiffCollTrack", true, "Reject tracks coming from different collisions (effective only for KFParticle w/o derived data)"}; + + Service ccdb; + o2::base::MatLayerCylSet* lut; + o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; + + o2::vertexing::DCAFitterN<3> df; + + int runNumber{0}; + float massXiPiPi{0.}; + float massXiPi0{0.}; + float massXiPi1{0.}; + double bz{0.}; + + using CascadesLinked = soa::Join; + using CascFull = soa::Join; + using KFCascadesLinked = soa::Join; + using KFCascFull = soa::Join; + + HistogramRegistry registry{"registry"}; + + void init(InitContext const&) + { + if ((doprocessXicplusWithDcaFitter + doprocessXicplusWithKFParticle) != 1) { + LOGP(fatal, "Only one process function can be enabled at a time."); + } + + // add histograms to registry + if (fillHistograms) { + registry.add("hMass3", "3-prong candidates;inv. mass (#Xi #pi #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{500, 2.3, 2.7}}}); + registry.add("hCovPVXX", "3-prong candidates;XX element of cov. matrix of prim. vtx. position (cm^{2});entries", {HistType::kTH1D, {{100, 0., 1.e-4}}}); + registry.add("hCovSVXX", "3-prong candidates;XX element of cov. matrix of sec. vtx. position (cm^{2});entries", {HistType::kTH1D, {{100, 0., 0.2}}}); + registry.add("hCovPVYY", "3-prong candidates;YY element of cov. matrix of prim. vtx. position (cm^{2});entries", {HistType::kTH1D, {{100, 0., 1.e-4}}}); + registry.add("hCovSVYY", "3-prong candidates;YY element of cov. matrix of sec. vtx. position (cm^{2});entries", {HistType::kTH1D, {{100, 0., 0.2}}}); + registry.add("hCovPVXZ", "3-prong candidates;XZ element of cov. matrix of prim. vtx. position (cm^{2});entries", {HistType::kTH1D, {{100, -1.e-4, 1.e-4}}}); + registry.add("hCovSVXZ", "3-prong candidates;XZ element of cov. matrix of sec. vtx. position (cm^{2});entries", {HistType::kTH1D, {{100, -1.e-4, 0.2}}}); + registry.add("hCovPVZZ", "3-prong candidates;ZZ element of cov. matrix of prim. vtx. position (cm^{2});entries", {HistType::kTH1D, {{100, 0., 1.e-4}}}); + registry.add("hCovSVZZ", "3-prong candidates;ZZ element of cov. matrix of sec. vtx. position (cm^{2});entries", {HistType::kTH1D, {{100, 0., 0.2}}}); + registry.add("hVertexerType", "Use KF or DCAFitterN;Vertexer type;entries", {HistType::kTH1F, {{2, -0.5, 1.5}}}); // See o2::aod::hf_cand::VertexerType + registry.add("hDcaXYProngs", "DCAxy of 3-prong candidates;#it{p}_{T} (GeV/#it{c};#it{d}_{xy}) (#mum);entries", {HistType::kTH2D, {{100, 0., 20.}, {200, -500., 500.}}}); + registry.add("hDcaZProngs", "DCAz of 3-prong candidates;#it{p}_{T} (GeV/#it{c};#it{d}_{z}) (#mum);entries", {HistType::kTH2D, {{100, 0., 20.}, {200, -500., 500.}}}); + } + + // fill hVertexerType histogram + if (doprocessXicplusWithDcaFitter && fillHistograms) { + registry.fill(HIST("hVertexerType"), aod::hf_cand::VertexerType::DCAFitter); + } + if (doprocessXicplusWithKFParticle && fillHistograms) { + registry.fill(HIST("hVertexerType"), aod::hf_cand::VertexerType::KfParticle); + } + + // initialize CCDB + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(ccdbPathLut)); + runNumber = 0; + + // initialize 3-prong vertex fitter + df.setPropagateToPCA(propagateToPCA); + df.setMaxR(maxR); + df.setMaxDZIni(maxDZIni); + df.setMinParamChange(minParamChange); + df.setMinRelChi2Change(minRelChi2Change); + df.setUseAbsDCA(useAbsDCA); + df.setWeightedFinalPCA(useWeightedFinalPCA); + } + + void processXicplusWithDcaFitter(aod::Collisions const&, + aod::HfCascLf3Prongs const& rowsTrackIndexXicPlus, + CascadesLinked const&, + CascFull const&, + aod::TracksWCovDca const&, + aod::BCsWithTimestamps const&) + { + // loop over triplets of track indices + for (const auto& rowTrackIndexXicPlus : rowsTrackIndexXicPlus) { + auto cascAodElement = rowTrackIndexXicPlus.cascade_as(); + if (!cascAodElement.has_cascData()) { + continue; + } + auto casc = cascAodElement.cascData_as(); + auto trackCharmBachelor0 = rowTrackIndexXicPlus.prong0_as(); + auto trackCharmBachelor1 = rowTrackIndexXicPlus.prong1_as(); + auto collision = rowTrackIndexXicPlus.collision(); + + // preselect cascade candidates + if (doCascadePreselection) { + if (std::abs(casc.dcaXYCascToPV()) > dcaXYToPVCascadeMax) { + continue; + } + if (std::abs(casc.mXi() - MassXiMinus) > massToleranceCascade) { + continue; + } + } + + //----------------------Set the magnetic field from ccdb--------------------------------------- + /// The static instance of the propagator was already modified in the HFTrackIndexSkimCreator, + /// but this is not true when running on Run2 data/MC already converted into AO2Ds. + auto bc = collision.bc_as(); + if (runNumber != bc.runNumber()) { + LOG(info) << ">>>>>>>>>>>> Current run number: " << runNumber; + initCCDB(bc, runNumber, ccdb, isRun2 ? ccdbPathGrp : ccdbPathGrpMag, lut, isRun2); + bz = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << ">>>>>>>>>>>> Magnetic field: " << bz; + } + df.setBz(bz); + + //----------------accessing particles in the decay chain------------- + auto trackPionFromXi = casc.bachelor_as(); // pion <- xi track from TracksWCovDca table + + //--------------------------info of V0 and cascades track from LF-tables--------------------------- + std::array vertexV0 = {casc.xlambda(), casc.ylambda(), casc.zlambda()}; + std::array pVecV0 = {casc.pxlambda(), casc.pylambda(), casc.pzlambda()}; + std::array vertexCasc = {casc.x(), casc.y(), casc.z()}; + std::array pVecCasc = {casc.px(), casc.py(), casc.pz()}; + std::array covCasc = {0.}; + + //----------------create cascade track------------------------------------------------------------ + constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int i = 0; i < 6; i++) { + covCasc[MomInd[i]] = casc.momentumCovMat()[i]; + covCasc[i] = casc.positionCovMat()[i]; + } + // create cascade track + o2::track::TrackParCov trackCasc; + if (trackPionFromXi.sign() > 0) { + trackCasc = o2::track::TrackParCov(vertexCasc, pVecCasc, covCasc, 1, true); + } else if (trackPionFromXi.sign() < 0) { + trackCasc = o2::track::TrackParCov(vertexCasc, pVecCasc, covCasc, -1, true); + } else { + continue; + } + trackCasc.setAbsCharge(1); + trackCasc.setPID(o2::track::PID::XiMinus); + + //----------------------------fit SV and create XicPlus track------------------ + auto trackParCovCharmBachelor0 = getTrackParCov(trackCharmBachelor0); + auto trackParCovCharmBachelor1 = getTrackParCov(trackCharmBachelor1); + + // reconstruct the 3-prong secondary vertex + if (df.process(trackCasc, trackParCovCharmBachelor0, trackParCovCharmBachelor1) == 0) { + continue; + } + + //----------------------------calculate physical properties----------------------- + // Charge of charm baryon + int signXic = casc.sign() < 0 ? +1 : -1; + + // get SV properties + const auto& secondaryVertex = df.getPCACandidate(); + auto chi2SV = df.getChi2AtPCACandidate(); + auto covMatrixSV = df.calcPCACovMatrixFlat(); + + // get track momenta + trackCasc = df.getTrack(0); + trackParCovCharmBachelor0 = df.getTrack(1); + trackParCovCharmBachelor1 = df.getTrack(2); + std::array pVecXi; + std::array pVecPi0; + std::array pVecPi1; + trackCasc.getPxPyPzGlo(pVecXi); + trackParCovCharmBachelor0.getPxPyPzGlo(pVecPi0); + trackParCovCharmBachelor1.getPxPyPzGlo(pVecPi1); + + // get invariant mass of Xic candidate + auto arrayMomenta = std::array{pVecXi, pVecPi0, pVecPi1}; + massXiPiPi = RecoDecay::m(std::move(arrayMomenta), std::array{MassXiMinus, MassPiPlus, MassPiPlus}); + + // get track impact parameters + // This modifies track momenta! + auto primaryVertex = getPrimaryVertex(collision); + auto covMatrixPV = primaryVertex.getCov(); + // calculate impact parameter + o2::dataformats::DCA impactParameterCasc; + o2::dataformats::DCA impactParameter0; + o2::dataformats::DCA impactParameter1; + trackCasc.propagateToDCA(primaryVertex, bz, &impactParameterCasc); + trackParCovCharmBachelor0.propagateToDCA(primaryVertex, bz, &impactParameter0); + trackParCovCharmBachelor1.propagateToDCA(primaryVertex, bz, &impactParameter1); + + // calculate cosine of pointing angle + std::array pvCoord = {collision.posX(), collision.posY(), collision.posZ()}; + double cpaLambda = casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + double cpaXYLambda = RecoDecay::cpaXY(pvCoord, vertexV0, pVecV0); + double cpaXi = casc.casccosPA(collision.posX(), collision.posY(), collision.posZ()); + double cpaXYXi = RecoDecay::cpaXY(pvCoord, vertexCasc, pVecCasc); + + // get invariant mass of Xi-pi pairs + auto arrayMomentaXiPi0 = std::array{pVecXi, pVecPi0}; + massXiPi0 = RecoDecay::m(std::move(arrayMomentaXiPi0), std::array{MassXiMinus, MassPiPlus}); + auto arrayMomentaXiPi1 = std::array{pVecXi, pVecPi1}; + massXiPi1 = RecoDecay::m(std::move(arrayMomentaXiPi1), std::array{MassXiMinus, MassPiPlus}); + + // get uncertainty of the decay length + double 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.)); + + //--------------------------------------------fill histograms---------------------------------------------------------------- + if (fillHistograms) { + // invariant mass + registry.fill(HIST("hMass3"), massXiPiPi); + // covariance matrix elements of PV + registry.fill(HIST("hCovPVXX"), covMatrixPV[0]); + registry.fill(HIST("hCovPVYY"), covMatrixPV[2]); + registry.fill(HIST("hCovPVXZ"), covMatrixPV[3]); + registry.fill(HIST("hCovPVZZ"), covMatrixPV[5]); + // covariance matrix elements of SV + registry.fill(HIST("hCovSVXX"), covMatrixSV[0]); + registry.fill(HIST("hCovSVYY"), covMatrixSV[2]); + registry.fill(HIST("hCovSVXZ"), covMatrixSV[3]); + registry.fill(HIST("hCovSVZZ"), covMatrixSV[5]); + // DCAs of prongs + registry.fill(HIST("hDcaXYProngs"), trackCasc.getPt(), impactParameterCasc.getY()); + registry.fill(HIST("hDcaXYProngs"), trackCharmBachelor0.pt(), impactParameter0.getY()); + registry.fill(HIST("hDcaXYProngs"), trackCharmBachelor1.pt(), impactParameter1.getY()); + registry.fill(HIST("hDcaZProngs"), trackCasc.getPt(), impactParameterCasc.getZ()); + registry.fill(HIST("hDcaZProngs"), trackCharmBachelor0.pt(), impactParameter0.getZ()); + registry.fill(HIST("hDcaZProngs"), trackCharmBachelor1.pt(), impactParameter1.getZ()); + } + + //---------------------------------fill candidate table rows------------------------------------------------------------------------------------------- + rowCandidateBase(collision.globalIndex(), + primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(), + covMatrixPV[0], covMatrixPV[2], covMatrixPV[5], + /*3-prong specific columns*/ + rowTrackIndexXicPlus.cascadeId(), rowTrackIndexXicPlus.prong0Id(), rowTrackIndexXicPlus.prong1Id(), + casc.bachelorId(), casc.posTrackId(), casc.negTrackId(), + secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], + covMatrixSV[0], covMatrixSV[2], covMatrixSV[5], + errorDecayLength, errorDecayLengthXY, + chi2SV, massXiPiPi, signXic, + pVecXi[0], pVecXi[1], pVecXi[2], + pVecPi0[0], pVecPi0[1], pVecPi0[2], + pVecPi1[0], pVecPi1[1], pVecPi1[2], + impactParameterCasc.getY(), impactParameter0.getY(), impactParameter1.getY(), + std::sqrt(impactParameterCasc.getSigmaY2()), std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()), + /*cascade specific columns*/ + vertexCasc[0], vertexCasc[1], vertexCasc[2], + vertexV0[0], vertexV0[1], vertexV0[2], + cpaXi, cpaXYXi, cpaLambda, cpaXYLambda, + massXiPi0, massXiPi1); + } // loop over track triplets + } + PROCESS_SWITCH(HfCandidateCreatorXicToXiPiPi, processXicplusWithDcaFitter, "Run candidate creator with DCAFitter.", true); + + void processXicplusWithKFParticle(aod::Collisions const&, + aod::HfCascLf3Prongs const& rowsTrackIndexXicPlus, + KFCascadesLinked const&, + KFCascFull const&, + aod::TracksWCovExtra const&, + aod::BCsWithTimestamps const&) + { + // loop over triplets of track indices + for (const auto& rowTrackIndexXicPlus : rowsTrackIndexXicPlus) { + auto cascAodElement = rowTrackIndexXicPlus.cascade_as(); + if (!cascAodElement.has_kfCascData()) { + continue; + } + auto casc = cascAodElement.kfCascData_as(); + auto trackCharmBachelor0 = rowTrackIndexXicPlus.prong0_as(); + auto trackCharmBachelor1 = rowTrackIndexXicPlus.prong1_as(); + auto collision = rowTrackIndexXicPlus.collision(); + + //-------------------preselect cascade candidates-------------------------------------- + if (doCascadePreselection) { + if (std::abs(casc.dcaXYCascToPV()) > dcaXYToPVCascadeMax) { + continue; + } + if (std::abs(casc.mXi() - MassXiMinus) > massToleranceCascade) { + continue; + } + } + + //----------------------Set the magnetic field from ccdb----------------------------- + /// The static instance of the propagator was already modified in the HFTrackIndexSkimCreator, + /// but this is not true when running on Run2 data/MC already converted into AO2Ds. + auto bc = collision.bc_as(); + if (runNumber != bc.runNumber()) { + LOG(info) << ">>>>>>>>>>>> Current run number: " << runNumber; + initCCDB(bc, runNumber, ccdb, isRun2 ? ccdbPathGrp : ccdbPathGrpMag, lut, isRun2); + bz = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << ">>>>>>>>>>>> Magnetic field: " << bz; + } + KFParticle::SetField(bz); + + //----------------------info of V0 and cascade tracks from LF-table------------------ + std::array vertexV0 = {casc.xlambda(), casc.ylambda(), casc.zlambda()}; + std::array pVecV0 = {casc.pxlambda(), casc.pylambda(), casc.pzlambda()}; + std::array vertexCasc = {casc.x(), casc.y(), casc.z()}; + std::array pVecCasc = {casc.px(), casc.py(), casc.pz()}; + + //----------------------Create XicPlus as KFParticle object------------------------------------------- + // initialize primary vertex + KFPVertex kfpVertex = createKFPVertexFromCollision(collision); + float covMatrixPV[6]; + kfpVertex.GetCovarianceMatrix(covMatrixPV); + KFParticle KFPV(kfpVertex); // for calculation of DCAs to PV + + // convert pion tracks into KFParticle object + KFPTrack kfpTrackCharmBachelor0 = createKFPTrackFromTrack(trackCharmBachelor0); + KFPTrack kfpTrackCharmBachelor1 = createKFPTrackFromTrack(trackCharmBachelor1); + KFParticle kfCharmBachelor0(kfpTrackCharmBachelor0, kPiPlus); + KFParticle kfCharmBachelor1(kfpTrackCharmBachelor1, kPiPlus); + + // create Xi as KFParticle object + // read {X,Y,Z,Px,Py,Pz} and corresponding covariance matrix from KF cascade Tables + std::array xyzpxpypz = {casc.x(), casc.y(), casc.z(), casc.px(), casc.py(), casc.pz()}; + float parPosMom[6]; + for (int i{0}; i < 6; ++i) { + parPosMom[i] = xyzpxpypz[i]; + } + // create KFParticle + KFParticle kfXi; + kfXi.Create(parPosMom, casc.kfTrackCovMat(), casc.sign(), casc.mXi()); + + // create XicPlus as KFParticle object + KFParticle kfXicPlus; + const KFParticle* kfDaughtersXicPlus[3] = {&kfCharmBachelor0, &kfCharmBachelor1, &kfXi}; + kfXicPlus.SetConstructMethod(kfConstructMethod); + try { + kfXicPlus.Construct(kfDaughtersXicPlus, 3); + } catch (std::runtime_error& e) { + LOG(debug) << "Failed to construct XicPlus : " << e.what(); + continue; + } + + // topological constraint + if (constrainXicPlusToPv) { + kfXicPlus.SetProductionVertex(KFPV); + } + auto covMatrixXicPlus = kfXicPlus.CovarianceMatrix(); + + // transport daughter particles to XicPlus decay vertex + kfCharmBachelor0.TransportToParticle(kfXicPlus); + kfCharmBachelor1.TransportToParticle(kfXicPlus); + kfXi.TransportToParticle(kfXicPlus); + + //---------------------calculate physical parameters of XicPlus candidate---------------------- + // sign of charm baryon + int signXic = casc.sign() < 0 ? +1 : -1; + + // get impact parameters of XicPlus daughters + float impactParameterPi0XY = 0., errImpactParameterPi0XY = 0.; + float impactParameterPi1XY = 0., errImpactParameterPi1XY = 0.; + float impactParameterXiXY = 0., errImpactParameterXiXY = 0.; + kfCharmBachelor0.GetDistanceFromVertexXY(KFPV, impactParameterPi0XY, errImpactParameterPi0XY); + kfCharmBachelor1.GetDistanceFromVertexXY(KFPV, impactParameterPi1XY, errImpactParameterPi1XY); + kfXi.GetDistanceFromVertexXY(KFPV, impactParameterXiXY, errImpactParameterXiXY); + + // calculate cosine of pointing angle + std::array pvCoord = {collision.posX(), collision.posY(), collision.posZ()}; + double cpaLambda = casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + double cpaXYLambda = RecoDecay::cpaXY(pvCoord, vertexV0, pVecV0); + double cpaXi = casc.casccosPA(collision.posX(), collision.posY(), collision.posZ()); + double cpaXYXi = RecoDecay::cpaXY(pvCoord, vertexCasc, pVecCasc); + + // get DCAs of Pi0-Pi1, Pi0-Xi, Pi1-Xi + float dcaXYPi0Pi1 = kfCharmBachelor0.GetDistanceFromParticleXY(kfCharmBachelor1); + float dcaXYPi0Xi = kfCharmBachelor0.GetDistanceFromParticleXY(kfXi); + float dcaXYPi1Xi = kfCharmBachelor1.GetDistanceFromParticleXY(kfXi); + + // mass of Xi-Pi0 pair + KFParticle kfXiPi0; + const KFParticle* kfXiResonanceDaughtersPi0[2] = {&kfXi, &kfCharmBachelor0}; + kfXiPi0.SetConstructMethod(kfConstructMethod); + try { + kfXiPi0.Construct(kfXiResonanceDaughtersPi0, 2); + massXiPi0 = kfXiPi0.GetMass(); + } catch (...) { + LOG(info) << "Failed to construct Xi(1530) with Pi 0"; + } + + // mass of Xi-Pi1 pair + KFParticle kfXiPi1; + const KFParticle* kfXiResonanceDaughtersPi1[2] = {&kfXi, &kfCharmBachelor1}; + kfXiPi1.SetConstructMethod(kfConstructMethod); + try { + kfXiPi1.Construct(kfXiResonanceDaughtersPi1, 2); + massXiPi1 = kfXiPi1.GetMass(); + } catch (...) { + LOG(info) << "Failed to construct Xi(1530) with Pi 1"; + } + + //-------------------------------fill histograms-------------------------------------------- + if (fillHistograms) { + // invariant mass + registry.fill(HIST("hMass3"), kfXicPlus.GetMass()); + // covariance matrix elements of PV + registry.fill(HIST("hCovPVXX"), covMatrixPV[0]); + registry.fill(HIST("hCovPVYY"), covMatrixPV[2]); + registry.fill(HIST("hCovPVXZ"), covMatrixPV[3]); + registry.fill(HIST("hCovPVZZ"), covMatrixPV[5]); + // covariance matrix elements of SV + registry.fill(HIST("hCovSVXX"), covMatrixXicPlus[0]); + registry.fill(HIST("hCovSVYY"), covMatrixXicPlus[2]); + registry.fill(HIST("hCovSVXZ"), covMatrixXicPlus[3]); + registry.fill(HIST("hCovSVZZ"), covMatrixXicPlus[5]); + // DCAs of prongs + registry.fill(HIST("hDcaXYProngs"), kfXi.GetPt(), impactParameterXiXY); + registry.fill(HIST("hDcaXYProngs"), kfCharmBachelor0.GetPt(), impactParameterPi0XY); + registry.fill(HIST("hDcaXYProngs"), kfCharmBachelor1.GetPt(), impactParameterPi1XY); + } + + //------------------------------fill candidate table rows-------------------------------------- + rowCandidateBase(collision.globalIndex(), + KFPV.GetX(), KFPV.GetY(), KFPV.GetZ(), + covMatrixPV[0], covMatrixPV[2], covMatrixPV[5], + /*3-prong specific columns*/ + rowTrackIndexXicPlus.cascadeId(), rowTrackIndexXicPlus.prong0Id(), rowTrackIndexXicPlus.prong1Id(), + casc.bachelorId(), casc.posTrackId(), casc.negTrackId(), + kfXicPlus.GetX(), kfXicPlus.GetY(), kfXicPlus.GetZ(), + kfXicPlus.GetErrX(), kfXicPlus.GetErrY(), kfXicPlus.GetErrZ(), + kfXicPlus.GetErrDecayLength(), kfXicPlus.GetErrDecayLengthXY(), + kfXicPlus.GetChi2(), kfXicPlus.GetMass(), signXic, + kfXi.GetPx(), kfXi.GetPy(), kfXi.GetPz(), + kfCharmBachelor0.GetPx(), kfCharmBachelor0.GetPy(), kfCharmBachelor0.GetPz(), + kfCharmBachelor1.GetPx(), kfCharmBachelor1.GetPy(), kfCharmBachelor1.GetPz(), + impactParameterXiXY, impactParameterPi0XY, impactParameterPi1XY, + errImpactParameterXiXY, errImpactParameterPi0XY, errImpactParameterPi1XY, + /*cascade specific columns*/ + casc.x(), casc.y(), casc.z(), + casc.xlambda(), casc.ylambda(), casc.zlambda(), + cpaXi, cpaXYXi, cpaLambda, cpaXYLambda, + massXiPi0, massXiPi1); + rowCandidateKF(casc.kfCascadeChi2(), casc.kfV0Chi2(), + dcaXYPi0Pi1, dcaXYPi0Xi, dcaXYPi1Xi); + } // loop over track triplets + } + PROCESS_SWITCH(HfCandidateCreatorXicToXiPiPi, processXicplusWithKFParticle, "Run candidate creator with KFParticle using derived data from HfTrackIndexSkimCreatorLfCascades.", false); +}; // struct + +/// Performs MC matching. +struct HfCandidateCreatorXicToXiPiPiExpressions { + Spawns rowCandidateXic; + Produces rowMcMatchRec; + Produces rowMcMatchGen; + + void init(InitContext const&) {} + + void processMc(aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles) + { + rowCandidateXic->bindExternalIndices(&tracks); + + int indexRec = -1; + int indexRecXicPlus = -1; + int8_t sign = 0; + int8_t flag = 0; + int8_t origin = 0; + int8_t debug = 0; + // for resonance matching: + std::vector arrDaughIndex; + std::array arrPDGDaugh; + std::array arrXiResonance = {3324, kPiPlus}; // 3324: Ξ(1530) + + // Match reconstructed candidates. + for (const auto& candidate : *rowCandidateXic) { + flag = 0; + sign = 0; + origin = RecoDecay::OriginType::None; + debug = 0; + arrDaughIndex.clear(); + + auto arrayDaughters = std::array{candidate.pi0_as(), // pi <- Xic + candidate.pi1_as(), // pi <- Xic + candidate.bachelor_as(), // pi <- cascade + candidate.posTrack_as(), // p <- lambda + candidate.negTrack_as()}; // pi <- lambda + auto arrayDaughtersCasc = std::array{candidate.bachelor_as(), + candidate.posTrack_as(), + candidate.negTrack_as()}; + auto arrayDaughtersV0 = std::array{candidate.posTrack_as(), + candidate.negTrack_as()}; + + // Xic → pi pi pi pi p + indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, Pdg::kXiCPlus, std::array{+kPiPlus, +kPiPlus, +kPiMinus, +kProton, +kPiMinus}, true, &sign, 4); + indexRecXicPlus = indexRec; + if (indexRec == -1) { + debug = 1; + } + if (indexRec > -1) { + // Xi- → pi pi p + indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughtersCasc, +kXiMinus, std::array{+kPiMinus, +kProton, +kPiMinus}, true, &sign, 2); + if (indexRec == -1) { + debug = 2; + } + if (indexRec > -1) { + // Lambda → p pi + indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughtersV0, +kLambda0, std::array{+kProton, +kPiMinus}, true, &sign, 1); + if (indexRec == -1) { + debug = 3; + } + if (indexRec > -1) { + RecoDecay::getDaughters(mcParticles.rawIteratorAt(indexRecXicPlus), &arrDaughIndex, std::array{0}, 1); + if (arrDaughIndex.size() == 2) { + for (auto iProng = 0u; iProng < arrDaughIndex.size(); ++iProng) { + auto daughI = mcParticles.rawIteratorAt(arrDaughIndex[iProng]); + arrPDGDaugh[iProng] = std::abs(daughI.pdgCode()); + } + if ((arrPDGDaugh[0] == arrXiResonance[0] && arrPDGDaugh[1] == arrXiResonance[1]) || (arrPDGDaugh[0] == arrXiResonance[1] && arrPDGDaugh[1] == arrXiResonance[0])) { + flag = sign * (1 << aod::hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiResPiToXiPiPi); + } else { + debug = 4; + } + } else { + flag = sign * (1 << aod::hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiPiPi); + } + } + } + } + + // Check whether the charm baryon is non-prompt (from a b quark). + if (flag != 0) { + auto particle = mcParticles.rawIteratorAt(indexRecXicPlus); + origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, true); + } + + rowMcMatchRec(flag, debug, origin); + } // close loop over candidates + + // Match generated particles. + for (const auto& particle : mcParticles) { + flag = 0; + sign = 0; + debug = 0; + origin = RecoDecay::OriginType::None; + arrDaughIndex.clear(); + + // Xic → Xi pi pi + if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kXiCPlus, std::array{+kXiMinus, +kPiPlus, +kPiPlus}, true, &sign, 2)) { + debug = 1; + // Xi- -> Lambda pi + auto cascMC = mcParticles.rawIteratorAt(particle.daughtersIds().front()); + if (RecoDecay::isMatchedMCGen(mcParticles, cascMC, +kXiMinus, std::array{+kLambda0, +kPiMinus}, true)) { + debug = 2; + // Lambda -> p pi + auto v0MC = mcParticles.rawIteratorAt(cascMC.daughtersIds().front()); + if (RecoDecay::isMatchedMCGen(mcParticles, v0MC, +kLambda0, std::array{+kProton, +kPiMinus}, true)) { + debug = 3; + + RecoDecay::getDaughters(particle, &arrDaughIndex, std::array{0}, 1); + if (arrDaughIndex.size() == 2) { + for (auto iProng = 0u; iProng < arrDaughIndex.size(); ++iProng) { + auto daughI = mcParticles.rawIteratorAt(arrDaughIndex[iProng]); + arrPDGDaugh[iProng] = std::abs(daughI.pdgCode()); + } + if ((arrPDGDaugh[0] == arrXiResonance[0] && arrPDGDaugh[1] == arrXiResonance[1]) || (arrPDGDaugh[0] == arrXiResonance[1] && arrPDGDaugh[1] == arrXiResonance[0])) { + flag = sign * (1 << aod::hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiResPiToXiPiPi); + } else { + debug = 4; + } + } else { + flag = sign * (1 << aod::hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiPiPi); + } + } + } + } + + // Check whether the charm baryon is non-prompt (from a b quark). + if (flag != 0) { + origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, true); + } + + rowMcMatchGen(flag, debug, origin); + } // close loop over generated particles + } // close process + PROCESS_SWITCH(HfCandidateCreatorXicToXiPiPiExpressions, processMc, "Process MC", false); +}; // close struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/TableProducer/candidateSelectorXicToXiPiPi.cxx b/PWGHF/TableProducer/candidateSelectorXicToXiPiPi.cxx new file mode 100644 index 00000000000..4928ab194fc --- /dev/null +++ b/PWGHF/TableProducer/candidateSelectorXicToXiPiPi.cxx @@ -0,0 +1,252 @@ +// 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 candidateSelectorXicToXiPiPi.cxx +/// \brief Ξc± → Ξ∓ π± π± candidate selector +/// +/// \author Phil Lennart Stahlhut , CERN + +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + +#include "Common/Core/TrackSelectorPID.h" + +#include "PWGHF/Core/SelectorCuts.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/Utils/utilsAnalysis.h" // findBin function + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::analysis; + +struct HfCandidateSelectorXicToXiPiPi { + Produces hfSelXicToXiPiPiCandidate; + + Configurable ptCandMin{"ptCandMin", 0., "Lower bound of candidate pT"}; + Configurable ptCandMax{"ptCandMax", 36., "Upper bound of candidate pT"}; + // topological cuts + Configurable> binsPt{"binsPt", std::vector{hf_cuts_xic_to_xi_pi_pi::vecBinsPt}, "pT bin limits"}; + Configurable> cuts{"cuts", {hf_cuts_xic_to_xi_pi_pi::cuts[0], hf_cuts_xic_to_xi_pi_pi::nBinsPt, hf_cuts_xic_to_xi_pi_pi::nCutVars, hf_cuts_xic_to_xi_pi_pi::labelsPt, hf_cuts_xic_to_xi_pi_pi::labelsCutVar}, "Xicplus candidate selection per pT bin"}; + // QA switch + Configurable activateQA{"activateQA", false, "Flag to enable QA histogram"}; + // Enable PID + Configurable usePid{"usePid", true, "Switch for PID selection at track level"}; + Configurable acceptPIDNotApplicable{"acceptPIDNotApplicable", true, "Switch to accept Status::NotApplicable [(NotApplicable for one detector) and (NotApplicable or Conditional for the other)] in PID selection"}; + // TPC PID + Configurable ptPidTpcMin{"ptPidTpcMin", 0.15, "Lower bound of track pT for TPC PID"}; + Configurable ptPidTpcMax{"ptPidTpcMax", 20., "Upper bound of track pT for TPC PID"}; + Configurable nSigmaTpcMax{"nSigmaTpcMax", 5., "Nsigma cut on TPC only"}; + Configurable nSigmaTpcCombinedMax{"nSigmaTpcCombinedMax", 5., "Nsigma cut on TPC combined with TOF"}; + // TOF PID + Configurable ptPidTofMin{"ptPidTofMin", 0.15, "Lower bound of track pT for TOF PID"}; + Configurable ptPidTofMax{"ptPidTofMax", 20., "Upper bound of track pT for TOF PID"}; + Configurable nSigmaTofMax{"nSigmaTofMax", 5., "Nsigma cut on TOF only"}; + Configurable nSigmaTofCombinedMax{"nSigmaTofCombinedMax", 5., "Nsigma cut on TOF combined with TPC"}; + + TrackSelectorPi selectorPion; + TrackSelectorPr selectorProton; + + using TracksPidWithSel = soa::Join; + + HistogramRegistry registry{"registry"}; + + void init(InitContext& initContext) + { + if (usePid) { + // pion + selectorPion.setRangePtTpc(ptPidTpcMin, ptPidTpcMax); + selectorPion.setRangeNSigmaTpc(-nSigmaTpcMax, nSigmaTpcMax); + selectorPion.setRangeNSigmaTpcCondTof(-nSigmaTpcCombinedMax, nSigmaTpcCombinedMax); + selectorPion.setRangePtTof(ptPidTofMin, ptPidTofMax); + selectorPion.setRangeNSigmaTof(-nSigmaTofMax, nSigmaTofMax); + selectorPion.setRangeNSigmaTofCondTpc(-nSigmaTofCombinedMax, nSigmaTofCombinedMax); + // proton + selectorProton.setRangePtTpc(ptPidTpcMin, ptPidTpcMax); + selectorProton.setRangeNSigmaTpc(-nSigmaTpcMax, nSigmaTpcMax); + selectorProton.setRangeNSigmaTpcCondTof(-nSigmaTpcCombinedMax, nSigmaTpcCombinedMax); + selectorProton.setRangePtTof(ptPidTofMin, ptPidTofMax); + selectorProton.setRangeNSigmaTof(-nSigmaTofMax, nSigmaTofMax); + selectorProton.setRangeNSigmaTofCondTpc(-nSigmaTofCombinedMax, nSigmaTofCombinedMax); + } + + if (activateQA) { + constexpr int kNBinsSelections = 1 + SelectionStep::NSelectionSteps; + std::string labels[kNBinsSelections]; + labels[0] = "No selection"; + labels[1 + SelectionStep::RecoSkims] = "Skims selection"; + labels[1 + SelectionStep::RecoTopol] = "Skims & Topological selections"; + labels[1 + SelectionStep::RecoPID] = "Skims & Topological & PID selections"; + labels[1 + SelectionStep::RecoMl] = "Skims & Topological & PID & ML selections"; + static const AxisSpec axisSelections = {kNBinsSelections, 0.5, kNBinsSelections + 0.5, ""}; + registry.add("hSelections", "Selections;;#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {axisSelections, {(std::vector)binsPt, "#it{p}_{T} (GeV/#it{c})"}}}); + for (int iBin = 0; iBin < kNBinsSelections; ++iBin) { + registry.get(HIST("hSelections"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin].data()); + } + } + } + + /// Conjugate-independent topological cuts + /// \param candidate is candidate + /// \return true if candidate passes all cuts + template + bool selectionTopol(const T1& hfCandXic) + { + auto candpT = hfCandXic.pt(); + int pTBin = findBin(binsPt, candpT); + if (pTBin == -1) { + return false; + } + + // check that the candidate pT is within the analysis range + if (candpT < ptCandMin || candpT >= ptCandMax) { + return false; + } + + // check candidate mass is within a defined mass window + if (std::abs(hfCandXic.invMassXic() - o2::constants::physics::MassXiCPlus) > cuts->get(pTBin, "m")) { + return false; + } + + // cosine of pointing angle + if (hfCandXic.cpa() <= cuts->get(pTBin, "cos pointing angle")) { + return false; + } + + // cosine of pointing angle XY + if (hfCandXic.cpaXY() <= cuts->get(pTBin, "cos pointing angle XY")) { + return false; + } + + // candidate maximum decay length + if (hfCandXic.decayLength() > cuts->get(pTBin, "max decay length")) { + return false; + } + + // candidate maximum decay length XY + if (hfCandXic.decayLengthXY() > cuts->get(pTBin, "max decay length XY")) { + return false; + } + + // candidate chi2PC + if (hfCandXic.chi2PCA() > cuts->get(pTBin, "chi2PCA")) { + return false; + } + + // maximum DCA of daughters + if ((std::abs(hfCandXic.impactParameter0()) > cuts->get(pTBin, "max impParXY Xi")) || + (std::abs(hfCandXic.impactParameter1()) > cuts->get(pTBin, "max impParXY Pi0")) || + (std::abs(hfCandXic.impactParameter2()) > cuts->get(pTBin, "max impParXY Pi1"))) { + return false; + } + + // cut on daughter pT + if (hfCandXic.ptProng0() < cuts->get(pTBin, "pT Xi") || + hfCandXic.ptProng1() < cuts->get(pTBin, "pT Pi0") || + hfCandXic.ptProng2() < cuts->get(pTBin, "pT Pi1")) { + return false; + } + + return true; + } + + /// Apply PID selection + /// \param pidTrackPi0 PID status of trackPi0 (prong1 of Xic candidate) + /// \param pidTrackPi1 PID status of trackPi1 (prong2 of Xic candidate) + /// \param pidTrackPr PID status of trackPr (positive daughter of V0 candidate) + /// \param pidTrackPiLam PID status of trackPiLam (negative daughter of V0 candidate) + /// \param pidTrackPiXi PID status of trackPiXi (Bachelor of cascade candidate) + /// \param acceptPIDNotApplicable switch to accept Status::NotApplicable + /// \return true if prongs of Xic candidate pass all selections + bool selectionPid(TrackSelectorPID::Status const pidTrackPi0, + TrackSelectorPID::Status const pidTrackPi1, + TrackSelectorPID::Status const pidTrackPr, + TrackSelectorPID::Status const pidTrackPiLam, + TrackSelectorPID::Status const pidTrackPiXi, + bool const acceptPIDNotApplicable) + { + if (!acceptPIDNotApplicable && (pidTrackPi0 != TrackSelectorPID::Accepted || pidTrackPi1 != TrackSelectorPID::Accepted || pidTrackPr != TrackSelectorPID::Accepted || pidTrackPiLam != TrackSelectorPID::Accepted || pidTrackPiXi != TrackSelectorPID::Accepted)) { + return false; + } + if (acceptPIDNotApplicable && (pidTrackPi0 == TrackSelectorPID::Rejected || pidTrackPi1 == TrackSelectorPID::Rejected || pidTrackPr == TrackSelectorPID::Rejected || pidTrackPiLam == TrackSelectorPID::Rejected || pidTrackPiXi == TrackSelectorPID::Rejected)) { + return false; + } + return true; + } + + void process(aod::HfCandXic const& hfCandsXic, + TracksPidWithSel const&) + { + for (const auto& hfCandXic : hfCandsXic) { + int statusXicToXiPiPi = 0; + auto ptCandXic = hfCandXic.pt(); + + if (activateQA) { + registry.fill(HIST("hSelections"), 1, ptCandXic); + } + + // No hfFlag -> by default skim selected + SETBIT(statusXicToXiPiPi, SelectionStep::RecoSkims); // RecoSkims = 0 --> statusXicToXiPiPi = 1 + if (activateQA) { + registry.fill(HIST("hSelections"), 2 + SelectionStep::RecoSkims, ptCandXic); + } + + // topological cuts + if (!selectionTopol(hfCandXic)) { + hfSelXicToXiPiPiCandidate(statusXicToXiPiPi); + continue; + } + SETBIT(statusXicToXiPiPi, SelectionStep::RecoTopol); // RecoTopol = 1 --> statusXicToXiPiPi = 3 + if (activateQA) { + registry.fill(HIST("hSelections"), 2 + SelectionStep::RecoTopol, ptCandXic); + } + + // track-level PID selection + if (usePid) { + auto trackPi0 = hfCandXic.pi0_as(); + auto trackPi1 = hfCandXic.pi1_as(); + auto trackV0PosDau = hfCandXic.posTrack_as(); + auto trackV0NegDau = hfCandXic.negTrack_as(); + auto trackPiFromXi = hfCandXic.bachelor_as(); + // assign proton and pion hypothesis to V0 daughters + auto trackPr = trackV0PosDau; + auto trackPiFromLam = trackV0NegDau; + if (hfCandXic.sign() < 0) { + trackPr = trackV0NegDau; + trackPiFromLam = trackV0PosDau; + } + // PID info + TrackSelectorPID::Status pidTrackPi0 = selectorPion.statusTpcAndTof(trackPi0); + TrackSelectorPID::Status pidTrackPi1 = selectorPion.statusTpcAndTof(trackPi1); + TrackSelectorPID::Status pidTrackPr = selectorProton.statusTpcAndTof(trackPr); + TrackSelectorPID::Status pidTrackPiLam = selectorPion.statusTpcAndTof(trackPiFromLam); + TrackSelectorPID::Status pidTrackPiXi = selectorPion.statusTpcAndTof(trackPiFromXi); + + if (!selectionPid(pidTrackPi0, pidTrackPi1, pidTrackPr, pidTrackPiLam, pidTrackPiXi, acceptPIDNotApplicable.value)) { + hfSelXicToXiPiPiCandidate(statusXicToXiPiPi); + continue; + } + SETBIT(statusXicToXiPiPi, SelectionStep::RecoPID); // RecoPID = 2 --> statusXicToXiPiPi = 7 + if (activateQA) { + registry.fill(HIST("hSelections"), 2 + SelectionStep::RecoPID, ptCandXic); + } + } + + hfSelXicToXiPiPiCandidate(statusXicToXiPiPi); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/TableProducer/treeCreatorXicToXiPiPi.cxx b/PWGHF/TableProducer/treeCreatorXicToXiPiPi.cxx new file mode 100644 index 00000000000..c28709ffc75 --- /dev/null +++ b/PWGHF/TableProducer/treeCreatorXicToXiPiPi.cxx @@ -0,0 +1,792 @@ +// 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 treeCreatorXicToXiPiPi.cxx +/// \brief Writer of Ξc± → Ξ∓ π± π± candidates in the form of flat tables to be stored in TTrees. +/// +/// \author Phil Lennart Stahlhut , CERN + +#include "CommonConstants/PhysicsConstants.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +namespace o2::aod +{ +namespace full +{ +// track indices +DECLARE_SOA_COLUMN(CandidateSelFlag, candidateSelFlag, int); //! Selection flag of candidate (output of candidateSelector) +DECLARE_SOA_INDEX_COLUMN_FULL(Xi, xi, int, Tracks, "_pi0"); +DECLARE_SOA_INDEX_COLUMN_FULL(Pi0, pi0, int, Tracks, "_pi0"); +DECLARE_SOA_INDEX_COLUMN_FULL(Pi1, pi1, int, Tracks, "_pi1"); +// vertices +DECLARE_SOA_COLUMN(XPv, xPv, float); +DECLARE_SOA_COLUMN(YPv, yPv, float); +DECLARE_SOA_COLUMN(ZPv, zPv, float); +DECLARE_SOA_COLUMN(XPvErr, xPvErr, float); +DECLARE_SOA_COLUMN(YPvErr, yPvErr, float); +DECLARE_SOA_COLUMN(ZPvErr, zPvErr, float); +DECLARE_SOA_COLUMN(XSv, xSv, float); +DECLARE_SOA_COLUMN(YSv, ySv, float); +DECLARE_SOA_COLUMN(ZSv, zSv, float); +DECLARE_SOA_COLUMN(Chi2Sv, chi2Sv, float); +DECLARE_SOA_COLUMN(XSvErr, xSvErr, float); +DECLARE_SOA_COLUMN(YSvErr, ySvErr, float); +DECLARE_SOA_COLUMN(ZSvErr, zSvErr, float); +DECLARE_SOA_COLUMN(XDecVtxXi, xDecVtxXi, float); +DECLARE_SOA_COLUMN(YDecVtxXi, yDecVtxXi, float); +DECLARE_SOA_COLUMN(ZDecVtxXi, zDecVtxXi, float); +DECLARE_SOA_COLUMN(Chi2XiVtx, chi2XiVtx, float); +DECLARE_SOA_COLUMN(XDecVtxLam, xDecVtxLam, float); +DECLARE_SOA_COLUMN(YDecVtxLam, yDecVtxLam, float); +DECLARE_SOA_COLUMN(ZDecVtxLam, zDecVtxLam, float); +DECLARE_SOA_COLUMN(Chi2LamVtx, chi2LamVtx, float); +// properties of XicPlus +DECLARE_SOA_COLUMN(Sign, sign, float); +DECLARE_SOA_COLUMN(E, e, float); //! Energy of candidate (GeV) +DECLARE_SOA_COLUMN(M, m, float); //! Invariant mass of candidate (GeV/c2) +DECLARE_SOA_COLUMN(P, p, float); //! Momentum of candidate (GeV/c) +DECLARE_SOA_COLUMN(Pt, pt, float); //! Transverse momentum of candidate (GeV/c) +DECLARE_SOA_COLUMN(Y, y, float); //! Rapidity of candidate +DECLARE_SOA_COLUMN(Eta, eta, float); //! Pseudorapidity of candidate +DECLARE_SOA_COLUMN(Phi, phi, float); //! Azimuth angle of candidate +DECLARE_SOA_COLUMN(Ct, ct, float); //! Proper lifetime time ctau of candidate (cm) +DECLARE_SOA_COLUMN(DecayLength, decayLength, float); //! Decay length of candidate (cm) +DECLARE_SOA_COLUMN(DecayLengthXY, decayLengthXY, float); //! Transverse decay length of candidate (cm) +DECLARE_SOA_COLUMN(DecayLengthNormalised, decayLengthNormalised, float); //! Normalised decay length of candidate +DECLARE_SOA_COLUMN(DecayLengthXYNormalised, decayLengthXYNormalised, float); //! Normalised transverse decay length of candidate +DECLARE_SOA_COLUMN(Cpa, cpa, float); //! Cosine pointing angle of candidate +DECLARE_SOA_COLUMN(CpaXY, cpaXY, float); //! Cosine pointing angle of candidate in transverse plane +// properties of daughter tracks +DECLARE_SOA_COLUMN(PtXi, ptXi, float); //! Transverse momentum of Xi (prong0) (GeV/c) +DECLARE_SOA_COLUMN(ImpactParameterXi, impactParameterXi, float); //! Impact parameter of Xi (prong0) +DECLARE_SOA_COLUMN(ImpactParameterNormalisedXi, impactParameterNormalisedXi, float); //! Normalised impact parameter of Xi (prong0) +DECLARE_SOA_COLUMN(PtPi0, ptPi0, float); //! Transverse momentum of Pi0 (prong1) (GeV/c) +DECLARE_SOA_COLUMN(ImpactParameterPi0, impactParameterPi0, float); //! Impact parameter of Pi0 (prong1) +DECLARE_SOA_COLUMN(ImpactParameterNormalisedPi0, impactParameterNormalisedPi0, float); //! Normalised impact parameter of Pi0 (prong1) +DECLARE_SOA_COLUMN(PtPi1, ptPi1, float); //! Transverse momentum of Pi1 (prong2) (GeV/c) +DECLARE_SOA_COLUMN(ImpactParameterPi1, impactParameterPi1, float); //! Normalised impact parameter of Pi1 (prong2) +DECLARE_SOA_COLUMN(ImpactParameterNormalisedPi1, impactParameterNormalisedPi1, float); //! Normalised impact parameter of Pi1 (prong2) +DECLARE_SOA_COLUMN(MaxNormalisedDeltaIP, maxNormalisedDeltaIP, float); //! Maximum normalized difference between measured and expected impact parameter of candidate prongs +DECLARE_SOA_COLUMN(CpaXi, cpaXi, float); +DECLARE_SOA_COLUMN(CpaXYXi, cpaXYXi, float); +DECLARE_SOA_COLUMN(CpaLam, cpaLam, float); +DECLARE_SOA_COLUMN(CpaXYLam, cpaXYLam, float); +DECLARE_SOA_COLUMN(DcaPi0Pi1, dcaPi0Pi1, float); +DECLARE_SOA_COLUMN(DcaPi0Xi, dcaPi0Xi, float); +DECLARE_SOA_COLUMN(DcaPi1Xi, dcaPi1Xi, float); +DECLARE_SOA_COLUMN(InvMassXiPi0, invMassXiPi0, float); +DECLARE_SOA_COLUMN(InvMassXiPi1, invMassXiPi1, float); +} // namespace full + +DECLARE_SOA_TABLE(HfCandXicToXiPiPiLites, "AOD", "HFXICXI2PILITE", + full::CandidateSelFlag, + full::XPv, + full::YPv, + full::ZPv, + full::XSv, + full::YSv, + full::ZSv, + full::Chi2Sv, + full::Sign, + full::E, + full::M, + full::P, + full::Pt, + full::Y, + full::Eta, + full::Phi, + full::Ct, + full::DecayLength, + full::DecayLengthXY, + full::Cpa, + full::CpaXY, + full::PtXi, + full::PtPi0, + full::PtPi1, + full::ImpactParameterXi, + full::ImpactParameterPi0, + full::ImpactParameterPi1, + full::CpaXi, + full::CpaXYXi, + full::CpaLam, + full::CpaXYLam, + hf_cand_xic_to_xi_pi_pi::FlagMcMatchRec); + +DECLARE_SOA_TABLE(HfCandXicToXiPiPiLiteKfs, "AOD", "HFXICXI2PILITKF", + full::CandidateSelFlag, + full::XPv, + full::YPv, + full::ZPv, + full::XSv, + full::YSv, + full::ZSv, + full::Chi2Sv, + full::Sign, + full::E, + full::M, + full::P, + full::Pt, + full::Y, + full::Eta, + full::Phi, + full::Ct, + full::DecayLength, + full::DecayLengthXY, + full::Cpa, + full::CpaXY, + full::PtXi, + full::PtPi0, + full::PtPi1, + full::ImpactParameterXi, + full::ImpactParameterPi0, + full::ImpactParameterPi1, + full::CpaXi, + full::CpaXYXi, + full::CpaLam, + full::CpaXYLam, + full::Chi2XiVtx, + full::Chi2LamVtx, + full::DcaPi0Pi1, + full::DcaPi0Xi, + full::DcaPi1Xi, + hf_cand_xic_to_xi_pi_pi::FlagMcMatchRec); + +DECLARE_SOA_TABLE(HfCandXicToXiPiPiFulls, "AOD", "HFXICXI2PIFULL", + full::CandidateSelFlag, + full::XPv, + full::YPv, + full::ZPv, + full::XPvErr, + full::YPvErr, + full::ZPvErr, + full::XSv, + full::YSv, + full::ZSv, + full::Chi2Sv, + full::XSvErr, + full::YSvErr, + full::ZSvErr, + full::XDecVtxXi, + full::YDecVtxXi, + full::ZDecVtxXi, + full::XDecVtxLam, + full::YDecVtxLam, + full::ZDecVtxLam, + full::Sign, + full::E, + full::M, + full::P, + full::Pt, + full::Y, + full::Eta, + full::Phi, + full::Ct, + full::DecayLength, + full::DecayLengthNormalised, + full::DecayLengthXY, + full::DecayLengthXYNormalised, + full::Cpa, + full::CpaXY, + full::PtXi, + full::PtPi0, + full::PtPi1, + full::ImpactParameterXi, + full::ImpactParameterNormalisedXi, + full::ImpactParameterPi0, + full::ImpactParameterNormalisedPi0, + full::ImpactParameterPi1, + full::ImpactParameterNormalisedPi1, + full::MaxNormalisedDeltaIP, + full::CpaXi, + full::CpaXYXi, + full::CpaLam, + full::CpaXYLam, + full::InvMassXiPi0, + full::InvMassXiPi1, + hf_cand_xic_to_xi_pi_pi::FlagMcMatchRec); + +DECLARE_SOA_TABLE(HfCandXicToXiPiPiFullKfs, "AOD", "HFXICXI2PIFULKF", + full::CandidateSelFlag, + full::XPv, + full::YPv, + full::ZPv, + full::XPvErr, + full::YPvErr, + full::ZPvErr, + full::XSv, + full::YSv, + full::ZSv, + full::Chi2Sv, + full::XSvErr, + full::YSvErr, + full::ZSvErr, + full::XDecVtxXi, + full::YDecVtxXi, + full::ZDecVtxXi, + full::XDecVtxLam, + full::YDecVtxLam, + full::ZDecVtxLam, + full::Sign, + full::E, + full::M, + full::P, + full::Pt, + full::Y, + full::Eta, + full::Phi, + full::Ct, + full::DecayLength, + full::DecayLengthNormalised, + full::DecayLengthXY, + full::DecayLengthXYNormalised, + full::Cpa, + full::CpaXY, + full::PtXi, + full::PtPi0, + full::PtPi1, + full::ImpactParameterXi, + full::ImpactParameterNormalisedXi, + full::ImpactParameterPi0, + full::ImpactParameterNormalisedPi0, + full::ImpactParameterPi1, + full::ImpactParameterNormalisedPi1, + full::MaxNormalisedDeltaIP, + full::CpaXi, + full::CpaXYXi, + full::CpaLam, + full::CpaXYLam, + full::InvMassXiPi0, + full::InvMassXiPi1, + full::Chi2XiVtx, + full::Chi2LamVtx, + full::DcaPi0Pi1, + full::DcaPi0Xi, + full::DcaPi1Xi, + hf_cand_xic_to_xi_pi_pi::FlagMcMatchRec); + +DECLARE_SOA_TABLE(HfCandXicToXiPiPiDauInds, "AOD", "HFXICXI2PIDAUIN", + full::XiId, + full::Pi0Id, + full::Pi1Id); + +DECLARE_SOA_TABLE(HfCandXicToXiPiPiFullEvs, "AOD", "HFXICXI2PIFULEV", + collision::BCId, + collision::NumContrib, + collision::PosX, + collision::PosY, + collision::PosZ); + +DECLARE_SOA_TABLE(HfCandXicToXiPiPiFullPs, "AOD", "HFXICXI2PIFULLP", + collision::BCId, + full::Pt, + full::Eta, + full::Phi, + full::Y, + hf_cand_xic_to_xi_pi_pi::FlagMcMatchGen); +} // namespace o2::aod + +/// Writes the full information in an output TTree +struct HfTreeCreatorXicToXiPiPi { + Produces rowCandidateLite; + Produces rowCandidateLiteKf; + Produces rowCandidateFull; + Produces rowCandidateFullKf; + Produces rowCandidateDauIndices; + Produces rowCandidateFullEvents; + Produces rowCandidateFullParticles; + + Configurable selectionFlagXic{"selectionXic", 1, "Selection Flag for Xic"}; + Configurable fillCandidateLiteTable{"fillCandidateLiteTable", false, "Switch to fill lite table with candidate properties"}; + Configurable fillCandidateDauIndexTable{"fillCandidateDauIndexTable", false, "Switch to fill table with Xic daughters track indices"}; + // parameters for production of training samples + Configurable fillOnlySignal{"fillOnlySignal", false, "Flag to fill derived tables with signal for ML trainings"}; + Configurable fillOnlyBackground{"fillOnlyBackground", false, "Flag to fill derived tables with background for ML trainings"}; + Configurable downSampleBkgFactor{"downSampleBkgFactor", 1., "Fraction of background candidates to keep for ML trainings"}; + Configurable ptMaxForDownSample{"ptMaxForDownSample", 10., "Maximum pt for the application of the downsampling factor"}; + + using SelectedCandidates = soa::Filtered>; + using SelectedCandidatesKf = soa::Filtered>; + using SelectedCandidatesMc = soa::Filtered>; + using SelectedCandidatesKfMc = soa::Filtered>; + using TracksWPid = soa::Join; + + Filter filterSelectCandidates = aod::hf_sel_candidate_xic::isSelXicToXiPiPi >= selectionFlagXic; + + Partition recSig = nabs(aod::hf_cand_xic_to_xi_pi_pi::flagMcMatchRec) != int8_t(0); + Partition recBg = nabs(aod::hf_cand_xic_to_xi_pi_pi::flagMcMatchRec) == int8_t(0); + Partition recSigKf = nabs(aod::hf_cand_xic_to_xi_pi_pi::flagMcMatchRec) != int8_t(0); + Partition recBgKf = nabs(aod::hf_cand_xic_to_xi_pi_pi::flagMcMatchRec) == int8_t(0); + + void init(InitContext const&) + { + } + + template + void fillEvent(const T& collision) + { + rowCandidateFullEvents( + collision.bcId(), + collision.numContrib(), + collision.posX(), + collision.posY(), + collision.posZ()); + } + + template + void fillIndexTable(const T& candidate) + { + rowCandidateDauIndices( + candidate.cascadeId(), + candidate.pi0Id(), + candidate.pi1Id()); + } + + template + void fillCandidateTable(const T& candidate) + { + int8_t flagMc = 0; + if constexpr (doMc) { + flagMc = candidate.flagMcMatchRec(); + } + if constexpr (!doKf) { + if (fillCandidateLiteTable) { + rowCandidateLite( + candidate.isSelXicToXiPiPi(), + candidate.posX(), + candidate.posY(), + candidate.posZ(), + candidate.xSecondaryVertex(), + candidate.ySecondaryVertex(), + candidate.zSecondaryVertex(), + candidate.chi2PCA(), + candidate.sign(), + candidate.e(o2::constants::physics::MassXiCPlus), + candidate.invMassXic(), + candidate.p(), + candidate.pt(), + candidate.y(o2::constants::physics::MassXiCPlus), + candidate.eta(), + candidate.phi(), + candidate.ct(o2::constants::physics::MassXiCPlus), + candidate.decayLength(), + candidate.decayLengthXY(), + candidate.cpa(), + candidate.cpaXY(), + candidate.ptProng0(), + candidate.ptProng1(), + candidate.ptProng2(), + candidate.impactParameter0(), + candidate.impactParameter1(), + candidate.impactParameter2(), + candidate.cosPaXi(), + candidate.cosPaXYXi(), + candidate.cosPaLambda(), + candidate.cosPaXYLambda(), + flagMc); + } else { + rowCandidateFull( + candidate.isSelXicToXiPiPi(), + candidate.posX(), + candidate.posY(), + candidate.posZ(), + candidate.xPvErr(), + candidate.yPvErr(), + candidate.zPvErr(), + candidate.xSecondaryVertex(), + candidate.ySecondaryVertex(), + candidate.zSecondaryVertex(), + candidate.chi2PCA(), + candidate.xSvErr(), + candidate.ySvErr(), + candidate.zSvErr(), + candidate.xDecayVtxXi(), + candidate.yDecayVtxXi(), + candidate.zDecayVtxXi(), + candidate.xDecayVtxLambda(), + candidate.yDecayVtxLambda(), + candidate.zDecayVtxLambda(), + candidate.sign(), + candidate.e(o2::constants::physics::MassXiCPlus), + candidate.invMassXic(), + candidate.p(), + candidate.pt(), + candidate.y(o2::constants::physics::MassXiCPlus), + candidate.eta(), + candidate.phi(), + candidate.ct(o2::constants::physics::MassXiCPlus), + candidate.decayLength(), + candidate.decayLengthNormalised(), + candidate.decayLengthXY(), + candidate.decayLengthXYNormalised(), + candidate.cpa(), + candidate.cpaXY(), + candidate.ptProng0(), + candidate.ptProng1(), + candidate.ptProng2(), + candidate.impactParameter0(), + candidate.impactParameterNormalised0(), + candidate.impactParameter1(), + candidate.impactParameterNormalised1(), + candidate.impactParameter2(), + candidate.impactParameterNormalised2(), + candidate.maxNormalisedDeltaIP(), + candidate.cosPaXi(), + candidate.cosPaXYXi(), + candidate.cosPaLambda(), + candidate.cosPaXYLambda(), + candidate.invMassXiPi0(), + candidate.invMassXiPi1(), + flagMc); + } + } else { + if (fillCandidateLiteTable) { + rowCandidateLiteKf( + candidate.isSelXicToXiPiPi(), + candidate.posX(), + candidate.posY(), + candidate.posZ(), + candidate.xSecondaryVertex(), + candidate.ySecondaryVertex(), + candidate.zSecondaryVertex(), + candidate.chi2PCA(), + candidate.sign(), + candidate.e(o2::constants::physics::MassXiCPlus), + candidate.invMassXic(), + candidate.p(), + candidate.pt(), + candidate.y(o2::constants::physics::MassXiCPlus), + candidate.eta(), + candidate.phi(), + candidate.ct(o2::constants::physics::MassXiCPlus), + candidate.decayLength(), + candidate.decayLengthXY(), + candidate.cpa(), + candidate.cpaXY(), + candidate.ptProng0(), + candidate.ptProng1(), + candidate.ptProng2(), + candidate.impactParameter0(), + candidate.impactParameter1(), + candidate.impactParameter2(), + candidate.cosPaXi(), + candidate.cosPaXYXi(), + candidate.cosPaLambda(), + candidate.cosPaXYLambda(), + candidate.kfCascadeChi2(), + candidate.kfV0Chi2(), + candidate.dcaPi0Pi1(), + candidate.dcaPi0Xi(), + candidate.dcaPi1Xi(), + flagMc); + } else { + rowCandidateFullKf( + candidate.isSelXicToXiPiPi(), + candidate.posX(), + candidate.posY(), + candidate.posZ(), + candidate.xPvErr(), + candidate.yPvErr(), + candidate.zPvErr(), + candidate.xSecondaryVertex(), + candidate.ySecondaryVertex(), + candidate.zSecondaryVertex(), + candidate.chi2PCA(), + candidate.xSvErr(), + candidate.ySvErr(), + candidate.zSvErr(), + candidate.xDecayVtxXi(), + candidate.yDecayVtxXi(), + candidate.zDecayVtxXi(), + candidate.xDecayVtxLambda(), + candidate.yDecayVtxLambda(), + candidate.zDecayVtxLambda(), + candidate.sign(), + candidate.e(o2::constants::physics::MassXiCPlus), + candidate.invMassXic(), + candidate.p(), + candidate.pt(), + candidate.y(o2::constants::physics::MassXiCPlus), + candidate.eta(), + candidate.phi(), + candidate.ct(o2::constants::physics::MassXiCPlus), + candidate.decayLength(), + candidate.decayLengthNormalised(), + candidate.decayLengthXY(), + candidate.decayLengthXYNormalised(), + candidate.cpa(), + candidate.cpaXY(), + candidate.ptProng0(), + candidate.ptProng1(), + candidate.ptProng2(), + candidate.impactParameter0(), + candidate.impactParameterNormalised0(), + candidate.impactParameter1(), + candidate.impactParameterNormalised1(), + candidate.impactParameter2(), + candidate.impactParameterNormalised2(), + candidate.maxNormalisedDeltaIP(), + candidate.cosPaXi(), + candidate.cosPaXYXi(), + candidate.cosPaLambda(), + candidate.cosPaXYLambda(), + candidate.invMassXiPi0(), + candidate.invMassXiPi1(), + candidate.kfCascadeChi2(), + candidate.kfV0Chi2(), + candidate.dcaPi0Pi1(), + candidate.dcaPi0Xi(), + candidate.dcaPi1Xi(), + flagMc); + } + } + } + + void processData(aod::Collisions const& collisions, + SelectedCandidates const& candidates, + TracksWPid const&) + { + // Filling event properties + rowCandidateFullEvents.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision); + } + + // Filling candidate properties + if (fillCandidateDauIndexTable) { + rowCandidateDauIndices.reserve(candidates.size()); + } + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(candidates.size()); + } else { + rowCandidateFull.reserve(candidates.size()); + } + for (const auto& candidate : candidates) { + if (fillOnlyBackground && downSampleBkgFactor < 1.) { + float pseudoRndm = candidate.ptProng1() * 1000. - (int64_t)(candidate.ptProng1() * 1000); + if (pseudoRndm >= downSampleBkgFactor && candidate.pt() < ptMaxForDownSample) { + continue; + } + } + fillCandidateTable(candidate); + if (fillCandidateDauIndexTable) { + fillIndexTable(candidate); + } + } + } + PROCESS_SWITCH(HfTreeCreatorXicToXiPiPi, processData, "Process data", true); + + void processDataKf(aod::Collisions const& collisions, + SelectedCandidatesKf const& candidates, + TracksWPid const&) + { + // Filling event properties + rowCandidateFullEvents.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision); + } + + // Filling candidate properties + if (fillCandidateDauIndexTable) { + rowCandidateDauIndices.reserve(candidates.size()); + } + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(candidates.size()); + } else { + rowCandidateFull.reserve(candidates.size()); + } + for (const auto& candidate : candidates) { + if (fillOnlyBackground && downSampleBkgFactor < 1.) { + float pseudoRndm = candidate.ptProng1() * 1000. - (int64_t)(candidate.ptProng1() * 1000); + if (pseudoRndm >= downSampleBkgFactor && candidate.pt() < ptMaxForDownSample) { + continue; + } + } + fillCandidateTable(candidate); + if (fillCandidateDauIndexTable) { + fillIndexTable(candidate); + } + } + } + PROCESS_SWITCH(HfTreeCreatorXicToXiPiPi, processDataKf, "Process data with KF Particle reconstruction", false); + + void processMc(aod::Collisions const& collisions, + aod::McCollisions const&, + SelectedCandidatesMc const& candidates, + soa::Join const& particles, + TracksWPid const&) + { + // Filling event properties + rowCandidateFullEvents.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision); + } + + // Filling candidate properties + if (fillOnlySignal) { + if (fillCandidateDauIndexTable) { + rowCandidateDauIndices.reserve(candidates.size()); + } + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(recSig.size()); + } else { + rowCandidateFull.reserve(recSig.size()); + } + for (const auto& candidate : recSig) { + fillCandidateTable(candidate); + if (fillCandidateDauIndexTable) { + fillIndexTable(candidate); + } + } + } else if (fillOnlyBackground) { + if (fillCandidateDauIndexTable) { + rowCandidateDauIndices.reserve(candidates.size()); + } + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(recBg.size()); + } else { + rowCandidateFull.reserve(recBg.size()); + } + for (const auto& candidate : recBg) { + float pseudoRndm = candidate.ptProng1() * 1000. - (int64_t)(candidate.ptProng1() * 1000); + if (candidate.pt() < ptMaxForDownSample && pseudoRndm >= downSampleBkgFactor) { + continue; + } + fillCandidateTable(candidate); + if (fillCandidateDauIndexTable) { + fillIndexTable(candidate); + } + } + } else { + if (fillCandidateDauIndexTable) { + rowCandidateDauIndices.reserve(candidates.size()); + } + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(candidates.size()); + } else { + rowCandidateFull.reserve(candidates.size()); + } + for (const auto& candidate : candidates) { + fillCandidateTable(candidate); + if (fillCandidateDauIndexTable) { + fillIndexTable(candidate); + } + } + } + + // Filling particle properties + rowCandidateFullParticles.reserve(particles.size()); + for (const auto& particle : particles) { + if (TESTBIT(std::abs(particle.flagMcMatchGen()), aod::hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiPiPi)) { + rowCandidateFullParticles( + particle.mcCollision().bcId(), + particle.pt(), + particle.eta(), + particle.phi(), + RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, o2::constants::physics::MassXiCPlus), + particle.flagMcMatchGen()); + } + } + } + PROCESS_SWITCH(HfTreeCreatorXicToXiPiPi, processMc, "Process MC", false); + + void processMcKf(aod::Collisions const& collisions, + aod::McCollisions const&, + SelectedCandidatesKfMc const& candidates, + soa::Join const& particles, + TracksWPid const&) + { + // Filling event properties + rowCandidateFullEvents.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision); + } + + // Filling candidate properties + if (fillOnlySignal) { + if (fillCandidateDauIndexTable) { + rowCandidateDauIndices.reserve(candidates.size()); + } + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(recSigKf.size()); + } else { + rowCandidateFull.reserve(recSigKf.size()); + } + for (const auto& candidate : recSigKf) { + fillCandidateTable(candidate); + if (fillCandidateDauIndexTable) { + fillIndexTable(candidate); + } + } + } else if (fillOnlyBackground) { + if (fillCandidateDauIndexTable) { + rowCandidateDauIndices.reserve(candidates.size()); + } + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(recBgKf.size()); + } else { + rowCandidateFull.reserve(recBgKf.size()); + } + for (const auto& candidate : recBgKf) { + float pseudoRndm = candidate.ptProng1() * 1000. - (int64_t)(candidate.ptProng1() * 1000); + if (candidate.pt() < ptMaxForDownSample && pseudoRndm >= downSampleBkgFactor) { + continue; + } + fillCandidateTable(candidate); + if (fillCandidateDauIndexTable) { + fillIndexTable(candidate); + } + } + } else { + if (fillCandidateDauIndexTable) { + rowCandidateDauIndices.reserve(candidates.size()); + } + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(candidates.size()); + } else { + rowCandidateFull.reserve(candidates.size()); + } + for (const auto& candidate : candidates) { + fillCandidateTable(candidate); + if (fillCandidateDauIndexTable) { + fillIndexTable(candidate); + } + } + } + + // Filling particle properties + rowCandidateFullParticles.reserve(particles.size()); + for (const auto& particle : particles) { + if (TESTBIT(std::abs(particle.flagMcMatchGen()), aod::hf_cand_xic_to_xi_pi_pi::DecayType::XicToXiPiPi)) { + rowCandidateFullParticles( + particle.mcCollision().bcId(), + particle.pt(), + particle.eta(), + particle.phi(), + RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, o2::constants::physics::MassXiCPlus), + particle.flagMcMatchGen()); + } + } + } + PROCESS_SWITCH(HfTreeCreatorXicToXiPiPi, processMcKf, "Process MC with KF Particle reconstruction", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}