diff --git a/Analysis/ALICE3/src/alice3-lutmaker.cxx b/Analysis/ALICE3/src/alice3-lutmaker.cxx
index fa512799f0338..3d8bd24cc2201 100644
--- a/Analysis/ALICE3/src/alice3-lutmaker.cxx
+++ b/Analysis/ALICE3/src/alice3-lutmaker.cxx
@@ -103,7 +103,7 @@ struct Alice3LutMaker {
       if (mcParticle.pdgCode() != pdg) {
         continue;
       }
-      if (selPrim.value && !MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) { // Requiring is physical primary
+      if (selPrim.value && !MC::isPhysicalPrimary(mcParticle)) { // Requiring is physical primary
         continue;
       }
 
@@ -133,7 +133,7 @@ struct Alice3LutMaker {
       if (mcParticle.pdgCode() != pdg) {
         continue;
       }
-      if (!MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) { // Requiring is physical primary
+      if (!MC::isPhysicalPrimary(mcParticle)) { // Requiring is physical primary
         continue;
       }
 
diff --git a/Analysis/ALICE3/src/alice3-qa-singleparticle.cxx b/Analysis/ALICE3/src/alice3-qa-singleparticle.cxx
index 268c4d7e35e3d..81938595e5064 100644
--- a/Analysis/ALICE3/src/alice3-qa-singleparticle.cxx
+++ b/Analysis/ALICE3/src/alice3-qa-singleparticle.cxx
@@ -118,13 +118,13 @@ struct Alice3SingleParticle {
           continue;
         }
         auto mother = mcParticle.mother0_as<aod::McParticles>();
-        if (MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) {
+        if (MC::isPhysicalPrimary(mcParticle)) {
           histos.get<TH1>(HIST("primaries"))->Fill(Form("%i", mother.pdgCode()), 1.f);
         } else {
           histos.get<TH1>(HIST("secondaries"))->Fill(Form("%i", mother.pdgCode()), 1.f);
         }
         if (doPrint) {
-          LOG(INFO) << "Track " << track.globalIndex() << " is a " << mcParticle.pdgCode() << " and comes from a " << mother.pdgCode() << " and is " << (MC::isPhysicalPrimary<aod::McParticles>(mcParticle) ? "" : "not") << " a primary";
+          LOG(INFO) << "Track " << track.globalIndex() << " is a " << mcParticle.pdgCode() << " and comes from a " << mother.pdgCode() << " and is " << (MC::isPhysicalPrimary(mcParticle) ? "" : "not") << " a primary";
         }
       }
     }
diff --git a/Analysis/ALICE3/src/pidFTOFqa.cxx b/Analysis/ALICE3/src/pidFTOFqa.cxx
index 4914de807c28a..7e0ee3f9335e5 100644
--- a/Analysis/ALICE3/src/pidFTOFqa.cxx
+++ b/Analysis/ALICE3/src/pidFTOFqa.cxx
@@ -152,7 +152,7 @@ struct ftofPidQaMC {
       if (pdgCode != 0 && abs(mcParticle.pdgCode()) != pdgCode) {
         continue;
       }
-      if (useOnlyPhysicsPrimary == 1 && !MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) { // Selecting primaries
+      if (useOnlyPhysicsPrimary == 1 && !MC::isPhysicalPrimary(mcParticle)) { // Selecting primaries
         histos.fill(HIST("p/Sec"), track.p());
         continue;
       }
diff --git a/Analysis/ALICE3/src/pidRICHqa.cxx b/Analysis/ALICE3/src/pidRICHqa.cxx
index 5cb700fd302ad..b89e315bbe514 100644
--- a/Analysis/ALICE3/src/pidRICHqa.cxx
+++ b/Analysis/ALICE3/src/pidRICHqa.cxx
@@ -167,7 +167,7 @@ struct richPidQaMc {
     if (abs(particle.pdgCode()) == PDGs[pidIndex]) {
       histos.fill(HIST(hnsigmaMC[pidIndex]), track.pt(), nsigma);
 
-      if (MC::isPhysicalPrimary<aod::McParticles>(particle)) { // Selecting primaries
+      if (MC::isPhysicalPrimary(particle)) { // Selecting primaries
         histos.fill(HIST(hnsigmaMCprm[pidIndex]), track.pt(), nsigma);
       } else {
         histos.fill(HIST(hnsigmaMCsec[pidIndex]), track.pt(), nsigma);
@@ -250,7 +250,7 @@ struct richPidQaMc {
       }
       histos.fill(HIST(hnsigma[pid_type]), track.pt(), nsigma);
       histos.fill(HIST(hdelta[pid_type]), track.p(), delta);
-      if (MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) { // Selecting primaries
+      if (MC::isPhysicalPrimary(mcParticle)) { // Selecting primaries
         histos.fill(HIST(hnsigmaprm[pid_type]), track.pt(), nsigma);
         histos.fill(HIST("p/Prim"), track.p());
       } else {
diff --git a/Analysis/Core/CMakeLists.txt b/Analysis/Core/CMakeLists.txt
index bc4a674b231b7..2b1266a3aa782 100644
--- a/Analysis/Core/CMakeLists.txt
+++ b/Analysis/Core/CMakeLists.txt
@@ -20,7 +20,6 @@ o2_target_root_dictionary(AnalysisCore
                       include/AnalysisCore/TrackSelection.h
                       include/AnalysisCore/TrackSelectionDefaults.h
                       include/AnalysisCore/TriggerAliases.h
-                      include/AnalysisCore/MC.h
               LINKDEF src/AnalysisCoreLinkDef.h)
 
 o2_add_executable(merger
diff --git a/Analysis/Core/include/AnalysisCore/MC.h b/Analysis/Core/include/AnalysisCore/MC.h
index 54e0c17cef813..5c5e725e5a33d 100644
--- a/Analysis/Core/include/AnalysisCore/MC.h
+++ b/Analysis/Core/include/AnalysisCore/MC.h
@@ -13,6 +13,7 @@
 #define MC_H
 
 #include "Framework/Logger.h"
+#include "Framework/AnalysisDataModel.h"
 
 #include "TPDGCode.h"
 
@@ -61,7 +62,7 @@ bool isStable(int pdg)
 
 // Ported from AliRoot AliStack::IsPhysicalPrimary
 template <typename Particles, typename Particle>
-bool isPhysicalPrimary(Particle& particle)
+bool isPhysicalPrimary(Particle const& particle)
 {
   // Test if a particle is a physical primary according to the following definition:
   // Particles produced in the collision including products of strong and
@@ -169,6 +170,12 @@ bool isPhysicalPrimary(Particle& particle)
   }
 }
 
+// Short hand for the standard type
+bool isPhysicalPrimary(o2::aod::McParticle const& particle)
+{
+  return isPhysicalPrimary<o2::aod::McParticles>(particle);
+}
+
 }; // namespace MC
 
 #endif
diff --git a/Analysis/Tasks/PID/qaTOFMC.cxx b/Analysis/Tasks/PID/qaTOFMC.cxx
index 8d4ba5314c95f..855922522470e 100644
--- a/Analysis/Tasks/PID/qaTOFMC.cxx
+++ b/Analysis/Tasks/PID/qaTOFMC.cxx
@@ -135,7 +135,7 @@ struct pidTOFTaskQA {
     if (abs(particle.pdgCode()) == PDGs[pidIndex]) {
       histos.fill(HIST(hnsigmaMC[pidIndex]), track.pt(), nsigma);
 
-      if (MC::isPhysicalPrimary<aod::McParticles>(particle)) { // Selecting primaries
+      if (MC::isPhysicalPrimary(particle)) { // Selecting primaries
         histos.fill(HIST(hnsigmaMCprm[pidIndex]), track.pt(), nsigma);
       } else {
         histos.fill(HIST(hnsigmaMCsec[pidIndex]), track.pt(), nsigma);
@@ -190,7 +190,7 @@ struct pidTOFTaskQA {
       histos.fill(HIST(hnsigma[pid_type]), t.pt(), nsigma);
       histos.fill(HIST("event/tofbeta"), t.p(), t.beta());
       const auto particle = t.mcParticle();
-      if (MC::isPhysicalPrimary<aod::McParticles>(particle)) { // Selecting primaries
+      if (MC::isPhysicalPrimary(particle)) { // Selecting primaries
         histos.fill(HIST(hnsigmaprm[pid_type]), t.pt(), nsigma);
         histos.fill(HIST("event/tofbetaPrm"), t.p(), t.beta());
       } else {
diff --git a/Analysis/Tasks/PWGLF/NucleiSpectraEfficiency.cxx b/Analysis/Tasks/PWGLF/NucleiSpectraEfficiency.cxx
index 19753e7b12643..cd6ef867b1765 100644
--- a/Analysis/Tasks/PWGLF/NucleiSpectraEfficiency.cxx
+++ b/Analysis/Tasks/PWGLF/NucleiSpectraEfficiency.cxx
@@ -70,7 +70,7 @@ struct NucleiSpectraEfficiencyGen {
     spectra.add("histGenPt", "generated particles", HistType::kTH1F, {ptAxis});
   }
 
-  void process(aod::McCollision const& mcCollision, aod::McParticles& mcParticles)
+  void process(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles)
   {
     //
     // loop over generated particles and fill generated particles
@@ -79,7 +79,7 @@ struct NucleiSpectraEfficiencyGen {
       if (mcParticleGen.pdgCode() != -1000020030) {
         continue;
       }
-      if (!MC::isPhysicalPrimary<aod::McParticles>(mcParticleGen)) {
+      if (!MC::isPhysicalPrimary(mcParticleGen)) {
         continue;
       }
       if (abs(mcParticleGen.y()) > 0.5) {
diff --git a/Analysis/Tasks/PWGLF/raacharged.cxx b/Analysis/Tasks/PWGLF/raacharged.cxx
index 373991a40c5d7..4ba0417463c27 100644
--- a/Analysis/Tasks/PWGLF/raacharged.cxx
+++ b/Analysis/Tasks/PWGLF/raacharged.cxx
@@ -152,7 +152,7 @@ struct raacharged {
 
   Configurable<bool> isMC{"isMC", 1, "0 - data, 1 - MC"};
 
-  void process(soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::Cents>::iterator const& collision, soa::Join<aod::Tracks, aod::McTrackLabels, aod::TrackSelection> const& tracks, aod::McParticles& mcParticles)
+  void process(soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::Cents>::iterator const& collision, soa::Join<aod::Tracks, aod::McTrackLabels, aod::TrackSelection> const& tracks, aod::McParticles const& mcParticles)
   {
     if (!collision.alias()[kINT7]) {
       return;
@@ -179,7 +179,7 @@ struct raacharged {
         continue;
       }
       const auto particle = track.mcParticle();
-      if (MC::isPhysicalPrimary<aod::McParticles>(particle)) {
+      if (MC::isPhysicalPrimary(particle)) {
         mcInfoVal = 0.0;
       } else {
         mcInfoVal = 1.0;
@@ -198,7 +198,7 @@ struct raacharged {
         if (abs(mcParticle.eta()) > 0.8) {
           continue;
         }
-        if (!MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) {
+        if (!MC::isPhysicalPrimary(mcParticle)) {
           continue;
         }
 
diff --git a/Analysis/Tasks/PWGLF/trackchecks.cxx b/Analysis/Tasks/PWGLF/trackchecks.cxx
index ba24fccd4259d..87b2394f6ef8c 100644
--- a/Analysis/Tasks/PWGLF/trackchecks.cxx
+++ b/Analysis/Tasks/PWGLF/trackchecks.cxx
@@ -90,7 +90,7 @@ struct TrackCheckTaskEvSel {
   //Filters
   Filter collfilter = nabs(aod::collision::posZ) < cfgCutVZ;
   void process(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels>>::iterator const& col,
-               soa::Join<aod::Tracks, aod::TracksExtra, aod::McTrackLabels>& tracks, aod::McParticles& mcParticles)
+               soa::Join<aod::Tracks, aod::TracksExtra, aod::McTrackLabels>& tracks, aod::McParticles const& mcParticles)
   {
 
     //event selection
@@ -110,7 +110,7 @@ struct TrackCheckTaskEvSel {
         const auto particle = track.mcParticle();
         int pdgcode = particle.pdgCode();
 
-        if (MC::isPhysicalPrimary<aod::McParticles>(particle)) { //is primary?
+        if (MC::isPhysicalPrimary(particle)) { //is primary?
           isPrimary = true;
         }
 
@@ -208,7 +208,7 @@ struct TrackCheckTaskEvSelTrackSel {
         const auto particle = track.mcParticle();
         int pdgcode = particle.pdgCode();
 
-        if (MC::isPhysicalPrimary<aod::McParticles>(particle)) { //is primary?
+        if (MC::isPhysicalPrimary(particle)) { //is primary?
           isPrimary = true;
         }
         //Calculate y
diff --git a/Analysis/Tasks/PWGPP/qaEfficiency.cxx b/Analysis/Tasks/PWGPP/qaEfficiency.cxx
index e33da8481e9d0..59208723852d3 100644
--- a/Analysis/Tasks/PWGPP/qaEfficiency.cxx
+++ b/Analysis/Tasks/PWGPP/qaEfficiency.cxx
@@ -229,7 +229,7 @@ struct QaTrackingEfficiency {
         return true;
       }
       histos.fill(h, 5);
-      if ((selPrim == 1) && (!MC::isPhysicalPrimary<o2::aod::McParticles>(p))) { // Requiring is physical primary
+      if ((selPrim == 1) && (!MC::isPhysicalPrimary(p))) { // Requiring is physical primary
         return true;
       }
       histos.fill(h, 6);
diff --git a/Analysis/Tasks/PWGPP/qaEventTrack.cxx b/Analysis/Tasks/PWGPP/qaEventTrack.cxx
index deb0d7252e3ab..55e71ba7cc60d 100644
--- a/Analysis/Tasks/PWGPP/qaEventTrack.cxx
+++ b/Analysis/Tasks/PWGPP/qaEventTrack.cxx
@@ -203,7 +203,7 @@ struct QaTrackingKine {
       if (pdgCodeSel != 0 && particle.pdgCode() != pdgCodeSel) { // Checking PDG code
         continue;
       }
-      if (MC::isPhysicalPrimary<o2::aod::McParticles>(particle)) {
+      if (MC::isPhysicalPrimary(particle)) {
         histos.fill(HIST("trackingPrm/pt"), t.pt());
         histos.fill(HIST("trackingPrm/eta"), t.eta());
         histos.fill(HIST("trackingPrm/phi"), t.phi());
@@ -223,7 +223,7 @@ struct QaTrackingKine {
       histos.fill(HIST("particle/pt"), particle.pt());
       histos.fill(HIST("particle/eta"), particle.eta());
       histos.fill(HIST("particle/phi"), particle.phi());
-      if (MC::isPhysicalPrimary<o2::aod::McParticles>(particle)) {
+      if (MC::isPhysicalPrimary(particle)) {
         histos.fill(HIST("particlePrm/pt"), particle.pt());
         histos.fill(HIST("particlePrm/eta"), particle.eta());
         histos.fill(HIST("particlePrm/phi"), particle.phi());
@@ -376,7 +376,7 @@ struct QaTrackingResolution {
       if (pdgCodeSel != 0 && particle.pdgCode() != pdgCodeSel) {
         continue;
       }
-      if (useOnlyPhysicsPrimary && !MC::isPhysicalPrimary<o2::aod::McParticles>(particle)) {
+      if (useOnlyPhysicsPrimary && !MC::isPhysicalPrimary(particle)) {
         continue;
       }
       const double deltaPt = track.pt() - particle.pt();
diff --git a/Analysis/Tutorials/src/mcHistograms.cxx b/Analysis/Tutorials/src/mcHistograms.cxx
index e5ce383a80672..740a202f89c03 100644
--- a/Analysis/Tutorials/src/mcHistograms.cxx
+++ b/Analysis/Tutorials/src/mcHistograms.cxx
@@ -38,14 +38,14 @@ struct AccessMcData {
   OutputObj<TH1F> etaH{TH1F("eta", "eta", 102, -2.01, 2.01)};
 
   // group according to McCollisions
-  void process(aod::McCollision const& mcCollision, aod::McParticles& mcParticles)
+  void process(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles)
   {
     // access MC truth information with mcCollision() and mcParticle() methods
     LOGF(info, "MC. vtx-z = %f", mcCollision.posZ());
     LOGF(info, "First: %d | Length: %d", mcParticles.begin().index(), mcParticles.size());
     int count = 0;
     for (auto& mcParticle : mcParticles) {
-      if (MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) {
+      if (MC::isPhysicalPrimary(mcParticle)) {
         phiH->Fill(mcParticle.phi());
         etaH->Fill(mcParticle.eta());
         count++;
@@ -61,7 +61,8 @@ struct AccessMcTruth {
   OutputObj<TH1F> phiDiff{TH1F("phiDiff", ";phi_{MC} - phi_{Rec}", 100, -M_PI, M_PI)};
 
   // group according to reconstructed Collisions
-  void process(soa::Join<aod::Collisions, aod::McCollisionLabels>::iterator const& collision, soa::Join<aod::Tracks, aod::McTrackLabels> const& tracks, aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions)
+  void process(soa::Join<aod::Collisions, aod::McCollisionLabels>::iterator const& collision, soa::Join<aod::Tracks, aod::McTrackLabels> const& tracks,
+               aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions)
   {
     // access MC truth information with mcCollision() and mcParticle() methods
     LOGF(info, "vtx-z (data) = %f | vtx-z (MC) = %f", collision.posZ(), collision.mcCollision().posZ());
@@ -84,11 +85,34 @@ struct AccessMcTruth {
   }
 };
 
+// Loop over MCColisions and get corresponding collisions (there can be more than one)
+// For each of them get the corresponding tracks
+struct LoopOverMcMatched {
+  OutputObj<TH1F> etaDiff{TH1F("etaDiff", ";eta_{MC} - eta_{Rec}", 100, -2, 2)};
+  void process(aod::McCollision const& mcCollision, soa::Join<aod::McCollisionLabels, aod::Collisions> const& collisions,
+               soa::Join<aod::Tracks, aod::McTrackLabels> const& tracks, aod::McParticles const& mcParticles)
+  {
+    // access MC truth information with mcCollision() and mcParticle() methods
+    LOGF(info, "MC collision at vtx-z = %f with %d mc particles and %d reconstructed collisions", mcCollision.posZ(), mcParticles.size(), collisions.size());
+    for (auto& collision : collisions) {
+      LOGF(info, "  Reconstructed collision at vtx-z = %f", collision.posZ());
+
+      // NOTE this will be replaced by a improved grouping in the future
+      auto groupedTracks = tracks.sliceBy(aod::track::collisionId, collision.globalIndex());
+      LOGF(info, "  which has %d tracks", groupedTracks.size());
+      for (auto& track : groupedTracks) {
+        etaDiff->Fill(track.mcParticle().eta() - track.eta());
+      }
+    }
+  }
+};
+
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {
   return WorkflowSpec{
     adaptAnalysisTask<VertexDistribution>(cfgc),
     adaptAnalysisTask<AccessMcData>(cfgc),
     adaptAnalysisTask<AccessMcTruth>(cfgc),
+    adaptAnalysisTask<LoopOverMcMatched>(cfgc),
   };
 }