diff --git a/Control/CalypsoExample/Generation/python/faser_parser.py b/Control/CalypsoExample/Generation/python/faser_parser.py
index 055dd7200201ce69823d798609cb991d85f4e794..03817efb8d9f8852a1dda650a78f7f4de8283d4e 100644
--- a/Control/CalypsoExample/Generation/python/faser_parser.py
+++ b/Control/CalypsoExample/Generation/python/faser_parser.py
@@ -44,6 +44,10 @@ def faser_pgparser():
                         help="Specify radius (in mm)")
     parser.add_argument("--angle", default=0.005, type=float_or_none,
                         help="Specify angular width (in Rad)")
+    parser.add_argument("--xpos", default=None, type=float,
+                        help="Specify x position of particles (in mm)")
+    parser.add_argument("--ypos", default=None, type=float,
+                        help="Specify y position of particles (in mm)")
     parser.add_argument("--zpos", default=None, type=float,
                         help="Specify z position of particles (in mm) (helpful to avoid FASERnu)")
 
diff --git a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
index 169badecf7cf2a392a95b098a858634c5a2523c6..7ec05e2aaf4df3d434577b3d7f423bcf56856722 100755
--- a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
+++ b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
@@ -148,9 +148,15 @@ if __name__ == '__main__':
 
     # -1000 is safely upstream of detector (to be checked)
     # Note zpos is in mm!
-    if args.zpos:
+    if args.zpos != None:
         sg_dict["z"] = args.zpos
 
+    if args.xpos != None:
+        sg_dict["x"] = args.xpos
+
+    if args.ypos != None:
+        sg_dict["y"] = args.ypos
+
     # Determine energy sampling
     if args.sampler == "lin":
         sg_dict["energy"] = PG.UniformSampler(args.minE*GeV, args.maxE*GeV)
@@ -159,11 +165,28 @@ if __name__ == '__main__':
     elif args.sampler == "const":
         sg_dict["energy"] = PG.ConstSampler(args.maxE*GeV)
     elif args.sampler == "hist":
-        fname, hname = args.hist_name.split(":")
-        sg_dict["energy"] = PG.TH1Sampler(fname, hname)
+        nargs = len(args.hist_name.split(":"))
+        if nargs == 2:
+            fname, hname = args.hist_name.split(":")
+            sg_dict["energy"] = PG.TH1Sampler(fname, hname)
+        elif nargs == 3:
+            fname, hname, scale = args.hist_name.split(":")
+            sg_dict["energy"] = PG.TH1Sampler(fname, hname, scale)
+        else:
+            print(f"Can't parse histogram {args.hist_name}!")
+            sys.exit(1)
+
     elif args.sampler == "hist2D":
-        fname, hname = args.hist_name.split(":")
-        sg_dict["energy"] = PG.TH2Sampler(fname, hname)        
+        nargs = len(args.hist_name.split(":"))
+        if nargs == 2:
+            fname, hname = args.hist_name.split(":")
+            sg_dict["energy"] = PG.TH2Sampler(fname, hname)
+        elif nargs == 4:
+            fname, hname, scalex, scaley = args.hist_name.split(":")
+            sg_dict["energy"] = PG.TH2Sampler(fname, hname, scalex, scaley)
+        else:
+            print(f"Can't parse histogram {args.hist_name}!")
+            sys.exit(1)
     else:
         print(f"Sampler {args.sampler} not known!")
         sys.exit(1)
@@ -255,8 +278,12 @@ if __name__ == '__main__':
 
     b = time.time()
     log.info("Run G4FaserAlg in " + str(b-a) + " seconds")
-#
-# Success should be 0
-#
-    sys.exit(not sc.isSuccess())
+
+# Signal errors
+if sc.isSuccess():
+    log.info("Execution succeeded")
+    sys.exit(0)
+else:
+    log.info("Execution failed, return 1")
+    sys.exit(1)
 
diff --git a/Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh b/Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh
index 9b8676752be083d66acaaf0d0fa4744be43050b5..441ef2923b397274b74f53b8530adf592dc3e76f 100755
--- a/Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh
+++ b/Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh
@@ -185,9 +185,12 @@ cd "${config_file_stem}-${seg_str}"
 # Run job
 if [[ -z "$tag" ]]; then
     faser_particlegun.py  "--conf=$config_path" "--segment=$seg_str"
+    gen_code=$?
 else
     faser_particlegun.py  "--conf=$config_path" "--segment=$seg_str" "--tag=$tag"
+    gen_code=$?
 fi
+echo "Return code: $gen_code"
 #
 # Print out ending time
 date
@@ -198,25 +201,51 @@ export EOS_MGM_URL=root://eospublic.cern.ch
 #
 if ! [ -z "$outdest" ]
 then
+    echo "Output directory:"
     ls -l
-    echo "copy *-HITS.root to $outdest"
-    mkdir -p $outdest
-    eos cp *-HITS.root ${outdest}/ || true
+    thefile=`ls *-HITS.root`
+    if [ $? -eq 0 ]; then
+	echo "copy $thefile to $outdest"
+	eos mkdir -p $outdest
+	eos cp $thefile ${outdest}/${thefile} || true
+
+	# Check that it worked
+	eos ls ${outdest}/${thefile} > /dev/null
+	if [ $? -eq 0 ]; then
+	    echo "file $thefile copied to $outdest"
+	    copy_code=0
+	else
+	    echo "didnt find $thefile in $outdest !"
+	    copy_code=1
+	fi
+    else
+	echo "ls *-xAOD.root returned nothing!"
+	copy_code=1
+    fi
+
 fi
 #
 # Also copy log file
 if ! [ -z "$logdest" ]
 then
     cd ..
+    echo "Working directory:"
     ls -l
     echo "copy $logfile to $logdest"
-    mkdir -p $logdest
+    eos mkdir -p $logdest
     eos cp $logfile $logdest/$logfile
 elif ! [ -z "$outdest" ]
 then 
     cd ..
     ls -l
     echo "copy $logfile to $outdest"
-    mkdir -p $outdest
+    eos mkdir -p $outdest
     eos cp $logfile $outdest/$logfile
 fi
+
+# Make sure to return an error is calypso failed
+if [ $gen_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+  exit 1
+else
+  exit 0
+fi
diff --git a/Generators/ParticleGun/python/histsampling.py b/Generators/ParticleGun/python/histsampling.py
index 9f541da7c972ac01170b8acfafc8001b6eea5a14..10549ef0b250b999c7ac66063bdb1d5ae7c90147 100644
--- a/Generators/ParticleGun/python/histsampling.py
+++ b/Generators/ParticleGun/python/histsampling.py
@@ -14,6 +14,7 @@ def load_hist(*args):
     """
     Load a histogram from a filename/TFile and histo name. If a single arg is
     provided, it has to be a histo object and will be cloned before return.
+
     """
     h = None
     if len(args) == 1 and issubclass(type(args[0]), ROOT.TH1):
@@ -25,6 +26,7 @@ def load_hist(*args):
             #f.Close()
         elif type(args[0]) is ROOT.TFile and type(args[1]) is str:
             h = args[0].Get(args[1]).Clone()
+
     if h is None:
         raise Exception("Error in histogram loading from " + args)
     else:
@@ -77,7 +79,7 @@ def get_random_bin(globalbins, cheights):
     raise Exception("Sample fell outside range of cumulative distribution?!?!")
 
 
-def get_random_x(h, globalbins, cheights, globalbin_to_axisbin):
+def get_random_x(h, globalbins, cheights, globalbin_to_axisbin, scale=1.):
     """
     Choose a random bin via get_random_bin, then pick a uniform random x
     point in that bin (without any attempt at estimating the in-bin distribution).
@@ -86,10 +88,10 @@ def get_random_x(h, globalbins, cheights, globalbin_to_axisbin):
     axisids = globalbin_to_axisbin.get(irand)
     assert axisids is not None
     xrand = random.uniform(h.GetXaxis().GetBinLowEdge(axisids[0]), h.GetXaxis().GetBinUpEdge(axisids[0]))
-    return xrand
+    return scale * xrand
 
 
-def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin):
+def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin, xscale=1., yscale=1.):
     """
     Choose a random bin via get_random_bin, then pick a uniform random x,y
     point in that bin (without any attempt at estimating the in-bin distribution).
@@ -99,34 +101,71 @@ def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin):
     assert axisids is not None
     xrand = random.uniform(h2.GetXaxis().GetBinLowEdge(axisids[0]), h2.GetXaxis().GetBinUpEdge(axisids[0]))
     yrand = random.uniform(h2.GetYaxis().GetBinLowEdge(axisids[1]), h2.GetYaxis().GetBinUpEdge(axisids[1]))
-    return xrand, yrand
+    return (xscale*xrand), (yscale*yrand)
 
 
 class TH1(object):
     "Minimal wrapper for ROOT TH1, for sampling consistency and easy loading"
 
     def __init__(self, *args):
-        self.th1 = load_hist(*args)
+        """ Args can be variable length, but these are now kwargs, so order matters
+
+        fname - filename
+        hname - histogram name
+        xscale - scaling factor for x axis variable
+        """
+        if len(args) > 2:
+            self.th1 = load_hist(*(args[0:2]))
+        else:
+            self.th1 = load_hist(*args)
+
         self.globalbins, self.globalbin_to_axisbin, self.cheights = None, None, None
 
+        if len(args) >= 3:
+            self.xscale = float(args[2])
+        else:
+            self.xscale = 1.
+
     def GetRandom(self):
         "A GetRandom that works for TH1s and uses Python random numbers"
         if self.globalbins is None or self.globalbin_to_axisbin is None or self.cheights is None:
             self.globalbins, self.globalbin_to_axisbin, self.cheights = get_sampling_vars(self.th1)
-        return get_random_x(self.th1, self.globalbins, self.cheights, self.globalbin_to_axisbin)
+        return get_random_x(self.th1, self.globalbins, self.cheights, self.globalbin_to_axisbin, self.xscale)
 
 
 class TH2(object):
     "Minimal wrapper for ROOT TH2, for easy loading and to allow 2D sampling"
 
     def __init__(self, *args):
-        self.th2 = load_hist(*args)
+        """ Args can be variable length, but these are now kwargs, so order matters
+
+        fname - filename
+        hname - histogram name
+        xscale - scaling factor for x axis variable
+        yscale - scaling factor for y axis variable
+        """
+        if len(args) > 2:
+            self.th2 = load_hist(*(args[0:2]))
+        else:
+            self.th2 = load_hist(*args)
+
         self.globalbins, self.globalbin_to_axisbin, self.cheights = None, None, None
 
+        if len(args) >= 3:
+            self.xscale = float(args[2])
+        else:
+            self.xscale = 1.
+
+        if len(args) >= 4:
+            self.yscale = float(args[3])
+        else:
+            self.yscale = 1.
+
+
     def GetRandom(self):
         "A GetRandom that works for TH2s"
         if self.globalbins is None or self.globalbin_to_axisbin is None or self.cheights is None:
             self.globalbins, self.globalbin_to_axisbin, self.cheights = get_sampling_vars(self.th2)
-        return get_random_xy(self.th2, self.globalbins, self.cheights, self.globalbin_to_axisbin)
+        return get_random_xy(self.th2, self.globalbins, self.cheights, self.globalbin_to_axisbin, self.xscale, self.yscale)
 
 
diff --git a/Generators/ParticleGun/python/samplers.py b/Generators/ParticleGun/python/samplers.py
index d6e540ea6cf7a747d9152202cd488c135fa27099..9126b607ca7246480e4129669f3933e236c2c126 100644
--- a/Generators/ParticleGun/python/samplers.py
+++ b/Generators/ParticleGun/python/samplers.py
@@ -633,14 +633,23 @@ class EThetaMPhiSampler(MomSampler):
         pt = p sin(theta)
         """
 
-        if self._theta is None:
-            e,theta = self.energy()
-        else:            
-            e = self.energy()
-            theta = self.theta()
-            
         m = self.mass()
-        p = math.sqrt( e**2 - m**2 )
+
+        count = 0
+        e = -1
+        while (e < m and count < 5):
+            count += 1
+            if self._theta is None:
+                e,theta = self.energy()
+            else:            
+                e = self.energy()
+                theta = self.theta()
+
+        try:            
+            p = math.sqrt( e**2 - m**2 )
+        except Exception:
+            raise Exception(f"Error generating E: {e} m: {m}!")
+
         pz = p * math.cos(theta)
         pt = p * math.sin(theta)
         phi = self.phi()
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8d8665ba924345e2fdcdb409384cdb66291446c3
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
@@ -0,0 +1,12 @@
+atlas_subdir(OverlayRDO)
+
+atlas_add_component(
+        OverlayRDO
+        src/OverlayRDOAlg.h src/SelectRDOAlg.h src/BinRDOAlg.h src/HistoRDOAlg.h
+        src/OverlayRDOAlg.cxx src/SelectRDOAlg.cxx src/BinRDOAlg.cxx src/HistoRDOAlg.cxx
+        src/component/OverlayRDO_entries.cxx
+        LINK_LIBRARIES AthenaBaseComps StoreGateLib TrackerRawData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack xAODFaserWaveform
+)
+
+atlas_install_python_modules(python/*.py)
+# atlas_install_scripts(test/*.py)
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..cfb2dfefe48927be61f6416bd206a78e0b754eab
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
@@ -0,0 +1,127 @@
+"""
+    Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+# from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+def HistoRDOAlgCfg(flags, **kwargs):
+
+    from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+    acc = FaserGeometryCfg(flags)
+    HistoRDOAlg = CompFactory.HistoRDOAlg("HistoRDOAlg",**kwargs)
+    acc.addEventAlgo(HistoRDOAlg)
+    HistoRDOAlg.OutputLevel = VERBOSE
+    # ItemList = []
+    # # ItemList += ["xAOD::EventInfo#" + chargePrefix + "EventInfo"]
+    # # ItemList += ["xAOD::EventAuxInfo#" + chargePrefix + "EventInfoAux."]
+    # ItemList += ["xAOD::EventInfo#EventInfo"]
+    # ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
+    # ItemList += ["TrackCollection#Orig_CKFTrackCollectionWithoutIFT"]
+    # ItemList += ["FaserSCT_RDO_Container#SCT_RDOs"]
+    # ItemList += ["FaserSCT_RDO_Container#SCT_EDGEMODE_RDOs"]
+    # # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
+
+    # acc.merge(OutputStreamCfg(ConfigFlags,"RDO", ItemList=ItemList, disableEventTag=True))
+
+    thistSvc = CompFactory.THistSvc()
+    thistSvc.Output += ["HIST2 DATAFILE='RDOtree.root' OPT='RECREATE'"]
+    acc.addService(thistSvc)
+
+    return acc
+
+if __name__ == "__main__":
+
+    import sys
+    from AthenaCommon.Logging import log, logging
+    from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+    from AthenaCommon.Configurable import Configurable
+    from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+    from AthenaConfiguration.TestDefaults import defaultTestFiles
+    from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+    from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+    # from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+    # Set up logging and new style config
+    log.setLevel(DEBUG)
+    Configurable.configurableRun3Behavior = True
+
+    # Configure
+    # ConfigFlags.Overlay.DataOverlay = False
+    ConfigFlags.Input.Files = [ 
+        '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin0RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin1RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin2RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin3RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin4RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin5RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin6RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin7RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin8RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin0RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin1RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin2RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin3RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin4RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin5RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin6RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin7RDO.pool.root',
+        '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin8RDO.pool.root',
+    ]
+    # ConfigFlags.Input.SecondaryFiles = [ 'Neg_RDO.pool.root' ]
+    # ConfigFlags.Input.Files = [ '/eos/experiment/faser/rec/2022/r0013/009171/Faser-Physics-009171-00006-r0013-xAOD.root']
+    # ConfigFlags.Input.SecondaryFiles = [ '/eos/experiment/faser/rec/2022/r0013/009166/Faser-Physics-009166-00485-r0013-xAOD.root' ]
+    # ConfigFlags.Output.RDOFileName = "Overlay.RDO.pool.root"
+    ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03"             # Always needed; must match FaserVersionS
+    ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3"               # Use data conditions for now
+    ConfigFlags.Input.ProjectName = "data21"                     # Needed to bypass autoconfig
+    ConfigFlags.Input.isMC = False                                # Needed to bypass autoconfig
+    ConfigFlags.GeoModel.FaserVersion     = "FASERNU-03"           # FASER geometry
+    ConfigFlags.Common.isOnline = False
+    ConfigFlags.GeoModel.Align.Dynamic = False
+    ConfigFlags.Beam.NumberOfCollisions = 0.
+
+    ConfigFlags.Detector.GeometryFaserSCT = True
+
+    ConfigFlags.lock()
+
+    # Core components
+    acc = MainServicesCfg(ConfigFlags)
+    acc.merge(PoolReadCfg(ConfigFlags))
+
+    # algorithm
+    acc.merge(HistoRDOAlgCfg(ConfigFlags))
+
+    # from SGComps.AddressRemappingConfig import AddressRemappingCfg
+    # acc.merge(AddressRemappingCfg([
+    #     "xAOD::EventInfo#EventInfo->" + ConfigFlags.Overlay.SigPrefix + "EventInfo",
+    #     "xAOD::EventAuxInfo#EventInfoAux.->" + ConfigFlags.Overlay.SigPrefix + "EventInfoAux.",
+    # ]))
+
+    # Hack to avoid problem with our use of MC databases when isMC = False
+    replicaSvc = acc.getService("DBReplicaSvc")
+    replicaSvc.COOLSQLiteVetoPattern = ""
+    replicaSvc.UseCOOLSQLite = True
+    replicaSvc.UseCOOLFrontier = False
+    replicaSvc.UseGeomSQLite = True
+
+    # Timing
+    #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+    # Dump config
+    # logging.getLogger('forcomps').setLevel(VERBOSE)
+    # acc.foreach_component("*").OutputLevel = VERBOSE
+    # acc.foreach_component("*ClassID*").OutputLevel = INFO
+    # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+    # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+    # acc.getService("StoreGateSvc").Dump = True
+    # acc.getService("ConditionStore").Dump = True
+    # acc.printConfig(withDetails=True)
+    # ConfigFlags.dump()
+
+    # Execute and finish
+    sc = acc.run(maxEvents=-1)
+
+    # Success should be 0
+    sys.exit(not sc.isSuccess())    
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..bce33065b965395662707d54a82606d1e7038b9e
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
@@ -0,0 +1,106 @@
+"""
+    Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+def OverlayRDOAlgCfg(flags, **kwargs):
+
+    from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+    acc = FaserGeometryCfg(flags)
+    OverlayRDOAlg = CompFactory.OverlayRDOAlg("OverlayRDOAlg",**kwargs)
+    acc.addEventAlgo(OverlayRDOAlg)
+
+    ItemList = []
+    # ItemList += ["xAOD::EventInfo#" + chargePrefix + "EventInfo"]
+    # ItemList += ["xAOD::EventAuxInfo#" + chargePrefix + "EventInfoAux."]
+    ItemList += ["xAOD::EventInfo#EventInfo"]
+    ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
+    ItemList += ["TrackCollection#Orig_CKFTrackCollectionWithoutIFT"]
+    ItemList += ["FaserSCT_RDO_Container#SCT_RDOs"]
+    ItemList += ["FaserSCT_RDO_Container#SCT_EDGEMODE_RDOs"]
+    # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,"RDO", ItemList=ItemList, disableEventTag=True))
+
+
+
+    return acc
+
+if __name__ == "__main__":
+
+    import sys
+    from AthenaCommon.Logging import log, logging
+    from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+    from AthenaCommon.Configurable import Configurable
+    from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+    from AthenaConfiguration.TestDefaults import defaultTestFiles
+    from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+    from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+    # from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+    # Set up logging and new style config
+    log.setLevel(DEBUG)
+    Configurable.configurableRun3Behavior = True
+
+    # Configure
+    ConfigFlags.Overlay.DataOverlay = False
+    ConfigFlags.Input.Files = [ 'Pos_RDO.pool.root']
+    ConfigFlags.Input.SecondaryFiles = [ 'Neg_RDO.pool.root' ]
+    # ConfigFlags.Input.Files = [ '/eos/experiment/faser/rec/2022/r0013/009171/Faser-Physics-009171-00006-r0013-xAOD.root']
+    # ConfigFlags.Input.SecondaryFiles = [ '/eos/experiment/faser/rec/2022/r0013/009166/Faser-Physics-009166-00485-r0013-xAOD.root' ]
+    ConfigFlags.Output.RDOFileName = "Overlay.RDO.pool.root"
+    ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03"             # Always needed; must match FaserVersionS
+    ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3"               # Use data conditions for now
+    ConfigFlags.Input.ProjectName = "data21"                     # Needed to bypass autoconfig
+    ConfigFlags.Input.isMC = False                                # Needed to bypass autoconfig
+    ConfigFlags.GeoModel.FaserVersion     = "FASERNU-03"           # FASER geometry
+    ConfigFlags.Common.isOnline = False
+    ConfigFlags.GeoModel.Align.Dynamic = False
+    ConfigFlags.Beam.NumberOfCollisions = 0.
+
+    ConfigFlags.Detector.GeometryFaserSCT = True
+
+    ConfigFlags.lock()
+
+    # Core components
+    acc = MainServicesCfg(ConfigFlags)
+    acc.merge(PoolReadCfg(ConfigFlags))
+
+    # algorithm
+    acc.merge(OverlayRDOAlgCfg(ConfigFlags))
+
+    # from SGComps.AddressRemappingConfig import AddressRemappingCfg
+    # acc.merge(AddressRemappingCfg([
+    #     "xAOD::EventInfo#EventInfo->" + ConfigFlags.Overlay.SigPrefix + "EventInfo",
+    #     "xAOD::EventAuxInfo#EventInfoAux.->" + ConfigFlags.Overlay.SigPrefix + "EventInfoAux.",
+    # ]))
+
+    # Hack to avoid problem with our use of MC databases when isMC = False
+    replicaSvc = acc.getService("DBReplicaSvc")
+    replicaSvc.COOLSQLiteVetoPattern = ""
+    replicaSvc.UseCOOLSQLite = True
+    replicaSvc.UseCOOLFrontier = False
+    replicaSvc.UseGeomSQLite = True
+
+    # Timing
+    #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+    # Dump config
+    # logging.getLogger('forcomps').setLevel(VERBOSE)
+    # acc.foreach_component("*").OutputLevel = VERBOSE
+    # acc.foreach_component("*ClassID*").OutputLevel = INFO
+    # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+    # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+    # acc.getService("StoreGateSvc").Dump = True
+    # acc.getService("ConditionStore").Dump = True
+    # acc.printConfig(withDetails=True)
+    # ConfigFlags.dump()
+
+    # Execute and finish
+    sc = acc.run(maxEvents=-1)
+
+    # Success should be 0
+    sys.exit(not sc.isSuccess())    
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..29534a9b3293511b430305770da134c073b10213
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
@@ -0,0 +1,156 @@
+"""
+    Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+chargePrefix = "Pos_"
+
+def SelectRDOAlgCfg(flags, **kwargs):
+
+    from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+    acc = FaserGeometryCfg(flags)
+    SelectRDOAlg = CompFactory.SelectRDOAlg("SelectRDOAlg",**kwargs)
+    SelectRDOAlg.AcceptPositive = ("Pos_" in chargePrefix)
+    acc.addEventAlgo(SelectRDOAlg)
+    acc.addEventAlgo(CompFactory.BinRDOAlg("Bin0RDO",Xmin=-66, Xmax=0,  Ymin=-66, Ymax=0))
+    acc.addEventAlgo(CompFactory.BinRDOAlg("Bin1RDO",Xmin=0,   Xmax=66, Ymin=-66, Ymax=0))
+    acc.addEventAlgo(CompFactory.BinRDOAlg("Bin2RDO",Xmin=-66, Xmax=0,  Ymin=0,   Ymax=66))
+    acc.addEventAlgo(CompFactory.BinRDOAlg("Bin3RDO",Xmin=0,   Xmax=66, Ymin=0,   Ymax=66))
+
+    acc.addEventAlgo(CompFactory.BinRDOAlg("Bin4RDO",Xmin=-33, Xmax=33,  Ymin=-33, Ymax=33))
+    acc.addEventAlgo(CompFactory.BinRDOAlg("Bin5RDO",Xmin=-33, Xmax=33,  Ymin=33,  Ymax=99))
+    acc.addEventAlgo(CompFactory.BinRDOAlg("Bin6RDO",Xmin=-33, Xmax=33,  Ymin=-99, Ymax=-33))
+    acc.addEventAlgo(CompFactory.BinRDOAlg("Bin7RDO",Xmin=-99, Xmax=-33, Ymin=-33, Ymax=33))
+    acc.addEventAlgo(CompFactory.BinRDOAlg("Bin8RDO",Xmin=33,  Xmax=99,  Ymin=-33, Ymax=33))
+
+    ItemList = []
+    # ItemList += ["xAOD::EventInfo#" + chargePrefix + "EventInfo"]
+    # ItemList += ["xAOD::EventAuxInfo#" + chargePrefix + "EventInfoAux."]
+    ItemList += ["xAOD::EventInfo#EventInfo"]
+    ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
+    ItemList += ["TrackCollection#" + chargePrefix + "CKFTrackCollectionWithoutIFT"]
+    ItemList += ["FaserSCT_RDO_Container#" + chargePrefix + "SCT_RDOs"]
+    ItemList += ["FaserSCT_RDO_Container#" + chargePrefix + "SCT_EDGEMODE_RDOs"]
+    # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin0RDO", ItemList=ItemList, disableEventTag=True))
+    osrdo0 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin0RDO")
+    osrdo0.RequireAlgs += ["SelectRDOAlg","Bin0RDO"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin1RDO", ItemList=ItemList, disableEventTag=True))
+    osrdo1 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin1RDO")
+    osrdo1.RequireAlgs += ["SelectRDOAlg","Bin1RDO"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin2RDO", ItemList=ItemList, disableEventTag=True))
+    osrdo2 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin2RDO")
+    osrdo2.RequireAlgs += ["SelectRDOAlg","Bin2RDO"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin3RDO", ItemList=ItemList, disableEventTag=True))
+    osrdo3 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin3RDO")
+    osrdo3.RequireAlgs += ["SelectRDOAlg","Bin3RDO"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin4RDO", ItemList=ItemList, disableEventTag=True))
+    osrdo4 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin4RDO")
+    osrdo4.RequireAlgs += ["SelectRDOAlg","Bin4RDO"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin5RDO", ItemList=ItemList, disableEventTag=True))
+    osrdo5 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin5RDO")
+    osrdo5.RequireAlgs += ["SelectRDOAlg","Bin5RDO"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin6RDO", ItemList=ItemList, disableEventTag=True))
+    osrdo6 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin6RDO")
+    osrdo6.RequireAlgs += ["SelectRDOAlg","Bin6RDO"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin7RDO", ItemList=ItemList, disableEventTag=True))
+    osrdo7 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin7RDO")
+    osrdo7.RequireAlgs += ["SelectRDOAlg","Bin7RDO"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin8RDO", ItemList=ItemList, disableEventTag=True))
+    osrdo8 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin8RDO")
+    osrdo8.RequireAlgs += ["SelectRDOAlg","Bin8RDO"]
+
+    return acc
+
+
+
+if __name__ == "__main__":
+
+    import sys
+    from AthenaCommon.Logging import log, logging
+    from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+    from AthenaCommon.Configurable import Configurable
+    from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+    from AthenaConfiguration.TestDefaults import defaultTestFiles
+    from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+    from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+
+    # Set up logging and new style config
+    log.setLevel(DEBUG)
+    Configurable.configurableRun3Behavior = True
+
+    # Configure
+    ConfigFlags.Input.Files = [ #'/eos/experiment/faser/rec/2022/r0013/009171/Faser-Physics-009171-00006-r0013-xAOD.root' 
+                                '/eos/experiment/faser/rec/2022/r0013/009166/Faser-Physics-009166-00485-r0013-xAOD.root'
+                              ]
+    # ConfigFlags.Input.Files = [ 'Faser-Physics-009171-00006-r0013-xAOD.root', 
+    #                             'Faser-Physics-009166-00485-r0013-xAOD.root']
+    # ConfigFlags.Output.RDOFileName = chargePrefix + "RDO.pool.root"
+    ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03"             # Always needed; must match FaserVersionS
+    ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3"               # Use data conditions for now
+    ConfigFlags.Input.ProjectName = "data21"                     # Needed to bypass autoconfig
+    ConfigFlags.Input.isMC = False                                # Needed to bypass autoconfig
+    ConfigFlags.GeoModel.FaserVersion     = "FASERNU-03"           # FASER geometry
+    ConfigFlags.Common.isOnline = False
+    ConfigFlags.GeoModel.Align.Dynamic = False
+    ConfigFlags.Beam.NumberOfCollisions = 0.
+
+    ConfigFlags.Detector.GeometryFaserSCT = True
+
+    ConfigFlags.lock()
+
+    # Core components
+    acc = MainServicesCfg(ConfigFlags)
+    acc.merge(PoolReadCfg(ConfigFlags))
+
+    # algorithm
+    acc.merge(SelectRDOAlgCfg(ConfigFlags, TrackCollection = "CKFTrackCollectionWithoutIFT", OutputTrackCollection = chargePrefix + "CKFTrackCollectionWithoutIFT"))
+
+    from SGComps.AddressRemappingConfig import AddressRemappingCfg
+    acc.merge(AddressRemappingCfg([
+        # "xAOD::EventInfo#EventInfo->" + chargePrefix + "EventInfo",
+        # "xAOD::EventAuxInfo#EventInfoAux.->" + chargePrefix + "EventInfoAux.",
+        # "TrackCollection#CKFTrackCollectionWithoutIFT->" + chargePrefix + "CKFTrackCollectionWithoutIFT",
+        "FaserSCT_RDO_Container#SCT_RDOs->" + chargePrefix + "SCT_RDOs",
+        "FaserSCT_RDO_Container#SCT_EDGEMODE_RDOs->" + chargePrefix + "SCT_EDGEMODE_RDOs",
+        # "Tracker::FaserSCT_ClusterContainer#SCT_ClusterContainer->" + chargePrefix + "SCT_ClusterContainer"
+    ]))
+
+    # Hack to avoid problem with our use of MC databases when isMC = False
+    replicaSvc = acc.getService("DBReplicaSvc")
+    replicaSvc.COOLSQLiteVetoPattern = ""
+    replicaSvc.UseCOOLSQLite = True
+    replicaSvc.UseCOOLFrontier = False
+    replicaSvc.UseGeomSQLite = True
+
+    # Timing
+    #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+    # Dump config
+    # logging.getLogger('forcomps').setLevel(VERBOSE)
+    # acc.foreach_component("*").OutputLevel = VERBOSE
+    # acc.foreach_component("*ClassID*").OutputLevel = INFO
+    # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+    # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+    # acc.getService("StoreGateSvc").Dump = True
+    # acc.getService("ConditionStore").Dump = True
+    # acc.printConfig(withDetails=True)
+    # ConfigFlags.dump()
+
+    # Execute and finish
+    sc = acc.run(maxEvents=-1)
+
+    # Success should be 0
+    sys.exit(not sc.isSuccess())    
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a917447dde480b9b82a316653fb8e76cf76a2156
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx
@@ -0,0 +1,55 @@
+#include "BinRDOAlg.h"
+#include "TrkTrack/Track.h"
+
+BinRDOAlg::BinRDOAlg(const std::string &name, 
+                           ISvcLocator *pSvcLocator)
+    : AthReentrantAlgorithm(name, pSvcLocator)
+{}
+
+
+StatusCode BinRDOAlg::initialize() 
+{
+  ATH_CHECK(m_trackCollection.initialize());
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode BinRDOAlg::execute(const EventContext &ctx) const 
+{
+
+  setFilterPassed(false, ctx);
+  SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx};
+  if (!trackCollection.isValid() || trackCollection->size() == 0) return StatusCode::SUCCESS;
+
+  // track cuts
+
+  double xUpstream {0.0};
+  double yUpstream {0.0};
+
+  for (const Trk::Track* track : *trackCollection)
+  {
+    if (track == nullptr) continue;
+    const Trk::TrackParameters* upstreamParameters {nullptr};
+    for (auto params : *(track->trackParameters()))
+    {
+        if (params->position().z() < 0) continue;  // Ignore IFT hits
+        if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+    }
+    xUpstream = upstreamParameters->position().x();
+    yUpstream = upstreamParameters->position().y();
+    break;
+  }
+
+
+  // set filter passed
+
+  setFilterPassed(((xUpstream >= m_xMin) && (xUpstream <= m_xMax) && (yUpstream >= m_yMin) && (yUpstream <= m_yMax)), ctx);
+  return StatusCode::SUCCESS;
+}
+
+
+StatusCode BinRDOAlg::finalize() 
+{
+  return StatusCode::SUCCESS;
+}
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..f8e86e844adf11e7a1a4021b9e8d7a11757ce58c
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h
@@ -0,0 +1,23 @@
+#pragma once
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "TrkTrack/TrackCollection.h"
+
+class BinRDOAlg : public AthReentrantAlgorithm {
+public:
+  BinRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+  virtual ~BinRDOAlg() = default;
+  virtual StatusCode initialize() override;
+  virtual StatusCode execute(const EventContext &ctx) const override;
+  virtual StatusCode finalize() override;
+
+private:
+
+  SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "CKFTrackCollectionWithoutIFT", "Input track collection name" };
+  DoubleProperty     m_xMin {this, "Xmin", -25.0, "Minimum x position to accept"};
+  DoubleProperty     m_xMax {this, "Xmax",  25.0, "Maximum x position to accept"};
+  DoubleProperty     m_yMin {this, "Ymin", -25.0, "Minimum y position to accept"};
+  DoubleProperty     m_yMax {this, "Ymax",  25.0, "Maximum y position to accept"};
+
+};
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..78096c60004a9e183bf83aa235ef72501ddee4c4
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
@@ -0,0 +1,247 @@
+#include "HistoRDOAlg.h"
+#include "TrkTrack/Track.h"
+#include <TTree.h>
+
+
+HistoRDOAlg::HistoRDOAlg(const std::string &name, 
+                         ISvcLocator *pSvcLocator)
+    : AthReentrantAlgorithm(name, pSvcLocator), 
+      AthHistogramming(name),
+      m_histSvc("THistSvc/THistSvc", name) {}
+
+StatusCode HistoRDOAlg::initialize()
+{
+    m_tree = new TTree("overlay", "Overlay tree");
+
+    m_tree->Branch("longTracks1", &m_longTracks1, "longTracks1/I");
+    m_tree->Branch("x1", &m_x1);
+    m_tree->Branch("y1", &m_y1);
+    m_tree->Branch("p1", &m_p1);
+    m_tree->Branch("theta1", &m_theta1);
+    m_tree->Branch("phi1", &m_phi1);
+    m_tree->Branch("charge1", &m_charge1);
+    m_tree->Branch("chiSquared1", &m_chiSquared1);
+    m_tree->Branch("nDoF1", &m_nDoF1);
+    m_tree->Branch("pPolar1", &m_pTotPolar1, "pPolar1/D");
+
+    m_tree->Branch("longTracks2", &m_longTracks2, "longTracks2/I");
+    m_tree->Branch("x2", &m_x2);
+    m_tree->Branch("y2", &m_y2);
+    m_tree->Branch("p2", &m_p2);
+    m_tree->Branch("theta2", &m_theta2);
+    m_tree->Branch("phi2", &m_phi2);
+    m_tree->Branch("charge2", &m_charge2);
+    m_tree->Branch("chiSquared2", &m_chiSquared2);
+    m_tree->Branch("nDoF2", &m_nDoF2);
+    m_tree->Branch("pPolar2", &m_pTotPolar2, "pPolar2/D");
+
+    m_tree->Branch("yAgreement", &m_yAgreement, "yAgreement/D");
+    m_tree->Branch("matchFraction", &m_matchFraction, "matchFraction/D");
+
+    ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree));
+
+    ATH_CHECK(m_trackCollection1.initialize());
+    ATH_CHECK(m_trackCollection2.initialize());
+
+    return StatusCode::SUCCESS;
+}
+
+StatusCode HistoRDOAlg::execute(const EventContext& ctx) const
+{
+
+    m_longTracks1 = 0;
+    m_x1.clear();
+    m_y1.clear();
+    m_p1.clear();
+    m_theta1.clear();
+    m_phi1.clear();
+    m_charge1.clear();
+    m_chiSquared1.clear();
+    m_nDoF1.clear();
+    m_pTotPolar1 = 0;
+
+    m_longTracks2 = 0;
+    m_x2.clear();
+    m_y2.clear();
+    m_p2.clear();
+    m_theta2.clear();
+    m_phi2.clear();
+    m_charge2.clear();
+    m_chiSquared2.clear();
+    m_nDoF2.clear();
+    m_pTotPolar2 = 0;
+
+    m_yAgreement = 0;
+    m_matchFraction = 0;
+
+    SG::ReadHandle<TrackCollection> trackCollection1 {m_trackCollection1, ctx};
+    // ATH_CHECK(trackCollection1.isValid());
+    SG::ReadHandle<TrackCollection> trackCollection2 {m_trackCollection2, ctx};
+    // ATH_CHECK(trackCollection2.isValid());
+
+    std::vector<const Trk::TrackParameters*> upstreamSingle {};
+    if (trackCollection1.isValid())
+    {
+        Amg::Vector3D pTot {0, 0, 0};
+        for (const Trk::Track* track : *trackCollection1)
+        {
+            if (track == nullptr) continue;
+            m_longTracks1++;
+            const Trk::TrackParameters* upstreamParameters {nullptr};
+            for (auto params : *(track->trackParameters()))
+            {
+                if (params->position().z() < 0) continue;  // Ignore IFT hits
+                if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+                ATH_MSG_VERBOSE("TrackCollection 1, Track " << m_longTracks1 << ", z = " << params->position().z() << ", y = " << params->position().y() << ", q/p(TeV) = " << params->charge()/(params->momentum().mag()/1e6));
+            }
+            upstreamSingle.push_back(upstreamParameters);
+            m_x1.push_back(upstreamParameters->position().x());
+            m_y1.push_back(upstreamParameters->position().y());
+            m_p1.push_back(upstreamParameters->momentum().mag());
+            m_charge1.push_back(upstreamParameters->charge());
+            m_theta1.push_back(asin(upstreamParameters->momentum().perp()/upstreamParameters->momentum().mag()));
+            m_phi1.push_back(atan2(upstreamParameters->momentum().y(), upstreamParameters->momentum().x()));
+            m_chiSquared1.push_back(track->fitQuality()->chiSquared());
+            m_nDoF1.push_back(track->fitQuality()->doubleNumberDoF());
+            pTot += upstreamParameters->momentum();
+        }
+        m_pTotPolar1 = asin(pTot.perp()/pTot.mag());
+    }
+
+    std::vector<const Trk::TrackParameters*> upstreamOverlay {};
+    if (trackCollection2.isValid())
+    {
+        Amg::Vector3D pTot {0, 0, 0};
+        for (const Trk::Track* track : *trackCollection2)
+        {
+            if (track == nullptr) continue;
+            m_longTracks2++;
+            const Trk::TrackParameters* upstreamParameters {nullptr};
+            for (auto params : *(track->trackParameters()))
+            {
+                if (params->position().z() < 0) continue;  // Ignore IFT hits
+                if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+                ATH_MSG_VERBOSE("TrackCollection 2, Track " << m_longTracks2 << ", z = " << params->position().z()<< ", y = " << params->position().y()<< ", q/p(TeV) = " << params->charge()/(params->momentum().mag()/1e6));
+            }
+            upstreamOverlay.push_back(upstreamParameters);
+            m_x2.push_back(upstreamParameters->position().x());
+            m_y2.push_back(upstreamParameters->position().y());
+            m_p2.push_back(upstreamParameters->momentum().mag());
+            m_charge2.push_back(upstreamParameters->charge());
+            m_theta2.push_back(asin(upstreamParameters->momentum().perp()/upstreamParameters->momentum().mag()));
+            m_phi2.push_back(atan2(upstreamParameters->momentum().y(), upstreamParameters->momentum().x()));
+            m_chiSquared2.push_back(track->fitQuality()->chiSquared());
+            m_nDoF2.push_back(track->fitQuality()->doubleNumberDoF());
+            pTot += upstreamParameters->momentum();
+        }
+        m_pTotPolar2 = asin(pTot.perp()/pTot.mag());
+    }
+
+    if (trackCollection1.isValid() && trackCollection2.isValid() && m_longTracks2 ==2)
+    {
+        double qOverPSingle1 = upstreamSingle[0]->charge()/upstreamSingle[0]->momentum().mag();
+        double qOverPCovSingle1 = (*upstreamSingle[0]->covariance())(4,4);
+        double qOverPSingle2 = upstreamSingle[1]->charge()/upstreamSingle[1]->momentum().mag();
+        double qOverPCovSingle2 = (*upstreamSingle[1]->covariance())(4,4);
+        double qOverPOverlay1 = upstreamOverlay[0]->charge()/upstreamOverlay[0]->momentum().mag();
+        double qOverPCovOverlay1 = (*upstreamOverlay[0]->covariance())(4,4);
+        double qOverPOverlay2 = upstreamOverlay[1]->charge()/upstreamOverlay[1]->momentum().mag();
+        double qOverPCovOverlay2 = (*upstreamOverlay[1]->covariance())(4,4);
+
+        double agreementDirect = pow(qOverPSingle1 - qOverPOverlay1, 2) / (qOverPCovSingle1 + qOverPCovOverlay1) +
+                                 pow(qOverPSingle2 - qOverPOverlay2, 2) / (qOverPCovSingle2 + qOverPCovOverlay2);
+        double agreementSwap   = pow(qOverPSingle1 - qOverPOverlay2, 2) / (qOverPCovSingle1 + qOverPCovOverlay2) +
+                                 pow(qOverPSingle2 - qOverPOverlay1, 2) / (qOverPCovSingle2 + qOverPCovOverlay1);
+
+        const Trk::Track* firstSingle {nullptr};
+        const Trk::Track* firstOverlay {nullptr};
+        const Trk::Track* secondSingle {nullptr};
+        const Trk::Track* secondOverlay {nullptr};
+        if (agreementDirect <= agreementSwap)
+        {
+            firstSingle = (*trackCollection1)[0];
+            firstOverlay = (*trackCollection2)[0];
+            secondSingle = (*trackCollection1)[1];
+            secondOverlay = (*trackCollection2)[1];
+            ATH_MSG_VERBOSE("Matched q/p = " << qOverPSingle1 << " with " << qOverPOverlay1 << ", and q/p = " << qOverPSingle2 << " with " << qOverPOverlay2);
+        }
+        else
+        {
+            firstSingle = (*trackCollection1)[0];
+            firstOverlay = (*trackCollection2)[1];
+            secondSingle = (*trackCollection1)[1];
+            secondOverlay = (*trackCollection2)[0];
+            ATH_MSG_VERBOSE("Matched q/p = " << qOverPSingle1 << " with " << qOverPOverlay2 << ", and q/p = " << qOverPSingle2 << " with " << qOverPOverlay1);
+        }
+        int nZMatch {0};
+        auto singleParameters = firstSingle->trackParameters();
+        auto overlayParameters = firstOverlay->trackParameters();
+        for (size_t iSingle = 0, iOverlay = 0; (iSingle < singleParameters->size()) && (iOverlay < overlayParameters->size());)
+        {
+            auto singleState = (*singleParameters)[iSingle];
+            auto overlayState = (*overlayParameters)[iOverlay];
+            auto singleCov = singleState->covariance();
+            auto overlayCov = overlayState->covariance();
+            if (abs(singleState->position().z() - overlayState->position().z()) < 0.03)
+            {
+                nZMatch++;
+                m_yAgreement += pow(singleState->position().y() - overlayState->position().y(), 2)/((*singleCov)(1,1) + (*overlayCov)(1,1));
+                iSingle++;
+                iOverlay++;
+                ATH_MSG_VERBOSE("Matched (1): " << singleState->position().z() << " and " << overlayState->position().z());
+            }
+            else if (singleState->position().z() < overlayState->position().z())
+            {
+                ATH_MSG_VERBOSE("z not matched (1): " << singleState->position().z() << " vs " << overlayState->position().z());
+                iSingle++;
+            }
+            else if (overlayState->position().z() < singleState->position().z())
+            {
+                ATH_MSG_VERBOSE("z not matched (1): " << singleState->position().z() << " vs " << overlayState->position().z());
+                iOverlay++;
+            }
+        }
+        double matchFract1 { ((double) nZMatch)/singleParameters->size() };
+        ATH_MSG_VERBOSE("MatchFract1: " << matchFract1);
+        nZMatch = 0;
+        singleParameters = secondSingle->trackParameters();
+        overlayParameters = secondOverlay->trackParameters();
+        for (size_t iSingle = 0, iOverlay = 0; (iSingle < singleParameters->size()) && (iOverlay < overlayParameters->size());)
+        {
+            auto singleState = (*singleParameters)[iSingle];
+            auto overlayState = (*overlayParameters)[iOverlay];
+            auto singleCov = singleState->covariance();
+            auto overlayCov = overlayState->covariance();
+            if (abs(singleState->position().z() - overlayState->position().z()) < 0.03)
+            {
+                nZMatch++;
+                m_yAgreement += pow(singleState->position().y() - overlayState->position().y(), 2)/((*singleCov)(1,1) + (*overlayCov)(1,1));
+                iSingle++;
+                iOverlay++;
+                ATH_MSG_VERBOSE("Matched (2): " << singleState->position().z() << " and " << overlayState->position().z());
+            }
+            else if (singleState->position().z() < overlayState->position().z())
+            {
+                ATH_MSG_VERBOSE("z not matched (2): " << singleState->position().z() << " vs " << overlayState->position().z());
+                iSingle++;
+            }
+            else if (overlayState->position().z() < singleState->position().z())
+            {
+                ATH_MSG_VERBOSE("z not matched (2): " << singleState->position().z() << " vs " << overlayState->position().z());
+                iOverlay++;
+            }
+        }
+        double matchFract2 { ((double) nZMatch)/singleParameters->size() };
+        ATH_MSG_VERBOSE("MatchFract2: " << matchFract2);
+        m_matchFraction = std::min(matchFract1, matchFract2);
+        ATH_MSG_VERBOSE("MatchFraction: " << m_matchFraction);
+    }
+
+    m_tree->Fill();
+    return StatusCode::SUCCESS;
+}
+
+StatusCode HistoRDOAlg::finalize()
+{
+    return StatusCode::SUCCESS;
+}
\ No newline at end of file
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..4fe18f7ad71e075b8b37b56518058da6b7be6f0e
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h
@@ -0,0 +1,58 @@
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "AthenaBaseComps/AthHistogramming.h"
+#include "TrkTrack/TrackCollection.h"
+
+#include <vector>
+
+class TTree;
+
+class HistoRDOAlg : public AthReentrantAlgorithm, AthHistogramming {
+public:
+  HistoRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+  virtual ~HistoRDOAlg() = default;
+  virtual StatusCode initialize() override;
+  virtual StatusCode execute(const EventContext &ctx) const override;
+  virtual StatusCode finalize() override;
+  const ServiceHandle <ITHistSvc> &histSvc() const;
+
+private:
+
+  ServiceHandle <ITHistSvc> m_histSvc;
+
+  SG::ReadHandleKey<TrackCollection> m_trackCollection1 { this, "TrackCollection1", "Pos_CKFTrackCollectionWithoutIFT", "Input track collection 1 name" };
+  SG::ReadHandleKey<TrackCollection> m_trackCollection2 { this, "TrackCollection2", "Neg_CKFTrackCollectionWithoutIFT", "Input track collection 2 name" };
+
+//   mutable TTree* m_tree;
+
+  mutable TTree* m_tree;
+
+  mutable int m_longTracks1;
+  mutable std::vector<double> m_p1;
+  mutable std::vector<double> m_theta1;
+  mutable std::vector<double> m_phi1;
+  mutable std::vector<double> m_x1;
+  mutable std::vector<double> m_y1;
+  mutable std::vector<double> m_charge1;
+  mutable std::vector<double> m_chiSquared1;
+  mutable std::vector<double> m_nDoF1;
+  mutable double m_pTotPolar1;
+
+  mutable int m_longTracks2;
+  mutable std::vector<double> m_p2;
+  mutable std::vector<double> m_theta2;
+  mutable std::vector<double> m_phi2;
+  mutable std::vector<double> m_x2;
+  mutable std::vector<double> m_y2;
+  mutable std::vector<double> m_charge2;
+  mutable std::vector<double> m_chiSquared2;
+  mutable std::vector<double> m_nDoF2;
+  mutable double m_pTotPolar2;
+
+  mutable double m_yAgreement;
+  mutable double m_matchFraction;
+
+};
+
+inline const ServiceHandle <ITHistSvc> &HistoRDOAlg::histSvc() const {
+  return m_histSvc;
+}
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9b88a96f9590462d5c03d6ffba601598185ec109
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
@@ -0,0 +1,186 @@
+#include "OverlayRDOAlg.h"
+#include "TrkTrack/Track.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "Identifier/Identifier.h"
+#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
+#include "TrackerReadoutGeometry/SiDetectorElement.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include <cmath>
+
+
+
+OverlayRDOAlg::OverlayRDOAlg(const std::string &name, 
+                                    ISvcLocator *pSvcLocator)
+    : AthReentrantAlgorithm(name, pSvcLocator)
+    {}
+
+StatusCode OverlayRDOAlg::initialize() 
+{
+  ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
+  ATH_CHECK(m_posTrackKey.initialize());
+  ATH_CHECK(m_negTrackKey.initialize());
+  ATH_CHECK(m_posRdoKey.initialize());
+  ATH_CHECK(m_negRdoKey.initialize());
+  ATH_CHECK(m_posEdgeModeRdoKey.initialize());
+  ATH_CHECK(m_negEdgeModeRdoKey.initialize());
+  ATH_CHECK(m_outRdoKey.initialize());
+  ATH_CHECK(m_outEdgeModeRdoKey.initialize());
+  ATH_CHECK(m_outputTrackCollection.initialize());
+
+  return StatusCode::SUCCESS;
+}
+
+
+StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const 
+{
+
+  SG::ReadHandle<TrackCollection> posTrackCollection {m_posTrackKey, ctx};
+  ATH_CHECK(posTrackCollection.isValid());
+  SG::ReadHandle<TrackCollection> negTrackCollection {m_negTrackKey, ctx};
+  ATH_CHECK(negTrackCollection.isValid());
+  SG::ReadHandle<FaserSCT_RDO_Container> posRdoContainer {m_posRdoKey, ctx};
+  ATH_CHECK(posRdoContainer.isValid());
+  SG::ReadHandle<FaserSCT_RDO_Container> negRdoContainer {m_negRdoKey, ctx};
+  ATH_CHECK(negRdoContainer.isValid());
+  SG::ReadHandle<FaserSCT_RDO_Container> posEdgeModeRdoContainer {m_posEdgeModeRdoKey, ctx};
+  ATH_CHECK(posEdgeModeRdoContainer.isValid());
+  SG::ReadHandle<FaserSCT_RDO_Container> negEdgeModeRdoContainer {m_negEdgeModeRdoKey, ctx};
+  ATH_CHECK(negEdgeModeRdoContainer.isValid());
+
+  // SCT_RDOs
+
+  std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> > rdoDB {};
+  processContainer(*posRdoContainer, rdoDB);
+  processContainer(*negRdoContainer, rdoDB);
+
+  SG::WriteHandle<FaserSCT_RDO_Container> outRdoContainer {m_outRdoKey, ctx};
+  ATH_CHECK(outRdoContainer.record(std::make_unique<FaserSCT_RDO_Container>(m_sctHelper->wafer_hash_max())));
+
+  processDB(rdoDB, outRdoContainer, [](int i) {return true;});
+
+  // SCT_EDGEMODE_RDOs
+
+  rdoDB.clear();
+  processContainer(*posEdgeModeRdoContainer, rdoDB);
+  processContainer(*negEdgeModeRdoContainer, rdoDB);
+
+  SG::WriteHandle<FaserSCT_RDO_Container> outEdgeModeRdoContainer {m_outEdgeModeRdoKey, ctx};
+  ATH_CHECK(outEdgeModeRdoContainer.record(std::make_unique<FaserSCT_RDO_Container>(m_sctHelper->wafer_hash_max())));
+
+  processDB(rdoDB, outEdgeModeRdoContainer, [](int hitPattern) { return (((hitPattern & 0x2) == 0 ) || ((hitPattern & 0x4) != 0) ) ? false : true;});
+
+  // Tracks
+
+  SG::WriteHandle<TrackCollection> outputTrackCollection {m_outputTrackCollection, ctx};
+  std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
+
+  for (auto theTrack : *posTrackCollection)
+  {
+    outputTracks->push_back(cloneTrack(theTrack) );
+  }
+
+  for (auto theTrack : *negTrackCollection)
+  {
+    outputTracks->push_back(cloneTrack(theTrack) );
+  }
+
+  ATH_CHECK(outputTrackCollection.record(std::move(outputTracks)));
+
+  return StatusCode::SUCCESS;
+}
+
+
+void 
+OverlayRDOAlg::processDB(const std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB, SG::WriteHandle<FaserSCT_RDO_Container>& outRdoContainer, std::function<bool(int)> lambda) const
+{
+  for (auto& entry : rdoDB)
+  {
+    IdentifierHash waferHash = entry.first;
+    Identifier id = m_sctHelper->wafer_id(waferHash);
+    auto vec = entry.second;
+    if (vec.size() == 0)
+    {
+      std::cout << "Unexpected zero-length collection vector for ID = " << id << " ..." << std::endl;
+      continue;
+    }
+    std::unique_ptr<FaserSCT_RDO_Collection> current_collection = std::make_unique<FaserSCT_RDO_Collection>(waferHash);
+    current_collection->setIdentifier(id);
+    std::map<int, unsigned int> stripMap;
+    for (auto& collection : vec)
+    {
+      for (auto rawData : *collection)
+      {
+        // const FaserSCT3_RawData* data = dynamic_cast<const FaserSCT3_RawData*>(rawData);
+        Identifier stripID = rawData->identify();
+        int stripNumber = m_sctHelper->strip(stripID);
+        int groupSize = rawData->getGroupSize();
+        unsigned int dataWord = rawData->getWord();
+        // int time = data->getTimeBin();
+        // std::cout << "Word: " << std::hex << dataWord << std::dec << " Time: " << time << std::endl;
+	for (int i = stripNumber; i < stripNumber + groupSize; i++)
+	  stripMap[i] |= dataWord;
+      }
+    }
+
+    for (auto stripEntry : stripMap)
+    {
+      Identifier rdoID {m_sctHelper->strip_id(id, stripEntry.first)};
+      if (!lambda(stripEntry.second>>22)) continue;
+      current_collection->emplace_back(new FaserSCT3_RawData(rdoID, stripEntry.second, std::vector<int> {}));
+    }
+
+    outRdoContainer->getWriteHandle(waferHash).addOrDelete(std::move(current_collection));
+  }
+}
+
+
+Trk::Track* OverlayRDOAlg::cloneTrack(const Trk::Track* theTrack) const
+{
+  Trk::TrackInfo i { theTrack->info() };
+  Trk::FitQuality* q = new Trk::FitQuality { *(theTrack->fitQuality()) };
+  DataVector<const Trk::TrackStateOnSurface>* s = new DataVector<const Trk::TrackStateOnSurface> {};  
+
+  for (auto oldState : (*theTrack->trackStateOnSurfaces()))
+  {
+    const Trk::CurvilinearParameters* oldParam = dynamic_cast<const Trk::CurvilinearParameters*>(oldState->trackParameters());
+    std::unique_ptr<Trk::TrackParameters> newParam { new Trk::CurvilinearParameters { *oldParam } };
+    Trk::TrackStateOnSurface* newState = new Trk::TrackStateOnSurface { nullptr, newParam.release() };
+    s->push_back(newState);
+  }
+  std::unique_ptr<DataVector<const Trk::TrackStateOnSurface>> sink {s};
+  return new Trk::Track {i, std::move(*sink) , q } ;
+}
+
+void OverlayRDOAlg::processContainer(const FaserSCT_RDO_Container& container, std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB) const
+{
+  size_t nLong{0};
+  size_t nData{0};
+  for( const auto& collection : container)
+  {
+      if (collection->empty()) continue;
+      Identifier id = collection->identify();
+      IdentifierHash hash = m_sctHelper->wafer_hash(id);
+      if (rdoDB.count(hash))
+      {
+        rdoDB[hash].push_back(collection);
+      }
+      else
+      {
+        rdoDB[hash] = std::vector<const FaserSCT_RDO_Collection*>{collection};
+      }
+      for(const auto& rawdata : *collection)
+      {
+        nData++;
+        if (rawdata->getGroupSize() > 1) nLong++;
+      }
+  }
+  std::cout << nLong << " long data out of " << nData << " total RDOs" << std::endl;
+}
+
+
+StatusCode OverlayRDOAlg::finalize() 
+{
+  return StatusCode::SUCCESS;
+}
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..eba2f90d7cde46873bbe9124e17876fef5059e5a
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
@@ -0,0 +1,38 @@
+#pragma once
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "TrkTrack/TrackCollection.h"
+#include "TrackerRawData/FaserSCT_RDO_Container.h"
+
+class TTree;
+class FaserSCT_ID;
+
+class OverlayRDOAlg : public AthReentrantAlgorithm {
+public:
+  OverlayRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+  virtual ~OverlayRDOAlg() = default;
+  virtual StatusCode initialize() override;
+  virtual StatusCode execute(const EventContext &ctx) const override;
+  virtual StatusCode finalize() override;
+
+private:
+
+  void processContainer(const FaserSCT_RDO_Container& container, std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB) const;
+  void processDB(const std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& db, SG::WriteHandle<FaserSCT_RDO_Container>& outContainer, std::function<bool(int)> lambda) const;
+  Trk::Track* cloneTrack(const Trk::Track* originalTrack) const;
+
+  SG::ReadHandleKey<TrackCollection> m_posTrackKey { this, "PosTrackCollection", "Pos_CKFTrackCollectionWithoutIFT", "Input track collection name" };
+  SG::ReadHandleKey<TrackCollection> m_negTrackKey { this, "NegTrackCollection", "Neg_CKFTrackCollectionWithoutIFT", "Input track collection name" };
+  SG::ReadHandleKey<FaserSCT_RDO_Container> m_posRdoKey { this, "PosRdoContainer", "Pos_SCT_RDOs"};
+  SG::ReadHandleKey<FaserSCT_RDO_Container> m_negRdoKey { this, "NegRdoContainer", "Neg_SCT_RDOs"};
+  SG::ReadHandleKey<FaserSCT_RDO_Container> m_posEdgeModeRdoKey { this, "PosEdgeModeRdoContainer", "Pos_SCT_EDGEMODE_RDOs"};
+  SG::ReadHandleKey<FaserSCT_RDO_Container> m_negEdgeModeRdoKey { this, "NegEdgeModeRdoContainer", "Neg_SCT_EDGEMODE_RDOs"};
+
+  SG::WriteHandleKey<FaserSCT_RDO_Container> m_outRdoKey{this, "OutputRDOObjectName", "SCT_RDOs", "Output RDO Object name"};
+  SG::WriteHandleKey<FaserSCT_RDO_Container> m_outEdgeModeRdoKey{this, "OutputEDGEMODEObjectName", "SCT_EDGEMODE_RDOs", "Output EDGEMODE Object name"};
+  SG::WriteHandleKey<TrackCollection> m_outputTrackCollection { this, "OutputTrackCollection", "Orig_CKFTrackCollectionWithoutIFT", "Output track collection name"};
+
+  const FaserSCT_ID* m_sctHelper;
+
+};
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b33bed9d18ac4d51f2e2dc7ae9fc52217cbffc67
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
@@ -0,0 +1,187 @@
+#include "SelectRDOAlg.h"
+#include "TrkTrack/Track.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterCollection.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "Identifier/Identifier.h"
+#include <cmath>
+
+SelectRDOAlg::SelectRDOAlg(const std::string &name, 
+                           ISvcLocator *pSvcLocator)
+    : AthReentrantAlgorithm(name, pSvcLocator)
+{}
+
+
+StatusCode SelectRDOAlg::initialize() 
+{
+  ATH_CHECK(m_trackCollection.initialize());
+  ATH_CHECK(m_segmentCollection.initialize());
+  ATH_CHECK(m_clusterContainer.initialize());
+  ATH_CHECK(m_triggerContainer.initialize());
+  ATH_CHECK(m_outputTrackCollection.initialize());
+
+  ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
+
+  return StatusCode::SUCCESS;
+}
+
+
+StatusCode SelectRDOAlg::execute(const EventContext &ctx) const 
+{
+
+  setFilterPassed(false, ctx);
+  SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx};
+  if (!trackCollection.isValid() || trackCollection->size() == 0) return StatusCode::SUCCESS;
+
+  SG::ReadHandle<TrackCollection> segmentCollection {m_segmentCollection, ctx};
+  if (!segmentCollection.isValid()) return StatusCode::SUCCESS;
+
+  SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> clusterContainer { m_clusterContainer, ctx };
+  if (!clusterContainer.isValid()) return StatusCode::SUCCESS;
+
+  SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
+  if (!triggerContainer.isValid()) return StatusCode::SUCCESS;
+
+  // segment cuts
+
+  std::vector<int> segmentCount {0, 0, 0, 0};
+  int nSegments {0};
+
+  for (const Trk::Track* segment : *segmentCollection)
+  {
+    if (segment == nullptr) continue;
+    nSegments++;
+    for (auto measurement : *(segment->measurementsOnTrack()))
+    {
+        const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
+        if (cluster != nullptr)
+        {
+            Identifier id = cluster->identify();
+            int station = m_sctHelper->station(id);
+            if (station >= 0 && station <= 3)
+            {
+              segmentCount[station]++;
+              break;
+            }
+        }
+    }
+  }
+  if ((nSegments > m_maxSegmentsTotal) || (segmentCount[1] > m_maxSegmentsStation) || (segmentCount[2] > m_maxSegmentsStation) || (segmentCount[3] > m_maxSegmentsStation))
+    return StatusCode::SUCCESS;
+
+  // track cuts
+
+  int nTracks {0};
+  int nLongTracks {0};
+  double chi2PerDoF {0};
+  double charge {0};
+  double maxRadius {0};
+  double momentum {0};
+  const Trk::Track* theTrack {nullptr};
+
+  for (const Trk::Track* track : *trackCollection)
+  {
+    if (track == nullptr) continue;
+    nTracks++;
+    std::set<int> stationMap;
+    std::set<std::pair<int, int>> layerMap;
+
+    // Check for hit in the three downstream stations
+    for (auto measurement : *(track->measurementsOnTrack()))
+    {
+        const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
+        if (cluster != nullptr)
+        {
+            Identifier id = cluster->identify();
+            int station = m_sctHelper->station(id);
+            int layer = m_sctHelper->layer(id);
+            stationMap.emplace(station);
+            layerMap.emplace(station, layer);
+        }
+    }
+    if (stationMap.count(1) == 0 || stationMap.count(2) == 0 || stationMap.count(3) == 0) continue;
+
+    int nLayers = std::count_if(layerMap.begin(), layerMap.end(), [](std::pair<int,int> p){return p.first != 0;});
+    if (nLayers < m_minLayers) continue;
+    nLongTracks++;
+    const Trk::TrackParameters* upstreamParameters {nullptr};
+    for (auto params : *(track->trackParameters()))
+    {
+        if (params->position().z() < 0) continue;  // Ignore IFT hits
+        double radius = sqrt(pow(params->position().x(), 2) + pow(params->position().y(), 2));
+        if (radius > maxRadius) maxRadius = radius;
+        if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+    }
+
+    momentum = upstreamParameters->momentum().mag()/1000;
+    chi2PerDoF = track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF();
+    charge = upstreamParameters->charge();
+    theTrack = track;
+  }
+
+  if ((nTracks != 1) || (nLongTracks != 1) || (momentum < m_minMomentum) || (momentum > m_maxMomentum) || (chi2PerDoF > m_maxChi2PerDoF) || ((charge>0) != m_acceptPositive) || (maxRadius > m_maxRadius))
+    return StatusCode::SUCCESS;
+
+  // cluster cuts
+
+  // std::vector<int> clusterCount {0, 0, 0, 0};
+  // for (auto collection : *clusterContainer)
+  // {
+  //   Identifier id = collection->identify();
+  //   int station = m_sctHelper->station(id);
+  //   if (station >= 0 && station <= 3) clusterCount[station] += collection->size();
+  // }
+  // if ((clusterCount[1] > m_maxClustersStation) || (clusterCount[2] > m_maxClustersStation) || (clusterCount[3] > m_maxClustersStation))
+  //   return StatusCode::SUCCESS;
+
+  // waveform charge cuts
+
+  std::vector<double> timingCharge {0, 0, 0, 0};
+  const int baseChannel = 8;
+
+  for (auto hit : *triggerContainer) {
+    if ((hit->hit_status()&2)==0) { // dont store secondary hits as they can overwrite the primary hit
+      int ch=hit->channel();
+      timingCharge[ch - baseChannel] += hit->integral()/50;
+    }
+  }
+  if ((timingCharge[0] + timingCharge[1] > m_maxTimingCharge) || (timingCharge[2] + timingCharge[3] > m_maxTimingCharge))
+    return StatusCode::SUCCESS;
+
+  // passes all cuts
+
+  // copy track without associated clusters
+  
+  SG::WriteHandle<TrackCollection> outputTrackCollection {m_outputTrackCollection, ctx};
+  std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
+
+  Trk::TrackInfo i { theTrack->info() };
+  Trk::FitQuality* q = new Trk::FitQuality { *(theTrack->fitQuality()) };
+  DataVector<const Trk::TrackStateOnSurface>* s = new DataVector<const Trk::TrackStateOnSurface> {};  
+
+  for (auto oldState : (*theTrack->trackStateOnSurfaces()))
+  {
+    const Trk::CurvilinearParameters* oldParam = dynamic_cast<const Trk::CurvilinearParameters*>(oldState->trackParameters());
+    std::unique_ptr<Trk::TrackParameters> newParam { new Trk::CurvilinearParameters { *oldParam } };
+    Trk::TrackStateOnSurface* newState = new Trk::TrackStateOnSurface { nullptr, newParam.release() };
+    s->push_back(newState);
+  }
+  std::unique_ptr<DataVector<const Trk::TrackStateOnSurface>> sink {s};
+  outputTracks->push_back(new Trk::Track {i, std::move(*sink) , q } );
+
+  ATH_CHECK(outputTrackCollection.record(std::move(outputTracks)));
+
+  // set filter passed
+
+  setFilterPassed(true, ctx);
+  return StatusCode::SUCCESS;
+}
+
+
+StatusCode SelectRDOAlg::finalize() 
+{
+  return StatusCode::SUCCESS;
+}
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..842282d534a23159625a1dabf59c3fbc0494c343
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
@@ -0,0 +1,43 @@
+#pragma once
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "TrkTrack/TrackCollection.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "xAODFaserWaveform/WaveformHitContainer.h"
+#include "xAODFaserWaveform/WaveformHit.h"
+#include <xAODEventInfo/EventInfo.h>
+
+
+class FaserSCT_ID;
+
+class SelectRDOAlg : public AthReentrantAlgorithm {
+public:
+  SelectRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+  virtual ~SelectRDOAlg() = default;
+  virtual StatusCode initialize() override;
+  virtual StatusCode execute(const EventContext &ctx) const override;
+  virtual StatusCode finalize() override;
+
+private:
+
+  const FaserSCT_ID* m_sctHelper;
+
+  SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "CKFTrackCollectionWithoutIFT", "Input track collection name" };
+  SG::ReadHandleKey<TrackCollection> m_segmentCollection { this, "SegmentCollection", "SegmentFit", "Track segment collection name" };
+  SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_clusterContainer { this, "ClusterContainer", "SCT_ClusterContainer" };
+  SG::ReadHandleKey<xAOD::WaveformHitContainer> m_triggerContainer { this, "TriggerContainer", "TriggerWaveformHits", "Trigger hit container name" };
+  SG::WriteHandleKey<TrackCollection> m_outputTrackCollection { this, "OutputTrackCollection", "Pos_CKFTrackCollectionWithoutIFT", "Output track collection name"};
+
+  IntegerProperty   m_minLayers {this, "MinLayers", 7, "Minimum hit layers to accept track"};
+  DoubleProperty    m_maxChi2PerDoF {this, "MaxChi2PerDoF", 25.0, "Maximum chi2 per degree of freedom to accept track"};
+  DoubleProperty    m_minMomentum {this, "MinMomentumGeV", 50.0, "Minimum momentum in GeV to accept track"};
+  DoubleProperty    m_maxMomentum {this, "MaxMomentumGeV", 5000.0, "Maximum momentum in GeV to accept track"};
+  DoubleProperty    m_maxRadius {this, "MaxRadiusMm", 95.0, "Maximum radius at first measurement to accept track"};
+  IntegerProperty   m_maxSegmentsTotal {this, "MaxSegmentsTotal", 4, "Maximum number of segments in three stations to accept track"};
+  IntegerProperty   m_maxSegmentsStation {this, "MaxSegmentsStation", 2, "Maximum number of segments in any single station to accept track"};
+  // IntegerProperty   m_maxClustersStation {this, "MaxClustersStation", 9, "Maximum number of clusters in any single station to accept track"};
+  DoubleProperty    m_maxTimingCharge {this, "MaxTimingChargePc", 70.0, "Maximum charge in pC recorded by upper or lower timing scintillator"};
+  BooleanProperty   m_acceptPositive {this, "AcceptPositive", true, "Accept positive (true) or negative (false) charged tracks"};
+
+};
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..13e61cb84654d0c3e4096ffa01b0f6efbeb6184c
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
@@ -0,0 +1,9 @@
+#include "../OverlayRDOAlg.h"
+#include "../SelectRDOAlg.h"
+#include "../BinRDOAlg.h"
+#include "../HistoRDOAlg.h"
+
+DECLARE_COMPONENT(OverlayRDOAlg)
+DECLARE_COMPONENT(SelectRDOAlg)
+DECLARE_COMPONENT(BinRDOAlg)
+DECLARE_COMPONENT(HistoRDOAlg)
diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py b/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
index f0de1582f2ccf47270203d4d1d6eb9a6bca341e9..eeaadd061bbc5b1c3c14af0a1e128e529d687cae 100644
--- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
+++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
@@ -15,9 +15,9 @@ def SegmentFitAlgBasicCfg(flags, **kwargs):
     Tracker__SegmentFitAlg=CompFactory.Tracker.SegmentFitAlg
     acc.addEventAlgo(Tracker__SegmentFitAlg(**kwargs))
     
-    thistSvc = CompFactory.THistSvc()
-    thistSvc.Output += ["HIST DATAFILE='SegmentFitHistograms.root' OPT='RECREATE'"]
-    acc.addService(thistSvc)
+    # thistSvc = CompFactory.THistSvc()
+    # thistSvc.Output += ["HIST DATAFILE='SegmentFitHistograms.root' OPT='RECREATE'"]
+    # acc.addService(thistSvc)
     return acc
 
 # with output defaults
diff --git a/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx b/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx
index 814fc31a5b82c9705472a283ba8abfd074ed067b..c8a7df00173d6503683ac47148baefd8bf54da4c 100644
--- a/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx
+++ b/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx
@@ -122,7 +122,7 @@ FaserSCT_SpacePoint* TrackerSpacePointMakerTool::makeSCT_SpacePoint(const Tracke
        {
          if (fabs(lambda1) > 1 + m_stripLengthTolerance)
          {
-           ATH_MSG_WARNING("Intersection lies outside the bounds of both strips");
+           ATH_MSG_DEBUG("Intersection lies outside the bounds of both strips");
            ok = false;
          }
        }
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index f990be9458e0ee2ac0a5bdef8b670703f1424c26..14616c6f42dd0e60df67c5b73358da599093130c 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -50,6 +50,7 @@ StatusCode CKF2::initialize() {
   ATH_CHECK(m_kalmanFitterTool1.retrieve());
   ATH_CHECK(m_createTrkTrackTool.retrieve());
   ATH_CHECK(m_trackCollection.initialize());
+  // ATH_CHECK(m_allTrackCollection.initialize());
   ATH_CHECK(m_eventInfoKey.initialize());
   if (m_performanceWriter && !m_noDiagnostics) {
     ATH_CHECK(m_performanceWriterTool.retrieve());
@@ -85,6 +86,9 @@ StatusCode CKF2::execute() {
   SG::WriteHandle trackContainer{m_trackCollection, ctx};
   std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
 
+  // SG::WriteHandle allTrackContainer{m_allTrackCollection, ctx};
+  // std::unique_ptr<TrackCollection> outputAllTracks = std::make_unique<TrackCollection>();
+
   std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry
       = m_trackingGeometryTool->trackingGeometry();
 
@@ -157,10 +161,19 @@ StatusCode CKF2::execute() {
       // TODO use status bits for different errors
       // result.error() == Acts::CombinatorialKalmanFilterError::NoTrackFound
       if (result.error() == Acts::PropagatorError::StepCountLimitReached ||
-          result.error() == Acts::CombinatorialKalmanFilterError::PropagationReachesMaxSteps) {
-        if (!eventInfo->updateErrorState(xAOD::EventInfo::SCT, xAOD::EventInfo::Error)) {
-          ATH_MSG_WARNING (" cannot set error state for SCT");
-        }
+          result.error() == Acts::CombinatorialKalmanFilterError::PropagationReachesMaxSteps) 
+      {
+          try
+          {
+            if (!eventInfo->updateErrorState(xAOD::EventInfo::SCT, xAOD::EventInfo::Error)) 
+            {
+              ATH_MSG_WARNING ("Cannot set error state for SCT.");
+            }
+          }
+          catch (...)
+          {
+              ATH_MSG_DEBUG ("SCT error state is locked.");
+          }
       }
       continue;
     }
@@ -183,6 +196,33 @@ StatusCode CKF2::execute() {
   // select all tracks with at least 13 heats and with 6 or less shared hits, starting from the best track
   // TODO use Gaudi parameters for the number of hits and shared hits
   // TODO allow shared hits only in the first station?
+  // std::vector<FaserActsRecMultiTrajectory> rawTrajectories {};
+  // for (auto raw : allTrajectories)
+  // {
+  //   rawTrajectories.push_back(raw.trajectory);
+  // }
+  // for (const FaserActsRecMultiTrajectory &traj : rawTrajectories) {
+  //   const auto params = traj.trackParameters(traj.tips().front());
+  //   ATH_MSG_DEBUG("Fitted parameters (raw)");
+  //   ATH_MSG_DEBUG("  params:   " << params.parameters().transpose());
+  //   ATH_MSG_DEBUG("  position: " << params.position(gctx).transpose());
+  //   ATH_MSG_DEBUG("  momentum: " << params.momentum().transpose());
+  //   ATH_MSG_DEBUG("  charge:   " << params.charge());
+  //   std::unique_ptr<Trk::Track> track = m_createTrkTrackTool->createTrack(gctx, traj);
+  //   if (track != nullptr) {
+  //     m_numberOfSelectedTracks++;
+  //     std::unique_ptr<Trk::Track> track2 = m_kalmanFitterTool1->fit(ctx, gctx, track.get(), Acts::BoundVector::Zero(), m_isMC, origin);
+  //     if (track2) {
+  //       outputAllTracks->push_back(std::move(track2));
+  //     } else {
+  //       outputAllTracks->push_back(std::move(track));
+  //       ATH_MSG_WARNING("Re-Fit failed.");
+  //     }
+  //   } else {
+  //     ATH_MSG_WARNING("CKF failed.");
+  //   }
+  // }
+
   std::vector<FaserActsRecMultiTrajectory> selectedTrajectories {};
   while (not allTrajectories.empty()) {
     TrajectoryInfo selected = allTrajectories.front();
@@ -224,6 +264,7 @@ StatusCode CKF2::execute() {
   if  (m_performanceWriter && !m_noDiagnostics) {
     ATH_CHECK(m_performanceWriterTool->write(gctx, selectedTrajectories));
   }
+  // ATH_CHECK(allTrackContainer.record(std::move(outputAllTracks)));
   ATH_CHECK(trackContainer.record(std::move(outputTracks)));
 
   return StatusCode::SUCCESS;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
index 170f0353904bef82095146c50bf98d5721ef2138..cd78bf72531931fc3ef7ef3b1b67edc854e46596 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
@@ -149,6 +149,7 @@ private:
   ToolHandle<CreateTrkTrackTool> m_createTrkTrackTool {this, "CreateTrkTrackTool", "CreateTrkTrackTool"};
   Gaudi::Property<bool> m_isMC {this, "isMC", false};
   SG::WriteHandleKey<TrackCollection> m_trackCollection { this, "OutputCollection", "CKFTrackCollection", "Output track collection name" };
+  // SG::WriteHandleKey<TrackCollection> m_allTrackCollection { this, "AllTrackCollection", "CKFAllTrackCollection", "Output all track collection name" };
   SG::WriteDecorHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", "EventInfo"};
 };
 
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
index 952f75287f594e5df31feb2b6a493bfe235a4086..bbae70e3b10be38ff5de9675796bbe509c8f2483 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
@@ -130,10 +130,16 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
   std::list<Seed> allSeeds;
   for (const Seed &seed : seeds) allSeeds.push_back(seed);
 
+  // allSeeds.sort([](const Seed &left, const Seed &right) {
+  //   if (left.size > right.size) return true;
+  //   if (left.size < right.size) return false;
+  //   if (left.chi2 < right.chi2) return true;
+  //   else return false;
+  // });
   allSeeds.sort([](const Seed &left, const Seed &right) {
-    if (left.size > right.size) return true;
-    if (left.size < right.size) return false;
-    if (left.chi2 < right.chi2) return true;
+    if (left.stations > right.stations) return true;
+    if (left.stations < right.stations) return false;
+    if (left.chi2/std::max(left.positions.size() + left.fakePositions.size() - left.constraints, 1UL) < right.chi2/std::max(right.positions.size() + right.fakePositions.size() - right.constraints, 1UL)) return true;
     else return false;
   });
 
@@ -141,6 +147,9 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
   while (not allSeeds.empty()) {
     Seed selected = allSeeds.front();
     selectedSeeds.push_back(selected);
+    // allSeeds.remove_if([&](const Seed &p) {
+    //   return ((p.clusterSet & selected.clusterSet).count() == p.clusterSet.count());
+    // });
     allSeeds.remove_if([&](const Seed &p) {
       return (p.size < 10) || ((p.clusterSet & selected.clusterSet).count() > 6);
     });
@@ -161,7 +170,9 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
   double origin = !selectedSeeds.empty() ? minSeed->minZ - 10 : 0;
   m_targetZPosition = origin;
   std::vector<Acts::CurvilinearTrackParameters> initParams {};
+  ATH_MSG_DEBUG("Sorted seed properties:");
   for (const Seed &seed : selectedSeeds) {
+    ATH_MSG_DEBUG("seed size: " << seed.size << ", chi2: " << seed.chi2);
     initParams.push_back(seed.get_params(origin, cov));
   }
 
@@ -264,13 +275,16 @@ CircleFitTrackSeedTool::Seed::Seed(const std::vector<Segment> &segments) :
 
   if (segments.size() > 1) {
     fakeFit();
+    constraints = 5;
   } else {
     momentum = 9999999.;
     charge = 1;
+    constraints = 2;
   }
 
   getChi2();
   size = clusters.size();
+  stations = segments.size();
 }
 
 
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
index 0467c92620eebd955f0a78f11cfdea3fd0cb0b1d..a4d53c629c24f53e1dfe4488ca7524755f117ad3 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
@@ -66,7 +66,7 @@ private:
 
     double c0, c1, cx, cy, r, chi2, momentum, charge, minZ;
     Acts::Vector3 direction;
-    size_t size;
+    size_t size, stations, constraints;
     Acts::CurvilinearTrackParameters get_params(double origin, Acts::BoundSymMatrix cov) const;
 
   private:
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
index 4928fc917e2b0bfe73e38848939690680fc3c19a..994689010203a269946c3c88392004a61a42adf4 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
@@ -61,9 +61,17 @@ CreateTrkTrackTool::createTrack(const Acts::GeometryContext &gctx, const FaserAc
           };
         }
       }
-      const Trk::TrackStateOnSurface *perState = new Trk::TrackStateOnSurface(clusterOnTrack, parm);
-      if (perState) {
-        finalTrajectory->insert(finalTrajectory->begin(), perState);
+      double nDoF = state.calibratedSize();
+      const Trk::FitQualityOnSurface* quality = new Trk::FitQualityOnSurface(state.chi2(), nDoF);
+      const Trk::TrackStateOnSurface* perState = new Trk::TrackStateOnSurface(
+									      clusterOnTrack, 
+									      parm,
+									      quality,
+									      nullptr,
+									      typePattern);
+      if (perState)
+      {
+	        finalTrajectory->insert(finalTrajectory->begin(), perState);
       }
     }
     return;