diff --git a/InnerDetector/InDetDigitization/PixelDigitization/python/ITkPixelDigitizationConfig.py b/InnerDetector/InDetDigitization/PixelDigitization/python/ITkPixelDigitizationConfig.py
index 47c6fd9d5647e42b2d6693638438d0663e6aded2..13ecf61d382b0639f4d94ae7d6dea663a0db1070 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/python/ITkPixelDigitizationConfig.py
+++ b/InnerDetector/InDetDigitization/PixelDigitization/python/ITkPixelDigitizationConfig.py
@@ -91,6 +91,7 @@ def ITkSensorSimPlanarToolCfg(flags, name="ITkSensorSimPlanarTool", **kwargs):
     kwargs.setdefault("PixelModuleData", "ITkPixelModuleData")
     SensorSimPlanarTool = CompFactory.SensorSimPlanarTool
     kwargs.setdefault("doRadDamage", flags.Digitization.DoPixelPlanarRadiationDamage)
+    kwargs.setdefault("doRadDamageTemplate", flags.Digitization.DoPixelPlanarRadiationDamageTemplate)
     if flags.Digitization.DoPixelPlanarRadiationDamage:
         # acc.merge(ITkPixelRadSimFluenceMapAlgCfg(flags))  # TODO: not supported yet
         pass
@@ -107,6 +108,7 @@ def ITkSensorSim3DToolCfg(flags, name="ITkSensorSim3DTool", **kwargs):
     kwargs.setdefault("PixelModuleData", "ITkPixelModuleData")
     SensorSim3DTool = CompFactory.SensorSim3DTool
     kwargs.setdefault("doRadDamage", flags.Digitization.DoPixel3DRadiationDamage)
+    kwargs.setdefault("doRadDamageTemplate", flags.Digitization.DoPixel3DRadiationDamageTemplate)
     if flags.Digitization.DoPixel3DRadiationDamage:
         # acc.merge(ITkPixelRadSimFluenceMapAlgCfg(flags))  # TODO: not supported yet
         pass
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py
index a201f33998dccd9306346b352433563e7847082a..46377a4513b7fa800da76ff14d6fb853aab1c878 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py
+++ b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfig.py
@@ -54,12 +54,14 @@ def SensorSimPlanarTool(name="SensorSimPlanarTool", **kwargs):
     kwargs.setdefault("SiPropertiesTool", ToolSvc.PixelSiPropertiesTool)
     kwargs.setdefault("LorentzAngleTool", ToolSvc.PixelLorentzAngleTool)
     kwargs.setdefault("doRadDamage", digitizationFlags.doPixelPlanarRadiationDamage.get_Value())
+    kwargs.setdefault("doRadDamageTemplate", digitizationFlags.doPixelPlanarRadiationDamageTemplate.get_Value())
     return CfgMgr.SensorSimPlanarTool(name, **kwargs)
 
 def SensorSim3DTool(name="SensorSim3DTool", **kwargs):
     from AthenaCommon.AppMgr import ToolSvc
     kwargs.setdefault("SiPropertiesTool", ToolSvc.PixelSiPropertiesTool)
     kwargs.setdefault("doRadDamage", digitizationFlags.doPixel3DRadiationDamage.get_Value())
+    kwargs.setdefault("doRadDamageTemplate", digitizationFlags.doPixel3DRadiationDamageTemplate.get_Value())
     return CfgMgr.SensorSim3DTool(name, **kwargs)
 
 def SensorSimTool(name="SensorSimTool", **kwargs):
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py
index 52044ce01d18b4e3db83a1fffff59f0144f3d424..11e09cd59128476addd00490fe69148d4f27da7c 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py
+++ b/InnerDetector/InDetDigitization/PixelDigitization/python/PixelDigitizationConfigNew.py
@@ -83,6 +83,7 @@ def SensorSimPlanarToolCfg(flags, name="SensorSimPlanarTool", **kwargs):
     kwargs.setdefault("LorentzAngleTool", LorentzTool)
     SensorSimPlanarTool = CompFactory.SensorSimPlanarTool
     kwargs.setdefault("doRadDamage", flags.Digitization.DoPixelPlanarRadiationDamage)
+    kwargs.setdefault("doRadDamageTemplate", flags.Digitization.DoPixelPlanarRadiationDamageTemplate)
     if flags.Digitization.DoPixelPlanarRadiationDamage:
         acc.merge(PixelRadSimFluenceMapAlgCfg(flags))
     acc.setPrivateTools(SensorSimPlanarTool(name, **kwargs))
@@ -97,6 +98,7 @@ def SensorSim3DToolCfg(flags, name="SensorSim3DTool", **kwargs):
     kwargs.setdefault("SiPropertiesTool", SiTool)
     SensorSim3DTool = CompFactory.SensorSim3DTool
     kwargs.setdefault("doRadDamage", flags.Digitization.DoPixel3DRadiationDamage)
+    kwargs.setdefault("doRadDamageTemplate", flags.Digitization.DoPixel3DRadiationDamageTemplate)
     if flags.Digitization.DoPixel3DRadiationDamage:
         acc.merge(PixelRadSimFluenceMapAlgCfg(flags))
     acc.setPrivateTools(SensorSim3DTool(name, **kwargs))
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
index acfccaa5c22be957faa7020dc3ce8238cd7fdff3..9a2a412feeb1f27e9883bfb7548e92e835ce80cb 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSim3DTool.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
 */
 
 /**
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
index 2cf68c9ef00d9b5417441bd8701e10030e7e7a9a..7b996bc87de8f9a9150edcd2b85cab3b58805034 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimPlanarTool.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "SensorSimPlanarTool.h"
@@ -7,6 +7,7 @@
 #include "PixelReadoutGeometry/PixelModuleDesign.h"
 #include "SiDigitization/SiSurfaceCharge.h"
 #include "InDetSimEvent/SiHit.h"
+#include "InDetSimEvent/SiTrackDistance.h"
 #include "InDetIdentifier/PixelID.h"
 #include "GeneratorObjects/HepMcParticleLink.h"
 #include "SiPropertiesTool/SiliconProperties.h"
@@ -231,6 +232,13 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
   double sensorThickness = Module.design().thickness();
   const InDet::SiliconProperties& siProperties = m_siPropertiesTool->getSiProperties(Module.identifyHash(), ctx);
 
+  // Prepare values that are needed for radiation damage corrections later
+  const HepGeom::Point3D<double>& startPosition = phit->localStartPosition(); 
+  const HepGeom::Point3D<double>& endPosition   = phit->localEndPosition(); 
+  const float zPosDiff = endPosition.Z - startPosition.Z;
+  const float tanTheta = (zPosDiff != 0.) ? (endPosition.Y - startPosition.Y) / zPosDiff : (endPosition.Y - startPosition.Y) * 9999.;
+  const float tanPhi   = (zPosDiff != 0.) ? (endPosition.X - startPosition.X) / zPosDiff : (endPosition.X - startPosition.X) * 9999.;
+
   int etaCells = p_design.columns();
   int phiCells = p_design.rows();
 
@@ -265,7 +273,7 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
   const double halfEtaPitch = 0.5*Module.etaPitch();
   const double halfPhiPitch = 0.5*Module.phiPitch();
 
-  if (m_doRadDamage && !(Module.isDBM()) && Module.isBarrel()) {
+  if (m_doRadDamage && !(Module.isDBM()) && Module.isBarrel() && !m_doRadDamageTemplate) {
     SG::ReadCondHandle<PixelRadiationDamageFluenceMapData> fluenceData(m_fluenceDataKey,ctx);
 
     std::pair<double, double> trappingTimes;
@@ -494,7 +502,7 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
     });
 
   } 
-  else { //If no radDamage, run original
+  else if (m_doRadDamage && m_doRadDamageTemplate && !(Module.isDBM()) && Module.isBarrel()){ // will run radiation damage but with the template method
     for (size_t i = 0; i < trfHitRecord.size(); i++) {
       std::pair<double, double> const& iHitRecord = trfHitRecord[i];
 
@@ -548,8 +556,10 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
           ed = energy_per_step * eleholePairEnergy;
         }
 
+        // prepare SiTrackDistance object needed
+        const SiTrackDistance trackDistance(tanTheta, tanPhi, startPosition.Z);
         //The following lines are adapted from SiDigitization's Inserter class
-        const SiSurfaceCharge scharge(chargePos, SiCharge(ed, pHitTime, SiCharge::track, particleLink));
+        const SiSurfaceCharge scharge(chargePos, SiCharge(ed, pHitTime, SiCharge::track, particleLink, trackDistance));
 
         const SiCellId& diode = Module.cellIdOfPosition(scharge.position());
 
@@ -559,7 +569,72 @@ StatusCode SensorSimPlanarTool::induceCharge(const TimedHitPtr<SiHit>& phit,
         }
       }//end cycle for charge
     }//trfHitRecord.size()
-  } //else: no radDamage, run original
+  } else { // run without radiation damage
+    for (size_t i = 0; i < trfHitRecord.size(); i++) {
+      std::pair<double, double> const& iHitRecord = trfHitRecord[i];
+
+      double eta_i = eta_0;
+      double phi_i = phi_0;
+      double depth_i = depth_0;
+      if (iTotalLength) {
+        eta_i += 1.0 * iHitRecord.first / iTotalLength * dEta;
+        phi_i += 1.0 * iHitRecord.first / iTotalLength * dPhi;
+        depth_i += 1.0 * iHitRecord.first / iTotalLength * dDepth;
+      }
+
+      // Distance between charge and readout side.  p_design->readoutSide() is
+      // +1 if readout side is in +ve depth axis direction and visa-versa.
+      double dist_electrode = 0.5 * sensorThickness - Module.design().readoutSide() * depth_i;
+      if (dist_electrode < 0) {
+        dist_electrode = 0;
+      }
+
+      // nonTrapping probability
+      double nontrappingProbability = 1.0;
+      if (Module.isDBM()) {
+        nontrappingProbability = exp(-dist_electrode / collectionDist);
+      }
+
+      for (int j = 0; j < ncharges; j++) {
+        // amount of energy to be converted into charges at current step
+        double energy_per_step = 1.0 * iHitRecord.second / 1.E+6 / ncharges;
+        // diffusion sigma
+        double rdif = this->m_diffusionConstant * std::sqrt(dist_electrode * coLorentz / 0.3);
+
+        // position at the surface
+        double phiRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+        double phi_drifted = phi_i + dist_electrode * tanLorentz + rdif * phiRand;
+        double etaRand = CLHEP::RandGaussZiggurat::shoot(rndmEngine);
+        double eta_drifted = eta_i + rdif * etaRand;
+
+        // Slim Edge for IBL planar sensors:
+        if (!(Module.isDBM()) && p_design.getReadoutTechnology() == InDetDD::PixelReadoutTechnology::FEI4) {
+          applyIBLSlimEdges(energy_per_step, eta_drifted);
+        }
+
+        // Get the charge position in Reconstruction local coordinates.
+        const SiLocalPosition& chargePos = Module.hitLocalToLocal(eta_drifted, phi_drifted);
+
+        // The parametrization of the sensor efficiency (if needed)
+        double ed = 0;
+        if (Module.isDBM()) {
+          ed = energy_per_step * eleholePairEnergy * nontrappingProbability * smearScale;
+        } else {
+          ed = energy_per_step * eleholePairEnergy;
+        }
+
+        //The following lines are adapted from SiDigitization's Inserter class
+        const SiSurfaceCharge scharge(chargePos, SiCharge(ed, pHitTime, SiCharge::track, particleLink));
+
+        const SiCellId& diode = Module.cellIdOfPosition(scharge.position());
+
+        if (diode.isValid()) {
+          const SiCharge& charge = scharge.charge();
+          chargedDiodes.add(diode, charge);
+        }
+      }//end cycle for charge
+    }
+  }
 
   return StatusCode::SUCCESS;
 }
diff --git a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimTool.h b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimTool.h
index 8496ca329c6668e017cac87ca64257a82ff69677..232801739b29de281b1eeaf7f7140950b9950c71 100644
--- a/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimTool.h
+++ b/InnerDetector/InDetDigitization/PixelDigitization/src/SensorSimTool.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
 */
 
 /**
@@ -70,7 +70,12 @@ protected:
 
   Gaudi::Property<bool> m_doRadDamage
   {
-    this, "doRadDamage", false, "doRadDmaage bool: should be flag"
+    this, "doRadDamage", false, "doRadDamage bool: should be flag"
+  };
+
+  Gaudi::Property<bool> m_doRadDamageTemplate
+  {
+    this, "doRadDamageTemplate", false, "doRadDamageTemplate bool: should be flag - tells the code to use the template version of the radiation damage"
   };
 
   SG::ReadCondHandleKey<PixelRadiationDamageFluenceMapData> m_fluenceDataKey
diff --git a/InnerDetector/InDetDigitization/SiDigitization/SiDigitization/SiChargedDiode.h b/InnerDetector/InDetDigitization/SiDigitization/SiDigitization/SiChargedDiode.h
index 41efa7d1e07d3294075561e9cb620f62984c3111..55ff0f0101855f5b9dff31ce88a23af8b1d94bfa 100755
--- a/InnerDetector/InDetDigitization/SiDigitization/SiDigitization/SiChargedDiode.h
+++ b/InnerDetector/InDetDigitization/SiDigitization/SiDigitization/SiChargedDiode.h
@@ -21,6 +21,7 @@
 
 // Data member classes
 #include "InDetSimEvent/SiTotalCharge.h"
+#include "InDetSimEvent/SiTrackDistance.h"
 #include "ReadoutGeometryBase/SiCellId.h"
 #include "ReadoutGeometryBase/SiReadoutCellId.h"
 
@@ -89,6 +90,8 @@ private:
   InDetDD::SiReadoutCellId m_readoutCell; //Readout cell associated to this diode
   int m_word;   // a flag for noise etc etc as in InDetSimData
   SiChargedDiode * m_nextInCluster; //the next strip to navigate to - allows traversing clusters since the SiChargedDiodeCollection is not guaranteed to be contiguous
+  float m_maxCharge{}; // maximum deposited charge
+  SiTrackDistance m_trackDistance{}; // track parameters for the maximum deposited charge
 };
 
 ///////////////////////////////////////////////////////////////////
@@ -122,11 +125,6 @@ inline SiChargedDiode * SiChargedDiode::nextInCluster()
   return m_nextInCluster;
 }
 
-inline void SiChargedDiode::add(const SiCharge &charge)
-{
-  m_totalCharge.add(charge);
-}
-
 inline void SiChargedDiode::add(const SiTotalCharge &totcharge)
 {
   m_totalCharge.add(totcharge);
diff --git a/InnerDetector/InDetDigitization/SiDigitization/src/SiChargedDiode.cxx b/InnerDetector/InDetDigitization/SiDigitization/src/SiChargedDiode.cxx
index 4b7b60e5e02b633a5f467a31abd7aef07b509b87..3e87bd3a79ddd884aa3a49a842889dd90febd0a2 100755
--- a/InnerDetector/InDetDigitization/SiDigitization/src/SiChargedDiode.cxx
+++ b/InnerDetector/InDetDigitization/SiDigitization/src/SiChargedDiode.cxx
@@ -23,6 +23,20 @@ SiChargedDiode::SiChargedDiode(const InDetDD::SiCellId & diode, const InDetDD::S
      m_nextInCluster(nextInCluster)
 {}
 
+// add another charge
+// If the current charge is large than max, update
+// If we update, replace the track parameters
+void SiChargedDiode::add(const SiCharge &charge) {
+  const float currentCharge = std::abs(charge.charge());
+  if (currentCharge > m_maxCharge) {
+    m_maxCharge = currentCharge;
+    // update the track distance parameters
+    m_trackDistance = charge.trackDistance();
+  }
+
+  // add the charge
+  m_totalCharge.add(charge);
+}
 
 ///////////////////////////////////////////////////////////////////
 // Input/Output stream functions:
diff --git a/InnerDetector/InDetSimEvent/InDetSimEvent/SiCharge.h b/InnerDetector/InDetSimEvent/InDetSimEvent/SiCharge.h
index 81365d2fc68908e1b002e4f1dbcb0febe080c68e..454333cb4cec8f8d1cc98ca62f07982030e63f1d 100755
--- a/InnerDetector/InDetSimEvent/InDetSimEvent/SiCharge.h
+++ b/InnerDetector/InDetSimEvent/InDetSimEvent/SiCharge.h
@@ -16,11 +16,14 @@
 #ifndef SITRACKEREVENT_SICHARGE_H
 #define SITRACKEREVENT_SICHARGE_H
 
+#include "InDetSimEvent/SiTrackDistance.h"
+
 #include <iostream>
 
 // Member classes
 #include "GeneratorObjects/HepMcParticleLink.h"
 
+
 class SiCharge {
 
 public:
@@ -45,6 +48,12 @@ public:
   SiCharge(const double& charge,const double& time,
 	   const Process& processType);
 
+  SiCharge(const double& charge,
+           const double& time,
+           const Process& processType,
+           const HepMcParticleLink& PL,
+           const SiTrackDistance& trackDistance);
+
   // Destructor:
   ~SiCharge();
 
@@ -70,6 +79,9 @@ public:
   // Particle Link of the particle generating the charge
   const HepMcParticleLink& particleLink() const;
 
+  // Get the SiTrackDistance parameters
+  const SiTrackDistance& trackDistance() const;
+
   ///////////////////////////////////////////////////////////////////
   // Non-const methods:
   ///////////////////////////////////////////////////////////////////
@@ -94,6 +106,7 @@ private:
   Process m_processType; // type of process which produced this charge
   //  int m_trackNumber; // track number in case of track process
   HepMcParticleLink m_partLink; //Replace the track number with a PL
+  SiTrackDistance m_trackDistance;
 };
 
 ///////////////////////////////////////////////////////////////////
@@ -127,6 +140,11 @@ inline const HepMcParticleLink& SiCharge::particleLink() const
   return m_partLink;
 }
 
+inline const SiTrackDistance& SiCharge::trackDistance() const
+{
+  return m_trackDistance;
+}
+
 ///////////////////////////////////////////////////////////////////
 // Input/Output stream functions:
 ///////////////////////////////////////////////////////////////////
diff --git a/InnerDetector/InDetSimEvent/InDetSimEvent/SiTrackDistance.h b/InnerDetector/InDetSimEvent/InDetSimEvent/SiTrackDistance.h
new file mode 100644
index 0000000000000000000000000000000000000000..c4a23e9ec8ac3d33167e57a19e1f20558c1f72de
--- /dev/null
+++ b/InnerDetector/InDetSimEvent/InDetSimEvent/SiTrackDistance.h
@@ -0,0 +1,90 @@
+/*
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef SITRACKEREVENT_SITRACKDISTANCE_H
+#define SITRACKEREVENT_SITRACKDISTANCE_H
+
+/**
+ * @brief Simple class that holds informattion about the angles (phi and theta)
+ * and distance (in x and y) to a truth track wrt a diode (pixel)
+ * 
+ */
+
+class SiTrackDistance {
+
+public:
+
+    /**
+     * @brief Construct a new Si Track Distance object with the default parameters
+     * 
+     */
+    SiTrackDistance();
+
+    /**
+     * @brief Construct a new Si Track Distance object
+     * 
+     * @param tanTheta tangent of the theta angle of the track 
+     * @param tanPhi tangent of the phi angle of the track
+     * @param zInit z position of the charge entrance
+     */
+    SiTrackDistance(float tanTheta, float tanPhi, float zInit);
+
+    /**
+     * @brief Destroy the Si Track Distance object
+     * 
+     */
+    ~SiTrackDistance() = default;
+
+    /**
+     * @brief Default copy/move/assignment
+     * 
+     */
+    SiTrackDistance(const SiTrackDistance&) = default; 
+    SiTrackDistance& operator=(const SiTrackDistance&) = default; 
+    SiTrackDistance(SiTrackDistance&&) = default; 
+    SiTrackDistance& operator=(SiTrackDistance&&) = default; 
+
+    /**
+     * @brief Get the tangent of the Theta angle of the track
+     * 
+     * @return float 
+     */
+    inline float tanTheta() const {return m_tanTheta;}
+    
+    /**
+     * @brief Get the tangent Phi angle of the track
+     * 
+     * @return float 
+     */
+    inline float tanPhi() const {return m_tanPhi;}
+    
+    /**
+     * @brief Get the initial Z position of the track
+     * 
+     * @return float 
+     */
+    inline float zInit() const {return m_zInit;}
+    
+private:
+
+    /**
+     * @brief tangent of the angle theta of the track
+     * 
+     */
+    float m_tanTheta{};
+
+    /**
+     * @brief tangent of the angle phi of the track
+     * 
+     */
+    float m_tanPhi{};
+
+    /**
+     * @brief Initial Z position to the track
+     * 
+     */
+    float m_zInit{};
+};
+
+#endif
diff --git a/InnerDetector/InDetSimEvent/src/SiCharge.cxx b/InnerDetector/InDetSimEvent/src/SiCharge.cxx
index 293445c94fcadc254bdb846654bf855475cdb3b6..4740367bf8c1404b0ea083c02509a69b15ef7dac 100755
--- a/InnerDetector/InDetSimEvent/src/SiCharge.cxx
+++ b/InnerDetector/InDetSimEvent/src/SiCharge.cxx
@@ -18,7 +18,8 @@ SiCharge::SiCharge(const SiCharge &charge) :
   m_charge(charge.m_charge),
   m_time(charge.m_time),
   m_processType(charge.m_processType),
-  m_partLink(charge.m_partLink)
+  m_partLink(charge.m_partLink),
+  m_trackDistance(charge.m_trackDistance)
 {}
 
 // Constructor with parameters:
@@ -27,7 +28,8 @@ SiCharge::SiCharge(const double& charge,const double& time,
   m_charge(charge),
   m_time(time),
   m_processType(processType),
-  m_partLink(PL)
+  m_partLink(PL),
+  m_trackDistance(SiTrackDistance())
 {}
 
 SiCharge::SiCharge(const double& charge,const double& time,
@@ -35,7 +37,21 @@ SiCharge::SiCharge(const double& charge,const double& time,
   m_charge(charge),
   m_time(time),
   m_processType(processType),
-  m_partLink(HepMcParticleLink())
+  m_partLink(HepMcParticleLink()),
+  m_trackDistance(SiTrackDistance())
+{}
+
+// Constructor with parameters:
+SiCharge::SiCharge(const double& charge,
+                   const double& time,
+                   const Process& processType,
+                   const HepMcParticleLink& PL,
+                   const SiTrackDistance& trackDistance) :
+  m_charge(charge),
+  m_time(time),
+  m_processType(processType),
+  m_partLink(PL),
+  m_trackDistance(trackDistance)
 {}
 
 
@@ -47,6 +63,7 @@ SiCharge &SiCharge::operator=(const SiCharge &charge)
     m_time=charge.m_time;
     m_processType=charge.m_processType;
     m_partLink=charge.m_partLink;
+    m_trackDistance=charge.m_trackDistance;
   } else {}
   return *this;
 }
diff --git a/InnerDetector/InDetSimEvent/src/SiTrackDistance.cxx b/InnerDetector/InDetSimEvent/src/SiTrackDistance.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..718cd6b138618a5dd8a6b3da8647b3ee3755d032
--- /dev/null
+++ b/InnerDetector/InDetSimEvent/src/SiTrackDistance.cxx
@@ -0,0 +1,16 @@
+/*
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "InDetSimEvent/SiTrackDistance.h"
+
+SiTrackDistance::SiTrackDistance()
+{   
+}
+
+SiTrackDistance::SiTrackDistance(float tanTheta, float tanPhi, float zInit) :
+    m_tanTheta(tanTheta),
+    m_tanPhi(tanPhi),
+    m_zInit(zInit)
+{   
+}
diff --git a/Simulation/Digitization/python/DigitizationConfigFlags.py b/Simulation/Digitization/python/DigitizationConfigFlags.py
index 6edc6d249e0c0aac8389bec35588599417d60ed1..716f69c056f32c695d9fa9a610ae21751296b479 100644
--- a/Simulation/Digitization/python/DigitizationConfigFlags.py
+++ b/Simulation/Digitization/python/DigitizationConfigFlags.py
@@ -73,8 +73,12 @@ def createDigitizationCfgFlags():
 
     # Run radiation damage simulation for pixel planar sensors
     flags.addFlag("Digitization.DoPixelPlanarRadiationDamage", False)
+    # Run the template version of the radiation damage
+    flags.addFlag("Digitization.DoPixelPlanarRadiationDamageTemplate", False)
     # Run radiation damage simulation for pixel 3D sensors
     flags.addFlag("Digitization.DoPixel3DRadiationDamage", False)
+    # Run the template version of the radiation damage
+    flags.addFlag("Digitization.DoPixel3DRadiationDamageTemplate", False)
 
     # for PileUp digitization
     # Bunch structure configuration
diff --git a/Simulation/Digitization/python/DigitizationFlags.py b/Simulation/Digitization/python/DigitizationFlags.py
index d13401b814402a522430a7ebfcd5fe9452e76466..d02dec24d8d58a2ec3b8696a21ca5497b25b20fb 100755
--- a/Simulation/Digitization/python/DigitizationFlags.py
+++ b/Simulation/Digitization/python/DigitizationFlags.py
@@ -179,6 +179,14 @@ class doPixelPlanarRadiationDamage(JobProperty):
     allowedTypes=['bool']
     StoredValue=False
 
+#
+class doPixelPlanarRadiationDamageTemplate(JobProperty):
+    """ Include radiation damage simulation using the templates for pixel planar sensors where possible?
+    """
+    statusOn=True
+    allowedTypes=['bool']
+    StoredValue=False
+
 #
 class doPixel3DRadiationDamage(JobProperty):
     """ Include radiation damage simulation for pixel planar sensors where possible?
@@ -187,6 +195,14 @@ class doPixel3DRadiationDamage(JobProperty):
     allowedTypes=['bool']
     StoredValue=False
 
+#
+class doPixel3DRadiationDamageTemplate(JobProperty):
+    """ Include radiation damage simulation using the templates for pixel planar sensors where possible?
+    """
+    statusOn=True
+    allowedTypes=['bool']
+    StoredValue=False
+
 #
 class overrideMetadata(JobProperty):
     """ If digi config differs from that stored in sim metadata use digi values.
@@ -838,7 +854,8 @@ jobproperties.add_Container(Digitization)
 
 
 # We want always the following flags in the container
-list_jobproperties=[doInDetNoise,doCaloNoise,doMuonNoise,doFwdNoise,doPixelPlanarRadiationDamage,doPixel3DRadiationDamage,\
+list_jobproperties=[doInDetNoise,doCaloNoise,doMuonNoise,doFwdNoise,\
+                    doPixelPlanarRadiationDamage,doPixel3DRadiationDamage,doPixelPlanarRadiationDamageTemplate,doPixel3DRadiationDamageTemplate,\
                     rndmSvc,rndmSeedList,rndmSeedOffset1,rndmSeedOffset2,readSeedsFromFile,\
                     rndmSeedInputFile,physicsList,overrideMetadata,doBichselSimulation,\
                     IOVDbGlobalTag,SimG4VersionUsed,numberOfCollisions,\