diff --git a/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernel.cxx b/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernel.cxx
index 02bcb6db259a99d48d21732e903468217ddcec1a..691265b50f2d0e49727b0605a16e1ec7ec56683a 100644
--- a/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernel.cxx
+++ b/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernel.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 // ISF_Algs includes
@@ -26,6 +26,8 @@
 #include "GeneratorObjects/McEventCollection.h"
 #include "GeneratorObjects/HepMcParticleLink.h"
 
+#undef ISFDEBUG
+
 ///////////////////////////////////////////////////////////////////
 // Public methods:
 ///////////////////////////////////////////////////////////////////
@@ -372,11 +374,17 @@ StatusCode ISF::SimKernel::execute()
 
   // -----------------------------------------------------------------------------------------------
   // Step 3: ISimulation KERNEL : loop over particle stack, until empty
+  unsigned int loopCounter{0};
   ATH_MSG_DEBUG( "Starting simulation loop, initial particle stack size: " << m_particleBroker->numParticles());
-  while ( true ) {
+  while ( m_particleBroker->numParticles() ) {
+    ++loopCounter;
+    ATH_MSG_VERBOSE("Main Loop pass no. " << loopCounter);
+    ATH_MSG_VERBOSE("Queue starts with " << m_particleBroker->numParticles() << " particles.");
     // get next vector of particles for simulation
     const ISF::ConstISFParticleVector &particles = m_particleBroker->popVector(m_maxParticleVectorSize);
+    const unsigned int numParticlesLeftInBroker = m_particleBroker->numParticles();
     int numParticles = particles.size();
+
     // particle vector empty -> end simulation
     if (numParticles==0) break;
 
@@ -391,6 +399,16 @@ StatusCode ISF::SimKernel::execute()
     ATH_MSG_DEBUG  ( "Took " << numParticles << " particles from queue (remaining: " << m_particleBroker->numParticles() << ")" );
     ATH_MSG_VERBOSE( " -> All particles will be sent to '" << m_simSvcNames[simID] << "' simulator (SimSvcID=" << simID << ")"  );
 
+    #ifdef ISFDEBUG
+    if (loopCounter>100 && numParticles<3) {
+      ATH_MSG_INFO("Main Loop pass no. " << loopCounter);
+      ATH_MSG_INFO("Selected " << numParticles << " particles to be processed by " << m_simSvcNames[simID]);
+      for ( const ISFParticle *particle : particles ) {
+        ATH_MSG_INFO(*particle);
+      }
+    }
+    #endif // ISFDEBUG
+
     // ensure that all particles in the vector have the same SimID
     for ( const ISFParticle *particle : particles ) {
       if ( particle->nextSimID() != simID ) {
@@ -409,13 +427,15 @@ StatusCode ISF::SimKernel::execute()
       // NB Passing only the hard-scatter McEventCollection is not
       // correct if Geant4 simulation were to be used for pile-up Hits
       // in Fast Chain.
+      ATH_MSG_VERBOSE("Selected " << particles.size() << " particles to be processed by " <<  m_simSvcNames[simID]);
       if (m_simSvcs[simID]->simulateVector(particles, m_outputHardScatterTruth.ptr()).isFailure()) {
         ATH_MSG_WARNING( "Simulation of particles failed in Simulator: " << m_simSvcNames[simID]);
       }
+      ATH_MSG_VERBOSE(m_simSvcNames[simID] << " returned " << m_particleBroker->numParticles()-numParticlesLeftInBroker << " new particles to be added to the queue." );
     }
 
   }
-  ATH_MSG_VERBOSE( "Ending simulation loop, no more particles in the stack" );
+  ATH_MSG_VERBOSE("Final status: queue contains " << m_particleBroker->numParticles() << " particles.");
   // -----------------------------------------------------------------------------------------------
 
 
diff --git a/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernelMT.cxx b/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernelMT.cxx
index 58697b51e4946d739c41f82bd7dc68a4130a4aaa..a95ea80a4f5c23747dbadd8f4bf61b1fae405333 100644
--- a/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernelMT.cxx
+++ b/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernelMT.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 /**
@@ -18,7 +18,9 @@
 
 // STL
 #include <queue>
+#include <utility>
 
+#undef ISFDEBUG
 
 ISF::SimKernelMT::SimKernelMT( const std::string& name, ISvcLocator* pSvcLocator ) :
     ::AthAlgorithm( name, pSvcLocator ),
@@ -30,6 +32,8 @@ ISF::SimKernelMT::SimKernelMT( const std::string& name, ISvcLocator* pSvcLocator
     declareProperty("CaloSimulationSelectors", m_simSelectors[AtlasDetDescr::fAtlasCalo] );
     declareProperty("MSSimulationSelectors", m_simSelectors[AtlasDetDescr::fAtlasMS] );
     declareProperty("CavernSimulationSelectors", m_simSelectors[AtlasDetDescr::fAtlasCavern] );
+    // tuning parameters
+    declareProperty("MaximumParticleVectorSize"  , m_maxParticleVectorSize             );
 }
 
 
@@ -126,7 +130,7 @@ StatusCode ISF::SimKernelMT::initialize() {
 
 StatusCode ISF::SimKernelMT::execute() {
 
-  // Release the event from all simulators (TODO: make the tools do this)
+  // Call setupEvent for all simulators (TODO: make the tools do this)
   for (auto& curSimTool: m_simulationTools) {
     if ( curSimTool ) {
       ATH_CHECK(curSimTool->setupEvent());
@@ -183,13 +187,17 @@ StatusCode ISF::SimKernelMT::execute() {
     particleQueue.push( particle );
   }
 
+  unsigned int loopCounter{0};
   // loop until there are no more particles to simulate
   ISF::ConstISFParticleVector particles{};
   const ISimulatorTool* lastSimulator{};
   ISFParticleContainer newSecondaries{};
   while ( particleQueue.size() ) {
-
+    ++loopCounter;
+    ATH_MSG_VERBOSE("Main Loop pass no. " << loopCounter);
+    ATH_MSG_VERBOSE("Queue starts with " << particleQueue.size() << " particles.");
     // Create a vector of particles with the same simulator
+    ISFParticleOrderedQueue tempQueue;
     while ( particleQueue.size() ) {
       auto particlePtr = particleQueue.top();
       ISFParticle& curParticle( *particlePtr );
@@ -207,19 +215,29 @@ StatusCode ISF::SimKernelMT::execute() {
         particles.push_back(particlePtr);
         lastSimulator=&simTool;
       }
-      else if (&simTool==lastSimulator) {
-        particles.push_back(particlePtr);
+      else if (&simTool!=lastSimulator || particles.size() >= m_maxParticleVectorSize ) {
+        // Change of simulator, end the current vector
+        tempQueue.push(particlePtr);
       }
       else {
-        // Change of simulator, end the current vector
-        particleQueue.push(particlePtr);
-        break;
+        particles.push_back(particlePtr);
+      }
+    }
+    particleQueue = std::move(tempQueue);
+    #ifdef ISFDEBUG
+    if (loopCounter>100 && particles.size()<3) {
+      ATH_MSG_INFO("Main Loop pass no. " << loopCounter);
+      ATH_MSG_INFO("Selected " << particles.size() << " particles to be processed by " << lastSimulator->name());
+      for ( const ISFParticle *particle : particles ) {
+        ATH_MSG_INFO(*particle);
       }
     }
+    #endif // ISFDEBUG
 
+    ATH_MSG_VERBOSE("Selected " << particles.size() << " particles to be processed by " << lastSimulator->name());
     // Run the simulation
     ATH_CHECK( lastSimulator->simulateVector( particles, newSecondaries, outputTruth.ptr() ) );
-
+    ATH_MSG_VERBOSE(lastSimulator->name() << " returned " << newSecondaries.size() << " new particles to be added to the queue." );
     // Register returned particles with the entry layer tool, set their order and enqueue them
     for ( auto* secondary : newSecondaries ) {
       m_entryLayerTool->registerParticle( *secondary );
@@ -242,6 +260,7 @@ StatusCode ISF::SimKernelMT::execute() {
     }
     particles.clear();
   }
+  ATH_MSG_VERBOSE("Final status: queue contains " << particleQueue.size() << " particles.");
 
   // Release the event from all simulators (TODO: make the tools do this)
   for (auto& curSimTool: m_simulationTools) {
diff --git a/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernelMT.h b/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernelMT.h
index ad492636dd1830b663b0aaf05eebc50846116d72..61b832d9ac1a85d7d60be75c0d900593ab3e5cd7 100644
--- a/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernelMT.h
+++ b/Simulation/ISF/ISF_Core/ISF_Algorithms/src/SimKernelMT.h
@@ -127,6 +127,8 @@ private:
   /// Map of the simulation flavours used in this job to the corresponding Simulation Services
   std::map<ISF::SimulationFlavor, ISimulatorTool*> m_simToolMap;
 
+  /// Number of particles simultaneously sent to simulator
+  size_t m_maxParticleVectorSize{10240};
 };
 
 } // namespace ISF
diff --git a/Simulation/ISF/ISF_Core/ISF_Services/src/ParticleBrokerDynamicOnReadIn.cxx b/Simulation/ISF/ISF_Core/ISF_Services/src/ParticleBrokerDynamicOnReadIn.cxx
index 652fe01d7b29e44743c4fb0b21e8cff6dd3ab90b..33b0d73167d2691fdd27f28a3d02e3c9c6d7d711 100644
--- a/Simulation/ISF/ISF_Core/ISF_Services/src/ParticleBrokerDynamicOnReadIn.cxx
+++ b/Simulation/ISF/ISF_Core/ISF_Services/src/ParticleBrokerDynamicOnReadIn.cxx
@@ -30,6 +30,8 @@
 // ROOT includes
 #include "TTree.h"
 
+#include <utility>
+
 // C include
 #include <assert.h>
 
@@ -485,7 +487,7 @@ const ISF::ConstISFParticleVector& ISF::ParticleBrokerDynamicOnReadIn::popVector
   if ( m_particles.size() ) {
 
     SimSvcID returnID = m_particles.top()->nextSimID();
-
+    ISFParticleOrderedQueue tempQueue;
     // loop as long as we have particles in the m_particles queue
     do {
       // get the next particle from the ordered queue
@@ -494,14 +496,17 @@ const ISF::ConstISFParticleVector& ISF::ParticleBrokerDynamicOnReadIn::popVector
 
       // if this particle has a different SimID, or the maximum size of the return vector is reached
       //   -> don't add any more particles to the m_popParticles std::vector
-      if ( curID != returnID || m_popParticles.size() >= maxVectorSize ) break;
-
-      // add this particle to the, later returned, m_popParticles std::vector
-      m_popParticles.push_back( curParticle );
+      if ( curID != returnID || m_popParticles.size() >= maxVectorSize ) {
+        tempQueue.push(m_particles.top()); //break;
+      }
+      else {
+        // add this particle to the, later returned, m_popParticles std::vector
+        m_popParticles.push_back( curParticle );
+      }
       // remove this particle from the ordered queue
       m_particles.pop();
     } while ( m_particles.size() ) ;
-
+    m_particles = std::move(tempQueue);
   }
   // return the popParticles vector
   return m_popParticles;