diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/python/PixelConditionsConfig.py b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/python/PixelConditionsConfig.py
index 3ab63159c016aa58e473287a7e1cb4ac07a31311..f26f7beab959df653de0e9843790f9620c82cdb7 100644
--- a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/python/PixelConditionsConfig.py
+++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/python/PixelConditionsConfig.py
@@ -53,16 +53,23 @@ def PixelConfigCondAlgCfg(flags, name="PixelConfigCondAlg", **kwargs):
         #====================================================================================
         # Run-dependent SIMULATION(digitization) parameters:
         #====================================================================================
-        # RUN2 2015/2016
+        # RUN2 2015/2016 (mc16a)
+        # The pixel conditions are matched with 2016 data (mc16a) at L=17.3fb-1 (run#303638).
         CondArgs.update(
-            BarrelToTThreshold2016       = [   -1,    5,    5,    5],
-            FEI3BarrelLatency2016        = [    0,  151,  256,  256],
-            FEI3BarrelHitDuplication2016 = [False,False,False,False],
-            FEI3BarrelSmallHitToT2016    = [   -1,   -1,   -1,   -1],
-            FEI3BarrelTimingSimTune2016  = [   -1, 2015, 2015, 2015],
-            BarrelCrossTalk2016          = [ 0.30, 0.06, 0.06, 0.06],
-            BarrelNoiseOccupancy2016     = [ 5e-8, 5e-8, 5e-8, 5e-8],
-            BarrelDisableProbability2016 = [ 9e-3, 9e-3, 9e-3, 9e-3],
+            BarrelToTThreshold2016       = [     -1,      5,      5,      5],
+            FEI3BarrelLatency2016        = [      0,    151,    256,    256],
+            FEI3BarrelHitDuplication2016 = [  False,  False,  False,  False],
+            FEI3BarrelSmallHitToT2016    = [     -1,     -1,     -1,     -1],
+            FEI3BarrelTimingSimTune2016  = [     -1,   2015,   2015,   2015],
+            BarrelCrossTalk2016          = [   0.30,   0.06,   0.06,   0.06],
+            BarrelNoiseOccupancy2016     = [   5e-8,   5e-8,   5e-8,   5e-8],
+            BarrelDisableProbability2016 = [   9e-3,   9e-3,   9e-3,   9e-3],
+            DefaultBarrelBiasVoltage2016 = [   80.0,  350.0,  200.0,  150.0],
+            BarrelFluence2016            = [0.80e14,1.61e14,0.71e14,0.48e14],
+            BarrelFluenceMap2016 = ["PixelDigitization/maps_IBL_PL_80V_fl0_8e14.root",
+                                    "PixelDigitization/maps_PIX_350V_fl1_61e14.root",
+                                    "PixelDigitization/maps_PIX_200V_fl0_71e14.root",
+                                    "PixelDigitization/maps_PIX_150V_fl0_48e14.root"],
             EndcapToTThreshold2016       = [    5,    5,    5],
             FEI3EndcapLatency2016        = [  256,  256,  256],
             FEI3EndcapHitDuplication2016 = [False,False,False],
@@ -85,16 +92,23 @@ def PixelConfigCondAlgCfg(flags, name="PixelConfigCondAlg", **kwargs):
             #                  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1735, 0.3380, 0.4733, 0.5829, 0.6730, 0.7516, 0.8234, 0.8916, 0.9595, 1.0]]
         )
         #====================================================================================
-        # RUN2 2017
+        # RUN2 2017 (mc16d)
+        # The pixel conditions are matched with 2017 data (mc16d) at L=69.0fb-1 (run#336506).
         CondArgs.update(
-            BarrelToTThreshold2017       = [   -1,    5,    5,    5],
-            FEI3BarrelLatency2017        = [    0,  151,  256,  256],
-            FEI3BarrelHitDuplication2017 = [False,False,False,False],
-            FEI3BarrelSmallHitToT2017    = [   -1,   -1,   -1,   -1],
-            FEI3BarrelTimingSimTune2017  = [   -1, 2018, 2018, 2018],
-            BarrelCrossTalk2017          = [ 0.30, 0.06, 0.06, 0.06],
-            BarrelNoiseOccupancy2017     = [ 5e-8, 5e-8, 5e-8, 5e-8],
-            BarrelDisableProbability2017 = [ 9e-3, 9e-3, 9e-3, 9e-3],
+            BarrelToTThreshold2017       = [     -1,      5,      5,      5],
+            FEI3BarrelLatency2017        = [      0,    151,    256,    256],
+            FEI3BarrelHitDuplication2017 = [  False,  False,  False,  False],
+            FEI3BarrelSmallHitToT2017    = [     -1,     -1,     -1,     -1],
+            FEI3BarrelTimingSimTune2017  = [     -1,   2018,   2018,   2018],
+            BarrelCrossTalk2017          = [   0.30,   0.06,   0.06,   0.06],
+            BarrelNoiseOccupancy2017     = [   5e-8,   5e-8,   5e-8,   5e-8],
+            BarrelDisableProbability2017 = [   9e-3,   9e-3,   9e-3,   9e-3],
+            DefaultBarrelBiasVoltage2017 = [  350.0,  350.0,  200.0,  150.0],
+            BarrelFluence2017            = [3.18e14,3.42e14,1.50e14,1.01e14],
+            BarrelFluenceMap2017 = ["PixelDigitization/maps_IBL_PL_350V_fl3_18e14.root",
+                                    "PixelDigitization/maps_PIX_350V_fl3_42e14.root",
+                                    "PixelDigitization/maps_PIX_200V_fl1_5e14.root",
+                                    "PixelDigitization/maps_PIX_150V_fl1_01e14.root"],
             EndcapToTThreshold2017       = [    5,    5,    5],
             FEI3EndcapLatency2017        = [  256,  256,  256],
             FEI3EndcapHitDuplication2017 = [False,False,False],
@@ -112,16 +126,23 @@ def PixelConfigCondAlgCfg(flags, name="PixelConfigCondAlg", **kwargs):
             PixelNoiseShape2017  = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2418, 0.4397, 0.5858, 0.6949, 0.7737, 0.8414, 0.8959, 0.9414, 0.9828, 1.0],
         )
         #====================================================================================
-        # RUN2 2018
+        # RUN2 2018 (mc16e)
+        # The pixel conditions are matched with 2018 data (mc16e) at L=119.4fb-1 (run#357193).
         CondArgs.update(
-            BarrelToTThreshold2018       = [   -1,    3,    5,    5],
-            FEI3BarrelLatency2018        = [    0,  151,  256,  256],
-            FEI3BarrelHitDuplication2018 = [False,False,False,False],
-            FEI3BarrelSmallHitToT2018    = [   -1,   -1,   -1,   -1],
-            FEI3BarrelTimingSimTune2018  = [   -1, 2018, 2018, 2018],
-            BarrelCrossTalk2018          = [ 0.30, 0.06, 0.06, 0.06],
-            BarrelNoiseOccupancy2018     = [ 5e-8, 5e-8, 5e-8, 5e-8],
-            BarrelDisableProbability2018 = [ 9e-3, 9e-3, 9e-3, 9e-3],
+            BarrelToTThreshold2018       = [     -1,      3,      5,      5],
+            FEI3BarrelLatency2018        = [      0,    151,    256,    256],
+            FEI3BarrelHitDuplication2018 = [  False,  False,  False,  False],
+            FEI3BarrelSmallHitToT2018    = [     -1,     -1,     -1,     -1],
+            FEI3BarrelTimingSimTune2018  = [     -1,   2018,   2018,   2018],
+            BarrelCrossTalk2018          = [   0.30,   0.06,   0.06,   0.06],
+            BarrelNoiseOccupancy2018     = [   5e-8,   5e-8,   5e-8,   5e-8],
+            BarrelDisableProbability2018 = [   9e-3,   9e-3,   9e-3,   9e-3],
+            DefaultBarrelBiasVoltage2018 = [  400.0,  400.0,  250.0,  250.0],
+            BarrelFluence2018            = [5.50e14,5.19e14,2.28e14,1.53e14],
+            BarrelFluenceMap2018 = ["PixelDigitization/maps_IBL_PL_400V_fl5_5e14.root",
+                                    "PixelDigitization/maps_PIX_400V_fl5_19e14.root",
+                                    "PixelDigitization/maps_PIX_250V_fl2_28e14.root",
+                                    "PixelDigitization/maps_PIX_250V_fl1_53e14.root"],
             EndcapToTThreshold2018       = [    5,    5,    5],
             FEI3EndcapLatency2018        = [  256,  256,  256],
             FEI3EndcapHitDuplication2018 = [False,False,False],
@@ -141,14 +162,19 @@ def PixelConfigCondAlgCfg(flags, name="PixelConfigCondAlg", **kwargs):
         #====================================================================================
         # RUN1
         CondArgs.update(
-            BarrelToTThresholdRUN1       = [    3,    3,    3],
-            FEI3BarrelLatencyRUN1        = [  256,  256,  256],
-            FEI3BarrelHitDuplicationRUN1 = [ True, True, True],
-            FEI3BarrelSmallHitToTRUN1    = [    7,    7,    7],
-            FEI3BarrelTimingSimTuneRUN1  = [ 2009, 2009, 2009],
-            BarrelCrossTalkRUN1          = [ 0.06, 0.06, 0.06],
-            BarrelNoiseOccupancyRUN1     = [ 5e-8, 5e-8, 5e-8],
-            BarrelDisableProbabilityRUN1 = [ 9e-3, 9e-3, 9e-3],
+            BarrelToTThresholdRUN1       = [      3,      3,      3],
+            FEI3BarrelLatencyRUN1        = [    256,    256,    256],
+            FEI3BarrelHitDuplicationRUN1 = [   True,   True,   True],
+            FEI3BarrelSmallHitToTRUN1    = [      7,      7,      7],
+            FEI3BarrelTimingSimTuneRUN1  = [   2009,   2009,   2009],
+            BarrelCrossTalkRUN1          = [   0.06,   0.06,   0.06],
+            BarrelNoiseOccupancyRUN1     = [   5e-8,   5e-8,   5e-8],
+            BarrelDisableProbabilityRUN1 = [   9e-3,   9e-3,   9e-3],
+            DefaultBarrelBiasVoltageRUN1 = [  150.0,  150.0,  150.0],
+            BarrelFluenceRUN1            = [1.01e14,0.44e14,0.30e14],
+            BarrelFluenceMapRUN1 = ["PixelDigitization/maps_PIX_250V_fl1_01e14.root",
+                                    "PixelDigitization/maps_PIX_150V_fl0_44e14.root",
+                                    "PixelDigitization/maps_PIX_150V_fl0_3e14.root"],
             EndcapToTThresholdRUN1       = [    3,    3,    3],
             FEI3EndcapLatencyRUN1        = [  256,  256,  256],
             FEI3EndcapHitDuplicationRUN1 = [ True, True, True],
@@ -163,10 +189,17 @@ def PixelConfigCondAlgCfg(flags, name="PixelConfigCondAlg", **kwargs):
         #====================================================================================
         # ITK
         CondArgs.update(
-            BarrelToTThresholdITK       = [    3,    3,    3,    3,    3],
-            BarrelCrossTalkITK          = [ 0.06, 0.06, 0.06, 0.06, 0.06],
-            BarrelNoiseOccupancyITK     = [ 5e-8, 5e-8, 5e-8, 5e-8, 5e-8],
-            BarrelDisableProbabilityITK = [ 9e-3, 9e-3, 9e-3, 9e-3, 9e-3],
+            BarrelToTThresholdITK       = [     3,     3,     3,     3,     3],
+            BarrelCrossTalkITK          = [  0.06,  0.06,  0.06,  0.06,  0.06],
+            BarrelNoiseOccupancyITK     = [  5e-8,  5e-8,  5e-8,  5e-8,  5e-8],
+            BarrelDisableProbabilityITK = [  9e-3,  9e-3,  9e-3,  9e-3,  9e-3],
+            DefaultBarrelBiasVoltageITK = [ 150.0, 150.0, 150.0, 150.0, 150.0],
+            BarrelFluenceITK            = [0.0e14,0.0e14,0.0e14,0.0e14,0.0e14],
+            BarrelFluenceMapITK = ["PixelDigitization/maps_IBL_PL_80V_fl0e14.root",
+                                   "PixelDigitization/maps_IBL_PL_80V_fl0e14.root",
+                                   "PixelDigitization/maps_IBL_PL_80V_fl0e14.root",
+                                   "PixelDigitization/maps_IBL_PL_80V_fl0e14.root",
+                                   "PixelDigitization/maps_IBL_PL_80V_fl0e14.root"],
             EndcapToTThresholdITK       = [    3,    3,    3,    3,    3,    3,    3,    3,    3,    3,    3,    3,    3,    3],
             EndcapCrossTalkITK          = [ 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06],
             EndcapNoiseOccupancyITK     = [ 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8],
@@ -181,8 +214,6 @@ def PixelConfigCondAlgCfg(flags, name="PixelConfigCondAlg", **kwargs):
         DefaultCalibrationParameterA=70.2,
         DefaultCalibrationParameterE=-3561.25,
         DefaultCalibrationParameterC=26000.0
-#        IBLChargeScale=1.0,
-#        IBLSpecificCorrection=False
     )
     # DCS parameters
     CondArgs.update(
diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
index 6a8c82ca2f37bfc59cb37f1318f75fad09b3efb3..db1962721e3ff97aad9fbb85b89adfb705091b3e 100644
--- a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
+++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.cxx
@@ -8,6 +8,10 @@
 #include <memory>
 #include <sstream>
 
+#include "PathResolver/PathResolver.h"
+#include "TFile.h"
+
+
 PixelConfigCondAlg::PixelConfigCondAlg(const std::string& name, ISvcLocator* pSvcLocator):
   ::AthReentrantAlgorithm(name, pSvcLocator)
 {
@@ -167,6 +171,9 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
   writeCdo -> setCablingMapToFile(m_cablingMapToFile);
   writeCdo -> setCablingMapFileName(m_cablingMapFileName);
 
+  // mapping files for radiation damage simulation
+  std::vector<std::string> mapsPath_list;
+
   int currentRunNumber = ctx.eventID().run_number();
   if (currentRunNumber<222222) {
     writeCdo -> setBarrelToTThreshold(m_BarrelToTThresholdRUN1);
@@ -177,6 +184,7 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     writeCdo -> setBarrelCrossTalk(m_BarrelCrossTalkRUN1);
     writeCdo -> setBarrelNoiseOccupancy(m_BarrelNoiseOccupancyRUN1);
     writeCdo -> setBarrelDisableProbability(m_BarrelDisableProbabilityRUN1);
+    writeCdo -> setDefaultBarrelBiasVoltage(m_BarrelBiasVoltageRUN1);
 
     writeCdo -> setEndcapToTThreshold(m_EndcapToTThresholdRUN1);
     writeCdo -> setFEI3EndcapLatency(m_FEI3EndcapLatencyRUN1);
@@ -186,6 +194,7 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     writeCdo -> setEndcapCrossTalk(m_EndcapCrossTalkRUN1);
     writeCdo -> setEndcapNoiseOccupancy(m_EndcapNoiseOccupancyRUN1);
     writeCdo -> setEndcapDisableProbability(m_EndcapDisableProbabilityRUN1);
+    writeCdo -> setDefaultEndcapBiasVoltage(m_EndcapBiasVoltageRUN1);
 
     // This is ad-hoc solution.
     for (size_t i=0; i<m_BLayerNoiseShapeRUN1.size(); i++) { writeCdo->setBarrelNoiseShape(0,m_BLayerNoiseShapeRUN1[i]); }
@@ -196,8 +205,14 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     for (size_t i=0; i<m_PixelNoiseShapeRUN1.size(); i++)  {
       for (size_t layer:{0,1,2}) { writeCdo->setEndcapNoiseShape(layer,m_PixelNoiseShapeRUN1[i]); }
     }
+
+    // Radiation damage simulation
+    writeCdo -> setFluenceLayer(m_BarrelFluenceRUN1);
+    for (size_t i=0; i<m_BarrelFluenceMapRUN1.size(); i++) {
+      mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMapRUN1[i]));
+    }
   }
-  else if (currentRunNumber<240000) {  // RUN2
+  else if (currentRunNumber<240000) {  // RUN2 (mc15)
     writeCdo -> setBarrelToTThreshold(m_BarrelToTThreshold2016);
     writeCdo -> setFEI3BarrelLatency(m_FEI3BarrelLatency2016);
     writeCdo -> setFEI3BarrelHitDuplication(m_FEI3BarrelHitDuplication2016);
@@ -206,6 +221,7 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     writeCdo -> setBarrelCrossTalk(m_BarrelCrossTalk2016);
     writeCdo -> setBarrelNoiseOccupancy(m_BarrelNoiseOccupancy2016);
     writeCdo -> setBarrelDisableProbability(m_BarrelDisableProbability2016);
+    writeCdo -> setDefaultBarrelBiasVoltage(m_BarrelBiasVoltage2016);
 
     writeCdo -> setEndcapToTThreshold(m_EndcapToTThreshold2016);
     writeCdo -> setFEI3EndcapLatency(m_FEI3EndcapLatency2016);
@@ -215,11 +231,13 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     writeCdo -> setEndcapCrossTalk(m_EndcapCrossTalk2016);
     writeCdo -> setEndcapNoiseOccupancy(m_EndcapNoiseOccupancy2016);
     writeCdo -> setEndcapDisableProbability(m_EndcapDisableProbability2016);
+    writeCdo -> setDefaultEndcapBiasVoltage(m_EndcapBiasVoltage2016);
 
     writeCdo -> setDBMToTThreshold(m_DBMToTThreshold2016);
     writeCdo -> setDBMCrossTalk(m_DBMCrossTalk2016);
     writeCdo -> setDBMNoiseOccupancy(m_DBMNoiseOccupancy2016);
     writeCdo -> setDBMDisableProbability(m_DBMDisableProbability2016);
+    writeCdo -> setDefaultDBMBiasVoltage(m_DBMBiasVoltage2016);
 
     /* 
        So far, Gaudi::Property does not support 2D vector.
@@ -244,17 +262,25 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     for (size_t i=0; i<m_IBLNoiseShape2016.size(); i++)    {
       for (size_t layer:{0,1,2}) { writeCdo->setDBMNoiseShape(layer,m_IBLNoiseShape2016[i]); }
     }
+
+    // Radiation damage simulation
+    writeCdo -> setFluenceLayer(m_BarrelFluence2016);
+    for (size_t i=0; i<m_BarrelFluenceMap2016.size(); i++) {
+      mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMap2016[i]));
+    }
   }
   else if (currentRunNumber<250000) {  // RUN4
     writeCdo -> setBarrelToTThreshold(m_BarrelToTThresholdITK);
     writeCdo -> setBarrelCrossTalk(m_BarrelCrossTalkITK);
     writeCdo -> setBarrelNoiseOccupancy(m_BarrelNoiseOccupancyITK);
     writeCdo -> setBarrelDisableProbability(m_BarrelDisableProbabilityITK);
+    writeCdo -> setDefaultBarrelBiasVoltage(m_BarrelBiasVoltageITK);
 
     writeCdo -> setEndcapToTThreshold(m_EndcapToTThresholdITK);
     writeCdo -> setEndcapCrossTalk(m_EndcapCrossTalkITK);
     writeCdo -> setEndcapNoiseOccupancy(m_EndcapNoiseOccupancyITK);
     writeCdo -> setEndcapDisableProbability(m_EndcapDisableProbabilityITK);
+    writeCdo -> setDefaultEndcapBiasVoltage(m_EndcapBiasVoltageITK);
 
     // This is ad-hoc solution.
     for (size_t i=0; i<m_InnermostNoiseShapeITK.size(); i++)     { writeCdo->setBarrelNoiseShape(0,m_InnermostNoiseShapeITK[i]); }
@@ -266,8 +292,15 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     for (size_t i=0; i<m_EndcapToTThresholdITK.size(); i++) {
       for (size_t j=0; j<m_PixelNoiseShapeITK.size(); j++)  { writeCdo->setEndcapNoiseShape(i,m_PixelNoiseShapeITK[j]); }
     }
+
+    // Radiation damage simulation
+    writeCdo -> setFluenceLayer(m_BarrelFluenceITK);
+    for (size_t i=0; i<m_BarrelFluenceMapITK.size(); i++) {
+      mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMapITK[i]));
+    }
+
   }
-  else if (currentRunNumber<300000) {  // RUN2 2015/2016
+  else if (currentRunNumber<300000) {  // RUN2 2015/2016 (mc16a)
     writeCdo -> setBarrelToTThreshold(m_BarrelToTThreshold2016);
     writeCdo -> setFEI3BarrelLatency(m_FEI3BarrelLatency2016);
     writeCdo -> setFEI3BarrelHitDuplication(m_FEI3BarrelHitDuplication2016);
@@ -276,6 +309,7 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     writeCdo -> setBarrelCrossTalk(m_BarrelCrossTalk2016);
     writeCdo -> setBarrelNoiseOccupancy(m_BarrelNoiseOccupancy2016);
     writeCdo -> setBarrelDisableProbability(m_BarrelDisableProbability2016);
+    writeCdo -> setDefaultBarrelBiasVoltage(m_BarrelBiasVoltage2016);
 
     writeCdo -> setEndcapToTThreshold(m_EndcapToTThreshold2016);
     writeCdo -> setFEI3EndcapLatency(m_FEI3EndcapLatency2016);
@@ -285,11 +319,13 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     writeCdo -> setEndcapCrossTalk(m_EndcapCrossTalk2016);
     writeCdo -> setEndcapNoiseOccupancy(m_EndcapNoiseOccupancy2016);
     writeCdo -> setEndcapDisableProbability(m_EndcapDisableProbability2016);
+    writeCdo -> setDefaultEndcapBiasVoltage(m_EndcapBiasVoltage2016);
 
     writeCdo -> setDBMToTThreshold(m_DBMToTThreshold2016);
     writeCdo -> setDBMCrossTalk(m_DBMCrossTalk2016);
     writeCdo -> setDBMNoiseOccupancy(m_DBMNoiseOccupancy2016);
     writeCdo -> setDBMDisableProbability(m_DBMDisableProbability2016);
+    writeCdo -> setDefaultDBMBiasVoltage(m_DBMBiasVoltage2016);
 
     // This is ad-hoc solution.
     for (size_t i=0; i<m_IBLNoiseShape2016.size(); i++)    { writeCdo->setBarrelNoiseShape(0,m_IBLNoiseShape2016[i]); }
@@ -305,8 +341,14 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     for (size_t i=0; i<m_IBLNoiseShape2016.size(); i++)    {
       for (size_t layer:{0,1,2}) { writeCdo->setDBMNoiseShape(layer,m_IBLNoiseShape2016[i]); }
     }
+
+    // Radiation damage simulation
+    writeCdo -> setFluenceLayer(m_BarrelFluence2016);
+    for (size_t i=0; i<m_BarrelFluenceMap2016.size(); i++) {
+      mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMap2016[i]));
+    }
   }
-  else if (currentRunNumber<310000) {  // RUN2 2017
+  else if (currentRunNumber<310000) {  // RUN2 2017 (mc16d)
     writeCdo -> setBarrelToTThreshold(m_BarrelToTThreshold2017);
     writeCdo -> setFEI3BarrelLatency(m_FEI3BarrelLatency2017);
     writeCdo -> setFEI3BarrelHitDuplication(m_FEI3BarrelHitDuplication2017);
@@ -315,6 +357,7 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     writeCdo -> setBarrelCrossTalk(m_BarrelCrossTalk2017);
     writeCdo -> setBarrelNoiseOccupancy(m_BarrelNoiseOccupancy2017);
     writeCdo -> setBarrelDisableProbability(m_BarrelDisableProbability2017);
+    writeCdo -> setDefaultBarrelBiasVoltage(m_BarrelBiasVoltage2017);
 
     writeCdo -> setEndcapToTThreshold(m_EndcapToTThreshold2017);
     writeCdo -> setFEI3EndcapLatency(m_FEI3EndcapLatency2017);
@@ -324,11 +367,13 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     writeCdo -> setEndcapCrossTalk(m_EndcapCrossTalk2017);
     writeCdo -> setEndcapNoiseOccupancy(m_EndcapNoiseOccupancy2017);
     writeCdo -> setEndcapDisableProbability(m_EndcapDisableProbability2017);
+    writeCdo -> setDefaultEndcapBiasVoltage(m_EndcapBiasVoltage2017);
 
     writeCdo -> setDBMToTThreshold(m_DBMToTThreshold2017);
     writeCdo -> setDBMCrossTalk(m_DBMCrossTalk2017);
     writeCdo -> setDBMNoiseOccupancy(m_DBMNoiseOccupancy2017);
     writeCdo -> setDBMDisableProbability(m_DBMDisableProbability2017);
+    writeCdo -> setDefaultDBMBiasVoltage(m_DBMBiasVoltage2017);
 
     // This is ad-hoc solution.
     for (size_t i=0; i<m_IBLNoiseShape2017.size(); i++)    { writeCdo->setBarrelNoiseShape(0,m_IBLNoiseShape2017[i]); }
@@ -344,8 +389,14 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     for (size_t i=0; i<m_IBLNoiseShape2017.size(); i++)    {
       for (size_t layer:{0,1,2}) { writeCdo->setDBMNoiseShape(layer,m_IBLNoiseShape2017[i]); }
     }
+
+    // Radiation damage simulation
+    writeCdo -> setFluenceLayer(m_BarrelFluence2017);
+    for (size_t i=0; i<m_BarrelFluenceMap2017.size(); i++) {
+      mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMap2017[i]));
+    }
   }
-  else {  // RUN2 2018
+  else {  // RUN2 2018 (mc16e)
     writeCdo -> setBarrelToTThreshold(m_BarrelToTThreshold2018);
     writeCdo -> setFEI3BarrelLatency(m_FEI3BarrelLatency2018);
     writeCdo -> setFEI3BarrelHitDuplication(m_FEI3BarrelHitDuplication2018);
@@ -354,6 +405,7 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     writeCdo -> setBarrelCrossTalk(m_BarrelCrossTalk2018);
     writeCdo -> setBarrelNoiseOccupancy(m_BarrelNoiseOccupancy2018);
     writeCdo -> setBarrelDisableProbability(m_BarrelDisableProbability2018);
+    writeCdo -> setDefaultBarrelBiasVoltage(m_BarrelBiasVoltage2018);
 
     writeCdo -> setEndcapToTThreshold(m_EndcapToTThreshold2018);
     writeCdo -> setFEI3EndcapLatency(m_FEI3EndcapLatency2018);
@@ -363,11 +415,13 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     writeCdo -> setEndcapCrossTalk(m_EndcapCrossTalk2018);
     writeCdo -> setEndcapNoiseOccupancy(m_EndcapNoiseOccupancy2018);
     writeCdo -> setEndcapDisableProbability(m_EndcapDisableProbability2018);
+    writeCdo -> setDefaultEndcapBiasVoltage(m_EndcapBiasVoltage2018);
 
     writeCdo -> setDBMToTThreshold(m_DBMToTThreshold2018);
     writeCdo -> setDBMCrossTalk(m_DBMCrossTalk2018);
     writeCdo -> setDBMNoiseOccupancy(m_DBMNoiseOccupancy2018);
     writeCdo -> setDBMDisableProbability(m_DBMDisableProbability2018);
+    writeCdo -> setDefaultDBMBiasVoltage(m_DBMBiasVoltage2018);
 
     // This is ad-hoc solution.
     for (size_t i=0; i<m_IBLNoiseShape2018.size(); i++)    { writeCdo->setBarrelNoiseShape(0,m_IBLNoiseShape2018[i]); }
@@ -383,8 +437,50 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
     for (size_t i=0; i<m_IBLNoiseShape2018.size(); i++)    {
       for (size_t layer:{0,1,2}) { writeCdo->setDBMNoiseShape(layer,m_IBLNoiseShape2018[i]); }
     }
+
+    // Radiation damage simulation
+    writeCdo -> setFluenceLayer(m_BarrelFluence2018);
+    for (size_t i=0; i<m_BarrelFluenceMap2018.size(); i++) {
+      mapsPath_list.push_back(PathResolverFindCalibFile(m_BarrelFluenceMap2018[i]));
+    }
   }
 
+  std::vector<TH3F*> ramoPotentialMap;
+  std::vector<TH2F*> lorentzMap_e;
+  std::vector<TH2F*> lorentzMap_h;
+  std::vector<TH2F*> distanceMap_e;
+  std::vector<TH2F*> 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);
+    std::unique_ptr<TFile> mapsFile(new TFile((mapsPath_list.at(i)).c_str())); //this is the ramo potential.
+
+    //Setup ramo weighting field map
+    TH3F* ramoPotentialMap_hold = 0;
+    ramoPotentialMap_hold = (TH3F*)mapsFile->Get("hramomap1");
+    if (ramoPotentialMap_hold==0) { ramoPotentialMap_hold=(TH3F*)mapsFile->Get("ramo3d"); }
+    if (ramoPotentialMap_hold==0) {
+      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);
+
+    std::unique_ptr<TH2F> lorentzMap_e_hold((TH2F*)mapsFile->Get("lorentz_map_e"));
+    std::unique_ptr<TH2F> lorentzMap_h_hold((TH2F*)mapsFile->Get("lorentz_map_h"));
+    std::unique_ptr<TH2F> distanceMap_h_hold((TH2F*)mapsFile->Get("edistance"));
+    std::unique_ptr<TH2F> distanceMap_e_hold((TH2F*)mapsFile->Get("hdistance"));
+
+    lorentzMap_e.push_back(lorentzMap_e_hold.get());
+    lorentzMap_h.push_back(lorentzMap_h_hold.get());
+    distanceMap_e.push_back(distanceMap_e_hold.get());
+    distanceMap_h.push_back(distanceMap_h_hold.get());
+  }
+  writeCdo -> setLorentzMap_e(lorentzMap_e);
+  writeCdo -> setLorentzMap_h(lorentzMap_h);
+  writeCdo -> setDistanceMap_e(distanceMap_e);
+  writeCdo -> setDistanceMap_h(distanceMap_h);
+  writeCdo -> setRamoPotentialMap(ramoPotentialMap);
+
   //=======================
   // Combine time interval
   //=======================
@@ -405,3 +501,4 @@ StatusCode PixelConfigCondAlg::execute(const EventContext& ctx) const {
   return StatusCode::SUCCESS;
 }
 
+
diff --git a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.h b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.h
index 2f22aa8a4a8fee954cae3e047d89088be4a6f435..458eb01203751236cf324404dc2997b48c893e1b 100644
--- a/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.h
+++ b/InnerDetector/InDetConditions/PixelConditionsAlgorithms/src/PixelConfigCondAlg.h
@@ -12,6 +12,7 @@
 #define PIXELCONFIGCONDALG
 
 #include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "GaudiKernel/ToolHandle.h"
 
 #include "StoreGate/ReadCondHandleKey.h"
 #include "AthenaPoolUtilities/CondAttrListCollection.h"
@@ -87,7 +88,9 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     //
     // MC Project:               RUN1            RUN2 mc16a            RUN2 mc16d            RUN2 mc16e 
     //  Year:                  - 2014             2015/2016                  2017                  2018
-    //  Run Number:           <222222                284500                300000                310000
+    // MC Run Number:         <222222                284500                300000                310000
+    // Reference run#:            ---                303638                336506                357193
+    // Luminosity(fb-1):                               17.3                  69.0                 119.4
     //                                 
     // Barrel:                             
     //  ToT:         [   3,   3,   3] [  -1,   5,   5,   5] [  -1,   5,   5,   5] [  -1,   3,   5,   5]
@@ -99,7 +102,9 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     //  NoiseOcc.:   [5e-8,5e-8,5e-8] [5e-8,5e-8,5e-8,5e-8] [5e-8,5e-8,5e-8,5e-8] [5e-8,5e-8,5e-8,5e-8]
     //  DisalbePix:  [9e-3,9e-3,9e-3] [9e-3,9e-3,9e-3,9e-3] [9e-3,9e-3,9e-3,9e-3] [9e-3,9e-3,9e-3,9e-3]
     //  NoiseShape:  [2018,2018,2018] [2018,2018,2018,2018] [2018,2018,2018,2018] [2018,2018,2018,2018]
-    //                                     
+    //  BiasVoltage: [ 150, 150, 150] [  80, 350, 200, 150] [ 350, 350, 200, 150] [ 400, 400, 250, 250]
+    //  Fluence(e14):[1.01,0.44,0.30] [0.80,1.61,0.71,0.48] [3.18,3.42,1.50,1.01] [5.50,5.19,2.28,1.53]
+    //
     // Endcap:                             
     //  ToT:         [   3,   3,   3]      [   5,   5,   5]      [   5,   5,   5]      [   5,   5,   5]
     //  Latency:     [ 256, 256, 256]      [ 256, 256, 256]      [ 256, 256, 256]      [ 256, 256, 256]
@@ -110,12 +115,15 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     //  NoiseOcc.:   [5e-8,5e-8,5e-8]      [5e-8,5e-8,5e-8]      [5e-8,5e-8,5e-8]      [5e-8,5e-8,5e-8]
     //  DisalbePix:  [9e-3,9e-3,9e-3]      [9e-3,9e-3,9e-3]      [9e-3,9e-3,9e-3]      [9e-3,9e-3,9e-3]
     //  NoiseShape:  [2018,2018,2018]      [2018,2018,2018]      [2018,2018,2018]      [2018,2018,2018]
-    //                                
+    //  BiasVoltage: [ 150, 150, 150]      [ 150, 150, 150]      [ 150, 150, 150]      [ 250, 250, 250]
+    //  Fluence(e14):[ n/a, n/a, n/a]      [ n/a, n/a, n/a]      [ n/a, n/a, n/a]      [ n/a, n/a, n/a]
+    //
     // DBM:                           
     //  ToT:                    [N/A]      [  -1,  -1,  -1]      [  -1,  -1,  -1]      [  -1,  -1,  -1]
     //  CrossTalk:              [N/A]      [0.06,0.06,0.06]      [0.06,0.06,0.06]      [0.06,0.06,0.06]
     //  NoiseOcc.:              [N/A]      [5e-8,5e-8,5e-8]      [5e-8,5e-8,5e-8]      [5e-8,5e-8,5e-8]
     //  DisalbePix:             [N/A]      [9e-3,9e-3,9e-3]      [9e-3,9e-3,9e-3]      [9e-3,9e-3,9e-3]
+    //  BiasVoltage:            [N/A]      [ 500, 500, 500]      [ 500, 500, 500]      [ 500, 500, 500]
     //
     // See  https://twiki.cern.ch/twiki/bin/view/Atlas/PixelConditionsRUN2
     // for further details.
@@ -133,6 +141,8 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     //  NoiseOcc.:   [5e-8,5e-8,5e-8,5e-8,5e-8]
     //  DisalbePix:  [9e-3,9e-3,9e-3,9e-3,9e-3]
     //  NoiseShape:  [2018,2018,2018,2018,2018]
+    //  BiasVoltage: [ 150, 150, 150, 150, 150]
+    //  Fluence(e14):[ 0.0, 0.0, 0.0, 0.0, 0.0]
     //                               
     // Endcap:                       
     //  ToT:         [   3,   3,   3,   3,   3,   3,   3,   3,   3,   3,   3,   3,   3,   3]
@@ -140,15 +150,29 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     //  NoiseOcc.:   [5e-8,5e-8,5e-8,5e-8,5e-8,5e-8,5e-8,5e-8,5e-8,5e-8,5e-8,5e-8,5e-8,5e-8]
     //  DisalbePix:  [9e-3,9e-3,9e-3,9e-3,9e-3,9e-3,9e-3,9e-3,9e-3,9e-3,9e-3,9e-3,9e-3,9e-3]
     //  NoiseShape:  [2018,2018,2018,2018,2018,2018,2018,2018,2018,2018,2018,2018,2018,2018]
+    //  BiasVoltage: [ 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150]
+    //  Fluence(e14):[ n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a, n/a]
     //
     //====================================================================================
 
-    // Variable definition is a bit ugly here...
     //====================================================================================
     // Barrel RUN2 2015/2016
     Gaudi::Property<std::vector<int>> m_BarrelToTThreshold2016
     {this, "BarrelToTThreshold2016", {-1, 5, 5, 5}, "ToT thresholds for barrel pixel layers in 2015/2016"};
 
+    Gaudi::Property<std::vector<float>> m_BarrelBiasVoltage2016
+    {this, "DefaultBarrelBiasVoltage2016", {80.0,350.0,200.0,150.0}, "Default barrel bias voltage in 2015/2016"};
+
+    Gaudi::Property<std::vector<double>> m_BarrelFluence2016
+    {this, "BarrelFluence2016", {0.80e14, 1.61e14, 0.71e14, 0.48e14}, "Barrel fluence for radiation damage in 2016"};
+
+    Gaudi::Property<std::vector<std::string>> m_BarrelFluenceMap2016
+    {this, "BarrelFluenceMap2016", {"PixelDigitization/maps_IBL_PL_80V_fl0_8e14.root",
+                                    "PixelDigitization/maps_PIX_350V_fl1_61e14.root",
+                                    "PixelDigitization/maps_PIX_200V_fl0_71e14.root",
+                                    "PixelDigitization/maps_PIX_150V_fl0_48e14.root"},
+                                    "Barrel fluence map for radiation damage in 2016"};
+
     Gaudi::Property<std::vector<int>> m_FEI3BarrelLatency2016
     {this, "FEI3BarrelLatency2016", {  0,151,256,256}, "FEI3 latency for barrel pixel layers in 2015/2016"};
 
@@ -188,6 +212,9 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_EndcapToTThreshold2016
     {this, "EndcapToTThreshold2016", { 5, 5, 5}, "ToT thresholds for endcap pixel layers in 2015/2016"};
 
+    Gaudi::Property<std::vector<float>> m_EndcapBiasVoltage2016
+    {this, "DefaultEndcapBiasVoltage2016", {150.0,150.0,150.0}, "Default endcap bias voltage in 2015/2016"};
+
     Gaudi::Property<std::vector<int>> m_FEI3EndcapLatency2016
     {this, "FEI3EndcapLatency2016", {256,256,256}, "FEI3 latency for endcap pixel layers in 2015/2016"};
 
@@ -213,6 +240,9 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_DBMToTThreshold2016
     {this, "DBMToTThreshold2016", {-1,-1,-1}, "ToT thresholds for DBM layers in 2015/2016"};
 
+    Gaudi::Property<std::vector<float>> m_DBMBiasVoltage2016
+    {this, "DefaultDBMBiasVoltage2016", {500.0,500.0,500.0}, "Default DBM bias voltage in 2015/2016"};
+
     Gaudi::Property<std::vector<double>> m_DBMCrossTalk2016
     {this, "DBMCrossTalk2016", {0.06,0.06,0.06}, "Cross-talk fraction of barrel DBM layers in 2015/2016"};
 
@@ -227,6 +257,19 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_BarrelToTThreshold2017
     {this, "BarrelToTThreshold2017", {-1, 5, 5, 5}, "ToT thresholds for barrel pixel layers in 2017"};
 
+    Gaudi::Property<std::vector<float>> m_BarrelBiasVoltage2017
+    {this, "DefaultBarrelBiasVoltage2017", {350.0,350.0,200.0,150.0}, "Default barrel bias voltage in 2017"};
+
+    Gaudi::Property<std::vector<double>> m_BarrelFluence2017
+    {this, "BarrelFluence2017", {3.18e14, 3.42e14, 1.50e14, 1.01e14}, "Barrel fluence for radiation damage in 2017"};
+
+    Gaudi::Property<std::vector<std::string>> m_BarrelFluenceMap2017
+    {this, "BarrelFluenceMap2017", {"PixelDigitization/maps_IBL_PL_350V_fl3_18e14.root",
+                                    "PixelDigitization/maps_PIX_350V_fl3_42e14.root",
+                                    "PixelDigitization/maps_PIX_200V_fl1_5e14.root",
+                                    "PixelDigitization/maps_PIX_150V_fl1_01e14.root"},
+                                    "Barrel fluence map for radiation damage in 2017"};
+
     Gaudi::Property<std::vector<int>> m_FEI3BarrelLatency2017
     {this, "FEI3BarrelLatency2017", {  0,151,256,256}, "FEI3 latency for barrel pixel layers in 2017"};
 
@@ -262,6 +305,9 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_EndcapToTThreshold2017
     {this, "EndcapToTThreshold2017", { 5, 5, 5}, "ToT thresholds for endcap pixel layers in 2017"};
 
+    Gaudi::Property<std::vector<float>> m_EndcapBiasVoltage2017
+    {this, "DefaultEndcapBiasVoltage2017", {150.0,150.0,150.0}, "Default endcap bias voltage in 2017"};
+
     Gaudi::Property<std::vector<int>> m_FEI3EndcapLatency2017
     {this, "FEI3EndcapLatency2017", {256,256,256}, "FEI3 latency for endcap pixel layers in 2017"};
 
@@ -287,6 +333,9 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_DBMToTThreshold2017
     {this, "DBMToTThreshold2017", {-1,-1,-1}, "ToT thresholds for DBM layers in 2017"};
 
+    Gaudi::Property<std::vector<float>> m_DBMBiasVoltage2017
+    {this, "DefaultDBMBiasVoltage2017", {500.0,500.0,500.0}, "Default DBM bias voltage in 2017"};
+
     Gaudi::Property<std::vector<double>> m_DBMCrossTalk2017
     {this, "DBMCrossTalk2017", {0.06,0.06,0.06}, "Cross-talk fraction of barrel DBM layers in 2017"};
 
@@ -301,6 +350,19 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_BarrelToTThreshold2018
     {this, "BarrelToTThreshold2018", {-1, 3, 5, 5}, "ToT thresholds for barrel pixel layers in 2018"};
 
+    Gaudi::Property<std::vector<float>> m_BarrelBiasVoltage2018
+    {this, "DefaultBarrelBiasVoltage2018", {400.0,400.0,250.0,250.0}, "Default barrel bias voltage in 2018"};
+
+    Gaudi::Property<std::vector<double>> m_BarrelFluence2018
+    {this, "BarrelFluence2018", {5.50e14, 5.19e14, 2.28e14, 1.53e14}, "Barrel fluence for radiation damage in 2018"};
+
+    Gaudi::Property<std::vector<std::string>> m_BarrelFluenceMap2018
+    {this, "BarrelFluenceMap2018", {"PixelDigitization/maps_IBL_PL_400V_fl5_5e14.root",
+                                    "PixelDigitization/maps_PIX_400V_fl5_19e14.root",
+                                    "PixelDigitization/maps_PIX_250V_fl2_28e14.root",
+                                    "PixelDigitization/maps_PIX_250V_fl1_53e14.root"},
+                                    "Barrel fluence map for radiation damage in 2018"};
+
     Gaudi::Property<std::vector<int>> m_FEI3BarrelLatency2018
     {this, "FEI3BarrelLatency2018", {  0,151,256,256}, "FEI3 latency for barrel pixel layers in 2018"};
 
@@ -336,6 +398,9 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_EndcapToTThreshold2018
     {this, "EndcapToTThreshold2018", { 5, 5, 5}, "ToT thresholds for endcap pixel layers in 2018"};
 
+    Gaudi::Property<std::vector<float>> m_EndcapBiasVoltage2018
+    {this, "DefaultEndcapBiasVoltage2018", {250.0,250.0,250.0}, "Default endcap bias voltage in 2018"};
+
     Gaudi::Property<std::vector<int>> m_FEI3EndcapLatency2018
     {this, "FEI3EndcapLatency2018", {256,256,256}, "FEI3 latency for endcap pixel layers in 2018"};
 
@@ -361,6 +426,9 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_DBMToTThreshold2018
     {this, "DBMToTThreshold2018", {-1,-1,-1}, "ToT thresholds for DBM layers in 2018"};
 
+    Gaudi::Property<std::vector<float>> m_DBMBiasVoltage2018
+    {this, "DefaultDBMBiasVoltage2018", {500.0,500.0,500.0}, "Default DBM bias voltage in 2018"};
+
     Gaudi::Property<std::vector<double>> m_DBMCrossTalk2018
     {this, "DBMCrossTalk2018", {0.06,0.06,0.06}, "Cross-talk fraction of barrel DBM layers in 2018"};
 
@@ -375,6 +443,18 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_BarrelToTThresholdRUN1
     {this, "BarrelToTThresholdRUN1", {3, 3, 3}, "ToT thresholds for barrel pixel layers in RUN1"};
 
+    Gaudi::Property<std::vector<float>> m_BarrelBiasVoltageRUN1
+    {this, "DefaultBarrelBiasVoltageRUN1", {150.0,150.0,150.0}, "Default barrel bias voltage in RUN1"};
+
+    Gaudi::Property<std::vector<double>> m_BarrelFluenceRUN1
+    {this, "BarrelFluenceRUN1", {1.01e14, 0.44e14, 0.30e14}, "Barrel fluence for radiation damage in RUN1"};
+
+    Gaudi::Property<std::vector<std::string>> m_BarrelFluenceMapRUN1
+    {this, "BarrelFluenceMapRUN1", {"PixelDigitization/maps_PIX_250V_fl1_01e14.root",   // this needs to be updated
+                                    "PixelDigitization/maps_PIX_150V_fl0_44e14.root",
+                                    "PixelDigitization/maps_PIX_150V_fl0_3e14.root"},
+                                    "Barrel fluence map for radiation damage in RUN1"};
+
     Gaudi::Property<std::vector<int>> m_FEI3BarrelLatencyRUN1
     {this, "FEI3BarrelLatencyRUN1", {256,256,256}, "FEI3 latency for barrel pixel layers in RUN1"};
 
@@ -407,6 +487,9 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_EndcapToTThresholdRUN1
     {this, "EndcapToTThresholdRUN1", { 3, 3, 3}, "ToT thresholds for endcap pixel layers in RUN1"};
 
+    Gaudi::Property<std::vector<float>> m_EndcapBiasVoltageRUN1
+    {this, "DefaultEndcapBiasVoltageRUN1", {150.0,150.0,150.0}, "Default endcap bias voltage in RUN1"};
+
     Gaudi::Property<std::vector<int>> m_FEI3EndcapLatencyRUN1
     {this, "FEI3EndcapLatencyRUN1", {256,256,256}, "FEI3 latency for endcap pixel layers in RUN1"};
 
@@ -433,6 +516,20 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_BarrelToTThresholdITK
     {this, "BarrelToTThresholdITK", {3,3,3,3,3}, "ToT thresholds for barrel pixel layers in ITK"};
 
+    Gaudi::Property<std::vector<float>> m_BarrelBiasVoltageITK
+    {this, "DefaultBarrelBiasVoltageITK", {150.0,150.0,150.0,150.0,150.0}, "Default barrel bias voltage in ITK"};
+
+    Gaudi::Property<std::vector<double>> m_BarrelFluenceITK
+    {this, "BarrelFluenceITK", {0.0e14, 0.0e14, 0.0e14, 0.0e14, 0.0e14}, "Barrel fluence for radiation damage in ITK"};
+
+    Gaudi::Property<std::vector<std::string>> m_BarrelFluenceMapITK
+    {this, "BarrelFluenceMapITK", {"PixelDigitization/maps_IBL_PL_80V_fl0e14.root",   // this needs to be updated
+                                   "PixelDigitization/maps_IBL_PL_80V_fl0e14.root",
+                                   "PixelDigitization/maps_IBL_PL_80V_fl0e14.root",
+                                   "PixelDigitization/maps_IBL_PL_80V_fl0e14.root",
+                                   "PixelDigitization/maps_IBL_PL_80V_fl0e14.root"},
+                                   "Barrel fluence map for radiation damage in ITK"};
+
     Gaudi::Property<std::vector<double>> m_BarrelCrossTalkITK
     {this, "BarrelCrossTalkITK", {0.06,0.06,0.06,0.06,0.06}, "Cross-talk fraction of barrel pixel layers in ITK"};
 
@@ -456,6 +553,9 @@ class PixelConfigCondAlg : public AthReentrantAlgorithm {
     Gaudi::Property<std::vector<int>> m_EndcapToTThresholdITK
     {this, "EndcapToTThresholdITK", {3,3,3,3,3,3,3,3,3,3,3,3,3,3}, "ToT thresholds for endcap pixel layers in ITK"};
 
+    Gaudi::Property<std::vector<float>> m_EndcapBiasVoltageITK
+    {this, "DefaultEndcapBiasVoltageITK", {150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0}, "Default endcap bias voltage in ITK"};
+
     Gaudi::Property<std::vector<double>> m_EndcapCrossTalkITK
     {this, "EndcapCrossTalkITK", {0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06}, "Cross-talk fraction of barrel endcap layers in ITK"};
 
diff --git a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h
index 06f32bd7addfc29daa2632f0f4ed559735e847a3..be5a24afc8a43c6c923beb43104c8d61b6e69426 100644
--- a/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h
+++ b/InnerDetector/InDetConditions/PixelConditionsData/PixelConditionsData/PixelModuleData.h
@@ -15,6 +15,9 @@
 #include <map>
 
 #include "AthenaKernel/CondCont.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
 
 class PixelModuleData {
   public:
@@ -140,6 +143,11 @@ class PixelModuleData {
     void setDefaultBiasVoltage(float biasVoltage);
     float getDefaultBiasVoltage() const;
 
+    void setDefaultBarrelBiasVoltage(std::vector<float> BarrelBiasVoltage);
+    void setDefaultEndcapBiasVoltage(std::vector<float> EndcapBiasVoltage);
+    void setDefaultDBMBiasVoltage(std::vector<float>    DBMBiasVoltage);
+    float getDefaultBiasVoltage(int bec, int layer) const;
+
     void setDefaultTemperature(float temperature);
     float getDefaultTemperature() const;
 
@@ -150,6 +158,23 @@ class PixelModuleData {
     void setCablingMapFileName(std::string cablingMapFileName);
     std::string getCablingMapFileName() const;
 
+    // Map for radiation damage simulation
+    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 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 setRamoPotentialMap(std::vector<TH3F*> ramoPotentialMap);
+    TH3F* getRamoPotentialMap(int layer) const;
+
     // Distortion parameters
     void setDistortionInputSource(int distortionInputSource);
     int getDistortionInputSource() const;
@@ -261,9 +286,20 @@ class PixelModuleData {
     float m_biasVoltage;
     float m_temperature;
 
+    std::vector<float> m_BarrelBiasVoltage;
+    std::vector<float> m_EndcapBiasVoltage;
+    std::vector<float> m_DBMBiasVoltage;
+
     bool        m_cablingMapToFile;
     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;
+
     int    m_distortionInputSource;
     int    m_distortionVersion;
     double m_distortionR1;
diff --git a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
index 0877f64ae306aa81dcbab238ddeafe389c192653..558bca2c312d1b7e1b3d39068040964fdbc1c16c 100644
--- a/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
+++ b/InnerDetector/InDetConditions/PixelConditionsData/src/PixelModuleData.cxx
@@ -284,6 +284,17 @@ float PixelModuleData::getDefaultQ2TotC() const { return m_paramC; }
 void PixelModuleData::setDefaultBiasVoltage(float biasVoltage) { m_biasVoltage=biasVoltage; }
 float PixelModuleData::getDefaultBiasVoltage() const { return m_biasVoltage; }
 
+void PixelModuleData::setDefaultBarrelBiasVoltage(std::vector<float> BarrelBiasVoltage) { m_BarrelBiasVoltage = BarrelBiasVoltage; }
+void PixelModuleData::setDefaultEndcapBiasVoltage(std::vector<float> EndcapBiasVoltage) { m_EndcapBiasVoltage = EndcapBiasVoltage; }
+void PixelModuleData::setDefaultDBMBiasVoltage(std::vector<float>    DBMBiasVoltage)    { m_DBMBiasVoltage = DBMBiasVoltage; }
+float PixelModuleData::getDefaultBiasVoltage(int bec, int layer) const {
+  float biasVoltage = 0.0;
+  if (std::abs(bec)==0 && layer<(int)m_BarrelBiasVoltage.size()) { biasVoltage=m_BarrelBiasVoltage.at(layer); }
+  if (std::abs(bec)==2 && layer<(int)m_EndcapBiasVoltage.size()) { biasVoltage=m_EndcapBiasVoltage.at(layer); }
+  if (std::abs(bec)==4 && layer<(int)m_DBMBiasVoltage.size())    { biasVoltage=m_DBMBiasVoltage.at(layer); }
+  return biasVoltage;
+}
+
 void PixelModuleData::setDefaultTemperature(float temperature) { m_temperature=temperature; }
 float PixelModuleData::getDefaultTemperature() const { return m_temperature; }
 
@@ -294,6 +305,23 @@ bool PixelModuleData::getCablingMapToFile() const { return m_cablingMapToFile; }
 void PixelModuleData::setCablingMapFileName(std::string cablingMapFileName) { m_cablingMapFileName = cablingMapFileName; }
 std::string PixelModuleData::getCablingMapFileName() const { return m_cablingMapFileName; }
 
+// Map for radiation damage simulation
+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::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); }
+
 // Distortion parameters
 void PixelModuleData::setDistortionInputSource(int distortionInputSource) { m_distortionInputSource = distortionInputSource; }
 int PixelModuleData::getDistortionInputSource() const { return m_distortionInputSource; }
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py
index fb8f71212c7585afccdd24cd908cef4492855855..7015d588204d9d6bb41189907f0d9674368a6627 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py
+++ b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py
@@ -54,15 +54,6 @@ def SensorSimPlanarTool(name="SensorSimPlanarTool", **kwargs):
     kwargs.setdefault("SiPropertiesTool", ToolSvc.PixelSiPropertiesTool)
     kwargs.setdefault("LorentzAngleTool", ToolSvc.PixelLorentzAngleTool)
     kwargs.setdefault("doRadDamage", digitizationFlags.doRadiationDamage.get_Value())
-    # TODO This is 2018 fluence setting. These parameters should be controlled by the conditions data.
-    kwargs.setdefault("fluence", 6.4)
-    kwargs.setdefault("fluenceB", 4.6)
-    kwargs.setdefault("fluence1", 2.1)
-    kwargs.setdefault("fluence2", 1.3)
-    kwargs.setdefault("voltage", 400)
-    kwargs.setdefault("voltageB", 400)
-    kwargs.setdefault("voltage1", 250)
-    kwargs.setdefault("voltage2", 250)
     return CfgMgr.SensorSimPlanarTool(name, **kwargs)
 
 def SensorSim3DTool(name="SensorSim3DTool", **kwargs):
@@ -204,15 +195,23 @@ def PixelConfigCondAlg_MC():
     #====================================================================================
     # Run-dependent SIMULATION(digitization) parameters:
     #====================================================================================
-    # RUN2 2015/2016
-    alg.BarrelToTThreshold2016       = [   -1,    5,    5,    5]
-    alg.FEI3BarrelLatency2016        = [    0,  151,  256,  256]
-    alg.FEI3BarrelHitDuplication2016 = [False,False,False,False]
-    alg.FEI3BarrelSmallHitToT2016    = [   -1,   -1,   -1,   -1]
-    alg.FEI3BarrelTimingSimTune2016  = [   -1, 2015, 2015, 2015]
-    alg.BarrelCrossTalk2016          = [ 0.30, 0.06, 0.06, 0.06]
-    alg.BarrelNoiseOccupancy2016     = [ 5e-8, 5e-8, 5e-8, 5e-8]
-    alg.BarrelDisableProbability2016 = [ 9e-3, 9e-3, 9e-3, 9e-3]
+    # RUN2 2015/2016 (mc16a)
+    # The pixel conditions are matched with 2016 data (mc16a) at L=17.3fb-1 (run#303638).
+    alg.BarrelToTThreshold2016       = [     -1,      5,      5,      5]
+    alg.FEI3BarrelLatency2016        = [      0,    151,    256,    256]
+    alg.FEI3BarrelHitDuplication2016 = [  False,  False,  False,  False]
+    alg.FEI3BarrelSmallHitToT2016    = [     -1,     -1,     -1,     -1]
+    alg.FEI3BarrelTimingSimTune2016  = [     -1,   2015,   2015,   2015]
+    alg.BarrelCrossTalk2016          = [   0.30,   0.06,   0.06,   0.06]
+    alg.BarrelNoiseOccupancy2016     = [   5e-8,   5e-8,   5e-8,   5e-8]
+    alg.BarrelDisableProbability2016 = [   9e-3,   9e-3,   9e-3,   9e-3]
+    alg.DefaultBarrelBiasVoltage2016 = [   80.0,  350.0,  200.0,  150.0]
+    alg.BarrelFluence2016            = [0.80e14,1.61e14,0.71e14,0.48e14]
+
+    alg.BarrelFluenceMap2016 = ["PixelDigitization/maps_IBL_PL_80V_fl0_8e14.root",
+                                "PixelDigitization/maps_PIX_350V_fl1_61e14.root",
+                                "PixelDigitization/maps_PIX_200V_fl0_71e14.root",
+                                "PixelDigitization/maps_PIX_150V_fl0_48e14.root"]
 
     alg.EndcapToTThreshold2016       = [    5,    5,    5]
     alg.FEI3EndcapLatency2016        = [  256,  256,  256]
@@ -222,11 +221,13 @@ def PixelConfigCondAlg_MC():
     alg.EndcapCrossTalk2016          = [ 0.06, 0.06, 0.06]
     alg.EndcapNoiseOccupancy2016     = [ 5e-8, 5e-8, 5e-8]
     alg.EndcapDisableProbability2016 = [ 9e-3, 9e-3, 9e-3]
+    alg.DefaultEndcapBiasVoltage2016 = [150.0,150.0,150.0]
 
     alg.DBMToTThreshold2016       = [   -1,   -1,   -1]
     alg.DBMCrossTalk2016          = [ 0.06, 0.06, 0.06]
     alg.DBMNoiseOccupancy2016     = [ 5e-8, 5e-8, 5e-8]
     alg.DBMDisableProbability2016 = [ 9e-3, 9e-3, 9e-3]
+    alg.DefaultDBMBiasVoltage2016 = [500.0,500.0,500.0]
 
     alg.IBLNoiseShape2016    = [0.0, 0.0330, 0.0, 0.3026, 0.5019, 0.6760, 0.8412, 0.9918, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
     alg.BLayerNoiseShape2016 = [0.0, 0.0, 0.0, 0.0, 0.2204, 0.5311, 0.7493, 0.8954, 0.9980, 1.0]
@@ -239,15 +240,23 @@ def PixelConfigCondAlg_MC():
     #                                     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1735, 0.3380, 0.4733, 0.5829, 0.6730, 0.7516, 0.8234, 0.8916, 0.9595, 1.0]]
 
     #====================================================================================
-    # RUN2 2017
-    alg.BarrelToTThreshold2017       = [   -1,    5,    5,    5]
-    alg.FEI3BarrelLatency2017        = [    0,  151,  256,  256]
-    alg.FEI3BarrelHitDuplication2017 = [False,False,False,False]
-    alg.FEI3BarrelSmallHitToT2017    = [   -1,   -1,   -1,   -1]
-    alg.FEI3BarrelTimingSimTune2017  = [   -1, 2018, 2018, 2018]
-    alg.BarrelCrossTalk2017          = [ 0.30, 0.06, 0.06, 0.06]
-    alg.BarrelNoiseOccupancy2017     = [ 5e-8, 5e-8, 5e-8, 5e-8]
-    alg.BarrelDisableProbability2017 = [ 9e-3, 9e-3, 9e-3, 9e-3]
+    # RUN2 2017 (mc16d)
+    # The pixel conditions are matched with 2017 data (mc16d) at L=69.0fb-1 (run#336506).
+    alg.BarrelToTThreshold2017       = [     -1,      5,      5,      5]
+    alg.FEI3BarrelLatency2017        = [      0,    151,    256,    256]
+    alg.FEI3BarrelHitDuplication2017 = [  False,  False,  False,  False]
+    alg.FEI3BarrelSmallHitToT2017    = [     -1,     -1,     -1,     -1]
+    alg.FEI3BarrelTimingSimTune2017  = [     -1,   2018,   2018,   2018]
+    alg.BarrelCrossTalk2017          = [   0.30,   0.06,   0.06,   0.06]
+    alg.BarrelNoiseOccupancy2017     = [   5e-8,   5e-8,   5e-8,   5e-8]
+    alg.BarrelDisableProbability2017 = [   9e-3,   9e-3,   9e-3,   9e-3]
+    alg.DefaultBarrelBiasVoltage2017 = [  350.0,  350.0,  200.0,  150.0]
+    alg.BarrelFluence2017            = [3.18e14,3.42e14,1.50e14,1.01e14]
+
+    alg.BarrelFluenceMap2017 = ["PixelDigitization/maps_IBL_PL_350V_fl3_18e14.root",
+                                "PixelDigitization/maps_PIX_350V_fl3_42e14.root",
+                                "PixelDigitization/maps_PIX_200V_fl1_5e14.root",
+                                "PixelDigitization/maps_PIX_150V_fl1_01e14.root"]
 
     alg.EndcapToTThreshold2017       = [    5,    5,    5]
     alg.FEI3EndcapLatency2017        = [  256,  256,  256]
@@ -257,26 +266,36 @@ def PixelConfigCondAlg_MC():
     alg.EndcapCrossTalk2017          = [ 0.06, 0.06, 0.06]
     alg.EndcapNoiseOccupancy2017     = [ 5e-8, 5e-8, 5e-8]
     alg.EndcapDisableProbability2017 = [ 9e-3, 9e-3, 9e-3]
+    alg.DefaultEndcapBiasVoltage2017 = [150.0,150.0,150.0]
 
     alg.DBMToTThreshold2017       = [   -1,   -1,   -1]
     alg.DBMCrossTalk2017          = [ 0.06, 0.06, 0.06]
     alg.DBMNoiseOccupancy2017     = [ 5e-8, 5e-8, 5e-8]
     alg.DBMDisableProbability2017 = [ 9e-3, 9e-3, 9e-3]
+    alg.DefaultDBMBiasVoltage2017 = [500.0,500.0,500.0]
 
     alg.IBLNoiseShape2017    = [0.0, 0.0330, 0.0, 0.3026, 0.5019, 0.6760, 0.8412, 0.9918, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
     alg.BLayerNoiseShape2017 = [0.0, 0.0, 0.0, 0.0, 0.2204, 0.5311, 0.7493, 0.8954, 0.9980, 1.0]
     alg.PixelNoiseShape2017  = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2418, 0.4397, 0.5858, 0.6949, 0.7737, 0.8414, 0.8959, 0.9414, 0.9828, 1.0]
 
     #====================================================================================
-    # RUN2 2018
-    alg.BarrelToTThreshold2018       = [   -1,    3,    5,    5]
-    alg.FEI3BarrelLatency2018        = [    0,  151,  256,  256]
-    alg.FEI3BarrelHitDuplication2018 = [False,False,False,False]
-    alg.FEI3BarrelSmallHitToT2018    = [   -1,   -1,   -1,   -1]
-    alg.FEI3BarrelTimingSimTune2018  = [   -1, 2018, 2018, 2018]
-    alg.BarrelCrossTalk2018          = [ 0.30, 0.06, 0.06, 0.06]
-    alg.BarrelNoiseOccupancy2018     = [ 5e-8, 5e-8, 5e-8, 5e-8]
-    alg.BarrelDisableProbability2018 = [ 9e-3, 9e-3, 9e-3, 9e-3]
+    # RUN2 2018 (mc16e)
+    # The pixel conditions are matched with 2018 data (mc16e) at L=119.4fb-1 (run#357193).
+    alg.BarrelToTThreshold2018       = [     -1,      3,      5,      5]
+    alg.FEI3BarrelLatency2018        = [      0,    151,    256,    256]
+    alg.FEI3BarrelHitDuplication2018 = [  False,  False,  False,  False]
+    alg.FEI3BarrelSmallHitToT2018    = [     -1,     -1,     -1,     -1]
+    alg.FEI3BarrelTimingSimTune2018  = [     -1,   2018,   2018,   2018]
+    alg.BarrelCrossTalk2018          = [   0.30,   0.06,   0.06,   0.06]
+    alg.BarrelNoiseOccupancy2018     = [   5e-8,   5e-8,   5e-8,   5e-8]
+    alg.BarrelDisableProbability2018 = [   9e-3,   9e-3,   9e-3,   9e-3]
+    alg.DefaultBarrelBiasVoltage2018 = [  400.0,  400.0,  250.0,  250.0]
+    alg.BarrelFluence2018            = [5.50e14,5.19e14,2.28e14,1.53e14]
+
+    alg.BarrelFluenceMap2018 = ["PixelDigitization/maps_IBL_PL_400V_fl5_5e14.root",
+                                "PixelDigitization/maps_PIX_400V_fl5_19e14.root",
+                                "PixelDigitization/maps_PIX_250V_fl2_28e14.root",
+                                "PixelDigitization/maps_PIX_250V_fl1_53e14.root"]
 
     alg.EndcapToTThreshold2018       = [    5,    5,    5]
     alg.FEI3EndcapLatency2018        = [  256,  256,  256]
@@ -286,11 +305,13 @@ def PixelConfigCondAlg_MC():
     alg.EndcapCrossTalk2018          = [ 0.06, 0.06, 0.06]
     alg.EndcapNoiseOccupancy2018     = [ 5e-8, 5e-8, 5e-8]
     alg.EndcapDisableProbability2018 = [ 9e-3, 9e-3, 9e-3]
+    alg.DefaultEndcapBiasVoltage2018 = [250.0,250.0,250.0]
 
     alg.DBMToTThreshold2018       = [   -1,   -1,   -1]
     alg.DBMCrossTalk2018          = [ 0.06, 0.06, 0.06]
     alg.DBMNoiseOccupancy2018     = [ 5e-8, 5e-8, 5e-8]
     alg.DBMDisableProbability2018 = [ 9e-3, 9e-3, 9e-3]
+    alg.DefaultDBMBiasVoltage2018 = [500.0,500.0,500.0]
 
     alg.IBLNoiseShape2018    = [0.0, 0.0330, 0.0, 0.3026, 0.5019, 0.6760, 0.8412, 0.9918, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
     alg.BLayerNoiseShape2018 = [0.0, 0.0, 0.0, 0.0, 0.2204, 0.5311, 0.7493, 0.8954, 0.9980, 1.0]
@@ -298,14 +319,20 @@ def PixelConfigCondAlg_MC():
 
     #====================================================================================
     # RUN1
-    alg.BarrelToTThresholdRUN1       = [    3,    3,    3]
-    alg.FEI3BarrelLatencyRUN1        = [  256,  256,  256]
-    alg.FEI3BarrelHitDuplicationRUN1 = [ True, True, True]
-    alg.FEI3BarrelSmallHitToTRUN1    = [    7,    7,    7]
-    alg.FEI3BarrelTimingSimTuneRUN1  = [ 2009, 2009, 2009]
-    alg.BarrelCrossTalkRUN1          = [ 0.06, 0.06, 0.06]
-    alg.BarrelNoiseOccupancyRUN1     = [ 5e-8, 5e-8, 5e-8]
-    alg.BarrelDisableProbabilityRUN1 = [ 9e-3, 9e-3, 9e-3]
+    alg.BarrelToTThresholdRUN1       = [      3,      3,      3]
+    alg.FEI3BarrelLatencyRUN1        = [    256,    256,    256]
+    alg.FEI3BarrelHitDuplicationRUN1 = [   True,   True,   True]
+    alg.FEI3BarrelSmallHitToTRUN1    = [      7,      7,      7]
+    alg.FEI3BarrelTimingSimTuneRUN1  = [   2009,   2009,   2009]
+    alg.BarrelCrossTalkRUN1          = [   0.06,   0.06,   0.06]
+    alg.BarrelNoiseOccupancyRUN1     = [   5e-8,   5e-8,   5e-8]
+    alg.BarrelDisableProbabilityRUN1 = [   9e-3,   9e-3,   9e-3]
+    alg.DefaultBarrelBiasVoltageRUN1 = [  150.0,  150.0,  150.0]
+    alg.BarrelFluenceRUN1            = [1.01e14,0.44e14,0.30e14]
+
+    alg.BarrelFluenceMapRUN1 = ["PixelDigitization/maps_PIX_250V_fl1_01e14.root",
+                                "PixelDigitization/maps_PIX_150V_fl0_44e14.root",
+                                "PixelDigitization/maps_PIX_150V_fl0_3e14.root"]
 
     alg.EndcapToTThresholdRUN1       = [    3,    3,    3]
     alg.FEI3EndcapLatencyRUN1        = [  256,  256,  256]
@@ -315,21 +342,31 @@ def PixelConfigCondAlg_MC():
     alg.EndcapCrossTalkRUN1          = [ 0.06, 0.06, 0.06]
     alg.EndcapNoiseOccupancyRUN1     = [ 5e-8, 5e-8, 5e-8]
     alg.EndcapDisableProbabilityRUN1 = [ 9e-3, 9e-3, 9e-3]
+    alg.DefaultEndcapBiasVoltageRUN1 = [150.0,150.0,150.0]
 
     alg.BLayerNoiseShapeRUN1 = [0.0, 0.0, 0.0, 0.0, 0.2204, 0.5311, 0.7493, 0.8954, 0.9980, 1.0]
     alg.PixelNoiseShapeRUN1  = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2418, 0.4397, 0.5858, 0.6949, 0.7737, 0.8414, 0.8959, 0.9414, 0.9828, 1.0]
 
     #====================================================================================
     # ITK
-    alg.BarrelToTThresholdITK       = [    3,    3,    3,    3,    3]
-    alg.BarrelCrossTalkITK          = [ 0.06, 0.06, 0.06, 0.06, 0.06]
-    alg.BarrelNoiseOccupancyITK     = [ 5e-8, 5e-8, 5e-8, 5e-8, 5e-8]
-    alg.BarrelDisableProbabilityITK = [ 9e-3, 9e-3, 9e-3, 9e-3, 9e-3]
+    alg.BarrelToTThresholdITK       = [     3,     3,     3,     3,     3]
+    alg.BarrelCrossTalkITK          = [  0.06,  0.06,  0.06,  0.06,  0.06]
+    alg.BarrelNoiseOccupancyITK     = [  5e-8,  5e-8,  5e-8,  5e-8,  5e-8]
+    alg.BarrelDisableProbabilityITK = [  9e-3,  9e-3,  9e-3,  9e-3,  9e-3]
+    alg.DefaultBarrelBiasVoltageITK = [ 150.0, 150.0, 150.0, 150.0, 150.0]
+    alg.BarrelFluenceITK            = [0.0e14,0.0e14,0.0e14,0.0e14,0.0e14]
+
+    alg.BarrelFluenceMapITK = ["PixelDigitization/maps_IBL_PL_80V_fl0e14.root",
+                               "PixelDigitization/maps_IBL_PL_80V_fl0e14.root",
+                               "PixelDigitization/maps_IBL_PL_80V_fl0e14.root",
+                               "PixelDigitization/maps_IBL_PL_80V_fl0e14.root",
+                               "PixelDigitization/maps_IBL_PL_80V_fl0e14.root"]
 
     alg.EndcapToTThresholdITK       = [    3,    3,    3,    3,    3,    3,    3,    3,    3,    3,    3,    3,    3,    3]
     alg.EndcapCrossTalkITK          = [ 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06]
     alg.EndcapNoiseOccupancyITK     = [ 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8, 5e-8]
     alg.EndcapDisableProbabilityITK = [ 9e-3, 9e-3, 9e-3, 9e-3, 9e-3, 9e-3, 9e-3, 9e-3, 9e-3, 9e-3, 9e-3, 9e-3, 9e-3, 9e-3]
+    alg.DefaultEndcapBiasVoltageITK = [150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0,150.0]
 
     alg.InnermostNoiseShapeITK     = [0.0, 1.0]
     alg.NextInnermostNoiseShapeITK = [0.0, 1.0]
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py
index 71563d2317a37231cf2d884d8c61a58b05437ed8..efc78d58aa4ff065d20539b90b0d667f19f2c5a7 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py
+++ b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py
@@ -84,15 +84,6 @@ def SensorSimPlanarToolCfg(flags, name="SensorSimPlanarTool", **kwargs):
     kwargs.setdefault("LorentzAngleTool", LorentzTool)
     SensorSimPlanarTool = CompFactory.SensorSimPlanarTool
     kwargs.setdefault("doRadDamage", flags.Digitization.DoRadiationDamage)
-    # TODO This is 2018 fluence setting. These parameters should be controlled by the conditions data.
-    kwargs.setdefault("fluence", 6.4)
-    kwargs.setdefault("fluenceB", 4.6)
-    kwargs.setdefault("fluence1", 2.1)
-    kwargs.setdefault("fluence2", 1.3)
-    kwargs.setdefault("voltage", 400)
-    kwargs.setdefault("voltageB", 400)
-    kwargs.setdefault("voltage1", 250)
-    kwargs.setdefault("voltage2", 250)
     acc.setPrivateTools(SensorSimPlanarTool(name, **kwargs))
     return acc
 
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
index df14c3d5aececbc9c238dc01ee0e660841e95da1..4444967edfbca94fb9f6f4be8adf5a6e473c649a 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
@@ -46,18 +46,7 @@ StatusCode SensorSimPlanarTool::initialize() {
 
   ATH_CHECK(m_lorentzAngleTool.retrieve());
 
-  //Calculate trapping times based on fluence (already includes check for fluence=0)
-  if (m_doRadDamage) {
-    std::pair<double,double> trappingTimes = m_radDamageUtil->getTrappingTimes( m_fluence );
-    m_trappingTimeElectrons = trappingTimes.first;
-    m_trappingTimeHoles = trappingTimes.second;
-  }
-  else {
-    m_fluence = 0;
-  }
-
   //If any fluence or voltage initialized negative use benchmark maps and not interpolation
-  bool doInterpolateEfield = (m_fluence > 0. && m_fluenceB > 0.  && m_fluence1 > 0. && m_fluence2 > 0. && m_voltage > 0. && m_voltageB > 0.  && m_voltage1 > 0. && m_voltage2 > 0.);
   std::vector<std::string> mapsPath_list;
   std::vector<std::string> TCADpath_list;
 
@@ -68,224 +57,54 @@ StatusCode SensorSimPlanarTool::initialize() {
   // For each layer one configuration
   TCADpath_list = {iblFiles, sensorFiles, sensorFiles, sensorFiles};           //IBL - 200um sensor depth, B layer - 20um, layer 1, layer 2
 
-  if (static_cast<int>(m_fluence)==1) {
-    ATH_MSG_INFO("Use benchmark point 1!");
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_IBL_PL_80V_fl0em10.root") );  //IBL  PL - Barrel
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_250V_fl7e13.root") );    //B-Layer - Barrel                                                                                                  
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_150V_fl3e13.root") );    //Layer-1 - Barrel
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_150V_fl2e13.root") );    //Layer-2 - Barrel
-
-    m_fluence_layers.push_back(1e-10);
-    m_fluence_layers.push_back(7e13);
-    m_fluence_layers.push_back(3e13);
-    m_fluence_layers.push_back(2e13);
-
-    m_voltage_layers.push_back(80);
-    m_voltage_layers.push_back(250);
-    m_voltage_layers.push_back(150);
-    m_voltage_layers.push_back(150);
-  }
-  else if (static_cast<int>(m_fluence)==2) {
-    ATH_MSG_INFO("Use benchmark point 2!");
-
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_IBL_PL_80V_fl1e14.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_350V_fl1_2e14.root") );                                                                                                            
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_200V_fl0_5e14.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_150V_fl3e13.root") );
-
-    m_fluence_layers.push_back(1e14);
-    m_fluence_layers.push_back(1.2e14);
-    m_fluence_layers.push_back(5e13);
-    m_fluence_layers.push_back(3e13);
-
-    m_voltage_layers.push_back(80);
-    m_voltage_layers.push_back(350);
-    m_voltage_layers.push_back(200);
-    m_voltage_layers.push_back(150);
-  }
-  else if (static_cast<int>(m_fluence)==3) {
-    ATH_MSG_INFO("Use benchmark point 3!");
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_IBL_PL_80V_fl2e14.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_350V_fl1_7e14.root") );                                                                                                            
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_200V_fl0_7e14.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_150V_fl4e13.root") );
-
-    m_fluence_layers.push_back(2e14);
-    m_fluence_layers.push_back(1.7e14);
-    m_fluence_layers.push_back(7e13);
-    m_fluence_layers.push_back(4e13);
-
-    m_voltage_layers.push_back(80);
-    m_voltage_layers.push_back(350);
-    m_voltage_layers.push_back(200);
-    m_voltage_layers.push_back(150);
-  }
-  else if (static_cast<int>(m_fluence)==4) {
-    ATH_MSG_INFO("Use benchmark point 4!");
-
-    mapsPath_list.push_back(  PathResolverFindCalibFile("PixelDigitization/maps_IBL_PL_150V_fl2e14.root") );
-    mapsPath_list.push_back(  PathResolverFindCalibFile("PixelDigitization/maps_PIX_350V_fl1_7e14.root") );                                                                                                            
-    mapsPath_list.push_back(  PathResolverFindCalibFile("PixelDigitization/maps_PIX_200V_fl0_7e14.root") );
-    mapsPath_list.push_back(  PathResolverFindCalibFile("PixelDigitization/maps_PIX_150V_fl4e13.root") );
-
-    m_fluence_layers.push_back(2e14);
-    m_fluence_layers.push_back(1.7e14);
-    m_fluence_layers.push_back(7e13);
-    m_fluence_layers.push_back(4e13);
-
-    m_voltage_layers.push_back(150);
-    m_voltage_layers.push_back(350);
-    m_voltage_layers.push_back(200);
-    m_voltage_layers.push_back(150);
-  }
-  else if (static_cast<int>(m_fluence)==5) {
-    ATH_MSG_INFO("Use benchmark point 5!");
-
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_IBL_PL_350V_fl5e14.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_350V_fl3_1e14.root") );                                                                                                            
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_200V_fl1_3e14.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_150V_fl8e13.root") );
-
-    m_fluence_layers.push_back(5e14);
-    m_fluence_layers.push_back(3.1e14);
-    m_fluence_layers.push_back(1.3e14);
-    m_fluence_layers.push_back(8e13);
-
-    m_voltage_layers.push_back(350);
-    m_voltage_layers.push_back(350);
-    m_voltage_layers.push_back(200);
-    m_voltage_layers.push_back(150);
-  }
-  else if (static_cast<int>(m_fluence)==6) {
-    ATH_MSG_INFO("Use benchmark point 6!");
-
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_IBL_400V_fl8_7e14.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_400V_fl4_6e14.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_250V_fl2_1e14.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_250V_fl1_3e14.root") );
-
-    m_fluence_layers.push_back(8.7e14);
-    m_fluence_layers.push_back(4.6e14);
-    m_fluence_layers.push_back(2.1e14);
-    m_fluence_layers.push_back(1.3e14);
-
-    m_voltage_layers.push_back(400);
-    m_voltage_layers.push_back(400);
-    m_voltage_layers.push_back(250);
-    m_voltage_layers.push_back(250);
-  }
-  else if (static_cast<int>(m_fluence)==7) {
-    ATH_MSG_INFO("Use benchmark point 7!");
-
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_IBL_endLHC.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_blayer_endLHC.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_L1_endLHC.root") );
-    mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_L2_endLHC.root") );
-
-    m_fluence_layers.push_back(2*8.7e14);
-    m_fluence_layers.push_back(2*4.6e14);
-    m_fluence_layers.push_back(2*2.1e14);
-    m_fluence_layers.push_back(2*1.3e14);
-
-    m_voltage_layers.push_back(400);
-    m_voltage_layers.push_back(400);
-    m_voltage_layers.push_back(250);
-    m_voltage_layers.push_back(250);
-  }
-
-  if (mapsPath_list.size()==0) {
-    if (doInterpolateEfield) {
-      ATH_MSG_INFO("No benchmark value set for fluence. Use interpolation.");
-      mapsPath_list.clear();
-      m_fluence_layers.clear();  
-      m_voltage_layers.clear();
-      //Set up default maps for ramoMap,
-      //but retrieve Efield from interpolation as well as Lorentz, time and distance map from E field
-      mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_IBL_PL_80V_fl0em10.root") );  //IBL  PL - Barrel
-      mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_150V_fl7e13.root") );    //B-Layer - Barrel 
-      mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_150V_fl3e13.root") );    //Layer-1 - Barrel 
-      mapsPath_list.push_back( PathResolverFindCalibFile("PixelDigitization/maps_PIX_150V_fl2e13.root") );    //Layer-2 - Barrel 
-      m_fluence_layers.push_back(m_fluence*1e14);
-      m_fluence_layers.push_back(m_fluenceB*1e14);
-      m_fluence_layers.push_back(m_fluence1*1e14);
-      m_fluence_layers.push_back(m_fluence2*1e14);
-
-      m_voltage_layers.push_back(m_voltage);
-      m_voltage_layers.push_back(m_voltageB);
-      m_voltage_layers.push_back(m_voltage1);
-      m_voltage_layers.push_back(m_voltage2);
-    }
-  }
-  else {
-    if (!doInterpolateEfield) {
-      ATH_MSG_INFO("m_fluence set to becnhmark value. Use becnhmark values for all layers.");
-    }
-    else {
-      ATH_MSG_WARNING("m_fluence set to benchmark value. Fluence and bias voltage for all layers overwritten accordingly.");
+  if (m_doInterpolateEfield) {
+    if (m_fluenceMap.size()==0 || m_fluenceLayer.size()==0 || m_voltageLayer.size()==0 || m_fluenceMap.size()!=m_fluenceLayer.size() || m_fluenceMap.size()!=m_voltageLayer.size()) {
+      ATH_MSG_INFO("Use interpolation, but the input map/fluence/valtage are not set.");
+      return StatusCode::FAILURE;
     }
-  }
-
-  if (m_doRadDamage) {
-    ATH_MSG_INFO("ibl       : Phi=" << m_fluence_layers.at(0) << "e14neq/cm2 U="<< m_voltage_layers.at(0) << "V");
-    ATH_MSG_INFO("B layer   : Phi=" << m_fluence_layers.at(1) << "e14neq/cm2 U="<< m_voltage_layers.at(1) << "V");
-    ATH_MSG_INFO("layer 1   : Phi=" << m_fluence_layers.at(2) << "e14neq/cm2 U="<< m_voltage_layers.at(2) << "V");
-    ATH_MSG_INFO("layer 2   : Phi=" << m_fluence_layers.at(3) << "e14neq/cm2 U="<< m_voltage_layers.at(3) << "V");
-  }
+    ATH_MSG_INFO("No benchmark value set for fluence. Use interpolation.");
 
-  // *****************************
-  // *** Setup Maps ****
-  // *****************************
-  //TODO This is only temporary until remotely stored maps and locally generated maps can be implemented 
-  //E field already implemented: needs fluence and bias voltage given as Property m_fluence, m_fluenceB, ...,  m_fluence1, ...
-  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);
-    if(doInterpolateEfield)ATH_MSG_INFO("Create E field via interpolation based on files from: " << TCADpath_list.at(i));
-    //std::unique_ptr<TFile>  mapsFile=std::make_unique<TFile>( (mapsPath_list.at(i)).c_str() ); //this is the ramo potential.
-    TFile* mapsFile=new TFile( (mapsPath_list.at(i)).c_str() ); //this is the ramo potential.
-
-    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=i;             //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2
-    //IBL Barrel doesn't exist. So the possible idexes should be: 0-0, 0-1, 0-2, 0-3, 1-1, 1-2, 1-3
-
-    //Setup ramo weighting field map
-    TH3F* ramoPotentialMap_hold;
-    ramoPotentialMap_hold=0;
-    ramoPotentialMap_hold=(TH3F*)mapsFile->Get("hramomap1");
-    if (ramoPotentialMap_hold==0) ramoPotentialMap_hold=(TH3F*)mapsFile->Get("ramo3d");
-    if (ramoPotentialMap_hold==0){
-      ATH_MSG_INFO("Did not find a Ramo potential map.  Will use an approximate form.");
-      ATH_MSG_WARNING("Not implemented yet - exit");
-      return StatusCode::FAILURE; //Obviously, remove this when gen. code is set up
-    }
-    m_ramoPotentialMap[Layer]=ramoPotentialMap_hold;
-    m_fluence_layersMaps[Layer]=m_fluence_layers.at(i);
-    //Now setup the E-field.
-    TH1F* eFieldMap_hold;
-    eFieldMap_hold= new TH1F() ;
-    if(doInterpolateEfield){
-      //ATH_MSG_INFO("Generating E field maps using interpolation.");
-      CHECK(m_radDamageUtil->generateEfieldMap( eFieldMap_hold, NULL, m_fluence_layers.at(i), m_voltage_layers.at(i), i, TCADpath_list.at(i), true));
-    }else{
-      //precomputed map
-      eFieldMap_hold=(TH1F*)mapsFile->Get("hEfield1D"); 
+    mapsPath_list.clear();
+    for (size_t i=0; i<m_fluenceMap.size(); i++) {
+      mapsPath_list.push_back(PathResolverFindCalibFile(m_fluenceMap[i]));
     }
 
-    if (eFieldMap_hold == 0){ 
-      ATH_MSG_INFO("Unable to load sensor e-field map.");
-      return StatusCode::FAILURE;
-    }
-    m_eFieldMap[Layer]=eFieldMap_hold;
+    // *****************************
+    // *** Setup Maps ****
+    // *****************************
+    //TODO This is only temporary until remotely stored maps and locally generated maps can be implemented 
+    //E field already implemented: needs fluence and bias voltage given as Property m_fluence, m_fluenceB, ...,  m_fluence1, ...
+    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);
+      ATH_MSG_INFO("Create E field via interpolation based on files from: " << TCADpath_list.at(i));
+      //std::unique_ptr<TFile>  mapsFile=std::make_unique<TFile>( (mapsPath_list.at(i)).c_str() ); //this is the ramo potential.
+      TFile* mapsFile=new TFile( (mapsPath_list.at(i)).c_str() ); //this is the ramo potential.
+
+      //Setup ramo weighting field map
+      TH3F* ramoPotentialMap_hold;
+      ramoPotentialMap_hold=0;
+      ramoPotentialMap_hold=(TH3F*)mapsFile->Get("hramomap1");
+      if (ramoPotentialMap_hold==0) ramoPotentialMap_hold=(TH3F*)mapsFile->Get("ramo3d");
+      if (ramoPotentialMap_hold==0){
+        ATH_MSG_INFO("Did not find a Ramo potential map.  Will use an approximate form.");
+        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);
+      //Now setup the E-field.
+      TH1F* eFieldMap_hold;
+      eFieldMap_hold= new TH1F() ;
+      //ATH_MSG_INFO("Generating E field maps using interpolation.");
+      CHECK(m_radDamageUtil->generateEfieldMap( eFieldMap_hold, NULL, m_fluenceLayer[i], m_voltageLayer[i], i, TCADpath_list.at(i), true));
 
-    TH2F* lorentzMap_e_hold    =new TH2F() ;
-    TH2F* lorentzMap_h_hold    =new TH2F() ;
-    TH2F* distanceMap_h_hold   =new TH2F() ;
-    TH2F* distanceMap_e_hold   =new TH2F() ;
-    TH1F* timeMap_e_hold       =new TH1F() ;
-    TH1F* timeMap_h_hold       =new TH1F() ;
+      TH2F* lorentzMap_e_hold  = new TH2F();
+      TH2F* lorentzMap_h_hold  = new TH2F();
+      TH2F* distanceMap_h_hold = new TH2F();
+      TH2F* distanceMap_e_hold = new TH2F();
+      TH1F* timeMap_e_hold     = new TH1F();
+      TH1F* timeMap_h_hold     = new TH1F();
 
-    if (doInterpolateEfield) {
       ATH_CHECK(m_radDamageUtil->generateDistanceTimeMap( distanceMap_e_hold, distanceMap_h_hold, timeMap_e_hold, timeMap_h_hold, lorentzMap_e_hold, lorentzMap_h_hold, eFieldMap_hold, NULL ));
       // For debugging and documentation: uncomment to save different maps which are based on the interpolated E field
       if(m_radDamageUtil->saveDebugMaps()){
@@ -304,27 +123,16 @@ StatusCode SensorSimPlanarTool::initialize() {
         prename.ReplaceAll( "_e","_h");
         lorentzMap_h_hold->SaveAs(prename);
       }
+      //Safetycheck
+      if (distanceMap_e_hold == 0 || distanceMap_h_hold == 0 || timeMap_e_hold == 0 || timeMap_h_hold == 0 || lorentzMap_e_hold == 0 || lorentzMap_h_hold == 0){
+        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);
     }
-    else {
-      //retrieve precomputed maps
-      lorentzMap_e_hold=(TH2F*)mapsFile->Get("lorentz_map_e");
-      lorentzMap_h_hold=(TH2F*)mapsFile->Get("lorentz_map_h");
-      distanceMap_h_hold=(TH2F*)mapsFile->Get("hdistance");
-      distanceMap_e_hold=(TH2F*)mapsFile->Get("edistance");
-      timeMap_e_hold=(TH1F*)mapsFile->Get("etimes");
-      timeMap_h_hold=(TH1F*)mapsFile->Get("htimes");
-    }
-    //Safetycheck
-    if (distanceMap_e_hold == 0 || distanceMap_h_hold == 0 || timeMap_e_hold == 0 || timeMap_h_hold == 0 || lorentzMap_e_hold == 0 || lorentzMap_h_hold == 0){
-      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_lorentzMap_e[Layer]=lorentzMap_e_hold;
-    m_lorentzMap_h[Layer]=lorentzMap_h_hold;
-    m_distanceMap_e[Layer]=distanceMap_e_hold;
-    m_distanceMap_h[Layer]=distanceMap_h_hold;
-    m_timeMap_e[Layer]=timeMap_e_hold;
-    m_timeMap_h[Layer]=timeMap_h_hold;
   }
   return StatusCode::SUCCESS;
 }
@@ -351,24 +159,16 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit> &phit, SiC
 	  }
   }
 
-  bool isBarrel = false;
-
   const PixelID* p_pixelId = static_cast<const PixelID *>(Module.getIdHelper());
   int layer=p_pixelId->layer_disk(Module.identify() );
-  int bec=1;
-  if (p_pixelId->is_barrel(Module.identify())) { bec=0; }
-
-  if (Module.isBarrel()) { isBarrel=true; }
 
-  std::pair<int, int> Layer;        // index for layer/end cap position
-  Layer.first=bec;                  //Barrel (0) or End Cap (1)   -    Maps only for Barrel at the moment. isBarrel will avoid zsh
-  Layer.second=layer;               //Layer: 0 = IBL Planar, 1=B-Layer, 2=Layer-1, 3=Layer-2
-  //IBL Barrel doesn't exist. So the possible idexes should be: 0-0, 0-1, 0-2, 0-3, 1-1, 1-2, 1-3
+  // retrieve conditions data
+  SG::ReadCondHandle<PixelModuleData> moduleData(m_moduleDataKey);
 
-  if (m_doRadDamage && isBarrel && m_fluence>0) {
-    std::pair<double,double> trappingTimes = m_radDamageUtil->getTrappingTimes( m_fluence_layersMaps[Layer] );
-    m_trappingTimeElectrons = trappingTimes.first;
-    m_trappingTimeHoles = trappingTimes.second;
+  std::pair<double,double> trappingTimes;
+  if (m_doRadDamage && Module.isBarrel()) {
+    if (m_doInterpolateEfield) { trappingTimes = m_radDamageUtil->getTrappingTimes(m_fluenceLayer[layer]); }
+    else                       { trappingTimes = m_radDamageUtil->getTrappingTimes(moduleData->getFluenceLayer(layer)); }
   }
 
   //Load values from energyDeposition
@@ -449,7 +249,7 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit> &phit, SiC
     int numBins_weightingPotential_y = 0;
     int numBins_weightingPotential_z = 0;
 
-    if (m_doRadDamage && m_fluence>0 && !(Module.isDBM()) && isBarrel) {
+    if (m_doRadDamage && !(Module.isDBM()) && Module.isBarrel()) {
       centreOfPixel_i = p_design.positionFromColumnRow(pixel_i.etaIndex(), pixel_i.phiIndex());
 
       //Find the displacment of the charge carriers from the centre of the pixel in +ve quadrant
@@ -464,12 +264,20 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit> &phit, SiC
       nnLoop_pixelPhiMin = std::max( -2, pixel_i.phiIndex() + 1 - phiCells );
 
       //Setup values to check for overflow when using maps
-      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();        
+      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
@@ -485,52 +293,66 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit> &phit, SiC
 
     for (int j=0; j<ncharges; j++) {
 
-      if (m_doRadDamage && m_fluence>0 && !(Module.isDBM()) && isBarrel) {
+      if (m_doRadDamage && !(Module.isDBM()) && Module.isBarrel()) {
         double u = CLHEP::RandFlat::shoot(0.,1.);
-        double drifttime_e = (-1.)*m_trappingTimeElectrons*TMath::Log(u); //ns
+        double drifttime_e = (-1.)*(trappingTimes.first)*TMath::Log(u); //ns
         u = CLHEP::RandFlat::shoot(0.,1.);
-        double drifttime_h = (-1.)*m_trappingTimeHoles*TMath::Log(u); //ns
+        double drifttime_h = (-1.)*(trappingTimes.second)*TMath::Log(u); //ns
 
         //Now, need the z-position at the trap.
-        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;
-        double depth_f_e = m_distanceMap_e[Layer]->GetBinContent( nbin_z_e_xbin,nbin_z_e_ybin );
-        double dz_e = fabs(dist_electrode - depth_f_e); 
-
+        double depth_f_e = 0.0;
+        double depth_f_h = 0.0;
+        double tanLorentz_e = 0.0;
+        double tanLorentz_h = 0.0;
         //TODO: the holes map does not currently extend for a drift time long enough that, any hole will reach
         //the corresponding electrode. This needs to be rectified by either (a) extrapolating the current map or
         //(b) making a new map with a y-axis (drift time) that extends to at least 18 ns so all charge carriers reach electrode.
         //However, if choose (b), will need to reduce granularity of map. 
-        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;
-        double depth_f_h = m_distanceMap_h[Layer]->GetBinContent( nbin_z_h_xbin,nbin_z_h_ybin );
+        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);       
+        }
+        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);       
+        }
+        double dz_e = fabs(dist_electrode - depth_f_e); 
         double dz_h = fabs(depth_f_h - dist_electrode);           
+        double coLorentz_e = std::sqrt(1.0+std::pow(tanLorentz_e,2));
 
         //Apply drift due to Lorentz force and diffusion
         double phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
 
-        int  nbin_Lorentz_e = m_lorentzMap_e[Layer]->FindBin(dist_electrode,depth_f_e);
-        tanLorentz = m_lorentzMap_e[Layer]->GetBinContent(nbin_Lorentz_e);       
-        coLorentz=sqrt(1+pow(tanLorentz,2));
-
         //Apply diffusion. rdif is teh max. diffusion
-        double rdif_e=this->m_diffusionConstant*sqrt( fabs(dist_electrode - depth_f_e)*coLorentz/0.3);
-        double phi_f_e=phi_i + dz_e*tanLorentz + rdif_e*phiRand;
+        double rdif_e=this->m_diffusionConstant*sqrt( fabs(dist_electrode - depth_f_e)*coLorentz_e/0.3);
+        double phi_f_e=phi_i + dz_e*tanLorentz_e + rdif_e*phiRand;
         double etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
         double eta_f_e=eta_i + rdif_e*etaRand;
 
         phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
-
-        int  nbin_Lorentz_h = m_lorentzMap_h[Layer]->FindBin(dist_electrode,depth_f_h);
-        tanLorentz = m_lorentzMap_h[Layer]->GetBinContent(nbin_Lorentz_h);       
-        coLorentz=sqrt(1+pow(tanLorentz,2));
-
-        double rdif_h=this->m_diffusionConstant*sqrt( fabs(dist_electrode - depth_f_h)*coLorentz/0.3);
-
-        double phi_f_h=phi_i + dz_h*tanLorentz + rdif_h*phiRand;
+        double coLorentz_h = std::sqrt(1.0+std::pow(tanLorentz_h,2));
+        double rdif_h=this->m_diffusionConstant*sqrt( fabs(dist_electrode - depth_f_h)*coLorentz_h/0.3);
+        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;
 
@@ -563,16 +385,6 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit> &phit, SiC
             double dEta_i_e = pixelEta_i - dEta_nn_centre;
             double dPhi_i_e = pixelPhi_i - dPhi_nn_centre;
 
-            int nbin_ramo_i_x = m_ramoPotentialMap[Layer]->GetXaxis()->FindBin( fabs( dPhi_i_e )*1000. );
-            int nbin_ramo_i_y = m_ramoPotentialMap[Layer]->GetYaxis()->FindBin( fabs( dEta_i_e )*1000. );
-            int nbin_ramo_i_z = m_ramoPotentialMap[Layer]->GetZaxis()->FindBin( dist_electrode*1000 );
-            //int nbin_ramo_i = m_ramoPotentialMap[0]->FindBin( fabs( dEta_i_e )*1000. , fabs( dPhi_i_e )*1000., dist_electrode*1000);
-
-            //Boundary check on maps
-            double ramo_i=0.;
-            if( nbin_ramo_i_x <= numBins_weightingPotential_x && nbin_ramo_i_y <= numBins_weightingPotential_y && nbin_ramo_i_z <=numBins_weightingPotential_z ){
-              ramo_i =m_ramoPotentialMap[Layer]->GetBinContent( nbin_ramo_i_x,nbin_ramo_i_y,nbin_ramo_i_z );
-            }     	
             //Find the displacment of the charge carriers from the centre of the pixel in +ve quadrant
             double pixelEta_f_e = eta_f_e - centreOfPixel_i.xEta() ;
             double pixelPhi_f_e = phi_f_e - centreOfPixel_i.xPhi() ;
@@ -581,45 +393,76 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit> &phit, SiC
             double pixelPhi_f_h = phi_f_h - centreOfPixel_i.xPhi() ;
 
             //Final position of charge carriers wrt nn centre
-            double dEta_f_e = pixelEta_f_e  - dEta_nn_centre ;
-            double dPhi_f_e = pixelPhi_f_e  - dPhi_nn_centre ;
-
-            int nbin_ramo_f_e_x = m_ramoPotentialMap[Layer]->GetXaxis()->FindBin( fabs( dPhi_f_e )*1000. );
-            int nbin_ramo_f_e_y = m_ramoPotentialMap[Layer]->GetYaxis()->FindBin( fabs( dEta_f_e )*1000. );
-            int nbin_ramo_f_e_z = m_ramoPotentialMap[Layer]->GetZaxis()->FindBin( depth_f_e*1000 );
-            //int nbin_ramo_f_e = m_ramoPotentialMap[0]->FindBin( fabs( dEta_f_e )*1000. , fabs( dPhi_f_e )*1000., depth_f_e*1000);
-            double ramo_f_e=0.;
-            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 ){
-              ramo_f_e =m_ramoPotentialMap[Layer]->GetBinContent( nbin_ramo_f_e_x,nbin_ramo_f_e_y,nbin_ramo_f_e_z );
-            }   
-
-            double dEta_f_h = pixelEta_f_h - dEta_nn_centre ;
-            double dPhi_f_h = pixelPhi_f_h - dPhi_nn_centre ;
-
-            int nbin_ramo_f_h_x = m_ramoPotentialMap[Layer]->GetXaxis()->FindBin( fabs( dPhi_f_h )*1000. );
-            int nbin_ramo_f_h_y = m_ramoPotentialMap[Layer]->GetYaxis()->FindBin( fabs( dEta_f_h )*1000. );
-            int nbin_ramo_f_h_z = m_ramoPotentialMap[Layer]->GetZaxis()->FindBin( depth_f_h*1000 );
-            //int nbin_ramo_f_h = m_ramoPotentialMap->FindBin( fabs( dEta_f_h )*1000. , fabs( dPhi_f_h )*1000., depth_f_h*1000);
+            double dEta_f_e = pixelEta_f_e - dEta_nn_centre;
+            double dPhi_f_e = pixelPhi_f_e - dPhi_nn_centre;
+
+            double dEta_f_h = pixelEta_f_h - dEta_nn_centre;
+            double dPhi_f_h = pixelPhi_f_h - dPhi_nn_centre;
+
             //Boundary check on maps
-            double ramo_f_h=0.;
-            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 ){
-              ramo_f_h =m_ramoPotentialMap[Layer]->GetBinContent( nbin_ramo_f_h_x,nbin_ramo_f_h_y,nbin_ramo_f_h_z );
-            } 
-
-            //Account for the imperfect binning that would cause charge to be double-counted
-            if(m_ramoPotentialMap[Layer]->GetZaxis()->FindBin(depth_f_h*1000) == m_ramoPotentialMap[Layer]->GetNbinsZ()+1) ramo_f_h=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.  || fabs(dPhi_f_e)>=Module.phiPitch()/2. ) ramo_f_e=0;
-              else if (fabs(dEta_f_e)<Module.etaPitch()/2.  && fabs(dPhi_f_e)<Module.phiPitch()/2.  ) ramo_f_e=1.;
+            double ramo_i   = 0.0;
+            double ramo_f_e = 0.0;
+            double ramo_f_h = 0.0;
+            if (m_doInterpolateEfield) {
+              int nbin_ramo_i_x = m_ramoPotentialMap[layer]->GetXaxis()->FindBin( fabs( dPhi_i_e )*1000. );
+              int nbin_ramo_i_y = m_ramoPotentialMap[layer]->GetYaxis()->FindBin( fabs( dEta_i_e )*1000. );
+              int nbin_ramo_i_z = m_ramoPotentialMap[layer]->GetZaxis()->FindBin( dist_electrode*1000 );
+              if (nbin_ramo_i_x<=numBins_weightingPotential_x && nbin_ramo_i_y<=numBins_weightingPotential_y && nbin_ramo_i_z<=numBins_weightingPotential_z) {
+                ramo_i = m_ramoPotentialMap[layer]->GetBinContent(nbin_ramo_i_x,nbin_ramo_i_y,nbin_ramo_i_z);
+              }
+              int nbin_ramo_f_e_x = m_ramoPotentialMap[layer]->GetXaxis()->FindBin(fabs(dPhi_f_e)*1000.0);
+              int nbin_ramo_f_e_y = m_ramoPotentialMap[layer]->GetYaxis()->FindBin(fabs(dEta_f_e)*1000.0);
+              int nbin_ramo_f_e_z = m_ramoPotentialMap[layer]->GetZaxis()->FindBin(depth_f_e*1000);
+              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) {
+                ramo_f_e = m_ramoPotentialMap[layer]->GetBinContent(nbin_ramo_f_e_x,nbin_ramo_f_e_y,nbin_ramo_f_e_z);
+              }
+              int nbin_ramo_f_h_x = m_ramoPotentialMap[layer]->GetXaxis()->FindBin( fabs( dPhi_f_h )*1000. );
+              int nbin_ramo_f_h_y = m_ramoPotentialMap[layer]->GetYaxis()->FindBin( fabs( dEta_f_h )*1000. );
+              int nbin_ramo_f_h_z = m_ramoPotentialMap[layer]->GetZaxis()->FindBin( depth_f_h*1000 );
+              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) {
+                ramo_f_h = m_ramoPotentialMap[layer]->GetBinContent(nbin_ramo_f_h_x,nbin_ramo_f_h_y,nbin_ramo_f_h_z);
+              }
+              //Account for the imperfect binning that would cause charge to be double-counted
+              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 { // use fluence value from conditions data
+              int nbin_ramo_i_x = moduleData->getRamoPotentialMap(layer)->GetXaxis()->FindBin( fabs( dPhi_i_e )*1000. );
+              int nbin_ramo_i_y = moduleData->getRamoPotentialMap(layer)->GetYaxis()->FindBin( fabs( dEta_i_e )*1000. );
+              int nbin_ramo_i_z = moduleData->getRamoPotentialMap(layer)->GetZaxis()->FindBin( dist_electrode*1000 );
+              if (nbin_ramo_i_x<=numBins_weightingPotential_x && nbin_ramo_i_y<=numBins_weightingPotential_y && nbin_ramo_i_z<=numBins_weightingPotential_z) {
+                ramo_i = moduleData->getRamoPotentialMap(layer)->GetBinContent(nbin_ramo_i_x,nbin_ramo_i_y,nbin_ramo_i_z);
+              }
+              int nbin_ramo_f_e_x = moduleData->getRamoPotentialMap(layer)->GetXaxis()->FindBin( fabs( dPhi_f_e )*1000. );
+              int nbin_ramo_f_e_y = moduleData->getRamoPotentialMap(layer)->GetYaxis()->FindBin( fabs( dEta_f_e )*1000. );
+              int nbin_ramo_f_e_z = moduleData->getRamoPotentialMap(layer)->GetZaxis()->FindBin( depth_f_e*1000 );
+              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) {
+                ramo_f_e = moduleData->getRamoPotentialMap(layer)->GetBinContent(nbin_ramo_f_e_x,nbin_ramo_f_e_y,nbin_ramo_f_e_z);
+              }
+              int nbin_ramo_f_h_x = moduleData->getRamoPotentialMap(layer)->GetXaxis()->FindBin( fabs( dPhi_f_h )*1000. );
+              int nbin_ramo_f_h_y = moduleData->getRamoPotentialMap(layer)->GetYaxis()->FindBin( fabs( dEta_f_h )*1000. );
+              int nbin_ramo_f_h_z = moduleData->getRamoPotentialMap(layer)->GetZaxis()->FindBin( depth_f_h*1000 );
+              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) {
+                ramo_f_h = moduleData->getRamoPotentialMap(layer)->GetBinContent(nbin_ramo_f_h_x,nbin_ramo_f_h_y,nbin_ramo_f_h_z);
+              }
+              //Account for the imperfect binning that would cause charge to be double-counted
+              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; }
+              }
             }
-
 
             //Given final position of charge carrier, find induced charge. The difference in Ramo weighting potential gives the fraction of charge induced.
             //The energy_per_step is transformed into charge with the eleholePair per Energy
             double induced_charge_e = (ramo_f_e - ramo_i) * energy_per_step * eleholePairEnergy;
             double induced_charge_h = -(ramo_f_h - ramo_i) * energy_per_step * eleholePairEnergy;
 
-
             //Collect charge in centre of each pixel, since location within pixel doesn't matter for record
             SiLocalPosition chargePos = Module.hitLocalToLocal( centreOfPixel_nn.xEta(), centreOfPixel_nn.xPhi() );
 
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h
index 3db9f5b5e22d13bca13f2be9bc14d169277d1857..f3f0fdda27e461b68f594385c5efe1ffaecf4e17 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.h
@@ -35,17 +35,11 @@ class SensorSimPlanarTool : public SensorSimTool {
     SensorSimPlanarTool();
 
     // Map for radiation damage simulation
-    std::map<std::pair<int, int>, TH3F*> m_ramoPotentialMap;
-    std::map<std::pair<int, int>, TH1F*> m_eFieldMap;
-    std::map<std::pair<int, int>, TH2F*> m_distanceMap_e;
-    std::map<std::pair<int, int>, TH2F*> m_distanceMap_h;
-    std::map<std::pair<int, int>, TH1F*> m_timeMap_e;
-    std::map<std::pair<int, int>, TH1F*> m_timeMap_h;
-    std::map<std::pair<int, int>, TH2F*> m_lorentzMap_e;
-    std::map<std::pair<int, int>, TH2F*> m_lorentzMap_h;
-
-    std::vector<double> m_fluence_layers,m_voltage_layers; //merging information from m_fluence* and m_voltage*
-    std::map<std::pair<int, int>, double> m_fluence_layersMaps;
+    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;
 
     Gaudi::Property<int> m_numberOfSteps
     {this, "numberOfSteps", 50, "Geant4:number of steps for PixelPlanar"};
@@ -56,35 +50,21 @@ class SensorSimPlanarTool : public SensorSimTool {
     Gaudi::Property<bool> m_doRadDamage
     {this, "doRadDamage", false, "doRadDmaage bool: should be flag"};
 
-    Gaudi::Property<double> m_trappingTimeElectrons
-    {this, "trappingTimeElectrons", 0.0, "Characteristic time till electron is trapped [ns]"};
+    Gaudi::Property<bool> m_doInterpolateEfield
+    {this, "doInterpolateEfield", false, "doInterpolateEfield bool: should be flag"};
 
-    Gaudi::Property<double> m_trappingTimeHoles
-    {this, "trappingTimeHoles", 0.0, "Characteristic time till hole is trapped [ns]"};
+    Gaudi::Property<std::vector<std::string>> m_fluenceMap
+    {this, "FluenceMap", {"PixelDigitization/maps_IBL_PL_400V_fl5_5e14.root",
+                          "PixelDigitization/maps_PIX_400V_fl5_19e14.root",
+                          "PixelDigitization/maps_PIX_250V_fl2_28e14.root",
+                          "PixelDigitization/maps_PIX_250V_fl1_53e14.root"},
+                          "Fluence map for radiation damage when interpolation method is activated"};
 
-    Gaudi::Property<double> m_fluence
-    {this, "fluence", 0, "this is the fluence benchmark, 0-6. 0 is unirradiated, 1 is start of Run 2, 5 is end of 2018 and 6 is projected end of 2018"};
+    Gaudi::Property<std::vector<double>> m_fluenceLayer
+    {this, "FluenceLayer", {5.50e14, 5.19e14, 2.28e14, 1.53e14}, "Fluence for radiation damage when interpolation method is activated"};
 
-    Gaudi::Property<double> m_fluenceB
-    {this, "fluenceB", -1, "fluence detector has recieved in neqcm2 at the B layer."};
-
-    Gaudi::Property<double> m_fluence1
-    {this, "fluence1", -1, "fluence detector has recieved in neqcm2 at the layer 1."};
-
-    Gaudi::Property<double> m_fluence2
-    {this, "fluence2", -1, "fluence detector has recieved in neqcm2 at the layer 2."};
-
-    Gaudi::Property<double> m_voltage
-    {this, "voltage", -1.0, "this is the bias voltage applied to the IBL - if not set use values from applied at benchmark points according to fluence"};
-
-    Gaudi::Property<double> m_voltageB
-    {this, "voltageB", -1.0, "bias voltage applied to the B layer."};
-
-    Gaudi::Property<double> m_voltage1
-    {this, "voltage1", -1.0, "bias voltage applied to the layer 1."};
-
-    Gaudi::Property<double> m_voltage2
-    {this, "voltage2", -1.0, "bias voltage applied to the layer 2."};
+    Gaudi::Property<std::vector<float>> m_voltageLayer
+    {this, "BiasVoltageLayer", {400.0,400.0,250.0,250.0}, "Bias voltage for radiation damage when interpolation method is activated"};
 
     ToolHandle<RadDamageUtil> m_radDamageUtil
     {this, "RadDamageUtil", "RadDamageUtil", "Rad Damage utility"};
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimTool.h
index fbf423714cb48e5b3bec9ea9c2407e923b8fba17..adb5b1a68d522f00a43125a20f7633f560b32876 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimTool.h
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimTool.h
@@ -23,6 +23,9 @@
 #include "PixelReadoutGeometry/PixelModuleDesign.h"
 #include "SiPropertiesTool/ISiPropertiesTool.h"
 
+#include "PixelConditionsData/PixelModuleData.h"
+#include "StoreGate/ReadCondHandleKey.h"
+
 static const InterfaceID IID_ISensorSimTool("SensorSimTool", 1, 0);
 
 class SensorSimTool:public AthAlgTool,virtual public IAlgTool {
@@ -39,6 +42,7 @@ class SensorSimTool:public AthAlgTool,virtual public IAlgTool {
     virtual StatusCode initialize() {
       ATH_CHECK(AthAlgTool::initialize()); 
       ATH_CHECK(m_siPropertiesTool.retrieve());
+      ATH_CHECK(m_moduleDataKey.initialize());
       return StatusCode::SUCCESS;
     }
 
@@ -52,7 +56,12 @@ class SensorSimTool:public AthAlgTool,virtual public IAlgTool {
     SensorSimTool();
 
   protected:
-    ToolHandle<ISiPropertiesTool>   m_siPropertiesTool{this, "SiPropertiesTool", "SiPropertiesTool", "Tool to retrieve SiProperties"};
+    ToolHandle<ISiPropertiesTool> m_siPropertiesTool
+    {this, "SiPropertiesTool", "SiPropertiesTool", "Tool to retrieve SiProperties"};
+
+    SG::ReadCondHandleKey<PixelModuleData> m_moduleDataKey
+    {this, "PixelModuleData", "PixelModuleData", "Pixel module data"};
+
 };
 
 #endif // PIXELDIGITIZATION_SensorSimTool_H