diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertexExt.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertexExt.h
index a228984f2ae5d..4fd36ed248677 100644
--- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertexExt.h
+++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertexExt.h
@@ -14,6 +14,7 @@
 
 #include "ReconstructionDataFormats/PrimaryVertex.h"
 #include "ReconstructionDataFormats/GlobalTrackID.h"
+#include "SimulationDataFormat/MCEventLabel.h"
 
 namespace o2
 {
@@ -27,6 +28,7 @@ struct PrimaryVertexExt : public PrimaryVertex {
   std::array<uint16_t, o2::dataformats::GlobalTrackID::Source::NSources> nSrc{};   // N contributors for each source type
   std::array<uint16_t, o2::dataformats::GlobalTrackID::Source::NSources> nSrcA{};  // N associated and passing cuts for each source type
   std::array<uint16_t, o2::dataformats::GlobalTrackID::Source::NSources> nSrcAU{}; // N ambgous associated and passing cuts for each source type
+  o2::MCEventLabel mcLb{};
   double FT0Time = -1.;                                                            // time of closest FT0 trigger
   float FT0A = -1;                                                               // amplitude of closest FT0 A side
   float FT0C = -1;                                                               // amplitude of closest FT0 C side
@@ -41,7 +43,7 @@ struct PrimaryVertexExt : public PrimaryVertex {
   std::string asString() const;
 #endif
 
-  ClassDefNV(PrimaryVertexExt, 6);
+  ClassDefNV(PrimaryVertexExt, 7);
 };
 
 #ifndef GPUCA_GPUCODE_DEVICE
diff --git a/Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx b/Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx
index 881ce9041ae04..38079d36a2e84 100644
--- a/Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx
+++ b/Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx
@@ -369,6 +369,9 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
     if (iv != nv - 1) {
       auto& pve = pveVec[iv];
       static_cast<o2::dataformats::PrimaryVertex&>(pve) = pvvec[iv];
+      if (mUseMC) {
+        pve.mcLb = recoData.getPrimaryVertexMCLabel(iv);
+      }
       // find best matching FT0 signal
       float bestTimeDiff = 1000, bestTime = -999;
       int bestFTID = -1;
@@ -611,13 +614,16 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
     vid[slot] = id;
   };
 
+  std::vector<o2::dataformats::PrimaryVertexExt> pveT; // neighbours in time
+  std::vector<o2::dataformats::PrimaryVertexExt> pveZ; // neighbours in Z
+  std::vector<int> idT(mMaxNeighbours), idZ(mMaxNeighbours);
+  std::vector<float> dT(mMaxNeighbours), dZ(mMaxNeighbours);
   for (int cnt = 0; cnt < nvtot; cnt++) {
+    pveT.clear(); // neighbours in time
+    pveZ.clear(); // neighbours in Z
+
     const auto& pve = pveVec[cnt];
     float tv = pve.getTimeStamp().getTimeStamp();
-    std::vector<o2::dataformats::PrimaryVertexExt> pveT(mMaxNeighbours); // neighbours in time
-    std::vector<o2::dataformats::PrimaryVertexExt> pveZ(mMaxNeighbours); // neighbours in Z
-    std::vector<int> idT(mMaxNeighbours), idZ(mMaxNeighbours);
-    std::vector<float> dT(mMaxNeighbours), dZ(mMaxNeighbours);
     for (int i = 0; i < mMaxNeighbours; i++) {
       idT[i] = idZ[i] = -1;
       dT[i] = mMaxVTTimeDiff;
@@ -666,14 +672,14 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
     }
     for (int i = 0; i < mMaxNeighbours; i++) {
       if (idT[i] != -1) {
-        pveT[i] = pveVec[idT[i]];
+        pveT.push_back(pveVec[idT[i]]);
       } else {
         break;
       }
     }
     for (int i = 0; i < mMaxNeighbours; i++) {
       if (idZ[i] != -1) {
-        pveZ[i] = pveVec[idZ[i]];
+        pveZ.push_back(pveVec[idZ[i]]);
       } else {
         break;
       }