From ddbcbba460320df80c9fae14da06ae3dd1d49e76 Mon Sep 17 00:00:00 2001 From: nzardosh Date: Wed, 22 Jan 2025 22:29:49 +0100 Subject: [PATCH] [PWGHF] Add the possibility to run MC Gen only for HF derived data (#9357) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Vít Kučera --- PWGHF/TableProducer/CMakeLists.txt | 5 + .../TableProducer/candidateCreator2Prong.cxx | 43 +--- .../TableProducer/candidateCreator3Prong.cxx | 92 +------ PWGHF/TableProducer/candidateCreatorB0.cxx | 16 +- PWGHF/TableProducer/candidateCreatorBplus.cxx | 28 +- PWGHF/TableProducer/candidateCreatorMcGen.cxx | 76 ++++++ .../derivedDataCreatorBplusToD0Pi.cxx | 7 + .../derivedDataCreatorD0ToKPi.cxx | 7 + .../derivedDataCreatorLcToPKPi.cxx | 7 + PWGHF/Utils/utilsMcGen.h | 239 ++++++++++++++++++ 10 files changed, 349 insertions(+), 171 deletions(-) create mode 100644 PWGHF/TableProducer/candidateCreatorMcGen.cxx create mode 100644 PWGHF/Utils/utilsMcGen.h diff --git a/PWGHF/TableProducer/CMakeLists.txt b/PWGHF/TableProducer/CMakeLists.txt index 0b1e485e837..550343621db 100644 --- a/PWGHF/TableProducer/CMakeLists.txt +++ b/PWGHF/TableProducer/CMakeLists.txt @@ -100,6 +100,11 @@ o2physics_add_dpl_workflow(candidate-creator-xicc PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(candidate-creator-mc-gen + SOURCES candidateCreatorMcGen.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + # Candidate selectors o2physics_add_dpl_workflow(candidate-selector-b0-to-d-pi diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index 7f148c31b2b..008a7132c8a 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -49,6 +49,7 @@ #include "PWGHF/Utils/utilsEvSelHf.h" #include "PWGHF/Utils/utilsPid.h" #include "PWGHF/Utils/utilsTrkCandHf.h" +#include "PWGHF/Utils/utilsMcGen.h" using namespace o2; using namespace o2::analysis; @@ -839,47 +840,7 @@ struct HfCandidateCreator2ProngExpressions { } continue; } - - // Match generated particles. - for (const auto& particle : mcParticlesPerMcColl) { - flag = 0; - origin = 0; - std::vector idxBhadMothers{}; - // Reject particles from background events - if (particle.fromBackgroundEvent() && rejectBackground) { - rowMcMatchGen(flag, origin, -1); - continue; - } - - // D0(bar) → π± K∓ - if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign)) { - flag = sign * (1 << DecayType::D0ToPiK); - } - - // J/ψ → e+ e− - if (flag == 0) { - if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kJPsi, std::array{+kElectron, -kElectron}, true)) { - flag = 1 << DecayType::JpsiToEE; - } - } - - // J/ψ → μ+ μ− - if (flag == 0) { - if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kJPsi, std::array{+kMuonPlus, -kMuonPlus}, true)) { - flag = 1 << DecayType::JpsiToMuMu; - } - } - - // Check whether the particle is non-prompt (from a b quark). - if (flag != 0) { - origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, false, &idxBhadMothers); - } - if (origin == RecoDecay::OriginType::NonPrompt) { - rowMcMatchGen(flag, origin, idxBhadMothers[0]); - } else { - rowMcMatchGen(flag, origin, -1); - } - } + hf_mc_gen::fillMcMatchGen2Prong(mcParticlesPerMcColl, rowMcMatchGen, rejectBackground); } } diff --git a/PWGHF/TableProducer/candidateCreator3Prong.cxx b/PWGHF/TableProducer/candidateCreator3Prong.cxx index 8a4ae2396de..5228b106098 100644 --- a/PWGHF/TableProducer/candidateCreator3Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator3Prong.cxx @@ -48,6 +48,7 @@ #include "PWGHF/Utils/utilsBfieldCCDB.h" #include "PWGHF/Utils/utilsEvSelHf.h" #include "PWGHF/Utils/utilsTrkCandHf.h" +#include "PWGHF/Utils/utilsMcGen.h" using namespace o2; using namespace o2::analysis; @@ -1048,96 +1049,7 @@ struct HfCandidateCreator3ProngExpressions { } continue; } - - // Match generated particles. - for (const auto& particle : mcParticlesPerMcColl) { - flag = 0; - origin = 0; - channel = 0; - arrDaughIndex.clear(); - std::vector idxBhadMothers{}; - // Reject particles from background events - if (particle.fromBackgroundEvent() && rejectBackground) { - rowMcMatchGen(flag, origin, channel, -1); - continue; - } - - // D± → π± K∓ π± - if (createDplus) { - if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kDPlus, std::array{+kPiPlus, -kKPlus, +kPiPlus}, true, &sign, 2)) { - flag = sign * (1 << DecayType::DplusToPiKPi); - } - } - - // Ds± → K± K∓ π± and D± → K± K∓ π± - if (flag == 0 && createDs) { - bool isDplus = false; - if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kDS, std::array{+kKPlus, -kKPlus, +kPiPlus}, true, &sign, 2)) { - // DecayType::DsToKKPi is used to flag both Ds± → K± K∓ π± and D± → K± K∓ π± - // TODO: move to different and explicit flags - flag = sign * (1 << DecayType::DsToKKPi); - } else if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kDPlus, std::array{+kKPlus, -kKPlus, +kPiPlus}, true, &sign, 2)) { - // DecayType::DsToKKPi is used to flag both Ds± → K± K∓ π± and D± → K± K∓ π± - // TODO: move to different and explicit flags - flag = sign * (1 << DecayType::DsToKKPi); - isDplus = true; - } - if (flag != 0) { - RecoDecay::getDaughters(particle, &arrDaughIndex, std::array{0}, 1); - if (arrDaughIndex.size() == 2) { - for (auto jProng = 0u; jProng < arrDaughIndex.size(); ++jProng) { - auto daughJ = mcParticles.rawIteratorAt(arrDaughIndex[jProng]); - arrPDGDaugh[jProng] = std::abs(daughJ.pdgCode()); - } - if ((arrPDGDaugh[0] == arrPDGResonantDPhiPi[0] && arrPDGDaugh[1] == arrPDGResonantDPhiPi[1]) || (arrPDGDaugh[0] == arrPDGResonantDPhiPi[1] && arrPDGDaugh[1] == arrPDGResonantDPhiPi[0])) { - channel = isDplus ? DecayChannelDToKKPi::DplusToPhiPi : DecayChannelDToKKPi::DsToPhiPi; - } else if ((arrPDGDaugh[0] == arrPDGResonantDKstarK[0] && arrPDGDaugh[1] == arrPDGResonantDKstarK[1]) || (arrPDGDaugh[0] == arrPDGResonantDKstarK[1] && arrPDGDaugh[1] == arrPDGResonantDKstarK[0])) { - channel = isDplus ? DecayChannelDToKKPi::DplusToK0starK : DecayChannelDToKKPi::DsToK0starK; - } - } - } - } - - // Λc± → p± K∓ π± - if (flag == 0 && createLc) { - if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kLambdaCPlus, std::array{+kProton, -kKPlus, +kPiPlus}, true, &sign, 2)) { - flag = sign * (1 << DecayType::LcToPKPi); - - // Flagging the different Λc± → p± K∓ π± decay channels - RecoDecay::getDaughters(particle, &arrDaughIndex, std::array{0}, 1); - if (arrDaughIndex.size() == 2) { - for (auto jProng = 0u; jProng < arrDaughIndex.size(); ++jProng) { - auto daughJ = mcParticles.rawIteratorAt(arrDaughIndex[jProng]); - arrPDGDaugh[jProng] = std::abs(daughJ.pdgCode()); - } - if ((arrPDGDaugh[0] == arrPDGResonant1[0] && arrPDGDaugh[1] == arrPDGResonant1[1]) || (arrPDGDaugh[0] == arrPDGResonant1[1] && arrPDGDaugh[1] == arrPDGResonant1[0])) { - channel = 1; - } else if ((arrPDGDaugh[0] == arrPDGResonant2[0] && arrPDGDaugh[1] == arrPDGResonant2[1]) || (arrPDGDaugh[0] == arrPDGResonant2[1] && arrPDGDaugh[1] == arrPDGResonant2[0])) { - channel = 2; - } else if ((arrPDGDaugh[0] == arrPDGResonant3[0] && arrPDGDaugh[1] == arrPDGResonant3[1]) || (arrPDGDaugh[0] == arrPDGResonant3[1] && arrPDGDaugh[1] == arrPDGResonant3[0])) { - channel = 3; - } - } - } - } - - // Ξc± → p± K∓ π± - if (flag == 0 && createXic) { - if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kXiCPlus, std::array{+kProton, -kKPlus, +kPiPlus}, true, &sign, 2)) { - flag = sign * (1 << DecayType::XicToPKPi); - } - } - - // Check whether the particle is non-prompt (from a b quark). - if (flag != 0) { - origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, false, &idxBhadMothers); - } - if (origin == RecoDecay::OriginType::NonPrompt) { - rowMcMatchGen(flag, origin, channel, idxBhadMothers[0]); - } else { - rowMcMatchGen(flag, origin, channel, -1); - } - } + hf_mc_gen::fillMcMatchGen3Prong(mcParticlesPerMcColl, rowMcMatchGen, rejectBackground, createDplus, createDs, createLc, createXic); } } diff --git a/PWGHF/TableProducer/candidateCreatorB0.cxx b/PWGHF/TableProducer/candidateCreatorB0.cxx index d82b1d31b38..27779e37ab5 100644 --- a/PWGHF/TableProducer/candidateCreatorB0.cxx +++ b/PWGHF/TableProducer/candidateCreatorB0.cxx @@ -34,6 +34,7 @@ #include "PWGHF/DataModel/CandidateSelectionTables.h" #include "PWGHF/Utils/utilsBfieldCCDB.h" #include "PWGHF/Utils/utilsTrkCandHf.h" +#include "PWGHF/Utils/utilsMcGen.h" using namespace o2; using namespace o2::analysis; @@ -452,20 +453,7 @@ struct HfCandidateCreatorB0Expressions { rowMcMatchRec(flag, origin, debug); } // rec - // Match generated particles. - for (const auto& particle : mcParticles) { - flag = 0; - origin = 0; - // B0 → D- π+ - if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kB0, std::array{-static_cast(Pdg::kDPlus), +kPiPlus}, true)) { - // D- → π- K+ π- - auto candDMC = mcParticles.rawIteratorAt(particle.daughtersIds().front()); - if (RecoDecay::isMatchedMCGen(mcParticles, candDMC, -static_cast(Pdg::kDPlus), std::array{-kPiPlus, +kKPlus, -kPiPlus}, true, &sign)) { - flag = sign * BIT(hf_cand_b0::DecayType::B0ToDPi); - } - } - rowMcMatchGen(flag, origin); - } // gen + hf_mc_gen::fillMcMatchGenB0(mcParticles, rowMcMatchGen); // gen } // processMc PROCESS_SWITCH(HfCandidateCreatorB0Expressions, processMc, "Process MC", false); }; // struct diff --git a/PWGHF/TableProducer/candidateCreatorBplus.cxx b/PWGHF/TableProducer/candidateCreatorBplus.cxx index 0d4753c6c98..4d98f65cdaa 100644 --- a/PWGHF/TableProducer/candidateCreatorBplus.cxx +++ b/PWGHF/TableProducer/candidateCreatorBplus.cxx @@ -37,6 +37,7 @@ #include "PWGHF/DataModel/CandidateSelectionTables.h" #include "PWGHF/Utils/utilsBfieldCCDB.h" #include "PWGHF/Utils/utilsTrkCandHf.h" +#include "PWGHF/Utils/utilsMcGen.h" using namespace o2; using namespace o2::analysis; @@ -365,7 +366,6 @@ struct HfCandidateCreatorBplusExpressions { int8_t signB = 0, signD0 = 0; int8_t flag = 0; int8_t origin = 0; - int kD0pdg = Pdg::kD0; // Match reconstructed candidates. // Spawned table can be used directly @@ -386,31 +386,7 @@ struct HfCandidateCreatorBplusExpressions { } rowMcMatchRec(flag, origin); } - - // Match generated particles. - for (const auto& particle : mcParticles) { - flag = 0; - origin = 0; - signB = 0; - signD0 = 0; - int indexGenD0 = -1; - - // B± → D0bar(D0) π± → (K± π∓) π± - std::vector arrayDaughterB; - if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kBPlus, std::array{-kD0pdg, +kPiPlus}, true, &signB, 1, &arrayDaughterB)) { - // D0(bar) → π± K∓ - for (auto iD : arrayDaughterB) { - auto candDaughterMC = mcParticles.rawIteratorAt(iD); - if (std::abs(candDaughterMC.pdgCode()) == kD0pdg) { - indexGenD0 = RecoDecay::isMatchedMCGen(mcParticles, candDaughterMC, Pdg::kD0, std::array{-kKPlus, +kPiPlus}, true, &signD0, 1); - } - } - if (indexGenD0 > -1) { - flag = signB * (1 << hf_cand_bplus::DecayType::BplusToD0Pi); - } - } - rowMcMatchGen(flag, origin); - } // B candidate + hf_mc_gen::fillMcMatchGenBplus(mcParticles, rowMcMatchGen); // gen } // process PROCESS_SWITCH(HfCandidateCreatorBplusExpressions, processMc, "Process MC", false); }; // struct diff --git a/PWGHF/TableProducer/candidateCreatorMcGen.cxx b/PWGHF/TableProducer/candidateCreatorMcGen.cxx new file mode 100644 index 00000000000..a1ccbc8e920 --- /dev/null +++ b/PWGHF/TableProducer/candidateCreatorMcGen.cxx @@ -0,0 +1,76 @@ +// 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 candidateCreatorMcGen.cxx +/// \brief McGen only selection of heavy-flavour particles +/// +/// \author Nima Zardoshti, nima.zardoshti@cern.ch, CERN + +#include +#include +#include + +#include + +#include "CommonConstants/PhysicsConstants.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" +#include "Framework/RunningWorkflowInfo.h" + +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/Utils/utilsMcGen.h" + +using namespace o2; +using namespace o2::analysis; +using namespace o2::framework; + +/// Reconstruction of heavy-flavour 2-prong decay candidates +struct HfCandidateCreatorMcGen { + + Produces rowMcMatchGen2Prong; + Produces rowMcMatchGen3Prong; + Produces rowMcMatchGenBplus; + Produces rowMcMatchGenB0; + Configurable fill2Prong{"fill2Prong", false, "fill table for 2 prong candidates"}; + Configurable fill3Prong{"fill3Prong", false, "fill table for 3 prong candidates"}; + Configurable fillBplus{"fillBplus", false, "fill table for for B+ candidates"}; + Configurable fillB0{"fillB0", false, "fill table for B0 candidates"}; + Configurable rejectBackground2Prong{"rejectBackground2Prong", false, "Reject particles from PbPb background for 2 prong candidates"}; + Configurable rejectBackground3Prong{"rejectBackground3Prong", false, "Reject particles from PbPb background for 3 prong candidates"}; + Configurable createDplus{"createDplus", false, "Create D+ in 3 prong"}; + Configurable createDs{"createDs", false, "Create Ds in 3 prong"}; + Configurable createLc{"createLc", false, "Create Lc in 3 prong"}; + Configurable createXic{"createXic", false, "Create Xic in 3 prong"}; + + void process(aod::McCollision const&, + aod::McParticles const& mcParticles) + { + if (fill2Prong) { + hf_mc_gen::fillMcMatchGen2Prong(mcParticles, rowMcMatchGen2Prong, rejectBackground2Prong); + } + if (fill3Prong) { + hf_mc_gen::fillMcMatchGen3Prong(mcParticles, rowMcMatchGen3Prong, rejectBackground3Prong, createDplus, createDs, createLc, createXic); + } + if (fillBplus) { + hf_mc_gen::fillMcMatchGenBplus(mcParticles, rowMcMatchGenBplus); + } + if (fillB0) { + hf_mc_gen::fillMcMatchGenB0(mcParticles, rowMcMatchGenB0); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/TableProducer/derivedDataCreatorBplusToD0Pi.cxx b/PWGHF/TableProducer/derivedDataCreatorBplusToD0Pi.cxx index c1fc262b0af..bedffe215b0 100644 --- a/PWGHF/TableProducer/derivedDataCreatorBplusToD0Pi.cxx +++ b/PWGHF/TableProducer/derivedDataCreatorBplusToD0Pi.cxx @@ -419,6 +419,13 @@ struct HfDerivedDataCreatorBplusToD0Pi { rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass); } PROCESS_SWITCH(HfDerivedDataCreatorBplusToD0Pi, processMcMlAll, "Process MC with ML", false); + + void processMcGenOnly(TypeMcCollisions const& mcCollisions, + MatchedGenCandidatesMc const& mcParticles) + { + rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass); + } + PROCESS_SWITCH(HfDerivedDataCreatorBplusToD0Pi, processMcGenOnly, "Process MC gen. only", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/TableProducer/derivedDataCreatorD0ToKPi.cxx b/PWGHF/TableProducer/derivedDataCreatorD0ToKPi.cxx index 92fdd9089b8..83f15ec5e80 100644 --- a/PWGHF/TableProducer/derivedDataCreatorD0ToKPi.cxx +++ b/PWGHF/TableProducer/derivedDataCreatorD0ToKPi.cxx @@ -516,6 +516,13 @@ struct HfDerivedDataCreatorD0ToKPi { rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass); } PROCESS_SWITCH(HfDerivedDataCreatorD0ToKPi, processMcWithKFParticleMlAll, "Process MC with KFParticle and ML", false); + + void processMcGenOnly(TypeMcCollisions const& mcCollisions, + MatchedGenCandidatesMc const& mcParticles) + { + rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass); + } + PROCESS_SWITCH(HfDerivedDataCreatorD0ToKPi, processMcGenOnly, "Process MC gen. only", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/TableProducer/derivedDataCreatorLcToPKPi.cxx b/PWGHF/TableProducer/derivedDataCreatorLcToPKPi.cxx index c574611228d..9c9becea741 100644 --- a/PWGHF/TableProducer/derivedDataCreatorLcToPKPi.cxx +++ b/PWGHF/TableProducer/derivedDataCreatorLcToPKPi.cxx @@ -387,6 +387,13 @@ struct HfDerivedDataCreatorLcToPKPi { rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass); } PROCESS_SWITCH(HfDerivedDataCreatorLcToPKPi, processMcMlAll, "Process MC with ML", false); + + void processMcGenOnly(TypeMcCollisions const& mcCollisions, + MatchedGenCandidatesMc const& mcParticles) + { + rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass); + } + PROCESS_SWITCH(HfDerivedDataCreatorLcToPKPi, processMcGenOnly, "Process MC gen. only", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/Utils/utilsMcGen.h b/PWGHF/Utils/utilsMcGen.h new file mode 100644 index 00000000000..0129bf6d840 --- /dev/null +++ b/PWGHF/Utils/utilsMcGen.h @@ -0,0 +1,239 @@ +// 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 utilsMcGen.h +/// \brief utility functions for HF Mc gen. workflows +/// +/// \author Nima Zardoshti, nima.zardoshti@cern.ch, CERN + +#ifndef PWGHF_UTILS_UTILSMCGEN_H_ +#define PWGHF_UTILS_UTILSMCGEN_H_ + +#include + +#include + +#include "CommonConstants/PhysicsConstants.h" + +#include "Common/Core/RecoDecay.h" + +#include "PWGHF/DataModel/CandidateReconstructionTables.h" + +namespace hf_mc_gen +{ + +template +void fillMcMatchGen2Prong(T const& mcParticles, U& rowMcMatchGen, bool rejectBackground) +{ + using namespace o2::constants::physics; + + // Match generated particles. + for (const auto& particle : mcParticles) { + int8_t flag = 0; + int8_t origin = 0; + int8_t sign = 0; + std::vector idxBhadMothers{}; + // Reject particles from background events + if (particle.fromBackgroundEvent() && rejectBackground) { + rowMcMatchGen(flag, origin, -1); + continue; + } + + // D0(bar) → π± K∓ + if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign)) { + flag = sign * (1 << o2::aod::hf_cand_2prong::DecayType::D0ToPiK); + } + + // J/ψ → e+ e− + if (flag == 0) { + if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kJPsi, std::array{+kElectron, -kElectron}, true)) { + flag = 1 << o2::aod::hf_cand_2prong::DecayType::JpsiToEE; + } + } + + // J/ψ → μ+ μ− + if (flag == 0) { + if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kJPsi, std::array{+kMuonPlus, -kMuonPlus}, true)) { + flag = 1 << o2::aod::hf_cand_2prong::DecayType::JpsiToMuMu; + } + } + + // Check whether the particle is non-prompt (from a b quark). + if (flag != 0) { + origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, false, &idxBhadMothers); + } + if (origin == RecoDecay::OriginType::NonPrompt) { + rowMcMatchGen(flag, origin, idxBhadMothers[0]); + } else { + rowMcMatchGen(flag, origin, -1); + } + } +} + +template +void fillMcMatchGen3Prong(T const& mcParticles, U& rowMcMatchGen, bool rejectBackground, bool createDplus, bool createDs, bool createLc, bool createXic) +{ + using namespace o2::constants::physics; + + // Match generated particles. + for (const auto& particle : mcParticles) { + int8_t flag = 0; + int8_t origin = 0; + int8_t channel = 0; + int8_t sign = 0; + std::vector arrDaughIndex; + std::vector idxBhadMothers{}; + std::array arrPDGDaugh; + std::array arrPDGResonant1 = {kProton, 313}; // Λc± → p± K* + std::array arrPDGResonant2 = {2224, kKPlus}; // Λc± → Δ(1232)±± K∓ + std::array arrPDGResonant3 = {102134, kPiPlus}; // Λc± → Λ(1520) π± + std::array arrPDGResonantDPhiPi = {333, kPiPlus}; // Ds± → Phi π± and D± → Phi π± + std::array arrPDGResonantDKstarK = {313, kKPlus}; // Ds± → K*(892)0bar K± and D± → K*(892)0bar K± + // Reject particles from background events + if (particle.fromBackgroundEvent() && rejectBackground) { + rowMcMatchGen(flag, origin, channel, -1); + continue; + } + + // D± → π± K∓ π± + if (createDplus) { + if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kDPlus, std::array{+kPiPlus, -kKPlus, +kPiPlus}, true, &sign, 2)) { + flag = sign * (1 << o2::aod::hf_cand_3prong::DecayType::DplusToPiKPi); + } + } + + // Ds± → K± K∓ π± and D± → K± K∓ π± + if (flag == 0 && createDs) { + bool isDplus = false; + if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kDS, std::array{+kKPlus, -kKPlus, +kPiPlus}, true, &sign, 2)) { + // DecayType::DsToKKPi is used to flag both Ds± → K± K∓ π± and D± → K± K∓ π± + // TODO: move to different and explicit flags + flag = sign * (1 << o2::aod::hf_cand_3prong::DecayType::DsToKKPi); + } else if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kDPlus, std::array{+kKPlus, -kKPlus, +kPiPlus}, true, &sign, 2)) { + // DecayType::DsToKKPi is used to flag both Ds± → K± K∓ π± and D± → K± K∓ π± + // TODO: move to different and explicit flags + flag = sign * (1 << o2::aod::hf_cand_3prong::DecayType::DsToKKPi); + isDplus = true; + } + if (flag != 0) { + RecoDecay::getDaughters(particle, &arrDaughIndex, std::array{0}, 1); + if (arrDaughIndex.size() == 2) { + for (auto jProng = 0u; jProng < arrDaughIndex.size(); ++jProng) { + auto daughJ = mcParticles.rawIteratorAt(arrDaughIndex[jProng]); + arrPDGDaugh[jProng] = std::abs(daughJ.pdgCode()); + } + if ((arrPDGDaugh[0] == arrPDGResonantDPhiPi[0] && arrPDGDaugh[1] == arrPDGResonantDPhiPi[1]) || (arrPDGDaugh[0] == arrPDGResonantDPhiPi[1] && arrPDGDaugh[1] == arrPDGResonantDPhiPi[0])) { + channel = isDplus ? o2::aod::hf_cand_3prong::DecayChannelDToKKPi::DplusToPhiPi : o2::aod::hf_cand_3prong::DecayChannelDToKKPi::DsToPhiPi; + } else if ((arrPDGDaugh[0] == arrPDGResonantDKstarK[0] && arrPDGDaugh[1] == arrPDGResonantDKstarK[1]) || (arrPDGDaugh[0] == arrPDGResonantDKstarK[1] && arrPDGDaugh[1] == arrPDGResonantDKstarK[0])) { + channel = isDplus ? o2::aod::hf_cand_3prong::DecayChannelDToKKPi::DplusToK0starK : o2::aod::hf_cand_3prong::DecayChannelDToKKPi::DsToK0starK; + } + } + } + } + + // Λc± → p± K∓ π± + if (flag == 0 && createLc) { + if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kLambdaCPlus, std::array{+kProton, -kKPlus, +kPiPlus}, true, &sign, 2)) { + flag = sign * (1 << o2::aod::hf_cand_3prong::DecayType::LcToPKPi); + + // Flagging the different Λc± → p± K∓ π± decay channels + RecoDecay::getDaughters(particle, &arrDaughIndex, std::array{0}, 1); + if (arrDaughIndex.size() == 2) { + for (auto jProng = 0u; jProng < arrDaughIndex.size(); ++jProng) { + auto daughJ = mcParticles.rawIteratorAt(arrDaughIndex[jProng]); + arrPDGDaugh[jProng] = std::abs(daughJ.pdgCode()); + } + if ((arrPDGDaugh[0] == arrPDGResonant1[0] && arrPDGDaugh[1] == arrPDGResonant1[1]) || (arrPDGDaugh[0] == arrPDGResonant1[1] && arrPDGDaugh[1] == arrPDGResonant1[0])) { + channel = 1; + } else if ((arrPDGDaugh[0] == arrPDGResonant2[0] && arrPDGDaugh[1] == arrPDGResonant2[1]) || (arrPDGDaugh[0] == arrPDGResonant2[1] && arrPDGDaugh[1] == arrPDGResonant2[0])) { + channel = 2; + } else if ((arrPDGDaugh[0] == arrPDGResonant3[0] && arrPDGDaugh[1] == arrPDGResonant3[1]) || (arrPDGDaugh[0] == arrPDGResonant3[1] && arrPDGDaugh[1] == arrPDGResonant3[0])) { + channel = 3; + } + } + } + } + + // Ξc± → p± K∓ π± + if (flag == 0 && createXic) { + if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kXiCPlus, std::array{+kProton, -kKPlus, +kPiPlus}, true, &sign, 2)) { + flag = sign * (1 << o2::aod::hf_cand_3prong::DecayType::XicToPKPi); + } + } + + // Check whether the particle is non-prompt (from a b quark). + if (flag != 0) { + origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, false, &idxBhadMothers); + } + if (origin == RecoDecay::OriginType::NonPrompt) { + rowMcMatchGen(flag, origin, channel, idxBhadMothers[0]); + } else { + rowMcMatchGen(flag, origin, channel, -1); + } + } +} + +template +void fillMcMatchGenBplus(T const& mcParticles, U& rowMcMatchGen) +{ + using namespace o2::constants::physics; + + // Match generated particles. + for (const auto& particle : mcParticles) { + int8_t flag = 0; + int8_t origin = 0; + int8_t signB = 0; + int8_t signD0 = 0; + int indexGenD0 = -1; + + // B± → D0bar(D0) π± → (K± π∓) π± + std::vector arrayDaughterB; + if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kBPlus, std::array{-Pdg::kD0, +kPiPlus}, true, &signB, 1, &arrayDaughterB)) { + // D0(bar) → π± K∓ + for (const auto iD : arrayDaughterB) { + auto candDaughterMC = mcParticles.rawIteratorAt(iD); + if (std::abs(candDaughterMC.pdgCode()) == Pdg::kD0) { + indexGenD0 = RecoDecay::isMatchedMCGen(mcParticles, candDaughterMC, Pdg::kD0, std::array{-kKPlus, +kPiPlus}, true, &signD0, 1); + } + } + if (indexGenD0 > -1) { + flag = signB * (1 << o2::aod::hf_cand_bplus::DecayType::BplusToD0Pi); + } + } + rowMcMatchGen(flag, origin); + } // B candidate +} + +template +void fillMcMatchGenB0(T const& mcParticles, U& rowMcMatchGen) +{ + using namespace o2::constants::physics; + + // Match generated particles. + for (const auto& particle : mcParticles) { + int8_t flag = 0; + int8_t origin = 0; + int8_t sign = 0; + // B0 → D- π+ + if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kB0, std::array{-static_cast(Pdg::kDPlus), +kPiPlus}, true)) { + // D- → π- K+ π- + auto candDMC = mcParticles.rawIteratorAt(particle.daughtersIds().front()); + if (RecoDecay::isMatchedMCGen(mcParticles, candDMC, -static_cast(Pdg::kDPlus), std::array{-kPiPlus, +kKPlus, -kPiPlus}, true, &sign)) { + flag = sign * BIT(o2::aod::hf_cand_b0::DecayType::B0ToDPi); + } + } + rowMcMatchGen(flag, origin); + } // gen +} + +} // namespace hf_mc_gen + +#endif // PWGHF_UTILS_UTILSMCGEN_H_