diff --git a/Generators/DIFGenerator/python/DIFSampler.py b/Generators/DIFGenerator/python/DIFSampler.py
index 007830b63c82ec06f19af6fbd3169cfc10ce8e2e..740fff10a80816485aa6b40e68e81aeb39326e96 100644
--- a/Generators/DIFGenerator/python/DIFSampler.py
+++ b/Generators/DIFGenerator/python/DIFSampler.py
@@ -6,7 +6,7 @@ import ParticleGun as PG
 from math import sqrt, sin, cos, acos
 from random import uniform
 from AthenaCommon.SystemOfUnits import TeV, GeV, MeV
-from AthenaCommon.PhysicalConstants import pi
+from AthenaCommon.PhysicalConstants import pi, c_light
 from ROOT import TLorentzVector
 
 class CylinderSampler(PG.Sampler):
@@ -14,11 +14,12 @@ class CylinderSampler(PG.Sampler):
     Sampler of position 4-vectors within a cylindrical volume.
     """
 
-    def __init__(self, rsq, phi, z, t=0):
+    def __init__(self, rsq, phi, z, t=0, axialTiming = True):
         self.rsq   = rsq
         self.phi = phi
         self.z   = z
         self.t   = t
+        self.axialTiming = axialTiming
 
     @property
     def rsq(self):
@@ -57,6 +58,8 @@ class CylinderSampler(PG.Sampler):
         phi = self.phi()
         z   = self.z()
         t   = self.t()
+        if self.axialTiming:
+            t += z/c_light
         #print "POS =", x, y, z, t
         return TLorentzVector(r * cos(phi), r * sin(phi), z, t)
 
diff --git a/Generators/FaserCosmicGenerator/python/cosmicSampler.py b/Generators/FaserCosmicGenerator/python/cosmicSampler.py
index 7a147e43638b7cc152bb18e2caff0cf30877b40d..dd0c841f0c739733b7e05a952f7802a07b56738f 100644
--- a/Generators/FaserCosmicGenerator/python/cosmicSampler.py
+++ b/Generators/FaserCosmicGenerator/python/cosmicSampler.py
@@ -5,6 +5,7 @@ from random import random
 from math import pi,sin,cos,acos,asin,sqrt
 import FaserCosmicGenerator.Range as r
 from numpy import array,add
+from AthenaCommon.PhysicalConstants import c_light
 import numpy as np
 
 
@@ -54,7 +55,8 @@ class CosmicSampler(Sampler):
         x, y, z    = CR.pos
         px, py, pz = CR.mom
 
-        self.genPosition.SetXYZT(x + self.x0, y + self.y0, z + self.z0, 0)
+        # impose vertical timing constraint
+        self.genPosition.SetXYZT(x + self.x0, y + self.y0, z + self.z0, -(y + self.y0)/(c_light * math.fabs(CR.costh)))
         self.genMomentum.SetPxPyPzE(px,py,pz,CR.Efinal)
 
         particles = []
diff --git a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
index d2e82a88d3fd6229b2e3ad33932ead7d5d650720..f75fabf921090eb0988b950189362a4d465cc1e0 100644
--- a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
+++ b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py
@@ -50,12 +50,14 @@ def FaserParticleGunSingleParticleCfg(ConfigFlags, **kwargs) :
                                            y = kwargs.setdefault("y", 0.0), 
                                            z = kwargs.setdefault("z", -3750.0),
                                            r = kwargs.setdefault("radius", 1.0),
-                                           t = kwargs.setdefault("t", 0.0) )
+                                           t = kwargs.setdefault("t", 0.0),
+                                           axialTiming = kwargs.setdefault("axialTiming", True)  )
     else:
         pg.sampler.pos = PG.PosSampler(x = kwargs.setdefault("x", [-5, 5]), 
                                        y = kwargs.setdefault("y", [-5, 5]), 
                                        z = kwargs.setdefault("z", -3750.0), 
-                                       t = kwargs.setdefault("t", 0.0) )
+                                       t = kwargs.setdefault("t", 0.0),
+                                       axialTiming = kwargs.setdefault("axialTiming", True) )
 
     return cfg
 
@@ -75,12 +77,14 @@ def FaserParticleGunSingleEcalParticleCfg(ConfigFlags, **kwargs) :
                                    y = kwargs.setdefault("y", 0.0), 
                                    z = kwargs.setdefault("z", 0.0),
                                    r = kwargs.setdefault("radius", 1.0),
-                                   t = kwargs.setdefault("t", 0.0) )
+                                   t = kwargs.setdefault("t", 0.0),
+                                   axialTiming = kwargs.setdefault("axialTiming", True) )
     else:
         pg.sampler.pos = PG.PosSampler(x = kwargs.setdefault("x", 0.0), 
                                    y = kwargs.setdefault("y", 0.0), 
                                    z = kwargs.setdefault("z", 0.0),
-                                   t = kwargs.setdefault("t", 0.0) )
+                                   t = kwargs.setdefault("t", 0.0),
+                                   axialTiming = kwargs.setdefault("axialTiming", True) )
 
     return cfg
 
@@ -98,7 +102,8 @@ def FaserParticleGunSingleEcalParticleCfg(ConfigFlags, **kwargs) :
     pg.sampler.pos = PG.PosSampler(x = kwargs.setdefault("x", [-5, 5]), 
                                    y = kwargs.setdefault("y", [-5, 5]), 
                                    z = kwargs.setdefault("z", 2730.0), 
-                                   t = kwargs.setdefault("t", 0.0) )
+                                   t = kwargs.setdefault("t", 0.0),
+                                   axialTiming = kwargs.setdefault("axialTiming", True) )
 
     return cfg
 '''
@@ -150,10 +155,11 @@ def FaserParticleGunDecayInFlightCfg(ConfigFlags, **kwargs) :
                                                   y = kwargs.setdefault("y", 0.0), 
                                                   z = kwargs.setdefault("z", -3750.0),
                                                   r = kwargs.setdefault("radius", 1.0),
-                                                  t = kwargs.setdefault("t", 0.0) )
+                                                  t = kwargs.setdefault("t", 0.0),
+                                                  axialTiming = kwargs.setdefault("axialTiming", True) )
     else:
         from DIFGenerator.DIFSampler import CylinderSampler
-        kwargs["mother_pos"]  = CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0)
+        kwargs["mother_pos"]  = CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0, True)
     
     if not "mother_mom" in kwargs:
         kwargs["mother_mom"] = PG.EThetaMPhiSampler(kwargs.setdefault("energy", ((1*TeV)**2 + (500*MeV)**2)**0.5),
diff --git a/Generators/FaserParticleGun/python/RadialPosSampler.py b/Generators/FaserParticleGun/python/RadialPosSampler.py
index 19b1db116991c082d0771699813df2fe94c0a1d9..8e1efdfd6437e83bf2e04490a1ed478340f833fd 100644
--- a/Generators/FaserParticleGun/python/RadialPosSampler.py
+++ b/Generators/FaserParticleGun/python/RadialPosSampler.py
@@ -1,6 +1,7 @@
 import random
 from math import pi, sin, cos, sqrt, log
 from ParticleGun.samplers import Sampler, mksampler
+from AthenaCommon.PhysicalConstants import c_light
 import ROOT
 
 class RadialPosSampler(Sampler):
@@ -8,12 +9,13 @@ class RadialPosSampler(Sampler):
     Sampler of Position 3-vectors, for modelling a beamspot, based on radius.
     """
 
-    def __init__(self, r, x, y, z, t = 0):
+    def __init__(self, r, x, y, z, t = 0, axialTiming = True):
         self.radius = r
         self.x = x
         self.y = y
         self.z = z
         self.t = t
+        self.axialTiming = axialTiming
 
     @property
     def z(self):
@@ -62,6 +64,8 @@ class RadialPosSampler(Sampler):
         y = self.y + r * sin(phi)
         z = self.z()
         t = self.t()
+        if self.axialTiming:
+            t += z/c_light
 
         return ROOT.TLorentzVector(x, y, z, t)
 
diff --git a/Generators/FlukaReader/python/FlukaReaderAlg.py b/Generators/FlukaReader/python/FlukaReaderAlg.py
index 4cc449e65e7a8285db742e1d05516ddd0915b206..ba4837c14ad0796ef7f357c94b59fa8ab7d7f9a8 100644
--- a/Generators/FlukaReader/python/FlukaReaderAlg.py
+++ b/Generators/FlukaReader/python/FlukaReaderAlg.py
@@ -2,6 +2,7 @@ from AthenaCommon.AppMgr import ServiceMgr as svcMgr
 from GeneratorModules.EvgenAlg import EvgenAlg
 from AthenaPython.PyAthena import StatusCode, EventInfo, EventID, EventType
 from AthenaCommon.SystemOfUnits import GeV, MeV, cm
+from AthenaCommon.PhysicalConstants import c_light
 from AthenaCommon.Constants import DEBUG
 
 from FaserCosmicGenerator import Range
@@ -266,7 +267,8 @@ class FlukaReader(EvgenAlg):
         ROOT.SetOwnership(mcEventInfo, False)
 
         # Create HepMC Vertex
-        pos = HepMC.FourVector(newentry["x"] * cm, newentry["y"] * cm, self.z, 0)
+        # Impose axial timing constraint
+        pos = HepMC.FourVector(newentry["x"] * cm, newentry["y"] * cm, self.z, self.z/c_light)
         gv = HepMC.GenVertex(pos)
 
         ROOT.SetOwnership(gv, False)
diff --git a/Generators/GenieReader/python/GenieReaderAlg.py b/Generators/GenieReader/python/GenieReaderAlg.py
index c97ea494b089a25ff88e4c2a1f2abdd4bf472c6b..373fa709c036b0f6ba28574f5768bcdef93d3de8 100644
--- a/Generators/GenieReader/python/GenieReaderAlg.py
+++ b/Generators/GenieReader/python/GenieReaderAlg.py
@@ -3,7 +3,8 @@
 from AthenaCommon.AppMgr import ServiceMgr as svcMgr
 from GeneratorModules.EvgenAlg import EvgenAlg
 from AthenaPython.PyAthena import StatusCode, EventInfo, EventID, EventType
-from AthenaCommon.SystemOfUnits import GeV, m
+from AthenaCommon.SystemOfUnits import GeV, m, nanosecond
+from AthenaCommon.PhysicalConstants import c_light
 import ROOT
 
 __author__ = "Dave Caser <dcasper@uci.edu>"
@@ -37,7 +38,8 @@ class GenieReader(EvgenAlg):
         ROOT.SetOwnership(mcEventId, False)
         ROOT.SetOwnership(mcEventInfo, False)
 
-        pos = HepMC.FourVector(self.evtStore["vx"]*m, self.evtStore["vy"]*m, self.evtStore["vz"]*m, 0)
+        # impose axial timing constraint
+        pos = HepMC.FourVector(self.evtStore["vx"]*m, self.evtStore["vy"]*m, self.evtStore["vz"]*m, (self.evtStore["vz"]*m/c_light))
         gv = HepMC.GenVertex(pos)
         ROOT.SetOwnership(gv, False)
         evt.add_vertex(gv)
diff --git a/Generators/ParticleGun/python/samplers.py b/Generators/ParticleGun/python/samplers.py
index 9126b607ca7246480e4129669f3933e236c2c126..b53848c0027a2fb832298e03a6e17a9d39a2773f 100644
--- a/Generators/ParticleGun/python/samplers.py
+++ b/Generators/ParticleGun/python/samplers.py
@@ -2,6 +2,7 @@
 
 import ROOT, math, random
 from ParticleGun.histsampling import TH1,TH2
+from AthenaCommon.PhysicalConstants import c_light
 
 ## For convenience
 PI = math.pi
@@ -309,11 +310,12 @@ class PosSampler(Sampler):
     Sampler of position 3-vectors, for modelling a beamspot.
     """
 
-    def __init__(self, x, y, z, t=0):
+    def __init__(self, x, y, z, t=0, axialTiming = False):
         self.x = x
         self.y = y
         self.z = z
         self.t = t
+        self.axialTiming = axialTiming
 
     @property
     def x(self):
@@ -352,6 +354,9 @@ class PosSampler(Sampler):
         y = self.y()
         z = self.z()
         t = self.t()
+        if self.axialTiming:
+            t += z/c_light
+
         #print "POS =", x, y, z, t
         return ROOT.TLorentzVector(x, y, z, t)