2121#include " Common/DataModel/EventSelection.h"
2222#include " Common/DataModel/TrackSelectionTables.h"
2323
24- #include " PWGJE/DataModel/JetHF .h"
24+ #include " PWGJE/DataModel/Jet .h"
2525#include " PWGJE/Core/JetUtilities.h"
2626#include " PWGHF/DataModel/CandidateReconstructionTables.h"
2727#include " PWGHF/DataModel/CandidateSelectionTables.h"
@@ -30,79 +30,93 @@ using namespace o2;
3030using namespace o2 ::framework;
3131using namespace o2 ::framework::expressions;
3232
33+ template <typename BaseJetCollection, typename TagJetCollection,
34+ typename BaseToTagMatchingTable, typename TagToBaseMatchingTable, typename HfCandidates>
3335struct JetMatchingHF {
36+ Configurable<float > maxMatchingDistance{" maxMatchingDistance" 0 .4f , " Max matching distance"
37+
38+ Produces<BaseToTagMatchingTable> jetsBaseToTag;
39+ Produces<TagToBaseMatchingTable> jetsTagToBase;
40+
41+ Preslice<TagJetCollection> mcJetsPerMcCollision = aod::jet::mcCollisionId;
42+
3443 using Collisions = soa::Join<aod::Collisions, aod::McCollisionLabels>;
3544 using Tracks = soa::Join<aod::Tracks, aod::McTrackLabels>;
36- using HfCandidates = soa::Join<aod::HfCand2Prong, aod::HfSelD0, aod::HfCand2ProngMcRec>;
3745 using McParticles = soa::Join<aod::McParticles, aod::HfCand2ProngMcGen>;
38- using DetectorLevelJets = soa::Join<aod::D0MCDJets, aod::D0MCDJetConstituents>;
39- using ParticleLevelJets = soa::Join<aod::D0MCPJets, aod::D0MCPJetConstituents>;
40-
41- Produces<aod::NewMatchedD0MCDJets> jetsPartToDetMatching;
42- Produces<aod::NewMatchedD0MCPJets> jetsDetToPartMatching;
43-
44- Preslice<ParticleLevelJets> ParticleLevelJetsPerMcCollision = aod::jet::mcCollisionId;
4546
4647 void init (InitContext const &)
4748 {
4849 }
4950
5051 void process (Collisions::iterator const & collision, aod::McCollisions const & mcCollisions,
51- ParticleLevelJets const & jetsMC, DetectorLevelJets const & jetsRec ,
52- HfCandidates const & hfcandidates, Tracks const & tracks ,
53- McParticles const & particlesMC )
52+ BaseJetCollection const & jetsBase, TagJetCollection const & jetsTag ,
53+ Tracks const & tracks, McParticles const & particlesMC ,
54+ HfCandidates const & hfcandidates )
5455 {
55- const auto jetsPL = jetsMC .sliceBy (ParticleLevelJetsPerMcCollision , collision.mcCollisionId ());
56+ const auto jetsPL = jetsTag .sliceBy (mcJetsPerMcCollision , collision.mcCollisionId ());
5657
57- // match rec to MC
58- for (const auto & jet : jetsRec) {
58+ // geometric matching
59+ std::vector<double > jetsBasePhi (jetsBase.size ());
60+ std::vector<double > jetsBaseEta (jetsBase.size ());
61+ for (auto jet : jetsBase) {
62+ jetsBasePhi.emplace_back (jet.phi ());
63+ jetsBaseEta.emplace_back (jet.eta ());
64+ }
65+ std::vector<double > jetsTagPhi (jetsTag.size ());
66+ std::vector<double > jetsTagEta (jetsTag.size ());
67+ for (auto & jet : jetsTag) {
68+ jetsTagPhi.emplace_back (jet.phi ());
69+ jetsTagEta.emplace_back (jet.eta ());
70+ }
71+ auto && [baseToTagIndexMap, tagToBaseIndexMap] = JetUtilities::MatchJetsGeometrically (jetsBasePhi, jetsBaseEta, jetsTagPhi, jetsTagEta, maxMatchingDistance);
72+
73+ // forward matching
74+ for (const auto & jet : jetsBase) {
5975 LOGF (info, " jet index: %d (coll %d, pt %g, phi %g) with %d tracks, %d HF candidates"
6076 jet.index (), jet.collisionId (), jet.pt (), jet.phi (), jet.tracks ().size (), jet.hfcandidates ().size ());
6177
62- const auto & cands = jet.hfcandidates_as <HfCandidates>();
78+ const auto & cands = jet.template hfcandidates_as <HfCandidates>();
6379 int matchedIdx = -1 ;
64- if ((cands.front ().flagMcMatchRec () & (1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) == 0 ) {
65- jetsDetToPartMatching (matchedIdx, -1 );
66- continue ;
67- }
68- for (const auto & cand : cands) {
69- const auto & daughter0 = cand.prong0_as <Tracks>();
70- const auto & daughter1 = cand.prong1_as <Tracks>();
71- if (!daughter0.has_mcParticle () || !daughter1.has_mcParticle ()) {
72- LOGF (warning, " Encountered candidate daughter (%d or %d) without MC particle" globalIndex (), daughter1.globalIndex ());
73- continue ;
74- }
75- const auto mother0Id = daughter0.mcParticle_as <McParticles>().mothers_as <McParticles>().front ().globalIndex ();
76- const auto mother1Id = daughter1.mcParticle_as <McParticles>().mothers_as <McParticles>().front ().globalIndex ();
77- LOGF (debug, " MC candidate %d with prongs: %d (MC %d), %d (MC %d)" globalIndex (),
78- daughter0.globalIndex (), daughter0.mcParticle_as <McParticles>().globalIndex (),
79- daughter1.globalIndex (), daughter1.mcParticle_as <McParticles>().globalIndex ());
80- LOGF (info, " MC ids of mothers: %d - %d"
81- if ((mother0Id == mother1Id) &&
82- std::abs (daughter0.mcParticle_as <McParticles>().mothers_as <McParticles>().front ().flagMcMatchGen ()) & (1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) {
83- LOGF (info, " D0 - looking for jet"
84- for (const auto & pjet : jetsPL) {
85- for (const auto & cand : pjet.hfcandidates_as <McParticles>()) {
86- if (mother0Id == cand.globalIndex ()) {
87- matchedIdx = pjet.globalIndex ();
88- LOGF (info, " Found match det to part: %d (pt %g) -> %d (pt %g)"
89- jet.globalIndex (), jet.pt (), matchedIdx, pjet.pt ());
80+ if ((cands.front ().flagMcMatchRec () & (1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) != 0 ) {
81+ for (const auto & cand : cands) {
82+ const auto & daughter0 = cand.template prong0_as <Tracks>();
83+ const auto & daughter1 = cand.template prong1_as <Tracks>();
84+ if (!daughter0.has_mcParticle () || !daughter1.has_mcParticle ()) {
85+ LOGF (warning, " Encountered candidate daughter (%d or %d) without MC particle" globalIndex (), daughter1.globalIndex ());
86+ continue ;
87+ }
88+ const auto mother0Id = daughter0.template mcParticle_as <McParticles>().template mothers_as <McParticles>().front ().globalIndex ();
89+ const auto mother1Id = daughter1.template mcParticle_as <McParticles>().template mothers_as <McParticles>().front ().globalIndex ();
90+ LOGF (debug, " MC candidate %d with prongs: %d (MC %d), %d (MC %d)" globalIndex (),
91+ daughter0.globalIndex (), daughter0.template mcParticle_as <McParticles>().globalIndex (),
92+ daughter1.globalIndex (), daughter1.template mcParticle_as <McParticles>().globalIndex ());
93+ LOGF (info, " MC ids of mothers: %d - %d"
94+ if ((mother0Id == mother1Id) &&
95+ std::abs (daughter0.template mcParticle_as <McParticles>().template mothers_as <McParticles>().front ().flagMcMatchGen ()) & (1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) {
96+ LOGF (info, " D0 - looking for jet"
97+ for (const auto & pjet : jetsPL) {
98+ for (const auto & cand : pjet.template hfcandidates_as <McParticles>()) {
99+ if (mother0Id == cand.globalIndex ()) {
100+ matchedIdx = pjet.globalIndex ();
101+ LOGF (info, " Found match det to part: %d (pt %g) -> %d (pt %g)"
102+ jet.globalIndex (), jet.pt (), matchedIdx, pjet.pt ());
103+ }
90104 }
91105 }
92106 }
93107 }
94108 }
95- jetsDetToPartMatching (matchedIdx, - 1 );
109+ jetsBaseToTag (matchedIdx, baseToTagIndexMap[jet. index ()]); // TODO: check usage of index
96110 }
97111
98- // match MC to rec
112+ // backward matching
99113 for (const auto & jet : jetsPL) {
100114 LOGF (info, " MC jet index: %d (coll %d) with %d tracks, %d HF candidates"
101115 jet.index (), jet.mcCollisionId (), jet.tracks ().size (), jet.hfcandidates ().size ());
102116
103117 int matchedIdx = -1 ;
104- for (const auto & cand : jet.hfcandidates_as <McParticles>()) {
105- const auto & daughters = cand.daughters_as <McParticles>();
118+ for (const auto & cand : jet.template hfcandidates_as <McParticles>()) {
119+ const auto & daughters = cand.template daughters_as <McParticles>();
106120 LOGF (info, " MC candidate %d with daughters %d, %d" globalIndex (), daughters.iteratorAt (0 ).globalIndex (), daughters.iteratorAt (1 ).globalIndex ());
107121 int index0 = -1 , index1 = -1 ;
108122 for (const auto & track : tracks) {
@@ -124,27 +138,30 @@ struct JetMatchingHF {
124138 int candIdx = 0 ;
125139 for (const auto & prong : hfcandidates) {
126140 LOGF (info, " checking prong %d with daughters %d-%d, %d-%d"
127- prong.globalIndex (), prong.prong0Id (), prong.prong0_as <Tracks>().globalIndex (), prong.prong1Id (), prong.prong1_as <Tracks>().globalIndex ());
128- if ((prong.prong0_as <Tracks>().globalIndex () == index0 && prong.prong1_as <Tracks>().globalIndex () == index1) ||
141+ prong.globalIndex (), prong.prong0Id (), prong.template prong0_as <Tracks>().globalIndex (), prong.prong1Id (), prong.template prong1_as <Tracks>().globalIndex ());
142+ if ((prong.template prong0_as <Tracks>().globalIndex () == index0 && prong.template prong1_as <Tracks>().globalIndex () == index1) ||
129143 (prong.prong0Id () == index1 && prong.prong1Id () == index0)) {
130144 candIdx = prong.globalIndex ();
131145 LOGF (info, " Found matching 2prong candidate: %d"
132146 }
133147 }
134- for (const auto & djet : jetsRec ) {
135- if (djet.hfcandidates_as <HfCandidates>().front ().globalIndex () == candIdx) {
148+ for (const auto & djet : jetsBase ) {
149+ if (djet.template hfcandidates_as <HfCandidates>().front ().globalIndex () == candIdx) {
136150 matchedIdx = djet.globalIndex ();
137151 LOGF (info, " Found match part to det: %d -> %d" globalIndex (), matchedIdx);
138152 }
139153 }
140154 }
141- jetsPartToDetMatching (matchedIdx, - 1 );
155+ jetsTagToBase (matchedIdx, tagToBaseIndexMap[jet. index ()]); // TODO: check usage of index
142156 }
143157 }
144158};
145159
146160WorkflowSpec defineDataProcessing (ConfigContext const & cfgc)
147161{
148162 return WorkflowSpec{
149- adaptAnalysisTask<JetMatchingHF>(cfgc, TaskName{" jet-matching-hf"
163+ adaptAnalysisTask<JetMatchingHF<soa::Join<aod::D0MCDJets, aod::D0MCDJetConstituents>,
164+ soa::Join<aod::D0MCPJets, aod::D0MCPJetConstituents>,
165+ aod::NewMatchedD0MCDJets, aod::NewMatchedD0MCPJets,
166+ soa::Join<aod::HfCand2Prong, aod::HfSelD0, aod::HfCand2ProngMcRec> > >(cfgc, TaskName{" jet-matching-hf"
150167}
0 commit comments