Skip to content

Commit a832ad8

Browse files
Marcellocostialibuild
Marcellocosti
and
alibuild
authored
[PWGHF] Implement gen and reco MC matching for correlated bkg decays (#11418)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent e62eb26 commit a832ad8

File tree

9 files changed

+865
-291
lines changed

9 files changed

+865
-291
lines changed

‎PWGHF/Core/DecayChannels.h‎

Lines changed: 17 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,13 @@ enum DecayChannelMain : int8_t {
7171
DsToPiPiPi, // π+ π− π+
7272
DsToPiPiPiPi0, // π+ π− π+ π0
7373
// D*+
74-
DstarToPiKPi, // π+ K− π+ (from [(D0 → π+ K−) π+])
74+
DstarToPiKPi, // π+ K− π+ (from [(D0 → π+ K−) π+])
75+
DstarToPiKPiPi0, // π+ K− π+ π0
76+
DstarToPiKPiPi0Pi0, // π+ K− π+ π0 π0
77+
DstarToPiKK, // π+ K− K+
78+
DstarToPiKKPi0, // π+ K− K+ π0
79+
DstarToPiPiPi, // π+ π− π+
80+
DstarToPiPiPiPi0, // π+ π− π+ π0
7581
// Λc+
7682
LcToPKPi, // p K− π+
7783
LcToPKPiPi0, // p K− π+ π0
@@ -102,6 +108,16 @@ enum DecayChannelResonant : int8_t {
102108
DsToF2_1270Pi, // f2(1270) π+
103109
DsToF0_1370K, // f0(1370) K+
104110
DsToEtaPi, // η π+
111+
// D*+
112+
DstarToD0ToRhoplusPi, // ρ+ π−
113+
DstarToD0ToRhoplusK, // ρ+ K−
114+
DstarToD0ToKstar0Pi0, // anti-K*0 π0
115+
DstarToD0ToKstarPi, // K*− π+
116+
DstarToDplusToPhiPi, // φ π+
117+
DstarToDplusToKstar0K, // anti-K*0 K+
118+
DstarToDplusToKstar1430_0K, // anti-K*0(1430) K+
119+
DstarToDplusToRho0Pi, // ρ0 π+
120+
DstarToDplusToF2_1270Pi, // f2(1270) π+
105121
// Λc+
106122
LcToPKstar0, // p K*0(892)
107123
LcToDeltaplusplusK, // Δ++ K−

‎PWGHF/DataModel/CandidateReconstructionTables.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -641,6 +641,9 @@ DECLARE_SOA_COLUMN(FlagMcMatchRec, flagMcMatchRec, int8_t); //! reconstruction l
641641
DECLARE_SOA_COLUMN(FlagMcMatchGen, flagMcMatchGen, int8_t); //! generator level
642642
DECLARE_SOA_COLUMN(OriginMcRec, originMcRec, int8_t); //! particle origin, reconstruction level
643643
DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int8_t); //! particle origin, generator level
644+
DECLARE_SOA_COLUMN(FlagMcDecayChanRec, flagMcDecayChanRec, int8_t); //! resonant decay channel flag, reconstruction level
645+
DECLARE_SOA_COLUMN(FlagMcDecayChanGen, flagMcDecayChanGen, int8_t); //! resonant decay channel flag, reconstruction level
646+
644647
// KF related properties
645648
DECLARE_SOA_COLUMN(KfGeoMassD0, kfGeoMassD0, float); //! mass of the D0 candidate from the KFParticle geometric fit
646649
DECLARE_SOA_COLUMN(KfGeoMassD0bar, kfGeoMassD0bar, float); //! mass of the D0bar candidate from the KFParticle geometric fit
@@ -754,6 +757,7 @@ DECLARE_SOA_TABLE(HfCand2ProngKF, "AOD", "HFCAND2PKF",
754757
DECLARE_SOA_TABLE(HfCand2ProngMcRec, "AOD", "HFCAND2PMCREC", //!
755758
hf_cand_2prong::FlagMcMatchRec,
756759
hf_cand_2prong::OriginMcRec,
760+
hf_cand_2prong::FlagMcDecayChanRec,
757761
hf_cand::PtBhadMotherPart,
758762
hf_cand::PdgBhadMotherPart,
759763
hf_cand::NTracksDecayed,
@@ -763,6 +767,7 @@ DECLARE_SOA_TABLE(HfCand2ProngMcRec, "AOD", "HFCAND2PMCREC", //!
763767
DECLARE_SOA_TABLE(HfCand2ProngMcGen, "AOD", "HFCAND2PMCGEN", //!
764768
hf_cand_2prong::FlagMcMatchGen,
765769
hf_cand_2prong::OriginMcGen,
770+
hf_cand_2prong::FlagMcDecayChanGen,
766771
hf_cand::IdxBhadMotherPart);
767772

768773
// cascade decay candidate table

‎PWGHF/TableProducer/candidateCreator2Prong.cxx‎

Lines changed: 122 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -20,45 +20,48 @@
2020
#define HomogeneousField // o2-linter: disable=name/macro (required by KFParticle)
2121
#endif
2222

23-
#include <memory>
24-
#include <string>
25-
#include <vector>
26-
27-
#include <KFParticleBase.h>
28-
#include <KFParticle.h>
29-
#include <KFPTrack.h>
30-
#include <KFPVertex.h>
31-
#include <KFVertex.h>
23+
#include "PWGHF/Core/CentralityEstimation.h"
24+
#include "PWGHF/Core/DecayChannels.h"
25+
#include "PWGHF/Core/SelectorCuts.h"
26+
#include "PWGHF/DataModel/CandidateReconstructionTables.h"
27+
#include "PWGHF/Utils/utilsBfieldCCDB.h"
28+
#include "PWGHF/Utils/utilsEvSelHf.h"
29+
#include "PWGHF/Utils/utilsMcGen.h"
30+
#include "PWGHF/Utils/utilsMcMatching.h"
31+
#include "PWGHF/Utils/utilsPid.h"
32+
#include "PWGHF/Utils/utilsTrkCandHf.h"
33+
#include "PWGLF/DataModel/mcCentrality.h"
3234

33-
#include <TPDGCode.h>
35+
#include "Common/Core/trackUtilities.h"
36+
#include "Tools/KFparticle/KFUtilities.h"
3437

3538
#include "CommonConstants/PhysicsConstants.h"
3639
#include "DCAFitter/DCAFitterN.h"
3740
#include "Framework/AnalysisTask.h"
3841
#include "Framework/HistogramRegistry.h"
39-
#include "Framework/runDataProcessing.h"
4042
#include "Framework/RunningWorkflowInfo.h"
43+
#include "Framework/runDataProcessing.h"
4144
#include "ReconstructionDataFormats/DCA.h"
4245

43-
#include "Common/Core/trackUtilities.h"
44-
#include "Tools/KFparticle/KFUtilities.h"
46+
#include <TPDGCode.h>
4547

46-
#include "PWGLF/DataModel/mcCentrality.h"
48+
#include <KFPTrack.h>
49+
#include <KFPVertex.h>
50+
#include <KFParticle.h>
51+
#include <KFParticleBase.h>
52+
#include <KFVertex.h>
4753

48-
#include "PWGHF/Core/CentralityEstimation.h"
49-
#include "PWGHF/Core/SelectorCuts.h"
50-
#include "PWGHF/DataModel/CandidateReconstructionTables.h"
51-
#include "PWGHF/Utils/utilsBfieldCCDB.h"
52-
#include "PWGHF/Utils/utilsEvSelHf.h"
53-
#include "PWGHF/Utils/utilsMcGen.h"
54-
#include "PWGHF/Utils/utilsPid.h"
55-
#include "PWGHF/Utils/utilsTrkCandHf.h"
54+
#include <memory>
55+
#include <string>
56+
#include <vector>
5657

5758
using namespace o2;
5859
using namespace o2::analysis;
5960
using namespace o2::hf_evsel;
6061
using namespace o2::hf_trkcandsel;
6162
using namespace o2::aod::hf_cand_2prong;
63+
using namespace o2::hf_decay;
64+
using namespace o2::hf_decay::hf_cand_2prong;
6265
using namespace o2::hf_centrality;
6366
using namespace o2::hf_occupancy;
6467
using namespace o2::constants::physics;
@@ -691,6 +694,7 @@ struct HfCandidateCreator2ProngExpressions {
691694
Configurable<bool> rejectBackground{"rejectBackground", true, "Reject particles from background events"};
692695
Configurable<bool> matchKinkedDecayTopology{"matchKinkedDecayTopology", false, "Match also candidates with tracks that decay with kinked topology"};
693696
Configurable<bool> matchInteractionsWithMaterial{"matchInteractionsWithMaterial", false, "Match also candidates with tracks that interact with material"};
697+
Configurable<bool> matchCorrelatedBackgrounds{"matchCorrelatedBackgrounds", false, "Match correlated background candidates"};
694698

695699
HfEventSelectionMc hfEvSelMc; // mc event selection and monitoring
696700

@@ -738,15 +742,18 @@ struct HfCandidateCreator2ProngExpressions {
738742
int indexRec = -1;
739743
int8_t sign = 0;
740744
int8_t flag = 0;
745+
int8_t channel = 0;
741746
int8_t origin = 0;
742747
int8_t nKinkedTracks = 0;
743748
int8_t nInteractionsWithMaterial = 0;
749+
constexpr std::size_t NDaughtersResonant{2u};
744750

745751
// Match reconstructed candidates.
746752
// Spawned table can be used directly
747753
for (const auto& candidate : *rowCandidateProng2) {
748754
flag = 0;
749755
origin = 0;
756+
channel = 0;
750757
auto arrayDaughters = std::array{candidate.prong0_as<aod::TracksWMc>(), candidate.prong1_as<aod::TracksWMc>()};
751758

752759
// Check whether the particle is from background events. If so, reject it.
@@ -762,47 +769,107 @@ struct HfCandidateCreator2ProngExpressions {
762769
}
763770
}
764771
if (fromBkg) {
765-
rowMcMatchRec(flag, origin, -1.f, 0, 0, 0);
772+
rowMcMatchRec(flag, origin, channel, -1.f, 0, 0, 0);
766773
continue;
767774
}
768775
}
769776
std::vector<int> idxBhadMothers{};
770777

771-
// D0(bar) → π± K∓
772-
if (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
773-
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(mcParticles, arrayDaughters, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign, 1, &nKinkedTracks, &nInteractionsWithMaterial);
774-
} else if (matchKinkedDecayTopology && !matchInteractionsWithMaterial) {
775-
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, false>(mcParticles, arrayDaughters, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign, 1, &nKinkedTracks);
776-
} else if (!matchKinkedDecayTopology && matchInteractionsWithMaterial) {
777-
indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, true>(mcParticles, arrayDaughters, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign, 1, nullptr, &nInteractionsWithMaterial);
778-
} else {
779-
indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign);
780-
}
781-
if (indexRec > -1) {
782-
flag = sign * (1 << DecayType::D0ToPiK);
783-
}
778+
if (matchCorrelatedBackgrounds) {
779+
indexRec = -1; // Index of the matched reconstructed candidate
780+
constexpr int FinalStateDepth = 2;
781+
constexpr int ResoDepth = 1;
782+
783+
// D0(bar) → π+ K−, π+ K− π0, π+ π−, π+ π− π0, K+ K−
784+
for (const auto& [chn, finalState] : hf_cand_2prong::daughtersD0Main) {
785+
std::array<int, 2> finalStateParts2Prong = std::array{finalState[0], finalState[1]};
786+
if (finalState.size() == 3) { // o2-linter: disable=magic-number (Partly Reco 3-prong decays)
787+
if (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
788+
indexRec = RecoDecay::getMatchedMCRec<false, false, true, true, true>(mcParticles, arrayDaughters, Pdg::kD0, finalStateParts2Prong, true, &sign, FinalStateDepth, &nKinkedTracks, &nInteractionsWithMaterial);
789+
} else if (matchKinkedDecayTopology && !matchInteractionsWithMaterial) {
790+
indexRec = RecoDecay::getMatchedMCRec<false, false, true, true, false>(mcParticles, arrayDaughters, Pdg::kD0, finalStateParts2Prong, true, &sign, FinalStateDepth, &nKinkedTracks);
791+
} else if (!matchKinkedDecayTopology && matchInteractionsWithMaterial) {
792+
indexRec = RecoDecay::getMatchedMCRec<false, false, true, false, true>(mcParticles, arrayDaughters, Pdg::kD0, finalStateParts2Prong, true, &sign, FinalStateDepth, nullptr, &nInteractionsWithMaterial);
793+
} else {
794+
indexRec = RecoDecay::getMatchedMCRec<false, false, true, false, false>(mcParticles, arrayDaughters, Pdg::kD0, finalStateParts2Prong, true, &sign, FinalStateDepth);
795+
}
784796

785-
// J/ψ → e+ e−
786-
if (flag == 0) {
787-
if (matchInteractionsWithMaterial) {
788-
indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, true>(mcParticles, arrayDaughters, Pdg::kJPsi, std::array{+kElectron, -kElectron}, true, &sign, 1, nullptr, &nInteractionsWithMaterial);
797+
if (indexRec > -1) {
798+
auto motherParticle = mcParticles.rawIteratorAt(indexRec);
799+
std::array<int, 3> finalStateParts2ProngAll = std::array{finalState[0], finalState[1], finalState[2]};
800+
changeFinalStatePdgSign(motherParticle.pdgCode(), +kPi0, finalStateParts2ProngAll);
801+
if (!RecoDecay::isMatchedMCGen(mcParticles, motherParticle, Pdg::kD0, finalStateParts2ProngAll, true, &sign, FinalStateDepth)) {
802+
indexRec = -1; // Reset indexRec if the generated decay does not match the reconstructed one does not match the reconstructed one
803+
}
804+
}
805+
} else if (finalState.size() == 2) { // o2-linter: disable=magic-number (Fully Reco 2-prong decays)
806+
if (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
807+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(mcParticles, arrayDaughters, Pdg::kD0, finalStateParts2Prong, true, &sign, FinalStateDepth, &nKinkedTracks, &nInteractionsWithMaterial);
808+
} else if (matchKinkedDecayTopology && !matchInteractionsWithMaterial) {
809+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, false>(mcParticles, arrayDaughters, Pdg::kD0, finalStateParts2Prong, true, &sign, FinalStateDepth, &nKinkedTracks);
810+
} else if (!matchKinkedDecayTopology && matchInteractionsWithMaterial) {
811+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, true>(mcParticles, arrayDaughters, Pdg::kD0, finalStateParts2Prong, true, &sign, FinalStateDepth, nullptr, &nInteractionsWithMaterial);
812+
} else {
813+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, false>(mcParticles, arrayDaughters, Pdg::kD0, finalStateParts2Prong, true, &sign, FinalStateDepth);
814+
}
815+
} else {
816+
LOG(fatal) << "Final state size not supported: " << finalState.size();
817+
continue;
818+
}
819+
if (indexRec > -1) {
820+
flag = sign * (1 << chn);
821+
822+
// Flag the resonant decay channel
823+
std::vector<int> arrResoDaughIndex = {};
824+
RecoDecay::getDaughters(mcParticles.rawIteratorAt(indexRec), &arrResoDaughIndex, std::array{0}, ResoDepth);
825+
std::array<int, NDaughtersResonant> arrPDGDaugh = {};
826+
if (arrResoDaughIndex.size() == NDaughtersResonant) {
827+
for (auto iProng = 0u; iProng < arrResoDaughIndex.size(); ++iProng) {
828+
auto daughI = mcParticles.rawIteratorAt(arrResoDaughIndex[iProng]);
829+
arrPDGDaugh[iProng] = daughI.pdgCode();
830+
}
831+
channel = flagResonantDecay(Pdg::kD0, arrPDGDaugh);
832+
}
833+
break;
834+
}
835+
}
836+
} else {
837+
// D0(bar) → π± K∓
838+
if (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
839+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(mcParticles, arrayDaughters, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign, 1, &nKinkedTracks, &nInteractionsWithMaterial);
840+
} else if (matchKinkedDecayTopology && !matchInteractionsWithMaterial) {
841+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, false>(mcParticles, arrayDaughters, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign, 1, &nKinkedTracks);
842+
} else if (!matchKinkedDecayTopology && matchInteractionsWithMaterial) {
843+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, true>(mcParticles, arrayDaughters, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign, 1, nullptr, &nInteractionsWithMaterial);
789844
} else {
790-
indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, Pdg::kJPsi, std::array{+kElectron, -kElectron}, true);
845+
indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign);
791846
}
792847
if (indexRec > -1) {
793-
flag = 1 << DecayType::JpsiToEE;
848+
flag = sign * (1 << DecayType::D0ToPiK);
794849
}
795-
}
796850

797-
// J/ψ → μ+ μ−
798-
if (flag == 0) {
799-
if (matchInteractionsWithMaterial) {
800-
indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, true>(mcParticles, arrayDaughters, Pdg::kJPsi, std::array{+kMuonPlus, -kMuonPlus}, true, &sign, 1, nullptr, &nInteractionsWithMaterial);
801-
} else {
802-
indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, Pdg::kJPsi, std::array{+kMuonPlus, -kMuonPlus}, true);
851+
// J/ψ → e+ e−
852+
if (flag == 0) {
853+
if (matchInteractionsWithMaterial) {
854+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, true>(mcParticles, arrayDaughters, Pdg::kJPsi, std::array{+kElectron, -kElectron}, true, &sign, 1, nullptr, &nInteractionsWithMaterial);
855+
} else {
856+
indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, Pdg::kJPsi, std::array{+kElectron, -kElectron}, true);
857+
}
858+
if (indexRec > -1) {
859+
flag = 1 << DecayType::JpsiToEE;
860+
}
803861
}
804-
if (indexRec > -1) {
805-
flag = 1 << DecayType::JpsiToMuMu;
862+
863+
// J/ψ → μ+ μ−
864+
if (flag == 0) {
865+
if (matchInteractionsWithMaterial) {
866+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, true>(mcParticles, arrayDaughters, Pdg::kJPsi, std::array{+kMuonPlus, -kMuonPlus}, true, &sign, 1, nullptr, &nInteractionsWithMaterial);
867+
} else {
868+
indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, Pdg::kJPsi, std::array{+kMuonPlus, -kMuonPlus}, true);
869+
}
870+
if (indexRec > -1) {
871+
flag = 1 << DecayType::JpsiToMuMu;
872+
}
806873
}
807874
}
808875

@@ -813,9 +880,9 @@ struct HfCandidateCreator2ProngExpressions {
813880
}
814881
if (origin == RecoDecay::OriginType::NonPrompt) {
815882
auto bHadMother = mcParticles.rawIteratorAt(idxBhadMothers[0]);
816-
rowMcMatchRec(flag, origin, bHadMother.pt(), bHadMother.pdgCode(), nKinkedTracks, nInteractionsWithMaterial);
883+
rowMcMatchRec(flag, origin, channel, bHadMother.pt(), bHadMother.pdgCode(), nKinkedTracks, nInteractionsWithMaterial);
817884
} else {
818-
rowMcMatchRec(flag, origin, -1.f, 0, nKinkedTracks, nInteractionsWithMaterial);
885+
rowMcMatchRec(flag, origin, channel, -1.f, 0, nKinkedTracks, nInteractionsWithMaterial);
819886
}
820887
}
821888

@@ -842,11 +909,11 @@ struct HfCandidateCreator2ProngExpressions {
842909
if (rejectionMask != 0) {
843910
// at least one event selection not satisfied --> reject all particles from this collision
844911
for (unsigned int i = 0; i < mcParticlesPerMcColl.size(); ++i) {
845-
rowMcMatchGen(0, 0, -1);
912+
rowMcMatchGen(0, 0, 0, -1);
846913
}
847914
continue;
848915
}
849-
hf_mc_gen::fillMcMatchGen2Prong(mcParticles, mcParticlesPerMcColl, rowMcMatchGen, rejectBackground);
916+
hf_mc_gen::fillMcMatchGen2Prong(mcParticles, mcParticlesPerMcColl, rowMcMatchGen, rejectBackground, matchCorrelatedBackgrounds);
850917
}
851918
}
852919

0 commit comments

Comments
 (0)