diff --git a/SimG4Core/Application/interface/RunManagerMT.h b/SimG4Core/Application/interface/RunManagerMT.h
index c07d6121e994b..3c77481355473 100644
--- a/SimG4Core/Application/interface/RunManagerMT.h
+++ b/SimG4Core/Application/interface/RunManagerMT.h
@@ -80,6 +80,10 @@ class RunManagerMT {
 private:
   void terminateRun();
 
+  void checkVoxels();
+
+  void setupVoxels();
+
   G4MTRunManagerKernel* m_kernel;
 
   CustomUIsession* m_UIsession;
diff --git a/SimG4Core/Application/python/g4SimHits_cfi.py b/SimG4Core/Application/python/g4SimHits_cfi.py
index ded467d761428..3a8ede10cf298 100644
--- a/SimG4Core/Application/python/g4SimHits_cfi.py
+++ b/SimG4Core/Application/python/g4SimHits_cfi.py
@@ -90,6 +90,9 @@
     UseParametrisedEMPhysics = cms.untracked.bool(True),
     ThresholdForGeometryExceptions = cms.double(0.1), ## in GeV
     TraceExceptions = cms.bool(False),
+    DefaultVoxelDensity = cms.double(2.0),
+    VoxelRegions = cms.vstring(),
+    VoxelDensityPerRegion = cms.vdouble(),
     CheckGeometry = cms.untracked.bool(False),
     OnlySDs = cms.vstring('ZdcSensitiveDetector', 'TotemT2ScintSensitiveDetector', 'TotemSensitiveDetector', 'RomanPotSensitiveDetector', 'PLTSensitiveDetector', 'MuonSensitiveDetector', 'MtdSensitiveDetector', 'BCM1FSensitiveDetector', 'EcalSensitiveDetector', 'CTPPSSensitiveDetector', 'BSCSensitiveDetector', 'CTPPSDiamondSensitiveDetector', 'FP420SensitiveDetector', 'BHMSensitiveDetector', 'CastorSensitiveDetector', 'CaloTrkProcessing', 'HcalSensitiveDetector', 'TkAccumulatingSensitiveDetector'),
     G4CheckOverlap = cms.untracked.PSet(
diff --git a/SimG4Core/Application/src/RunManagerMT.cc b/SimG4Core/Application/src/RunManagerMT.cc
index 222a88b2988a7..d1b779c4f56e9 100644
--- a/SimG4Core/Application/src/RunManagerMT.cc
+++ b/SimG4Core/Application/src/RunManagerMT.cc
@@ -202,6 +202,8 @@ void RunManagerMT::initG4(const DDCompactView* pDD,
     }
   }
 
+  setupVoxels();
+
   m_stateManager->SetNewState(G4State_Init);
   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT: G4State is Init";
   m_kernel->InitializePhysics();
@@ -216,6 +218,10 @@ void RunManagerMT::initG4(const DDCompactView* pDD,
     throw cms::Exception("LogicError") << "G4RunManagerKernel initialization failed!";
   }
 
+  if (m_check) {
+    checkVoxels();
+  }
+
   if (m_StorePhysicsTables) {
     std::ostringstream dir;
     dir << m_PhysicsTablesDir << '\0';
@@ -301,3 +307,55 @@ void RunManagerMT::terminateRun() {
   m_runTerminated = true;
   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT::terminateRun done";
 }
+
+void RunManagerMT::checkVoxels() {
+  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
+  int numLV = lvs->size();
+  edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT: nLV=" << numLV;
+  int nvox = 0;
+  int nslice = 0;
+  for (int i = 0; i < numLV; ++i) {
+    auto lv = (*lvs)[i];
+    auto nd = lv->GetNoDaughters();
+    auto vox = lv->GetVoxelHeader();
+    auto sma = lv->GetSmartless();
+    auto reg = lv->GetRegion();
+    size_t nsl = (nullptr == vox) ? 0 : vox->GetNoSlices();
+    if (0 < nsl) {
+      nslice += nsl;
+      std::string rname = (nullptr != reg) ? reg->GetName() : "";
+      edm::LogVerbatim("Voxels") << " " << i << ". Nd=" << nd << " Nsl=" << nsl << " Smartless=" << sma << " "
+                                 << lv->GetName() << " Region: " << rname;
+      ++nvox;
+    }
+  }
+  edm::LogVerbatim("SimG4CoreApplication")
+      << "RunManagerMT: nLV=" << numLV << " NlvVox=" << nvox << " Nslices=" << nslice;
+}
+
+void RunManagerMT::setupVoxels() {
+  double density = m_p.getParameter<double>("DefaultVoxelDensity");
+  std::vector<std::string> rnames = m_p.getParameter<std::vector<std::string> >("VoxelRegions");
+  std::vector<double> rdensities = m_p.getParameter<std::vector<double> >("VoxelDensityPerRegion");
+  int nr = 0;
+  std::size_t n = rnames.size();
+  if (n == rdensities.size()) {
+    nr = (int)n;
+  }
+  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
+  for (auto& lv : *lvs) {
+    double den = density;
+    if (0 < nr) {
+      std::string nam = lv->GetRegion()->GetName();
+      for (int i = 0; i < nr; ++i) {
+        if (nam == rnames[i]) {
+          den = rdensities[i];
+          break;
+        }
+      }
+    }
+    lv->SetSmartless(den);
+  }
+  edm::LogVerbatim("SimG4CoreApplication")
+      << "RunManagerMT: default voxel density=" << density << "; number of regions with special density " << nr;
+}