diff --git a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h
index 186509f54164c41c66954bb669a994d6ba16a162..5db0f2ace53d2377e2959f7a8cf3f635145254c6 100644
--- a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h
+++ b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h
@@ -83,7 +83,6 @@ class sTgcDigitMaker : public AthMessaging {
     double lowEdge; // low side of the interval in ns
     double kParameter;
     double thetaParameter;
-    double mostProbableTime;
   };
 
   /**
@@ -144,6 +143,8 @@ class sTgcDigitMaker : public AthMessaging {
 
   /** Find the gamma pdf parameters of a given distance */
   GammaParameter getGammaParameter(double distance) const;
+  /** Get the most probable time of arrival */
+  double getMostProbableArrivalTime(double distance) const;
 
   /** Energy threshold value for each chamber */
   double m_energyThreshold[N_STATIONNAME][N_STATIONETA][N_STATIONPHI][N_MULTIPLET][N_GASGAP][N_CHANNELTYPE]{};
@@ -166,6 +167,8 @@ class sTgcDigitMaker : public AthMessaging {
 
   // Parameters of the gamma pdf required for determining digit time
   std::vector<GammaParameter> m_gammaParameter;
+  // 4th-order polymonial describing the most probable time as function of the distance of closest approach
+  std::vector<double> m_mostProbableArrivalTime;
 
   // Time offset to add to Strip timing
   std::vector<double> m_timeOffsetStrip;
diff --git a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/share/sTGC_Digitization_timeArrival.dat b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/share/sTGC_Digitization_timeArrival.dat
index 21db68524e9ef95c4ffafa7f82e016001eba4485..598864e505bc87ffa336a77201f4eeae641e74b7 100644
--- a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/share/sTGC_Digitization_timeArrival.dat
+++ b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/share/sTGC_Digitization_timeArrival.dat
@@ -3,37 +3,47 @@
 # Time of arrival is obtained from Garfield simulation. It is parametrized as 
 # a function of the distance of closest approach, i.e. the shortest distance 
 # from a particle's trajectory to the nearest wire.
-# The approach is the following: dividing the Garfield data into 9 intervals 
-# of 0.1 mm (distance of closest approach), then fit the time distribution in
-# each interval.
+# The approach is the following: dividing the Garfield data into 11 intervals, 
+# then fit the time distribution in each interval.
 # 
 # In each interval, the time of arrival distribution can be described by a 
 # gamma distribution with three parameters: 
 #  - the usual k and theta parameters, and 
-#  - t0--the shift from zero.
+#  - t0, which is the earliest time of arrival for the interval and also the
+#    only parameter that varies with the distance of closest approach.
 # The gamma probability distribution function is:
 #   P(t) = 1 / (\Gamma(k) * {\theta}^k) * (t - t0)^{k-1} * exp(-(t - t0) / {\theta})
-# where t is the time of arrival to be selected and t0 is the earliest time of 
-# arrival allowed.
+# where t is the time of arrival to be selected.
 #
-# The current implementation in Athena use the parameters k and theta, and the 
-# most probable value (peak of the gamma pdf) instead of t0. First, a gamma 
-# random number generator is launched to choose a time, say a value val, then 
-# it is shifted by the best-fitted most probable value minus (k-1)*theta. The 
-# best-fitted most probable value (mpv) already includes t0, while the product 
-# (k-1)*theta gives the most probable value of the gamma pdf begining at origin.
-# In other word, t0 = mpv - (k-1)*theta and the resulting time of arrival is 
-# given by 
+# The current implementation in Athena uses the parameters k and theta, and the 
+# most probable value of time instead of t0. First, a gamma random number 
+# generator is launched to choose a time, say a value val, then it is shifted 
+# by the best-fitted most probable value minus (k-1)*theta. The  best-fitted 
+# most probable value (mpv) already includes t0, while the product (k-1)*theta 
+# gives the most probable value of a gamma pdf begining at origin.
+# In other word, t0 = mpv - (k-1)*theta for the intervals considered, except for
+# the 1st interval (distance less than 0.050mm) where both t0 and mpv are 
+# set to zero nanosecond. Also, mpv is parametrized as a 4th-order polymonial 
+# in order to get a continuous time of arrival as a function of distance.
+# So the resulting time of arrival is given by 
 #       time = val + (mpv - (k-1) * theta).
 #
 # Below are the parameters to be read. Line starting with '#' is not read.
-# keyword  low_edge[mm]  k  theta  mpv[ns]
-bin 0.0  1.000  0.833  0.0  
-bin 0.1  2.342  0.455  0.899
-bin 0.2  3.474  0.391  2.220
-bin 0.3  4.153  0.376  3.733
-bin 0.4  5.778  0.321  5.441
-bin 0.5  5.046  0.364  7.186
-bin 0.6  5.981  0.392  9.300
-bin 0.7  5.771  0.533  11.997
-bin 0.8  3.068  1.488  16.095
+# IMPORTANT: the intervals must be sorted in ascending order of low_edge.
+# keyword  low_edge[mm]  k  theta
+bin 0.00  0.500  1.460
+bin 0.05  3.236  0.335
+bin 0.15  7.270  0.219
+bin 0.25  8.406  0.236
+bin 0.35 10.809  0.223 
+bin 0.45  8.017  0.267
+bin 0.55  5.859  0.366
+bin 0.65  8.973  0.290
+bin 0.75  9.319  0.290
+bin 0.80  6.282  0.517
+bin 0.85  4.680  1.044
+
+# Most probable time is described using a 4th-order polymonial
+#   mpv(d) = p0 + p1 * d + p2 * d^2 + p3 * d^3 + p4 * d^4
+# keyword  p0  p1  p2  p3  p4
+mpv 0.106 -2.031 71.36 -128.2 87.16
diff --git a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx
index 0aabb85c363424da0b89d71a781d13d2a2571702..0313ba7796c7cd07c5567884f9633fc1ea27e4fc 100644
--- a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx
+++ b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx
@@ -231,33 +231,55 @@ std::unique_ptr<sTgcDigitCollection> sTgcDigitMaker::executeDigi(const sTGCSimHi
   // Distance should be in the range [0, 0.9] mm, unless particle passes through 
   // the wire plane near the edges
   double wire_pitch = detEl->wirePitch();
-  if ((dist_wire > -9.) && (std::abs(dist_wire) > (wire_pitch / 2))) {
-    ATH_MSG_DEBUG("Distance to the nearest wire (" << std::abs(dist_wire) << ") is greater than expected.");
+  // Absolute value of the distance
+  double abs_dist_wire = std::abs(dist_wire);
+  if ((dist_wire > -9.) && (abs_dist_wire > (wire_pitch / 2))) {
+    ATH_MSG_DEBUG("Distance to the nearest wire (" << abs_dist_wire << ") is greater than expected.");
+  }
+
+  // Do not digitize hits that are too far from the nearest wire
+  if (abs_dist_wire > wire_pitch) {
+    return nullptr;
   }
 
   // Get the gamma pdf parameters associated with the distance of closest approach.
-  GammaParameter gamma_par = getGammaParameter(std::abs(dist_wire));
+  //GammaParameter gamma_par = getGammaParameter(abs_dist_wire);
+  double par_kappa = (getGammaParameter(abs_dist_wire)).kParameter; 
+  double par_theta = (getGammaParameter(abs_dist_wire)).thetaParameter; 
+  double most_prob_time = getMostProbableArrivalTime(abs_dist_wire);
   // Compute the most probable value of the gamma pdf
-  double gamma_mpv = (gamma_par.kParameter - 1) * gamma_par.thetaParameter;
-  // If the most probable value is near zero, then ensure it is zero
-  if ((gamma_par.mostProbableTime) < 0.1) {gamma_mpv = 0.;}
-  double t0_par = gamma_par.mostProbableTime - gamma_mpv;
+  double gamma_mpv = (par_kappa - 1) * par_theta;
+  // If the most probable value is less than zero, then set it to zero
+  if (gamma_mpv < 0.) {gamma_mpv = 0.;}
+  double t0_par = most_prob_time - gamma_mpv;
 
   // Digit time follows a gamma distribution, so a value val is 
-  // chosen using a gamma random generator then shifted by t0
+  // chosen using a gamma random generator then is shifted by t0
   // to account for drift time.
   // Note: CLHEP::RandGamma takes the parameters k and lambda, 
   // where lambda = 1 / theta.
-  double digit_time = t0_par + CLHEP::RandGamma::shoot(rndmEngine, gamma_par.kParameter, 1/gamma_par.thetaParameter);
-  if (digit_time < 0.0) {
-    // Ensure the digit time is positive
-    digit_time = -1.0 * digit_time;
+  double digit_time = t0_par + CLHEP::RandGamma::shoot(rndmEngine, par_kappa, 1/par_theta);
+
+  // Sometimes, digit_time is negative because t0_par can be negative. 
+  // In such case, discard the negative value and shoot RandGamma for another value.
+  // However, if that has already been done many times then set digit_time to zero 
+  // in order to avoid runaway loop.
+  const int shoot_limit = 4;
+  int shoot_counter = 0;
+  while (digit_time < 0.) {
+    if (shoot_counter > shoot_limit) {
+      digit_time = 0.;
+      break;
+    }
+    digit_time = t0_par + CLHEP::RandGamma::shoot(rndmEngine, par_kappa, 1/par_theta);
+    ++shoot_counter;
   }
+
   ATH_MSG_DEBUG("sTgcDigitMaker distance = " << dist_wire 
                 << ", time = " << digit_time
-                << ", k parameter = " << gamma_par.kParameter
-                << ", theta parameter = " << gamma_par.thetaParameter
-                << ", most probable time = " << gamma_par.mostProbableTime);
+                << ", k parameter = " << par_kappa
+                << ", theta parameter = " << par_theta
+                << ", most probable time = " << most_prob_time);
 
   //// HV efficiency correction
   if (m_doEfficiencyCorrection){
@@ -1167,9 +1189,12 @@ StatusCode sTgcDigitMaker::readFileOfTimeArrival() {
     std::istringstream iss(line);
     iss >> key;
     if (key.compare("bin") == 0) {
-      iss >> param.lowEdge >> param.kParameter >> param.thetaParameter >> param.mostProbableTime;
+      iss >> param.lowEdge >> param.kParameter >> param.thetaParameter;
       m_gammaParameter.push_back(param);
-    } 
+    } else if (key.compare("mpv") == 0)  {
+      double mpt;
+      while (iss >> mpt) {m_mostProbableArrivalTime.push_back(mpt);}
+    }
   }
 
   // Close the file
@@ -1182,12 +1207,12 @@ sTgcDigitMaker::GammaParameter sTgcDigitMaker::getGammaParameter(double distance
 
   double d = distance;
   if (d < 0.) {
-    ATH_MSG_WARNING("getGammaParameter: expecting a positive distance, but got negative value: " << d
-                     << ". Proceed to the calculation with the absolute value.");
+    ATH_MSG_WARNING("getGammaParameter: expecting a positive distance, but got a negative value: " << d
+                     << ". Proceed to the calculation using its absolute value.");
     d = -1.0 * d;
   }
 
-  // Find the parameters assuming the container is sorted
+  // Find the parameters assuming the container is sorted in ascending order of 'lowEdge'
   int index{-1};
   for (auto& par: m_gammaParameter) {
     if (distance < par.lowEdge) {
@@ -1198,6 +1223,25 @@ sTgcDigitMaker::GammaParameter sTgcDigitMaker::getGammaParameter(double distance
   return m_gammaParameter.at(index);
 }
 
+//+++++++++++++++++++++++++++++++++++++++++++++++
+double sTgcDigitMaker::getMostProbableArrivalTime(double distance) const {
+
+  double d = distance;
+  if (d < 0.) {
+    ATH_MSG_WARNING("getMostProbableArrivalTime: expecting a positive distance, but got a negative value: " << d
+                     << ". Proceed to the calculation using its absolute value.");
+    d = -1.0 * d;
+  }
+
+  double mpt = m_mostProbableArrivalTime.at(0) 
+               + m_mostProbableArrivalTime.at(1) * d
+               + m_mostProbableArrivalTime.at(2) * d * d 
+               + m_mostProbableArrivalTime.at(3) * d * d * d
+               + m_mostProbableArrivalTime.at(4) * d * d * d * d;
+  return mpt;
+}
+
+//+++++++++++++++++++++++++++++++++++++++++++++++
 StatusCode sTgcDigitMaker::readFileOfTimeOffsetStrip() {
   // Verify the file sTGC_Digitization_timeOffsetStrip.dat exists
   const std::string file_name = "sTGC_Digitization_timeOffsetStrip.dat";
@@ -1242,6 +1286,7 @@ StatusCode sTgcDigitMaker::readFileOfTimeOffsetStrip() {
   return StatusCode::SUCCESS;
 }
 
+//+++++++++++++++++++++++++++++++++++++++++++++++
 double sTgcDigitMaker::getTimeOffsetStrip(int neighbor_index) const {
   if ((!m_timeOffsetStrip.empty()) && (neighbor_index >= 0)) {
     // Return the last element if out of range