diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
index 7024b858a1a6bd0a29b07e0465a0dfb83280c57a..2ba97d9cfcb5ad4a85ae66cddd3cfab82a24626d 100644
--- a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
+++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "PixelConfigCondAlg.h"
@@ -458,12 +458,11 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     }
   }
 
-  std::vector<TH3F*> ramoPotentialMap;
-  std::vector<TH2F*> lorentzMap_e;
-  std::vector<TH2F*> lorentzMap_h;
-  std::vector<TH2F*> distanceMap_e;
-  std::vector<TH2F*> distanceMap_h;
-
+  std::vector<PixelHistoConverter> ramoPotentialMap;
+  std::vector<PixelHistoConverter> lorentzMap_e;
+  std::vector<PixelHistoConverter> lorentzMap_h;
+  std::vector<PixelHistoConverter> distanceMap_e;
+  std::vector<PixelHistoConverter> distanceMap_h;
   for (unsigned int i=0; i<mapsPath_list.size(); i++) {
     ATH_MSG_INFO("Using maps located in: "<<mapsPath_list.at(i) << " for layer No." << i);
     TFile* mapsFile = new TFile((mapsPath_list.at(i)).c_str()); //this is the ramo potential.
@@ -476,12 +475,17 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
       ATH_MSG_FATAL("Did not find a Ramo potential map and an approximate form is available yet. Exit...");
       return StatusCode::FAILURE;
     }
-    ramoPotentialMap.push_back(ramoPotentialMap_hold);
 
-    lorentzMap_e.push_back((TH2F*)mapsFile->Get("lorentz_map_e"));
-    lorentzMap_h.push_back((TH2F*)mapsFile->Get("lorentz_map_h"));
-    distanceMap_e.push_back((TH2F*)mapsFile->Get("edistance"));
-    distanceMap_h.push_back((TH2F*)mapsFile->Get("hdistance"));
+    ramoPotentialMap.emplace_back();
+    ATH_CHECK(ramoPotentialMap.back().setHisto3D(ramoPotentialMap_hold));
+    lorentzMap_e.emplace_back();
+    lorentzMap_h.emplace_back();
+    distanceMap_e.emplace_back();
+    distanceMap_h.emplace_back();
+    ATH_CHECK(lorentzMap_e.back().setHisto2D((TH2F*)mapsFile->Get("lorentz_map_e")));
+    ATH_CHECK(lorentzMap_h.back().setHisto2D((TH2F*)mapsFile->Get("lorentz_map_h")));
+    ATH_CHECK(distanceMap_e.back().setHisto2D((TH2F*)mapsFile->Get("edistance")));
+    ATH_CHECK(distanceMap_h.back().setHisto2D((TH2F*)mapsFile->Get("hdistance")));
   }
   writeCdo -> setLorentzMap_e(lorentzMap_e);
   writeCdo -> setLorentzMap_h(lorentzMap_h);
diff --git a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelHistoConverter.h b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelHistoConverter.h
new file mode 100644
index 0000000000000000000000000000000000000000..13e03ce4813e87b0d7d8a9b0ef61c10622bb3a44
--- /dev/null
+++ b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelHistoConverter.h
@@ -0,0 +1,77 @@
+/*
+   Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+ */
+/**
+ * @file PixelDigitization/PixelHistoConverter.h
+ * @author Tomas Dado <tomas.dado@cern.ch>
+ * @date February, 2021
+ * @brief A helper class to conver ROOT TH* objects into simple vectors for speed up
+ */
+
+#ifndef PIXELCONDITIONSDATA_PIXELHISTOCONVRTER_H
+#define PIXELCONDITIONSDATA_PIXELHISTOCONVRTER_H
+
+#include "AthenaBaseComps/AthMessaging.h"
+#include "GaudiKernel/StatusCode.h"
+
+#include <vector>
+
+/// forward class declaration
+class TAxis;
+class TH1;
+class TH2;
+class TH3;
+
+class PixelHistoConverter{
+public:
+  PixelHistoConverter();
+  virtual ~PixelHistoConverter() = default;
+
+  StatusCode setHisto1D(const TH1* histo);
+  StatusCode setHisto2D(const TH2* histo);
+  StatusCode setHisto3D(const TH3* histo);
+
+  inline float getContent(const std::size_t& x) const {
+    return m_content[x];
+  }
+  inline float getContent(const std::size_t& x, const std::size_t& y) const {
+    return m_content[x + y*(m_xAxis.nBins)];
+  }
+  inline float getContent(const std::size_t& x, const std::size_t& y, const std::size_t& z) const {
+    return m_content[x + m_xAxis.nBins*(y + (m_yAxis.nBins * z))];
+  }
+
+  inline bool isOverflowZ(const float value) const {
+    return (value >= m_zAxis.max);
+  }
+  bool isFirstZ(const float value) const;
+  float getBinX(const float value) const;
+  float getBinY(const float value) const;
+  float getBinZ(const float value) const;
+
+private:
+  struct Axis {
+    std::size_t nBins;
+    float min;
+    float max;
+    float width;
+  };
+
+  Axis m_xAxis;
+  Axis m_yAxis;
+  Axis m_zAxis;
+
+  std::vector<float> m_content;
+
+  bool setAxis(Axis& axis, const TAxis* rootAxis);
+
+  inline std::size_t findBin(const Axis& axis, const float value) const {
+    if (value <= axis.min) return 0;
+    if (value >= axis.max) return (axis.nBins - 1);
+
+    return ((value - axis.min) * axis.width);
+  }
+
+};
+
+#endif
diff --git a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h
index e1efb342b00c5038223ed7eb4917f80083297610..4d9fa5c2fad7ed4734295c0193027ecb82a1a6ee 100644
--- a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h
+++ b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
 */
 /**
  * @file PixelConditionsData/PixelModuleData.h
@@ -15,6 +15,7 @@
 #include <map>
 
 #include "AthenaKernel/CondCont.h"
+#include "PixelConditionsData/PixelHistoConverter.h"
 #include "TH1.h"
 #include "TH2.h"
 #include "TH3.h"
@@ -167,18 +168,18 @@ class PixelModuleData {
     void setFluenceLayer(std::vector<double> fluenceLayer);
     double getFluenceLayer(int layer) const;
 
-    void setLorentzMap_e(std::vector<TH2F*> lorentzMap_e);
-    void setLorentzMap_h(std::vector<TH2F*> lorentzMap_h);
-    TH2F* getLorentzMap_e(int layer) const;
-    TH2F* getLorentzMap_h(int layer) const;
+    void setLorentzMap_e(std::vector<PixelHistoConverter> lorentzMap_e);
+    void setLorentzMap_h(std::vector<PixelHistoConverter> lorentzMap_h);
+    const PixelHistoConverter& getLorentzMap_e(int layer) const;
+    const PixelHistoConverter& getLorentzMap_h(int layer) const;
 
-    void setDistanceMap_e(std::vector<TH2F*> distanceMap_e);
-    void setDistanceMap_h(std::vector<TH2F*> distanceMap_h);
-    TH2F* getDistanceMap_e(int layer) const;
-    TH2F* getDistanceMap_h(int layer) const;
+    void setDistanceMap_e(std::vector<PixelHistoConverter> distanceMap_e);
+    void setDistanceMap_h(std::vector<PixelHistoConverter> distanceMap_h);
+    const PixelHistoConverter& getDistanceMap_e(int layer) const;
+    const PixelHistoConverter& getDistanceMap_h(int layer) const;
 
-    void setRamoPotentialMap(std::vector<TH3F*> ramoPotentialMap);
-    TH3F* getRamoPotentialMap(int layer) const;
+    void setRamoPotentialMap(std::vector<PixelHistoConverter> ramoPotentialMap);
+    const PixelHistoConverter& getRamoPotentialMap(int layer) const;
 
     // Distortion parameters
     void setDistortionInputSource(int distortionInputSource);
@@ -302,11 +303,11 @@ class PixelModuleData {
     std::string m_cablingMapFileName;
 
     std::vector<double> m_fluenceLayer;
-    std::vector<TH2F*> m_lorentzMap_e;
-    std::vector<TH2F*> m_lorentzMap_h;
-    std::vector<TH2F*> m_distanceMap_e;
-    std::vector<TH2F*> m_distanceMap_h;
-    std::vector<TH3F*> m_ramoPotentialMap;
+    std::vector<PixelHistoConverter> m_lorentzMap_e;
+    std::vector<PixelHistoConverter> m_lorentzMap_h;
+    std::vector<PixelHistoConverter> m_distanceMap_e;
+    std::vector<PixelHistoConverter> m_distanceMap_h;
+    std::vector<PixelHistoConverter> m_ramoPotentialMap;
 
     int    m_distortionInputSource;
     int    m_distortionVersion;
diff --git a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelHistoConverter.cxx b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelHistoConverter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..70ffe84b4d151b5dedc9bb8aa604047f71c5b075
--- /dev/null
+++ b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelHistoConverter.cxx
@@ -0,0 +1,145 @@
+/*
+   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+#include "PixelConditionsData/PixelHistoConverter.h"
+
+#include "StoreGate/StoreGateSvc.h"
+#include "GaudiKernel/ISvcLocator.h"
+
+#include "TAxis.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+
+#include <iostream>
+
+PixelHistoConverter::PixelHistoConverter()
+{
+}
+
+StatusCode PixelHistoConverter::setHisto1D(const TH1* histo) {
+
+  if (!histo) {
+    return StatusCode::FAILURE;
+  }
+
+  if (!setAxis(m_xAxis, histo->GetXaxis())) {
+    return StatusCode::FAILURE;
+  }
+
+  /// fill the content
+  const std::size_t xSize = m_xAxis.nBins;
+  m_content.resize(xSize);
+ 
+  for (std::size_t x = 0; x < xSize; ++x) {
+    m_content.at(x) = histo->GetBinContent(x+1);
+  }
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode PixelHistoConverter::setHisto2D(const TH2* histo) {
+  if (!histo) {
+    return StatusCode::FAILURE;
+  }
+
+  if (!setAxis(m_xAxis, histo->GetXaxis())) {
+    return StatusCode::FAILURE;
+  }
+  if (!setAxis(m_yAxis, histo->GetYaxis())) {
+    return StatusCode::FAILURE;
+  }
+
+  /// fill the content use linearized version for performance reasons
+  const std::size_t xSize = m_xAxis.nBins;
+  const std::size_t ySize = m_yAxis.nBins;
+  m_content.resize(xSize*ySize);
+ 
+  for (std::size_t x = 0; x < xSize; ++x) {
+    for (std::size_t y = 0; y < ySize; ++y) {
+      const std::size_t position = x + y*xSize;
+      m_content.at(position) = histo->GetBinContent(x+1,y+1);
+    }
+  }
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode PixelHistoConverter::setHisto3D(const TH3* histo) {
+  if (!histo) {
+    return StatusCode::FAILURE;
+  }
+
+  if (!setAxis(m_xAxis, histo->GetXaxis())) {
+    return StatusCode::FAILURE;
+  }
+  if (!setAxis(m_yAxis, histo->GetYaxis())) {
+    return StatusCode::FAILURE;
+  }
+  if (!setAxis(m_zAxis, histo->GetZaxis())) {
+    return StatusCode::FAILURE;
+  }
+
+  /// fill the content use linearized version for performance reasons
+  const std::size_t xSize = m_xAxis.nBins;
+  const std::size_t ySize = m_yAxis.nBins;
+  const std::size_t zSize = m_zAxis.nBins;
+  m_content.resize(xSize*ySize*zSize);
+ 
+  for (std::size_t x = 0; x < xSize; ++x) {
+    for (std::size_t y = 0; y < ySize; ++y) {
+      for (std::size_t z = 0; z < zSize; ++z) {
+        const std::size_t position = x + xSize*(y + (ySize * z));
+        m_content.at(position) = histo->GetBinContent(x+1,y+1,z+1);
+      }
+    }
+  }
+
+  return StatusCode::SUCCESS;
+}
+
+bool PixelHistoConverter::isFirstZ(const float value) const {
+  return (getBinZ(value) == 0);
+}
+
+float PixelHistoConverter::getBinX(const float value) const {
+  return findBin(m_xAxis, value);
+}
+  
+float PixelHistoConverter::getBinY(const float value) const {
+  return findBin(m_yAxis, value);
+}
+
+float PixelHistoConverter::getBinZ(const float value) const {
+  return findBin(m_zAxis, value);
+}
+
+bool PixelHistoConverter::setAxis(Axis& axis, const TAxis* rootAxis) {
+
+  if (!rootAxis) {
+    return false;
+  }
+
+  axis.nBins = rootAxis->GetNbins();
+  axis.min   = rootAxis->GetXmin();
+  axis.max   = rootAxis->GetXmax();
+
+  if (axis.nBins < 1) {
+    return false;
+  }
+
+  /// check if the histogram has equidistant bins
+  const float width = rootAxis->GetBinWidth(1);
+  for (std::size_t ibin = 2; ibin <= axis.nBins; ++ibin) {
+    /// use a threshold for imperfect binning
+    if (std::abs(rootAxis->GetBinWidth(ibin) - width) > 0.01*width) {
+      return false;
+    }
+  }
+
+  /// storing as 1/width to avoid (slow) division in retrieving
+  axis.width = 1.*axis.nBins/(axis.max - axis.min);
+
+  return true;
+}
diff --git a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
index 16da803924fcd45e2b570b9451e491b4fdcceaea..b815435be8486b168abe3aa5f2664ea6e0add5ac 100644
--- a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
+++ b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "PixelConditionsData/PixelModuleData.h"
@@ -319,19 +319,18 @@ std::string PixelModuleData::getCablingMapFileName() const { return m_cablingMap
 void PixelModuleData::setFluenceLayer(std::vector<double> fluenceLayer) { m_fluenceLayer = fluenceLayer; }
 double PixelModuleData::getFluenceLayer(int layer) const { return m_fluenceLayer.at(layer); }
 
-void PixelModuleData::setLorentzMap_e(std::vector<TH2F*> lorentzMap_e) { m_lorentzMap_e = lorentzMap_e; }
-void PixelModuleData::setLorentzMap_h(std::vector<TH2F*> lorentzMap_h) { m_lorentzMap_h = lorentzMap_h; }
-TH2F* PixelModuleData::getLorentzMap_e(int layer) const { return m_lorentzMap_e.at(layer); }
-TH2F* PixelModuleData::getLorentzMap_h(int layer) const { return m_lorentzMap_h.at(layer); }
+void PixelModuleData::setLorentzMap_e(std::vector<PixelHistoConverter> lorentzMap_e) { m_lorentzMap_e = lorentzMap_e; }
+void PixelModuleData::setLorentzMap_h(std::vector<PixelHistoConverter> lorentzMap_h) { m_lorentzMap_h = lorentzMap_h; }
+const PixelHistoConverter& PixelModuleData::getLorentzMap_e(int layer) const { return m_lorentzMap_e.at(layer); }
+const PixelHistoConverter& PixelModuleData::getLorentzMap_h(int layer) const { return m_lorentzMap_h.at(layer); }
 
-void PixelModuleData::setDistanceMap_e(std::vector<TH2F*> distanceMap_e) { m_distanceMap_e = distanceMap_e; }
-void PixelModuleData::setDistanceMap_h(std::vector<TH2F*> distanceMap_h) { m_distanceMap_h = distanceMap_h; }
-TH2F* PixelModuleData::getDistanceMap_e(int layer) const { return m_distanceMap_e.at(layer); }
-TH2F* PixelModuleData::getDistanceMap_h(int layer) const { return m_distanceMap_h.at(layer); }
-
-void PixelModuleData::setRamoPotentialMap(std::vector<TH3F*> ramoPotentialMap) { m_ramoPotentialMap = ramoPotentialMap; }
-TH3F* PixelModuleData::getRamoPotentialMap(int layer) const { return m_ramoPotentialMap.at(layer); }
+void PixelModuleData::setDistanceMap_e(std::vector<PixelHistoConverter> distanceMap_e) { m_distanceMap_e = distanceMap_e; }
+void PixelModuleData::setDistanceMap_h(std::vector<PixelHistoConverter> distanceMap_h) { m_distanceMap_h = distanceMap_h; }
+const PixelHistoConverter& PixelModuleData::getDistanceMap_e(int layer) const { return m_distanceMap_e.at(layer); }
+const PixelHistoConverter& PixelModuleData::getDistanceMap_h(int layer) const { return m_distanceMap_h.at(layer); }
 
+void PixelModuleData::setRamoPotentialMap(std::vector<PixelHistoConverter> ramoPotentialMap) { m_ramoPotentialMap = ramoPotentialMap; }
+const PixelHistoConverter& PixelModuleData::getRamoPotentialMap(int layer) const { return m_ramoPotentialMap.at(layer); }
 // Distortion parameters
 void PixelModuleData::setDistortionInputSource(int distortionInputSource) { m_distortionInputSource = distortionInputSource; }
 int PixelModuleData::getDistortionInputSource() const { return m_distortionInputSource; }
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx
index 7e1b009b4afd5a985d56568952a6424b070cc9b8..6e36ec275af73a69a9bb84c73302ce886dcc12aa 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+   Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "SensorSim3DTool.h"
@@ -105,8 +105,8 @@ StatusCode SensorSim3DTool::initialize() {
       return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up
     }
     //ramoPotentialMap.push_back(ramoPotentialMap_hold);
-    m_ramoPotentialMap[Layer] = ramoPotentialMap_hold;
-    m_fluence_layersMaps[Layer] = m_fluence_layers.at(i);
+    m_ramoPotentialMap.emplace_back();
+    ATH_CHECK(m_ramoPotentialMap.back().setHisto3D(ramoPotentialMap_hold));
 
     ATH_MSG_INFO("Using fluence " << m_fluence_layers.at(i));
 
@@ -118,7 +118,8 @@ StatusCode SensorSim3DTool::initialize() {
       ATH_MSG_INFO("Unable to load sensor e-field map, so generating one using approximations.");
       return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up
     }
-    m_eFieldMap[Layer] = eFieldMap_hold;
+    m_eFieldMap.emplace_back();
+    ATH_CHECK(m_eFieldMap.back().setHisto2D(eFieldMap_hold));
 
     TH3F* xPositionMap_e_hold;
     TH3F* xPositionMap_h_hold;
@@ -140,26 +141,22 @@ StatusCode SensorSim3DTool::initialize() {
     timeMap_e_hold = (TH2F* ) mapsFile->Get("etimes");
     timeMap_h_hold = (TH2F* ) mapsFile->Get("htimes");
     //Now, determine the time to reach the electrode and the trapping position.
-    if (xPositionMap_e_hold == 0 || xPositionMap_h_hold == 0 || yPositionMap_e_hold == 0 || yPositionMap_h_hold == 0 ||
-        timeMap_e_hold == 0 || timeMap_h_hold == 0) {
-      ATH_MSG_INFO("Unable to load at least one of teh distance/time maps, so generating all using approximations.");
-      return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up
-    }
-    m_xPositionMap_e[Layer] = xPositionMap_e_hold;
-    m_xPositionMap_h[Layer] = xPositionMap_h_hold;
-    m_yPositionMap_e[Layer] = yPositionMap_e_hold;
-    m_yPositionMap_h[Layer] = yPositionMap_h_hold;
-    m_timeMap_e[Layer] = timeMap_e_hold;
-    m_timeMap_h[Layer] = timeMap_h_hold;
+    m_xPositionMap_e.emplace_back();
+    m_xPositionMap_h.emplace_back();
+    m_yPositionMap_e.emplace_back();
+    m_yPositionMap_h.emplace_back();
+    ATH_CHECK(m_xPositionMap_e.back().setHisto3D(xPositionMap_e_hold));
+    ATH_CHECK(m_xPositionMap_h.back().setHisto3D(xPositionMap_h_hold));
+    ATH_CHECK(m_yPositionMap_e.back().setHisto3D(yPositionMap_e_hold));
+    ATH_CHECK(m_yPositionMap_h.back().setHisto3D(yPositionMap_h_hold));
+    m_timeMap_e.emplace_back();
+    m_timeMap_h.emplace_back();
+    ATH_CHECK(m_timeMap_e.back().setHisto2D(timeMap_e_hold));
+    ATH_CHECK(m_timeMap_h.back().setHisto2D(timeMap_h_hold));
 
     // Get average charge data (for charge chunk effect correction)
-    m_avgChargeMap_e = 0;
-    m_avgChargeMap_h = 0;
-    m_avgChargeMap_e = (TH2F*) mapsFile->Get("avgCharge_e");
-    m_avgChargeMap_h = (TH2F*) mapsFile->Get("avgCharge_h");
-    if (m_avgChargeMap_e == 0 || m_avgChargeMap_h == 0) {
-      ATH_MSG_INFO("Unsuccessful picking up histogram: m_avgChargeMap");
-    }
+    ATH_CHECK(m_avgChargeMap_e.setHisto2D((TH2F*) mapsFile->Get("avgCharge_e")));
+    ATH_CHECK(m_avgChargeMap_h.setHisto2D((TH2F*) mapsFile->Get("avgCharge_h")));
   }
 
   return StatusCode::SUCCESS;
@@ -359,13 +356,8 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
 
             ATH_MSG_DEBUG(" -- diffused position w.r.t. pixel edge = " << xposDiff << "  " << yposDiff);
 
-            double average_charge = m_avgChargeMap_e->GetBinContent(m_avgChargeMap_e->GetYaxis()->FindBin(
-                                                                      y_pix * 1000),
-                                                                    m_avgChargeMap_e->GetXaxis()->FindBin(x_pix *
-                                                                                                          1000));
-            if (isHole) average_charge =
-                m_avgChargeMap_h->GetBinContent(m_avgChargeMap_h->GetYaxis()->FindBin(
-                                                  y_pix * 1000), m_avgChargeMap_h->GetXaxis()->FindBin(x_pix * 1000));
+            float average_charge = isHole ? m_avgChargeMap_h.getContent(m_avgChargeMap_h.getBinY(1e3*y_pix), m_avgChargeMap_h.getBinX(1e3*x_pix)) :
+                                            m_avgChargeMap_e.getContent(m_avgChargeMap_e.getBinY(1e3*y_pix), m_avgChargeMap_e.getBinX(1e3*x_pix));
 
             ATH_MSG_DEBUG(" -- driftTime, timeToElectrode = " << driftTime << "  " << timeToElectrode);
 
@@ -379,6 +371,9 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
             for (int i = -1; i <= 1; i++) {
               double xNeighbor = i * pixel_size_x;
               // -- loop in the y-coordinate
+              const std::size_t index = 0;
+              const std::size_t ramo_init_bin_y  = m_ramoPotentialMap[index].getBinY(1000*(x_pix + pixel_size_x * 3 - xNeighbor));
+              const std::size_t ramo_final_bin_y = m_ramoPotentialMap[index].getBinY(1000*(xposFinal + pixel_size_x * 3 - xNeighbor));
               for (int j = -1; j <= 1; j++) {
                 double yNeighbor = j * pixel_size_y;
 
@@ -391,10 +386,8 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
 
                 //Ramo map over 500umx350um pixel area
                 //Ramo init different if charge diffused into neighboring pixel -> change primary pixel!!
-                double ramoInit = getRamoPotential(y_pix + pixel_size_y * 1 / 2 - yNeighbor,
-                                                   x_pix + pixel_size_x * 3 - xNeighbor);
-                double ramoFinal = getRamoPotential(yposFinal + pixel_size_y * 1 / 2 - yNeighbor,
-                                                    xposFinal + pixel_size_x * 3 - xNeighbor);
+                float ramoInit  = m_ramoPotentialMap[index].getContent(m_ramoPotentialMap[index].getBinX(1000*(y_pix + 0.5*pixel_size_y - yNeighbor)), ramo_init_bin_y);
+                float ramoFinal = m_ramoPotentialMap[index].getContent(m_ramoPotentialMap[index].getBinX(1000*(yposFinal + 0.5*pixel_size_y - yNeighbor)), ramo_final_bin_y);
 
                 // Record deposit
                 double eHitRamo = (1 - 2 * isHole) * eHit * (ramoFinal - ramoInit);
@@ -508,8 +501,6 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
             double ed = es_current * eleholePairEnergy * ccprob_neighbor;
 
             // -- pixel coordinates --> module coordinates
-            //double x_mod = x_neighbor - half_pixel_size_x + pixel_size_x*nPixX -module_size_x/2.;
-            //double y_mod = y_neighbor - half_pixel_size_y + pixel_size_y*nPixY -module_size_y/2.;
             double x_mod = x_neighbor + pixel_size_x / 2 + pixel_size_x * nPixX - module_size_x / 2.;
             double y_mod = y_neighbor + pixel_size_y / 2 + pixel_size_y * nPixY - module_size_y / 2.;
             SiLocalPosition chargePos = Module.hitLocalToLocal(y_mod, x_mod);
@@ -531,7 +522,7 @@ StatusCode SensorSim3DTool::induceCharge(const TimedHitPtr<SiHit>& phit,
 }
 
 // read the Charge Collection Prob Map from text file
-StatusCode SensorSim3DTool::readProbMap(std::string fileE) {
+StatusCode SensorSim3DTool::readProbMap(const std::string& fileE) {
   std::string line;
   const std::string fileName = fileE;
   std::string inputFile = PathResolverFindCalibFile(fileName);
@@ -567,20 +558,18 @@ StatusCode SensorSim3DTool::readProbMap(std::string fileE) {
 }
 
 // -- Print out the Charge Collection Probability map (full map)
-StatusCode SensorSim3DTool::printProbMap(std::string readout) {
+StatusCode SensorSim3DTool::printProbMap(const std::string& readout) const {
   if (readout == "FEI4") {
-    for (std::multimap<std::pair<int, int>, double >::iterator it = m_probMapFEI4.begin(); it != m_probMapFEI4.end();
-         ++it) {
+    for (const auto& it : m_probMapFEI4) {
       ATH_MSG_DEBUG(
-        "read full probMap  FEI4 --- bin x " << it->first.first << "   bin y  " << it->first.second << "    prob  " <<
-          it->second);
+        "read full probMap  FEI4 --- bin x " << it.first.first << "   bin y  " << it.first.second << "    prob  " <<
+          it.second);
     }
   } else if (readout == "FEI3") {
-    for (std::multimap<std::pair<int, int>, double >::iterator it = m_probMapFEI3.begin(); it != m_probMapFEI3.end();
-         ++it) {
+    for (const auto& it : m_probMapFEI3) {
       ATH_MSG_DEBUG(
-        "read full probMap  FEI3 --- bin x " << it->first.first << "   bin y  " << it->first.second << "    prob  " <<
-          it->second);
+        "read full probMap  FEI3 --- bin x " << it.first.first << "   bin y  " << it.first.second << "    prob  " <<
+          it.second);
     }
   } else {
     ATH_MSG_ERROR("Error in printout Charge Coll Prob Maps! (readout should contain FEI3 or FEI4 strings) ");
@@ -590,7 +579,7 @@ StatusCode SensorSim3DTool::printProbMap(std::string readout) {
 }
 
 // -- Returns the Charge Collection Probability at a given point (bin_x,bin_y)
-double SensorSim3DTool::getProbMapEntry(std::string readout, int binx, int biny) {
+double SensorSim3DTool::getProbMapEntry(const std::string& readout, int binx, int biny) const {
   std::pair<int, int> doublekey(binx, biny);
   double echarge;
   if (readout == "FEI4") {
@@ -607,15 +596,8 @@ double SensorSim3DTool::getProbMapEntry(std::string readout, int binx, int biny)
 }
 
 double SensorSim3DTool::getElectricField(double x, double y) {
-  std::pair < int, int > Layer; // index for layer/end cap position
-  Layer.first = 0; //Barrel (0) or End Cap (1)   -    Now only for Barrel. If we want to add End Caps, put them at
-                   // Layer.first=1
-  Layer.second = 0; //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2
-
-  int n_binx = m_eFieldMap[Layer]->GetXaxis()->FindBin(x * 1000); //position coordinates in um to return electric field
-                                                                  // in V/cm
-  int n_biny = m_eFieldMap[Layer]->GetYaxis()->FindBin(y * 1000);
-  double electricField = m_eFieldMap[Layer]->GetBinContent(n_binx, n_biny);
+  std::size_t index = 0;
+  double electricField = m_eFieldMap[index].getContent(m_eFieldMap[index].getBinX(1e3*x), m_eFieldMap[index].getBinY(1e3*y));
   return electricField * 1.0E-7; //return efield in MV/mm (for mobility calculation)
 }
 
@@ -656,79 +638,36 @@ double SensorSim3DTool::getDriftTime(bool isHoleBit) {
 }
 
 double SensorSim3DTool::getTimeToElectrode(double x, double y, bool isHoleBit) {
-  std::pair < int, int > Layer; // index for layer/end cap position
-  Layer.first = 0; //Barrel (0) or End Cap (1)   -    Now only for Barrel. If we want to add End Caps, put them at
-                   // Layer.first=1
-  Layer.second = 0; //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2
-
-  // Uses (x,y) position in um to return time to electrode in ns
+  std::size_t index = 0;
   double timeToElectrode = 0;
   if (!isHoleBit) {
-    int n_binx = m_timeMap_e[Layer]->GetXaxis()->FindBin(x * 1000); //convert from mm to um
-    int n_biny = m_timeMap_e[Layer]->GetYaxis()->FindBin(y * 1000);
-    timeToElectrode = m_timeMap_e[Layer]->GetBinContent(n_binx, n_biny);
-  }
-  if (isHoleBit) {
-    int n_binx = m_timeMap_h[Layer]->GetXaxis()->FindBin(x * 1000);
-    int n_biny = m_timeMap_h[Layer]->GetYaxis()->FindBin(y * 1000);
-    timeToElectrode = m_timeMap_h[Layer]->GetBinContent(n_binx, n_biny);
+    timeToElectrode = m_timeMap_e[index].getContent(m_timeMap_e[index].getBinX(1e3*x), m_timeMap_e[index].getBinY(1e3*y));
+  } else {
+    timeToElectrode = m_timeMap_h[index].getContent(m_timeMap_h[index].getBinX(1e3*x), m_timeMap_h[index].getBinY(1e3*y));
   }
   return timeToElectrode; //[ns]
 }
 
 double SensorSim3DTool::getTrappingPositionX(double initX, double initY, double driftTime, bool isHoleBit) {
-  std::pair < int, int > Layer; // index for layer/end cap position
-  Layer.first = 0; //Barrel (0) or End Cap (1)   -    Now only for Barrel. If we want to add End Caps, put them at
-                   // Layer.first=1
-  Layer.second = 0; //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2
-
+  std::size_t index = 0;
   double finalX = initX;
   if (!isHoleBit) {
-    int n_binx = m_xPositionMap_e[Layer]->GetXaxis()->FindBin(initX * 1000);
-    int n_biny = m_xPositionMap_e[Layer]->GetYaxis()->FindBin(initY * 1000);
-    int n_bint = m_xPositionMap_e[Layer]->GetZaxis()->FindBin(driftTime);
-    finalX = m_xPositionMap_e[Layer]->GetBinContent(n_binx, n_biny, n_bint);
+    finalX = m_xPositionMap_e[index].getContent(m_xPositionMap_e[index].getBinX(1e3*initX), m_xPositionMap_e[index].getBinY(1e3*initY), m_xPositionMap_e[index].getBinZ(driftTime));
   } else {
-    int n_binx = m_xPositionMap_h[Layer]->GetXaxis()->FindBin(initX * 1000);
-    int n_biny = m_xPositionMap_h[Layer]->GetYaxis()->FindBin(initY * 1000);
-    int n_bint = m_xPositionMap_h[Layer]->GetZaxis()->FindBin(driftTime);
-    finalX = m_xPositionMap_h[Layer]->GetBinContent(n_binx, n_biny, n_bint);
+    finalX = m_xPositionMap_h[index].getContent(m_xPositionMap_h[index].getBinX(1e3*initX), m_xPositionMap_h[index].getBinY(1e3*initY), m_xPositionMap_h[index].getBinZ(driftTime));
   }
 
   return finalX * 1e-3; //[mm]
 }
 
 double SensorSim3DTool::getTrappingPositionY(double initX, double initY, double driftTime, bool isHoleBit) {
-  std::pair < int, int > Layer; // index for layer/end cap position
-  Layer.first = 0; //Barrel (0) or End Cap (1)   -    Now only for Barrel. If we want to add End Caps, put them at
-                   // Layer.first=1
-  Layer.second = 0; //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2
-
+  std::size_t index = 0;
   double finalY = initY;
   if (!isHoleBit) {
-    int n_binx = m_yPositionMap_e[Layer]->GetXaxis()->FindBin(initX * 1000);
-    int n_biny = m_yPositionMap_e[Layer]->GetYaxis()->FindBin(initY * 1000);
-    int n_bint = m_yPositionMap_e[Layer]->GetZaxis()->FindBin(driftTime);
-    finalY = m_yPositionMap_e[Layer]->GetBinContent(n_binx, n_biny, n_bint);
+    finalY = m_yPositionMap_e[index].getContent(m_yPositionMap_e[index].getBinX(1e3*initX), m_yPositionMap_e[index].getBinY(1e3*initY), m_yPositionMap_e[index].getBinZ(driftTime));
   } else {
-    int n_binx = m_yPositionMap_h[Layer]->GetXaxis()->FindBin(initX * 1000);
-    int n_biny = m_yPositionMap_h[Layer]->GetYaxis()->FindBin(initY * 1000);
-    int n_bint = m_yPositionMap_h[Layer]->GetZaxis()->FindBin(driftTime);
-    finalY = m_yPositionMap_h[Layer]->GetBinContent(n_binx, n_biny, n_bint);
+    finalY = m_yPositionMap_h[index].getContent(m_yPositionMap_h[index].getBinX(1e3*initX), m_yPositionMap_h[index].getBinY(1e3*initY), m_yPositionMap_h[index].getBinZ(driftTime));
   }
 
   return finalY * 1e-3; //[mm]
 }
-
-double SensorSim3DTool::getRamoPotential(double x, double y) {
-  std::pair < int, int > Layer; // index for layer/end cap position
-  Layer.first = 0; //Barrel (0) or End Cap (1)   -    Now only for Barrel. If we want to add End Caps, put them at
-                   // Layer.first=1
-  Layer.second = 0; //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2
-
-
-  int n_binx = m_ramoPotentialMap[Layer]->GetXaxis()->FindBin(x * 1000);
-  int n_biny = m_ramoPotentialMap[Layer]->GetYaxis()->FindBin(y * 1000);
-  double ramoPotential = m_ramoPotentialMap[Layer]->GetBinContent(n_binx, n_biny);
-  return ramoPotential;
-}
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
index 803e99eff0c7b17708923731e033f9b703d408b3..ea37734ebf756fcdfd8c4b0c37f0c7fc8b5a0ba7 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
@@ -17,6 +17,7 @@
 
 #include "GaudiKernel/ToolHandle.h"
 #include "RadDamageUtil.h"
+#include "PixelConditionsData/PixelHistoConverter.h"
 
 
 class SensorSim3DTool: public SensorSimTool {
@@ -38,9 +39,9 @@ public:
 
 
   // 3D sensor simulation using probability density map (used in RUN-2 (no radiation damage)
-  StatusCode readProbMap(std::string);
-  StatusCode printProbMap(std::string);
-  double getProbMapEntry(std::string, int, int);
+  StatusCode readProbMap(const std::string&);
+  StatusCode printProbMap(const std::string&) const;
+  double getProbMapEntry(const std::string&, int, int) const;
 
   double getElectricField(double x, double y);
   double getMobility(double electricField, bool isHoleBit);
@@ -48,7 +49,6 @@ public:
   double getTimeToElectrode(double x, double y, bool isHoleBit);
   double getTrappingPositionX(double initX, double initY, double driftTime, bool isHoleBit);
   double getTrappingPositionY(double initX, double initY, double driftTime, bool isHoleBit);
-  double getRamoPotential(double x, double y);
 private:
   SensorSim3DTool();
 
@@ -57,19 +57,18 @@ private:
   std::multimap<std::pair<int, int>, double> m_probMapFEI3;
 
   // Map for radiation damage simulation
-  std::map<std::pair<int, int>, TH3F*> m_ramoPotentialMap;
-  std::map<std::pair<int, int>, TH2F*> m_eFieldMap;
-  std::map<std::pair<int, int>, TH3F*> m_xPositionMap_e;
-  std::map<std::pair<int, int>, TH3F*> m_xPositionMap_h;
-  std::map<std::pair<int, int>, TH3F*> m_yPositionMap_e;
-  std::map<std::pair<int, int>, TH3F*> m_yPositionMap_h;
-  std::map<std::pair<int, int>, TH2F*> m_timeMap_e;
-  std::map<std::pair<int, int>, TH2F*> m_timeMap_h;
-  TH2F* m_avgChargeMap_e;
-  TH2F* m_avgChargeMap_h;
+  std::vector<PixelHistoConverter> m_ramoPotentialMap;
+  std::vector<PixelHistoConverter> m_eFieldMap;
+  std::vector<PixelHistoConverter> m_xPositionMap_e;
+  std::vector<PixelHistoConverter> m_xPositionMap_h;
+  std::vector<PixelHistoConverter> m_yPositionMap_e;
+  std::vector<PixelHistoConverter> m_yPositionMap_h;
+  std::vector<PixelHistoConverter> m_timeMap_e;
+  std::vector<PixelHistoConverter> m_timeMap_h;
+  PixelHistoConverter m_avgChargeMap_e;
+  PixelHistoConverter m_avgChargeMap_h;
 
   std::vector<double> m_fluence_layers;
-  std::map<std::pair<int, int>, double> m_fluence_layersMaps;
 
   Gaudi::Property<std::string> m_cc_prob_file_fei3
   {
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
index b6b379f47437aede0817ade2d01b1dce1c9933ac..8622bf6ec9f0bec78d10776b965c10a98b6eebd7 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
@@ -94,7 +94,8 @@ StatusCode SensorSimPlanarTool::initialize() {
         ATH_MSG_WARNING("Not implemented yet - exit");
         return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up
       }
-      m_ramoPotentialMap.push_back(ramoPotentialMap_hold);
+      m_ramoPotentialMap.emplace_back();
+      ATH_CHECK(m_ramoPotentialMap.back().setHisto3D(ramoPotentialMap_hold));
       //Now setup the E-field.
       TH1F* eFieldMap_hold;
       eFieldMap_hold = new TH1F();
@@ -135,13 +136,15 @@ StatusCode SensorSimPlanarTool::initialize() {
         ATH_MSG_INFO("Unable to load at least one of the distance/time/Lorentz angle maps.");
         return StatusCode::FAILURE;//Obviously, remove this when gen. code is set up
       }
-      m_distanceMap_e.push_back(distanceMap_e_hold);
-      m_distanceMap_h.push_back(distanceMap_h_hold);
-      m_lorentzMap_e.push_back(lorentzMap_e_hold);
-      m_lorentzMap_h.push_back(lorentzMap_h_hold);
+      m_distanceMap_e.emplace_back();
+      m_distanceMap_h.emplace_back();
+      ATH_CHECK(m_distanceMap_e.back().setHisto2D(distanceMap_e_hold));
+      ATH_CHECK(m_distanceMap_h.back().setHisto2D(distanceMap_h_hold));
+      m_lorentzMap_e.emplace_back();
+      m_lorentzMap_h.emplace_back();
+      ATH_CHECK(m_lorentzMap_e.back().setHisto2D(lorentzMap_e_hold));
+      ATH_CHECK(m_lorentzMap_h.back().setHisto2D(lorentzMap_h_hold));
     }
-    // Define necessary variables to describe hists stored in ramoPotentialMap
-    // to skip calling FindBin during induceCharge for speed reasons
   }
   return StatusCode::SUCCESS;
 }
@@ -248,16 +251,10 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
 
     SiLocalPosition centreOfPixel_i;
 
-    int nnLoop_pixelEtaMax = 0;
-    int nnLoop_pixelEtaMin = 0;
-    int nnLoop_pixelPhiMax = 0;
-    int nnLoop_pixelPhiMin = 0;
-
-    int numBins_driftTime_e = 0;
-    int numBins_driftTime_h = 0;
-    int numBins_weightingPotential_x = 0;
-    int numBins_weightingPotential_y = 0;
-    int numBins_weightingPotential_z = 0;
+    int nnLoop_pixelEtaMax(0);
+    int nnLoop_pixelEtaMin(0);
+    int nnLoop_pixelPhiMax(0);
+    int nnLoop_pixelPhiMin(0);
 
     if (m_doRadDamage && !(Module.isDBM()) && Module.isBarrel()) {
       centreOfPixel_i = p_design.positionFromColumnRow(pixel_i.etaIndex(), pixel_i.phiIndex());
@@ -268,21 +265,6 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
 
       nnLoop_pixelPhiMax = std::min(1, pixel_i.phiIndex());
       nnLoop_pixelPhiMin = std::max(-1, pixel_i.phiIndex() + 1 - phiCells);
-
-      //Setup values to check for overflow when using maps
-      if (m_doInterpolateEfield) {
-        numBins_driftTime_e = m_distanceMap_e[layer]->GetNbinsY(); //Returns nBins = totalBins - underflow - overflow
-        numBins_driftTime_h = m_distanceMap_h[layer]->GetNbinsY();
-        numBins_weightingPotential_x = m_ramoPotentialMap[layer]->GetNbinsX();
-        numBins_weightingPotential_y = m_ramoPotentialMap[layer]->GetNbinsY();
-        numBins_weightingPotential_z = m_ramoPotentialMap[layer]->GetNbinsZ();
-      } else { // use fluence value from conditions data
-        numBins_driftTime_e = moduleData->getDistanceMap_e(layer)->GetNbinsY();
-        numBins_driftTime_h = moduleData->getDistanceMap_h(layer)->GetNbinsY();
-        numBins_weightingPotential_x = moduleData->getRamoPotentialMap(layer)->GetNbinsX();
-        numBins_weightingPotential_y = moduleData->getRamoPotentialMap(layer)->GetNbinsY();
-        numBins_weightingPotential_z = moduleData->getRamoPotentialMap(layer)->GetNbinsZ();
-      }
     }
 
     // Distance between charge and readout side.  p_design->readoutSide() is
@@ -298,33 +280,10 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
       nontrappingProbability = exp(-dist_electrode / collectionDist);
     }
 
-    if (m_doInterpolateEfield) {
-      m_ramo_x_binMap = 1000. *
-                        (m_ramoPotentialMap[layer]->GetNbinsX() /
-                         (m_ramoPotentialMap[layer]->GetXaxis()->GetXmax() -
-                          m_ramoPotentialMap[layer]->GetXaxis()->GetXmin()));
-      m_ramo_y_binMap = 1000. *
-                        (m_ramoPotentialMap[layer]->GetNbinsY() /
-                         (m_ramoPotentialMap[layer]->GetYaxis()->GetXmax() -
-                          m_ramoPotentialMap[layer]->GetYaxis()->GetXmin()));
-      m_ramo_z_binMap = 1000. *
-                        (m_ramoPotentialMap[layer]->GetNbinsZ() /
-                         (m_ramoPotentialMap[layer]->GetZaxis()->GetXmax() -
-                          m_ramoPotentialMap[layer]->GetZaxis()->GetXmin()));
-    } else {
-      m_ramo_x_binMap = 1000. *
-                        (moduleData->getRamoPotentialMap(layer)->GetNbinsX() /
-                         (moduleData->getRamoPotentialMap(layer)->GetXaxis()->GetXmax() -
-                          moduleData->getRamoPotentialMap(layer)->GetXaxis()->GetXmin()));
-      m_ramo_y_binMap = 1000. *
-                        (moduleData->getRamoPotentialMap(layer)->GetNbinsY() /
-                         (moduleData->getRamoPotentialMap(layer)->GetYaxis()->GetXmax() -
-                          moduleData->getRamoPotentialMap(layer)->GetYaxis()->GetXmin()));
-      m_ramo_z_binMap = 1000. *
-                        (moduleData->getRamoPotentialMap(layer)->GetNbinsZ() /
-                         (moduleData->getRamoPotentialMap(layer)->GetZaxis()->GetXmax() -
-                          moduleData->getRamoPotentialMap(layer)->GetZaxis()->GetXmin()));
-    }
+    const std::size_t distance_f_e_bin_x = m_doInterpolateEfield ? m_distanceMap_e[layer].getBinX(dist_electrode) : moduleData->getDistanceMap_e(layer).getBinX(dist_electrode);
+    const std::size_t distance_f_h_bin_x = m_doInterpolateEfield ? m_distanceMap_h[layer].getBinX(dist_electrode) : moduleData->getDistanceMap_h(layer).getBinX(dist_electrode);
+    const std::size_t tanLorentz_e_bin_x = m_doInterpolateEfield ? m_lorentzMap_e[layer].getBinX(dist_electrode) : moduleData->getLorentzMap_e(layer).getBinX(dist_electrode);
+    const std::size_t tanLorentz_h_bin_x = m_doInterpolateEfield ? m_lorentzMap_h[layer].getBinX(dist_electrode) : moduleData->getLorentzMap_h(layer).getBinX(dist_electrode);
 
     for (int j = 0; j < ncharges; j++) {
       if (m_doRadDamage && !(Module.isDBM()) && Module.isBarrel()) {
@@ -344,39 +303,15 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
         // electrode.
         //However, if choose (b), will need to reduce granularity of map.
         if (m_doInterpolateEfield) {
-          int nbin_z_e_xbin = m_distanceMap_e[layer]->GetXaxis()->FindBin(dist_electrode);
-          int nbin_z_e_ybin = m_distanceMap_e[layer]->GetYaxis()->FindBin(drifttime_e);
-          if (nbin_z_e_ybin > numBins_driftTime_e) {
-            nbin_z_e_ybin = numBins_driftTime_e;
-          }
-          depth_f_e = m_distanceMap_e[layer]->GetBinContent(nbin_z_e_xbin, nbin_z_e_ybin);
-          int nbin_z_h_xbin = m_distanceMap_h[layer]->GetXaxis()->FindBin(dist_electrode);
-          int nbin_z_h_ybin = m_distanceMap_h[layer]->GetYaxis()->FindBin(drifttime_h);
-          if (nbin_z_h_ybin > numBins_driftTime_h) {
-            nbin_z_h_ybin = numBins_driftTime_h;
-          }
-          depth_f_h = m_distanceMap_h[layer]->GetBinContent(nbin_z_h_xbin, nbin_z_h_ybin);
-          int nbin_Lorentz_e = m_lorentzMap_e[layer]->FindBin(dist_electrode, depth_f_e);
-          tanLorentz_e = m_lorentzMap_e[layer]->GetBinContent(nbin_Lorentz_e);
-          int nbin_Lorentz_h = m_lorentzMap_h[layer]->FindBin(dist_electrode, depth_f_h);
-          tanLorentz_h = m_lorentzMap_h[layer]->GetBinContent(nbin_Lorentz_h);
+          depth_f_e = m_distanceMap_e[layer].getContent(distance_f_e_bin_x, m_distanceMap_e[layer].getBinY(drifttime_e));
+          depth_f_h = m_distanceMap_h[layer].getContent(distance_f_h_bin_x, m_distanceMap_h[layer].getBinY(drifttime_h));
+          tanLorentz_e = m_lorentzMap_e[layer].getContent(tanLorentz_e_bin_x, m_lorentzMap_e[layer].getBinY(depth_f_e));
+          tanLorentz_h = m_lorentzMap_h[layer].getContent(tanLorentz_h_bin_x, m_lorentzMap_h[layer].getBinY(depth_f_h));
         } else { // use fluence value from conditions data
-          int nbin_z_e_xbin = moduleData->getDistanceMap_e(layer)->GetXaxis()->FindBin(dist_electrode);
-          int nbin_z_e_ybin = moduleData->getDistanceMap_e(layer)->GetYaxis()->FindBin(drifttime_e);
-          if (nbin_z_e_ybin > numBins_driftTime_e) {
-            nbin_z_e_ybin = numBins_driftTime_e;
-          }
-          depth_f_e = moduleData->getDistanceMap_e(layer)->GetBinContent(nbin_z_e_xbin, nbin_z_e_ybin);
-          int nbin_z_h_xbin = moduleData->getDistanceMap_h(layer)->GetXaxis()->FindBin(dist_electrode);
-          int nbin_z_h_ybin = moduleData->getDistanceMap_h(layer)->GetYaxis()->FindBin(drifttime_h);
-          if (nbin_z_h_ybin > numBins_driftTime_h) {
-            nbin_z_h_ybin = numBins_driftTime_h;
-          }
-          depth_f_h = moduleData->getDistanceMap_h(layer)->GetBinContent(nbin_z_h_xbin, nbin_z_h_ybin);
-          int nbin_Lorentz_e = moduleData->getLorentzMap_e(layer)->FindBin(dist_electrode, depth_f_e);
-          tanLorentz_e = moduleData->getLorentzMap_e(layer)->GetBinContent(nbin_Lorentz_e);
-          int nbin_Lorentz_h = moduleData->getLorentzMap_h(layer)->FindBin(dist_electrode, depth_f_h);
-          tanLorentz_h = moduleData->getLorentzMap_h(layer)->GetBinContent(nbin_Lorentz_h);
+          depth_f_e = moduleData->getDistanceMap_e(layer).getContent(distance_f_e_bin_x, moduleData->getDistanceMap_e(layer).getBinY(drifttime_e));
+          depth_f_h = moduleData->getDistanceMap_h(layer).getContent(distance_f_h_bin_x, moduleData->getDistanceMap_h(layer).getBinY(drifttime_h));
+          tanLorentz_e = moduleData->getLorentzMap_e(layer).getContent(tanLorentz_e_bin_x, moduleData->getLorentzMap_e(layer).getBinY(depth_f_e));
+          tanLorentz_h = moduleData->getLorentzMap_h(layer).getContent(tanLorentz_h_bin_x, moduleData->getLorentzMap_h(layer).getBinY(depth_f_h));
         }
         double dz_e = fabs(dist_electrode - depth_f_e);
         double dz_h = fabs(depth_f_h - dist_electrode);
@@ -397,7 +332,7 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
         double phi_f_h = phi_i + dz_h * tanLorentz_h + rdif_h * phiRand;
         etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
         double eta_f_h = eta_i + rdif_h * etaRand;
-
+        
         // amount of energy to be converted into charges at current step
         double energy_per_step = 1.0 * iHitRecord.second / 1.E+6 / ncharges;
 
@@ -407,21 +342,11 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
           ATH_CHECK(applyIBLSlimEdges(energy_per_step, eta_f_h));
         }
 
-        int nbin_ramo_f_e_z;
-        int nbin_ramo_f_h_z;
-        // distinction necessary because of min(z) = -0.5
-        float minz = 1;
-        if (layer != 0) minz = 1.5;
-        nbin_ramo_f_e_z = int( minz + depth_f_e * m_ramo_z_binMap );
-        nbin_ramo_f_h_z = int( minz + depth_f_h * m_ramo_z_binMap );
-
-        // Check for overflow in ramo hists in z-direction
-        if (nbin_ramo_f_h_z > numBins_weightingPotential_z) {
-          nbin_ramo_f_h_z = numBins_weightingPotential_z + 1;
-        }
-        if (nbin_ramo_f_e_z > numBins_weightingPotential_z) {
-          nbin_ramo_f_e_z = numBins_weightingPotential_z + 1;
-        }
+        const std::size_t ramo_f_e_bin_z = m_doInterpolateEfield ? m_ramoPotentialMap[layer].getBinZ(1e3*depth_f_e) : moduleData->getRamoPotentialMap(layer).getBinZ(1e3*depth_f_e);
+        const std::size_t ramo_f_h_bin_z = m_doInterpolateEfield ? m_ramoPotentialMap[layer].getBinZ(1e3*depth_f_h) : moduleData->getRamoPotentialMap(layer).getBinZ(1e3*depth_f_h);
+
+        const bool isFirstZ_e = m_doInterpolateEfield ? m_ramoPotentialMap[layer].isFirstZ(1e3*depth_f_e) : moduleData->getRamoPotentialMap(layer).isFirstZ(1e3*depth_f_e);
+        const bool isOverflowZ_h = m_doInterpolateEfield ? m_ramoPotentialMap[layer].isOverflowZ(1e3*depth_f_h) : moduleData->getRamoPotentialMap(layer).isOverflowZ(1e3*depth_f_h);
 
         //Loop over nearest neighbours in x and y
         //We assume that the lateral diffusion is minimal
@@ -472,74 +397,31 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
             double ramo_f_e = 0.0;
             double ramo_f_h = 0.0;
 
-            int nbin_ramo_f_e_x = int( 1 + std::abs(dPhi_f_e) * m_ramo_x_binMap );
-            int nbin_ramo_f_e_y = int( 1 + std::abs(dEta_f_e) * m_ramo_y_binMap );
-
-            // Check for overflow in ramo hists in x- and y-direction
-            if (nbin_ramo_f_e_x > numBins_weightingPotential_x) {
-              nbin_ramo_f_e_x = numBins_weightingPotential_x + 1;
-            }
-            if (nbin_ramo_f_e_y > numBins_weightingPotential_y) {
-              nbin_ramo_f_e_y = numBins_weightingPotential_y + 1;
-            }
-
-            if (nbin_ramo_f_e_x <= numBins_weightingPotential_x && nbin_ramo_f_e_y <= numBins_weightingPotential_y &&
-                nbin_ramo_f_e_z <= numBins_weightingPotential_z) {
+            if (!isFirstZ_e) {
               if (m_doInterpolateEfield) {
-                ramo_f_e = m_ramoPotentialMap[layer]->GetBinContent(nbin_ramo_f_e_x, nbin_ramo_f_e_y, nbin_ramo_f_e_z);
+                ramo_f_e = m_ramoPotentialMap[layer].getContent(m_ramoPotentialMap[layer].getBinX(1e3*std::abs(dPhi_f_e)), m_ramoPotentialMap[layer].getBinY(1e3*std::abs(dEta_f_e)), ramo_f_e_bin_z);
               } else {
-                ramo_f_e = moduleData->getRamoPotentialMap(layer)->GetBinContent(nbin_ramo_f_e_x, nbin_ramo_f_e_y,
-                                                                                 nbin_ramo_f_e_z);
+                ramo_f_e = moduleData->getRamoPotentialMap(layer).getContent(moduleData->getRamoPotentialMap(layer).getBinX(1e3*std::abs(dPhi_f_e)), moduleData->getRamoPotentialMap(layer).getBinY(1e3*std::abs(dEta_f_e)), ramo_f_e_bin_z);
               }
             }
 
-            int nbin_ramo_f_h_x = int( 1 + std::abs(dPhi_f_h) * m_ramo_x_binMap );
-            int nbin_ramo_f_h_y = int( 1 + std::abs(dEta_f_h) * m_ramo_y_binMap );
-
-            // Check for overflow in ramo hists in x- and y-direction
-            if (nbin_ramo_f_h_x > numBins_weightingPotential_x) {
-              nbin_ramo_f_h_x = numBins_weightingPotential_x + 1;
-            }
-            if (nbin_ramo_f_h_y > numBins_weightingPotential_y) {
-              nbin_ramo_f_h_y = numBins_weightingPotential_y + 1;
-            }
-
-            if (nbin_ramo_f_h_x <= numBins_weightingPotential_x && nbin_ramo_f_h_y <= numBins_weightingPotential_y &&
-                nbin_ramo_f_h_z <= numBins_weightingPotential_z) {
+            if (!isOverflowZ_h) {
               if (m_doInterpolateEfield) {
-                ramo_f_h = m_ramoPotentialMap[layer]->GetBinContent(nbin_ramo_f_h_x, nbin_ramo_f_h_y, nbin_ramo_f_h_z);
+                ramo_f_h = m_ramoPotentialMap[layer].getContent(m_ramoPotentialMap[layer].getBinX(1e3*std::abs(dPhi_f_h)), m_ramoPotentialMap[layer].getBinY(1e3*std::abs(dEta_f_h)), ramo_f_h_bin_z);
               } else {
-                ramo_f_h = moduleData->getRamoPotentialMap(layer)->GetBinContent(nbin_ramo_f_h_x, nbin_ramo_f_h_y,
-                                                                                 nbin_ramo_f_h_z);
+                ramo_f_h = moduleData->getRamoPotentialMap(layer).getContent(moduleData->getRamoPotentialMap(layer).getBinX(1e3*std::abs(dPhi_f_h)), moduleData->getRamoPotentialMap(layer).getBinY(1e3*std::abs(dEta_f_h)), ramo_f_h_bin_z);
               }
             }
             //Account for the imperfect binning that would cause charge to be double-counted
+            if (isOverflowZ_h) {
+              ramo_f_h = 0;
+            }
 
-            if (m_doInterpolateEfield) {
-              if (m_ramoPotentialMap[layer]->GetZaxis()->FindBin(depth_f_h * 1000) ==
-                  m_ramoPotentialMap[layer]->GetNbinsZ() + 1) {
-                ramo_f_h = 0.0;
-              } //this means the hole has reached the back end
-
-              if (m_ramoPotentialMap[layer]->GetZaxis()->FindBin(depth_f_e * 1000) == 1) {
-                if (fabs(dEta_f_e) >= Module.etaPitch() / 2.0 || fabs(dPhi_f_e) >= Module.phiPitch() / 2.0) {
-                  ramo_f_e = 0.0;
-                } else if (fabs(dEta_f_e) < Module.etaPitch() / 2.0 && fabs(dPhi_f_e) < Module.phiPitch() / 2.0) {
-                  ramo_f_e = 1.0;
-                }
-              }
-            } else {
-              if (moduleData->getRamoPotentialMap(layer)->GetZaxis()->FindBin(depth_f_h * 1000)
-                  == moduleData->getRamoPotentialMap(layer)->GetNbinsZ() + 1) {
-                ramo_f_h = 0;
-              } //this means the hole has reached the back end
-
-              if (moduleData->getRamoPotentialMap(layer)->GetZaxis()->FindBin(depth_f_e * 1000) == 1) {
-                if (fabs(dEta_f_e) >= Module.etaPitch() / 2.0 || fabs(dPhi_f_e) >= Module.phiPitch() / 2.0) {
-                  ramo_f_e = 0.0;
-                } else if (fabs(dEta_f_e) < Module.etaPitch() / 2.0 && fabs(dPhi_f_e) < Module.phiPitch() / 2.0) {
-                  ramo_f_e = 1.0;
-                }
+            if (isFirstZ_e) {
+              if (fabs(dEta_f_e) >= Module.etaPitch() / 2.0 || fabs(dPhi_f_e) >= Module.phiPitch() / 2.0) {
+                ramo_f_e = 0.0;
+              } else if (fabs(dEta_f_e) < Module.etaPitch() / 2.0 && fabs(dPhi_f_e) < Module.phiPitch() / 2.0) {
+                ramo_f_e = 1.0;
               }
             }
 
@@ -561,6 +443,7 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
             if (diode.isValid()) {
               chargedDiodes.add(diode, charge);
             } //IF
+
           } //For q
         } //for p
       } else { //If no radDamage, run original
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h
index 085e28f748554d1e9625c5cdb10f833300c8188d..8dd2757762e66cd44bb101716e65329d98c76e0f 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h
@@ -13,9 +13,10 @@
 #define PIXELDIGITIZATION_SensorSimPlanarTool_H
 
 #include "AthenaBaseComps/AthAlgTool.h"
-#include "SensorSimTool.h"
 #include "InDetCondTools/ISiLorentzAngleTool.h"
+#include "SensorSimTool.h"
 #include "RadDamageUtil.h"
+#include "PixelConditionsData/PixelHistoConverter.h"
 
 
 class SensorSimPlanarTool: public SensorSimTool {
@@ -41,16 +42,11 @@ private:
   SensorSimPlanarTool();
 
   // Map for radiation damage simulation
-  std::vector<TH3F*> m_ramoPotentialMap;
-  std::vector<TH2F*> m_distanceMap_e;
-  std::vector<TH2F*> m_distanceMap_h;
-  std::vector<TH2F*> m_lorentzMap_e;
-  std::vector<TH2F*> m_lorentzMap_h;
-
-  // maps to directly get factor to calculate bin instead of calling FindBin
-  double m_ramo_x_binMap;
-  double m_ramo_y_binMap;
-  double m_ramo_z_binMap;
+  std::vector<PixelHistoConverter> m_ramoPotentialMap;
+  std::vector<PixelHistoConverter> m_distanceMap_e;
+  std::vector<PixelHistoConverter> m_distanceMap_h;
+  std::vector<PixelHistoConverter> m_lorentzMap_e;
+  std::vector<PixelHistoConverter> m_lorentzMap_h;
 
   Gaudi::Property<int> m_numberOfSteps
   {