diff --git a/PWGJE/DataModel/JetSubstructure.h b/PWGJE/DataModel/JetSubstructure.h
index bca64ff1597..15ccb7d78d7 100644
--- a/PWGJE/DataModel/JetSubstructure.h
+++ b/PWGJE/DataModel/JetSubstructure.h
@@ -48,21 +48,37 @@ DECLARE_SOA_COLUMN(JetNConstituents, jetNConstituents, int);   //!
 } // namespace jetoutput
 
 // Defines the jet substrcuture table definition
-#define JETSUBSTRUCTURE_TABLE_DEF(_collision_type_, _jet_type_, _cand_type_, _name_, _description_)                                                                                                                                                                                  \
-  namespace _name_##substructure                                                                                                                                                                                                                                                     \
-  {                                                                                                                                                                                                                                                                                  \
-    DECLARE_SOA_INDEX_COLUMN(_jet_type_, jet);                                                                                                                                                                                                                                       \
-    DECLARE_SOA_INDEX_COLUMN_FULL(Candidate, candidate, int, _cand_type_, "_0");                                                                                                                                                                                                     \
-    DECLARE_SOA_DYNAMIC_COLUMN(Dummy##_jet_type_, dummy##_jet_type_, []() -> int { return 0; });                                                                                                                                                                                     \
-  }                                                                                                                                                                                                                                                                                  \
-  DECLARE_SOA_TABLE(_jet_type_##Substructures, "AOD", _description_ "SS", jetsubstructure::Zg, jetsubstructure::Rg, jetsubstructure::Nsd, _name_##substructure::Dummy##_jet_type_<>);                                                                                                \
-  DECLARE_SOA_TABLE(_jet_type_##Output, "AOD", _description_ "O", jetoutput::_collision_type_##Id, _name_##substructure::_jet_type_##Id, _name_##substructure::Candidate##Id, jetoutput::JetPt, jetoutput::JetPhi, jetoutput::JetEta, jetoutput::JetR, jetoutput::JetNConstituents); \
+#define JETSUBSTRUCTURE_TABLE_DEF(_collision_type_, _jet_type_, _matched_jet_type_, _cand_type_, _name_, _description_)                                                                                                                                                                                                                                                                                                     \
+  namespace _name_##substructure                                                                                                                                                                                                                                                                                                                                                                                            \
+  {                                                                                                                                                                                                                                                                                                                                                                                                                         \
+    DECLARE_SOA_INDEX_COLUMN(_jet_type_, jet);                                                                                                                                                                                                                                                                                                                                                                              \
+    DECLARE_SOA_INDEX_COLUMN_FULL(Candidate, candidate, int, _cand_type_, "_0");                                                                                                                                                                                                                                                                                                                                            \
+    DECLARE_SOA_DYNAMIC_COLUMN(Dummy##_jet_type_, dummy##_jet_type_, []() -> int { return 0; });                                                                                                                                                                                                                                                                                                                            \
+  }                                                                                                                                                                                                                                                                                                                                                                                                                         \
+                                                                                                                                                                                                                                                                                                                                                                                                                            \
+  namespace _name_##geomatched                                                                                                                                                                                                                                                                                                                                                                                              \
+  {                                                                                                                                                                                                                                                                                                                                                                                                                         \
+    DECLARE_SOA_ARRAY_INDEX_COLUMN_FULL(_matched_jet_type_, matchedJetGeo, int32_t, _matched_jet_type_##s, "_geo");                                                                                                                                                                                                                                                                                                         \
+  }                                                                                                                                                                                                                                                                                                                                                                                                                         \
+                                                                                                                                                                                                                                                                                                                                                                                                                            \
+  namespace _name_##ptmatched                                                                                                                                                                                                                                                                                                                                                                                               \
+  {                                                                                                                                                                                                                                                                                                                                                                                                                         \
+    DECLARE_SOA_ARRAY_INDEX_COLUMN_FULL(_matched_jet_type_, matchedJetPt, int32_t, _matched_jet_type_##s, "_pt");                                                                                                                                                                                                                                                                                                           \
+  }                                                                                                                                                                                                                                                                                                                                                                                                                         \
+                                                                                                                                                                                                                                                                                                                                                                                                                            \
+  namespace _name_##candmatched                                                                                                                                                                                                                                                                                                                                                                                             \
+  {                                                                                                                                                                                                                                                                                                                                                                                                                         \
+    DECLARE_SOA_ARRAY_INDEX_COLUMN_FULL(_matched_jet_type_, matchedJetCand, int32_t, _matched_jet_type_##s, "_hf");                                                                                                                                                                                                                                                                                                         \
+  }                                                                                                                                                                                                                                                                                                                                                                                                                         \
+                                                                                                                                                                                                                                                                                                                                                                                                                            \
+  DECLARE_SOA_TABLE(_jet_type_##Substructures, "AOD", _description_ "SS", jetsubstructure::Zg, jetsubstructure::Rg, jetsubstructure::Nsd, _name_##substructure::Dummy##_jet_type_<>);                                                                                                                                                                                                                                       \
+  DECLARE_SOA_TABLE(_jet_type_##Output, "AOD", _description_ "O", jetoutput::_collision_type_##Id, _name_##substructure::_jet_type_##Id, _name_##substructure::Candidate##Id, _name_##geomatched::_matched_jet_type_##Ids, _name_##ptmatched::_matched_jet_type_##Ids, _name_##candmatched::_matched_jet_type_##Ids, jetoutput::JetPt, jetoutput::JetPhi, jetoutput::JetEta, jetoutput::JetR, jetoutput::JetNConstituents); \
   DECLARE_SOA_TABLE(_jet_type_##SubstructureOutput, "AOD", _description_ "SSO", _name_##substructure::_jet_type_##Id, jetsubstructure::Zg, jetsubstructure::Rg, jetsubstructure::Nsd);
 
-#define JETSUBSTRUCTURE_TABLES_DEF(_jet_type_, _cand_type_, _description_)                                                               \
-  JETSUBSTRUCTURE_TABLE_DEF(Collision, _jet_type_##Jet, _cand_type_, _jet_type_##jet, _description_)                                     \
-  JETSUBSTRUCTURE_TABLE_DEF(Collision, _jet_type_##MCDetectorLevelJet, _cand_type_, _jet_type_##mcdetectorleveljet, _description_ "MCD") \
-  JETSUBSTRUCTURE_TABLE_DEF(McCollision, _jet_type_##MCParticleLevelJet, McParticles, _jet_type_##mcparticleleveljet, _description_ "MCP")
+#define JETSUBSTRUCTURE_TABLES_DEF(_jet_type_, _cand_type_, _description_)                                                                                               \
+  JETSUBSTRUCTURE_TABLE_DEF(Collision, _jet_type_##Jet, _jet_type_##Jet, _cand_type_, _jet_type_##jet, _description_)                                                    \
+  JETSUBSTRUCTURE_TABLE_DEF(Collision, _jet_type_##MCDetectorLevelJet, _jet_type_##MCParticleLevelJet, _cand_type_, _jet_type_##mcdetectorleveljet, _description_ "MCD") \
+  JETSUBSTRUCTURE_TABLE_DEF(McCollision, _jet_type_##MCParticleLevelJet, _jet_type_##MCDetectorLevelJet, McParticles, _jet_type_##mcparticleleveljet, _description_ "MCP")
 
 JETSUBSTRUCTURE_TABLES_DEF(Charged, HfCand2Prong, "C");
 JETSUBSTRUCTURE_TABLES_DEF(D0Charged, HfCand2Prong, "D0");
diff --git a/PWGJE/Tasks/jetsubstructurehfoutput.cxx b/PWGJE/Tasks/jetsubstructurehfoutput.cxx
index 1ddb9fba000..579ab6c0497 100644
--- a/PWGJE/Tasks/jetsubstructurehfoutput.cxx
+++ b/PWGJE/Tasks/jetsubstructurehfoutput.cxx
@@ -39,10 +39,14 @@ using namespace o2::framework::expressions;
 // NB: runDataProcessing.h must be included after customize!
 #include "Framework/runDataProcessing.h"
 
-template <typename CollisionTable, typename JetTable, typename CandidateTable, typename OutputTable, typename SubstructureOutputTable>
+template <typename CollisionTable, typename McCollisionTable, typename JetTableData, typename CandidateTableData, typename OutputTableData, typename SubstructureOutputTableData, typename JetTableMCD, typename CandidateTableMCD, typename OutputTableMCD, typename SubstructureOutputTableMCD, typename JetTableMCP, typename CandidateTableMCP, typename OutputTableMCP, typename SubstructureOutputTableMCP>
 struct JetSubstructureHFOutputTask {
-  Produces<OutputTable> jetOutputTable;
-  Produces<SubstructureOutputTable> jetSubstructureOutputTable;
+  Produces<OutputTableData> jetOutputTableData;
+  Produces<SubstructureOutputTableData> jetSubstructureOutputTableData;
+  Produces<OutputTableMCD> jetOutputTableMCD;
+  Produces<SubstructureOutputTableMCD> jetSubstructureOutputTableMCD;
+  Produces<OutputTableMCP> jetOutputTableMCP;
+  Produces<SubstructureOutputTableMCP> jetSubstructureOutputTableMCP;
 
   Configurable<float> jetPtMin{"jetPtMin", 0.0, "minimum jet pT cut"};
   Configurable<std::vector<double>> jetRadii{"jetRadii", std::vector<double>{0.4}, "jet resolution parameters"};
@@ -56,84 +60,100 @@ struct JetSubstructureHFOutputTask {
 
   Filter jetSelection = aod::jet::pt >= jetPtMin;
 
-  template <typename T, typename U, typename V>
-  void fillTables(T const& collision, U const& jet, V const& cands)
+  template <typename T, typename U, typename V, typename M, typename N>
+  void fillTables(T const& collision, U const& jet, V const& cand, M& jetOutputTable, N& jetSubstructureOutputTable, std::vector<int> geoMatching, std::vector<int> ptMatching, std::vector<int> candMatching)
   {
-    auto cand = cands[0];
-    jetOutputTable(collision.globalIndex(), jet.globalIndex(), cand.globalIndex(), jet.pt(), jet.phi(), jet.eta(), jet.r(), jet.tracks().size());
+    jetOutputTable(collision.globalIndex(), jet.globalIndex(), cand.globalIndex(), geoMatching, ptMatching, candMatching, jet.pt(), jet.phi(), jet.eta(), jet.r(), jet.tracks().size() + jet.hfcandidates().size());
     jetSubstructureOutputTable(jet.globalIndex(), jet.zg(), jet.rg(), jet.nsd());
   }
 
-  void processDummy(aod::Collision const& collision)
-  {
-  }
+  void processDummy(typename CollisionTable::iterator const& collision) {}
   PROCESS_SWITCH(JetSubstructureHFOutputTask, processDummy, "Dummy process function turned on by default", true);
 
-  void processOutput(typename CollisionTable::iterator const& collision, JetTable const& jets, // add template back
-                     CandidateTable const& candidates, aod::Tracks const& tracks)
+  void processOutputData(typename CollisionTable::iterator const& collision,
+                         JetTableData const& jets)
   {
+    std::vector<int> geoMatching{-1};
+    std::vector<int> ptMatching{-1};
+    std::vector<int> candMatching{-1};
     for (const auto& jet : jets) {
       for (const auto& jetRadiiValue : jetRadiiValues) {
         if (jet.r() == round(jetRadiiValue * 100.0f)) {
-          auto cands = jet.template hfcandidates_as<CandidateTable>();
-          fillTables(collision, jet, cands);
+          auto cands = jet.template hfcandidates_as<CandidateTableData>();
+          auto cand = cands[0];
+          fillTables(collision, jet, cand, jetOutputTableData, jetSubstructureOutputTableData, geoMatching, ptMatching, candMatching);
         }
       }
     }
   }
-  PROCESS_SWITCH(JetSubstructureHFOutputTask, processOutput, "HF jet substructure output on data", false);
-};
-using JetSubstructureOutputDataD0 = JetSubstructureHFOutputTask<aod::Collisions, soa::Filtered<soa::Join<aod::D0ChargedJets, aod::D0ChargedJetConstituents, aod::D0ChargedJetSubstructures>>, soa::Join<aod::HfCand2Prong, aod::HfSelD0>, aod::D0ChargedJetOutput, aod::D0ChargedJetSubstructureOutput>;
-using JetSubstructureOutputMCDetectorLevelD0 = JetSubstructureHFOutputTask<aod::Collisions, soa::Filtered<soa::Join<aod::D0ChargedMCDetectorLevelJets, aod::D0ChargedMCDetectorLevelJetConstituents, aod::D0ChargedMCDetectorLevelJetSubstructures>>, soa::Join<aod::HfCand2Prong, aod::HfSelD0, aod::HfCand2ProngMcRec>, aod::D0ChargedMCDetectorLevelJetOutput, aod::D0ChargedMCDetectorLevelJetSubstructureOutput>;
-using JetSubstructureOutputMCParticleLevelD0 = JetSubstructureHFOutputTask<aod::McCollisions, soa::Filtered<soa::Join<aod::D0ChargedMCParticleLevelJets, aod::D0ChargedMCParticleLevelJetConstituents, aod::D0ChargedMCParticleLevelJetSubstructures>>, aod::McParticles, aod::D0ChargedMCParticleLevelJetOutput, aod::D0ChargedMCParticleLevelJetSubstructureOutput>;
+  PROCESS_SWITCH(JetSubstructureHFOutputTask, processOutputData, "hf jet substructure output Data", false);
 
-// using JetSubstructureOutputDataLc = JetSubstructureHFOutputTask<aod::Collisions, soa::Filtered<soa::Join<aod::LcChargedJets, aod::LcChargedJetConstituents, aod::LcChargedJetSubstructures>>, soa::Join<aod::HfCand3Prong, aod::HfSelLc>, aod::LcChargedJetOutput, aod::LcChargedJetSubstructureOutput>;
-//  using JetSubstructureOutputMCDetectorLevelLc = JetSubstructureHFOutputTask<aod::Collisions, soa::Filtered<soa::Join<aod::LcChargedMCDetectorLevelJets, aod::LcChargedMCDetectorLevelJetConstituents, aod::LcChargedMCDetectorLevelJetSubstructures>>, soa::Join<aod::HfCand3Prong, aod::HfSelLc, aod::HfCand3ProngMcRec>,  aod::LcChargedMCDetectorLevelJetOutput, aod::LcChargedMCDetectorLevelJetSubstructureOutput>;
-//  using JetSubstructureOutputMCParticleLevelLc = JetSubstructureHFOutputTask<aod::McCollisions, soa::Filtered<soa::Join<aod::LcChargedMCParticleLevelJets, aod::LcChargedMCParticleLevelJetConstituents, aod::LcChargedMCParticleLevelJetSubstructures>>, aod::McParticles, aod::LcChargedMCParticleLevelJetOutput, aod::LcChargedMCParticleLevelJetSubstructureOutput>;
+  void processOutputMCD(typename CollisionTable::iterator const& collision,
+                        JetTableMCD const& mcdjets,
+                        JetTableMCP const& mcpjets)
+  {
+    for (const auto& mcdjet : mcdjets) {
+      for (const auto& jetRadiiValue : jetRadiiValues) {
+        if (mcdjet.r() == round(jetRadiiValue * 100.0f)) {
+          auto cands = mcdjet.template hfcandidates_as<CandidateTableMCD>();
+          auto cand = cands[0];
+          std::vector<int> geoMatching;
+          std::vector<int> ptMatching;
+          std::vector<int> candMatching;
+          for (auto& mcpjet : mcdjet.template matchedJetGeo_as<JetTableMCP>()) {
+            geoMatching.push_back(mcpjet.globalIndex());
+          }
+          for (auto& mcpjet : mcdjet.template matchedJetPt_as<JetTableMCP>()) {
+            ptMatching.push_back(mcpjet.globalIndex());
+          }
+          for (auto& mcpjet : mcdjet.template matchedJetCand_as<JetTableMCP>()) {
+            candMatching.push_back(mcpjet.globalIndex());
+          }
+          fillTables(collision, mcdjet, cand, jetOutputTableMCD, jetSubstructureOutputTableMCD, geoMatching, ptMatching, candMatching);
+        }
+      }
+    }
+  }
+  PROCESS_SWITCH(JetSubstructureHFOutputTask, processOutputMCD, "hf jet substructure output MCD", false);
 
-// using JetSubstructureOutputDataBplus = JetSubstructureHFOutputTask<aod::Collisions, soa::Filtered<soa::Join<aod::BplusChargedJets, aod::BplusChargedJetConstituents, aod::BplusChargedJetSubstructures>> , soa::Join<aod::HfCandBplus, aod::HfSelBplusToD0Pi>, aod::BplusChargedJetOutput, aod::BplusChargedJetSubstructureOutput>;
-// using JetSubstructureOutputMCDetectorLevelBplus = JetSubstructureHFOutputTask<aod::Collisions, soa::Filtered<soa::Join<aod::BplusChargedMCDetectorLevelJets, aod::BplusChargedMCDetectorLevelJetConstituents, aod::BplusChargedMCDetectorLevelJetSubstructures>>, soa::Join<aod::HfCandBplus, aod::HfSelBplusToD0Pi, aod::HfCandBplusMcRec>,  aod::BplusChargedMCDetectorLevelJetOutput, aod::BplusChargedMCDetectorLevelJetSubstructureOutput>;
-// using JetSubstructureOutputMCParticleLevelBplus = JetSubstructureHFOutputTask<aod::McCollisions, soa::Filtered<soa::Join<aod::BplusChargedMCParticleLevelJets, aod::BplusChargedMCParticleLevelJetConstituents, aod::BplusChargedMCParticleLevelJetSubstructures>>, aod::McParticles, aod::BplusChargedMCParticleLevelJetOutput, aod::BplusChargedMCParticleLevelJetSubstructureOutput>;
+  void processOutputMCP(typename McCollisionTable::iterator const& collision,
+                        JetTableMCP const& mcpjets,
+                        JetTableMCD const& mcdjets)
+  {
+    for (const auto& mcpjet : mcpjets) {
+      for (const auto& jetRadiiValue : jetRadiiValues) {
+        if (mcpjet.r() == round(jetRadiiValue * 100.0f)) {
+          auto cands = mcpjet.template hfcandidates_as<CandidateTableMCP>();
+          auto cand = cands[0];
+          std::vector<int> geoMatching;
+          std::vector<int> ptMatching;
+          std::vector<int> candMatching;
+          for (auto& mcdjet : mcpjet.template matchedJetGeo_as<JetTableMCD>()) {
+            geoMatching.push_back(mcdjet.globalIndex());
+          }
+          for (auto& mcdjet : mcpjet.template matchedJetPt_as<JetTableMCD>()) {
+            ptMatching.push_back(mcdjet.globalIndex());
+          }
+          for (auto& mcdjet : mcpjet.template matchedJetCand_as<JetTableMCD>()) {
+            candMatching.push_back(mcdjet.globalIndex());
+          }
+          fillTables(collision, mcpjet, cand, jetOutputTableMCP, jetSubstructureOutputTableMCP, geoMatching, ptMatching, candMatching);
+        }
+      }
+    }
+  }
+  PROCESS_SWITCH(JetSubstructureHFOutputTask, processOutputMCP, "hf jet substructure output MCP", false);
+};
+using JetSubstructureOutputD0 = JetSubstructureHFOutputTask<aod::Collisions, aod::McCollisions, soa::Filtered<soa::Join<aod::D0ChargedJets, aod::D0ChargedJetConstituents, aod::D0ChargedJetSubstructures>>, soa::Join<aod::HfCand2Prong, aod::HfSelD0>, aod::D0ChargedJetOutput, aod::D0ChargedJetSubstructureOutput, soa::Filtered<soa::Join<aod::D0ChargedMCDetectorLevelJets, aod::D0ChargedMCDetectorLevelJetConstituents, aod::D0ChargedMCDetectorLevelJetSubstructures, aod::D0ChargedMCDetectorLevelJetsMatchedToD0ChargedMCParticleLevelJets>>, soa::Join<aod::HfCand2Prong, aod::HfSelD0, aod::HfCand2ProngMcRec>, aod::D0ChargedMCDetectorLevelJetOutput, aod::D0ChargedMCDetectorLevelJetSubstructureOutput, soa::Filtered<soa::Join<aod::D0ChargedMCParticleLevelJets, aod::D0ChargedMCParticleLevelJetConstituents, aod::D0ChargedMCParticleLevelJetSubstructures, aod::D0ChargedMCParticleLevelJetsMatchedToD0ChargedMCDetectorLevelJets>>, aod::McParticles, aod::D0ChargedMCParticleLevelJetOutput, aod::D0ChargedMCParticleLevelJetSubstructureOutput>;
+// using JetSubstructureOutputLc = JetSubstructureHFOutputTask<aod::Collisions, aod::McCollisions, soa::Filtered<soa::Join<aod::LcChargedJets, aod::LcChargedJetConstituents, aod::LcChargedJetSubstructures>>, soa::Join<aod::HfCand3Prong, aod::HfSelLc>, aod::LcChargedJetOutput, aod::LcChargedJetSubstructureOutput,soa::Filtered<soa::Join<aod::LcChargedMCDetectorLevelJets, aod::LcChargedMCDetectorLevelJetConstituents, aod::LcChargedMCDetectorLevelJetSubstructures,aod::LcChargedMCDetectorLevelJetsMatchedToLcChargedMCParticleLevelJets>>, soa::Join<aod::HfCand3Prong, aod::HfSelLc, aod::HfCand3ProngMcRec>, aod::LcChargedMCDetectorLevelJetOutput, aod::LcChargedMCDetectorLevelJetSubstructureOutput, soa::Filtered<soa::Join<aod::LcChargedMCParticleLevelJets, aod::LcChargedMCParticleLevelJetConstituents, aod::LcChargedMCParticleLevelJetSubstructures,aod::LcChargedMCParticleLevelJetsMatchedToLcChargedMCDetectorLevelJets>>, aod::McParticles, aod::LcChargedMCParticleLevelJetOutput, aod::LcChargedMCParticleLevelJetSubstructureOutput>;
+// using JetSubstructureOutputBplus = JetSubstructureHFOutputTask<aod::Collisions, aod::McCollisions, soa::Filtered<soa::Join<aod::BplusChargedJets, aod::BplusChargedJetConstituents, aod::BplusChargedJetSubstructures>>, soa::Join<aod::HfCandBplus, aod::HfSelBplusToD0Pi>, aod::BplusChargedJetOutput, aod::BplusChargedJetSubstructureOutput,soa::Filtered<soa::Join<aod::BplusChargedMCDetectorLevelJets, aod::BplusChargedMCDetectorLevelJetConstituents, aod::BplusChargedMCDetectorLevelJetSubstructures,aod::BplusChargedMCDetectorLevelJetsMatchedToBplusChargedMCParticleLevelJets>>, soa::Join<aod::HfCandBplus, aod::HfSelBplusToD0Pi, aod::HfCandBplusMcRec>, aod::BplusChargedMCDetectorLevelJetOutput, aod::BplusChargedMCDetectorLevelJetSubstructureOutput, soa::Filtered<soa::Join<aod::BplusChargedMCParticleLevelJets, aod::BplusChargedMCParticleLevelJetConstituents, aod::BplusChargedMCParticleLevelJetSubstructures,aod::BplusChargedMCParticleLevelJetsMatchedToBplusChargedMCDetectorLevelJets>>, aod::McParticles, aod::BplusChargedMCParticleLevelJetOutput, aod::BplusChargedMCParticleLevelJetSubstructureOutput>;
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {
   std::vector<o2::framework::DataProcessorSpec> tasks;
 
-  tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputDataD0>(cfgc,
-                                                                    SetDefaultProcesses{},
-                                                                    TaskName{"jet-substructure-output-D0-data"}));
-
-  tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputMCDetectorLevelD0>(cfgc,
-                                                                               SetDefaultProcesses{},
-                                                                               TaskName{"jet-substructure-output-D0-mcd"}));
-
-  tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputMCParticleLevelD0>(cfgc,
-                                                                               SetDefaultProcesses{},
-                                                                               TaskName{"jet-substructure-output-D0-mcp"}));
-  /*
-   tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputDataLc>(cfgc,
-                                                                     SetDefaultProcesses{},
-                                                                     TaskName{"jet-substructure-output-Lc-data"}));
-
-     tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputMCDetectorLevelLc>(cfgc,
-                                                                                  SetDefaultProcesses{},
-                                                                                  TaskName{"jet-substructure-output-Lc-mcd"}));
-
-     tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputMCParticleLevelLc>(cfgc,
-                                                                                  SetDefaultProcesses{},
-                                                                                  TaskName{"jet-substructure-output-Lc-mcp"}));
-
-     tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputDataBplus>(cfgc,
-                                                                          SetDefaultProcesses{},
-                                                                          TaskName{"jet-substructure-output-Bplus-data"}));
-
-     tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputMCDetectorLevelBplus>(cfgc,
-                                                                                     SetDefaultProcesses{},
-                                                                                     TaskName{"jet-substructure-output-Bplus-mcd"}));
-
-     tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputMCParticleLevelBplus>(cfgc,
-                                                                                     SetDefaultProcesses{},
-                                                                                     TaskName{"jet-substructure-output-Bplus-mcp"}));
-   */
+  tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputD0>(cfgc, SetDefaultProcesses{}, TaskName{"jet-substructured0-output"}));
+  // tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputLc>(cfgc, SetDefaultProcesses{}, TaskName{"jet-substructurelc-output"}));
+  // tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputBplus>(cfgc, SetDefaultProcesses{}, TaskName{"jet-substructurebplus-output"}));
   return WorkflowSpec{tasks};
 }
diff --git a/PWGJE/Tasks/jetsubstructureoutput.cxx b/PWGJE/Tasks/jetsubstructureoutput.cxx
index 2844eeca383..ce55a16cbce 100644
--- a/PWGJE/Tasks/jetsubstructureoutput.cxx
+++ b/PWGJE/Tasks/jetsubstructureoutput.cxx
@@ -36,10 +36,14 @@ using namespace o2::framework::expressions;
 // NB: runDataProcessing.h must be included after customize!
 #include "Framework/runDataProcessing.h"
 
-template <typename CollisionTable, typename JetTable, typename OutputTable, typename SubstructureOutputTable>
+template <typename CollisionTable, typename McCollisionTable, typename JetTableData, typename OutputTableData, typename SubstructureOutputTableData, typename JetTableMCD, typename OutputTableMCD, typename SubstructureOutputTableMCD, typename JetTableMCP, typename OutputTableMCP, typename SubstructureOutputTableMCP>
 struct JetSubstructureOutputTask {
-  Produces<OutputTable> jetOutputTable;
-  Produces<SubstructureOutputTable> jetSubstructureOutputTable;
+  Produces<OutputTableData> jetOutputTableData;
+  Produces<SubstructureOutputTableData> jetSubstructureOutputTableData;
+  Produces<OutputTableMCD> jetOutputTableMCD;
+  Produces<SubstructureOutputTableMCD> jetSubstructureOutputTableMCD;
+  Produces<OutputTableMCP> jetOutputTableMCP;
+  Produces<SubstructureOutputTableMCP> jetSubstructureOutputTableMCP;
 
   Configurable<float> jetPtMin{"jetPtMin", 0.0, "minimum jet pT cut"};
   Configurable<std::vector<double>> jetRadii{"jetRadii", std::vector<double>{0.4}, "jet resolution parameters"};
@@ -53,42 +57,85 @@ struct JetSubstructureOutputTask {
 
   Filter jetSelection = aod::jet::pt >= jetPtMin;
 
-  template <typename T, typename U>
-  void fillTables(T const& collision, U const& jets)
+  template <typename T, typename U, typename V, typename M>
+  void fillTables(T const& collision, U const& jet, V& jetOutputTable, M& jetSubstructureOutputTable, std::vector<int> geoMatching, std::vector<int> ptMatching, std::vector<int> candMatching)
   {
+    jetOutputTable(collision.globalIndex(), jet.globalIndex(), -1, geoMatching, ptMatching, candMatching, jet.pt(), jet.phi(), jet.eta(), jet.r(), jet.tracks().size());
+    jetSubstructureOutputTable(jet.globalIndex(), jet.zg(), jet.rg(), jet.nsd());
+  }
+
+  void processDummy(typename CollisionTable::iterator const& collision) {}
+  PROCESS_SWITCH(JetSubstructureOutputTask, processDummy, "Dummy process function turned on by default", true);
+
+  void processOutputData(typename CollisionTable::iterator const& collision,
+                         JetTableData const& jets)
+  {
+    std::vector<int> geoMatching{-1};
+    std::vector<int> ptMatching{-1};
+    std::vector<int> candMatching{-1};
     for (const auto& jet : jets) {
       for (const auto& jetRadiiValue : jetRadiiValues) {
         if (jet.r() == round(jetRadiiValue * 100.0f)) {
-          jetOutputTable(collision.globalIndex(), jet.globalIndex(), -1, jet.pt(), jet.phi(), jet.eta(), jet.r(), jet.tracks().size());
-          jetSubstructureOutputTable(jet.globalIndex(), jet.zg(), jet.rg(), jet.nsd());
+          fillTables(collision, jet, jetOutputTableData, jetSubstructureOutputTableData, geoMatching, ptMatching, candMatching);
         }
       }
     }
   }
+  PROCESS_SWITCH(JetSubstructureOutputTask, processOutputData, "jet substructure output Data", false);
 
-  void processDummy(typename CollisionTable::iterator const& collision) {}
-  PROCESS_SWITCH(JetSubstructureOutputTask, processDummy, "Dummy process function turned on by default", true);
+  void processOutputMCD(typename CollisionTable::iterator const& collision,
+                        JetTableMCD const& mcdjets,
+                        JetTableMCP const& mcpjets)
+  {
+    std::vector<int> candMatching{-1};
+    for (const auto& mcdjet : mcdjets) {
+      for (const auto& jetRadiiValue : jetRadiiValues) {
+        if (mcdjet.r() == round(jetRadiiValue * 100.0f)) {
+          std::vector<int> geoMatching;
+          std::vector<int> ptMatching;
+          for (auto& mcpjet : mcdjet.template matchedJetGeo_as<JetTableMCP>()) {
+            geoMatching.push_back(mcpjet.globalIndex());
+          }
+          for (auto& mcpjet : mcdjet.template matchedJetPt_as<JetTableMCP>()) {
+            ptMatching.push_back(mcpjet.globalIndex());
+          }
+          fillTables(collision, mcdjet, jetOutputTableMCD, jetSubstructureOutputTableMCD, geoMatching, ptMatching, candMatching);
+        }
+      }
+    }
+  }
+  PROCESS_SWITCH(JetSubstructureOutputTask, processOutputMCD, "jet substructure output MCD", false);
 
-  void processOutput(typename CollisionTable::iterator const& collision,
-                     JetTable const& jets)
+  void processOutputMCP(typename McCollisionTable::iterator const& collision,
+                        JetTableMCP const& mcpjets,
+                        JetTableMCD const& mcdjets)
   {
-    fillTables(collision, jets);
+    std::vector<int> candMatching{-1};
+    for (const auto& mcpjet : mcpjets) {
+      for (const auto& jetRadiiValue : jetRadiiValues) {
+        if (mcpjet.r() == round(jetRadiiValue * 100.0f)) {
+          std::vector<int> geoMatching;
+          std::vector<int> ptMatching;
+          for (auto& mcdjet : mcpjet.template matchedJetGeo_as<JetTableMCD>()) {
+            geoMatching.push_back(mcdjet.globalIndex());
+          }
+          for (auto& mcdjet : mcpjet.template matchedJetPt_as<JetTableMCD>()) {
+            ptMatching.push_back(mcdjet.globalIndex());
+          }
+          fillTables(collision, mcpjet, jetOutputTableMCP, jetSubstructureOutputTableMCP, geoMatching, ptMatching, candMatching);
+        }
+      }
+    }
   }
-  PROCESS_SWITCH(JetSubstructureOutputTask, processOutput, "jet substructure output", false);
+  PROCESS_SWITCH(JetSubstructureOutputTask, processOutputMCP, "jet substructure output MCP", false);
 };
-using JetSubstructureOutputData = JetSubstructureOutputTask<aod::Collisions, soa::Filtered<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents, aod::ChargedJetSubstructures>>, aod::ChargedJetOutput, aod::ChargedJetSubstructureOutput>;
-using JetSubstructureOutputMCDetectorLevel = JetSubstructureOutputTask<aod::Collisions, soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetSubstructures>>, aod::ChargedMCDetectorLevelJetOutput, aod::ChargedMCDetectorLevelJetSubstructureOutput>;
-using JetSubstructureOutputMCParticleLevel = JetSubstructureOutputTask<aod::McCollisions, soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents, aod::ChargedMCParticleLevelJetSubstructures>>, aod::ChargedMCParticleLevelJetOutput, aod::ChargedMCParticleLevelJetSubstructureOutput>;
+using JetSubstructureOutput = JetSubstructureOutputTask<aod::Collisions, aod::McCollisions, soa::Filtered<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents, aod::ChargedJetSubstructures>>, aod::ChargedJetOutput, aod::ChargedJetSubstructureOutput, soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetSubstructures, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>>, aod::ChargedMCDetectorLevelJetOutput, aod::ChargedMCDetectorLevelJetSubstructureOutput, soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents, aod::ChargedMCParticleLevelJetSubstructures, aod::ChargedMCParticleLevelJetsMatchedToChargedMCDetectorLevelJets>>, aod::ChargedMCParticleLevelJetOutput, aod::ChargedMCParticleLevelJetSubstructureOutput>;
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {
   std::vector<o2::framework::DataProcessorSpec> tasks;
 
-  tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputData>(cfgc, SetDefaultProcesses{}, TaskName{"jet-substructure-output-data"}));
-
-  tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputMCDetectorLevel>(cfgc, SetDefaultProcesses{}, TaskName{"jet-substructure-output-mcd"}));
-
-  tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutputMCParticleLevel>(cfgc, SetDefaultProcesses{}, TaskName{"jet-substructure-output-mcp"}));
+  tasks.emplace_back(adaptAnalysisTask<JetSubstructureOutput>(cfgc, SetDefaultProcesses{}, TaskName{"jet-substructure-output"}));
 
   return WorkflowSpec{tasks};
 }