diff --git a/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py
index 793f73401b54c91724bdb6283fa4c6099811628c..1b8a0355ebe97a514954cd5282b2bef462a0f0c4 100644
--- a/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py
+++ b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py
@@ -9,14 +9,14 @@ from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
 from WaveformConditionsTools.WaveformCableMappingConfig import WaveformCableMappingCfg
 
 # One stop shopping for normal FASER data
-def CaloWaveformDigitizationCfg(flags):
+def CaloWaveformDigitizationCfg(flags, **kwargs):
     """ Return all algorithms and tools for Waveform digitization """
     acc = ComponentAccumulator()
 
     if not flags.Input.isMC:
         return acc
 
-    acc.merge(CaloWaveformDigiCfg(flags, "CaloWaveformDigiAlg"))
+    acc.merge(CaloWaveformDigiCfg(flags, "CaloWaveformDigiAlg", **kwargs))
     acc.merge(CaloWaveformDigitizationOutputCfg(flags))
     acc.merge(WaveformCableMappingCfg(flags))
 
@@ -27,23 +27,24 @@ def CaloWaveformDigiCfg(flags, name="CaloWaveformDigiAlg", **kwargs):
 
     acc = ComponentAccumulator()
 
-    tool = CompFactory.WaveformDigitisationTool(name="CaloWaveformDigtisationTool", **kwargs)
+    tool = CompFactory.WaveformDigitisationTool(name="CaloWaveformDigtisationTool")
+    kwargs.setdefault("WaveformDigitisationTool", tool)
     
     kwargs.setdefault("CaloHitContainerKey", "EcalHits")
     kwargs.setdefault("WaveformContainerKey", "CaloWaveforms")
 
-    digiAlg = CompFactory.CaloWaveformDigiAlg(name, **kwargs)
-    kwargs.setdefault("WaveformDigitisationTool", tool)
-    
-    digiAlg.CB_alpha = -0.9
-    digiAlg.CB_n = 10
-    digiAlg.CB_sigma = 4
-    digiAlg.CB_mean = 820
-    digiAlg.CB_norm = 2
+    kwargs.setdefault("CB_alpha", -0.9)
+    kwargs.setdefault("CB_n", 10)
+    kwargs.setdefault("CB_sigma", 4)
+    kwargs.setdefault("CB_mean", 820) # Time in ns
+    kwargs.setdefault("CB_norm", 4)  # Low gain default, use 20 for high gain
 
-    digiAlg.base_mean = 15000
-    digiAlg.base_rms = 3
+    kwargs.setdefault("base_mean", 15000)
+    kwargs.setdefault("base_rms", 3)
 
+    digiAlg = CompFactory.CaloWaveformDigiAlg(name, **kwargs)
+    print(f"CaloWaveformDigiAlg normalization: {digiAlg.CB_norm}")
+    print(f"CaloWaveformDigiAlg mean time: {digiAlg.CB_mean}")
 
     acc.addEventAlgo(digiAlg)
 
diff --git a/Control/CalypsoExample/Digitization/scripts/faserMDC_digi.py b/Control/CalypsoExample/Digitization/scripts/faserMDC_digi.py
index 0d94cf0c344b4010627ad63826c4d7aa4732a580..14a46818f6c334123f8f62f35212f72c57e25740 100755
--- a/Control/CalypsoExample/Digitization/scripts/faserMDC_digi.py
+++ b/Control/CalypsoExample/Digitization/scripts/faserMDC_digi.py
@@ -23,6 +23,8 @@ parser.add_argument("run_type", nargs="?", default="",
                     help="Specify run type (if it can't be parsed from path)")
 parser.add_argument("-t", "--tag", default="",
                     help="Specify digi tag (to append to output filename)")
+parser.add_argument("--highCaloGain", action='store_true',
+                    help="Use high gain settings for calo PMTs")
 parser.add_argument("-n", "--nevts", type=int, default=-1,
                     help="Specify number of events to process (default: all)")
 parser.add_argument("-v", "--verbose", action='store_true', 
@@ -130,15 +132,19 @@ from FaserSCT_Digitization.FaserSCT_DigitizationConfigNew import FaserSCT_Digiti
 acc.merge(FaserSCT_DigitizationCfg(ConfigFlags))
 
 from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg
-acc.merge(CaloWaveformDigitizationCfg(ConfigFlags))
+if args.highCaloGain:
+    calo_norm = 20.
+else:
+    calo_norm = 4.
+acc.merge(CaloWaveformDigitizationCfg(ConfigFlags, CB_norm=calo_norm))
 
 from ScintDigiAlgs.ScintDigiAlgsConfig import ScintWaveformDigitizationCfg
 acc.merge(ScintWaveformDigitizationCfg(ConfigFlags))
 
 # Configure verbosity    
-# ConfigFlags.dump()
 if args.verbose:
     acc.foreach_component("*").OutputLevel = VERBOSE
+    ConfigFlags.dump()
 
 else:
     acc.foreach_component("*").OutputLevel = INFO
diff --git a/Control/CalypsoExample/Digitization/scripts/faserMDC_digi_merge.py b/Control/CalypsoExample/Digitization/scripts/faserMDC_digi_merge.py
index 9c656de660dc7957aa222b691e50af457a870bf0..cef68cb5668923644c36080cc7ed1e9e1d85950f 100755
--- a/Control/CalypsoExample/Digitization/scripts/faserMDC_digi_merge.py
+++ b/Control/CalypsoExample/Digitization/scripts/faserMDC_digi_merge.py
@@ -27,6 +27,8 @@ parser.add_argument("-f", "--files", type=int, default=5,
                     help="Specify number of input files to run in one batch")
 parser.add_argument("-t", "--tag", default="",
                     help="Specify digi tag (to append to output filename)")
+parser.add_argument("--highCaloGain", action='store_true',
+                    help="Use high gain settings for calo PMTs")
 parser.add_argument("-n", "--nevts", type=int, default=-1,
                     help="Specify number of events to process (default: all)")
 parser.add_argument("-v", "--verbose", action='store_true', 
@@ -190,15 +192,19 @@ from FaserSCT_Digitization.FaserSCT_DigitizationConfigNew import FaserSCT_Digiti
 acc.merge(FaserSCT_DigitizationCfg(ConfigFlags))
 
 from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg
-acc.merge(CaloWaveformDigitizationCfg(ConfigFlags))
+if args.highCaloGain:
+    calo_norm = 20.
+else:
+    calo_norm = 4.
+acc.merge(CaloWaveformDigitizationCfg(ConfigFlags, CB_norm=calo_norm))
 
 from ScintDigiAlgs.ScintDigiAlgsConfig import ScintWaveformDigitizationCfg
 acc.merge(ScintWaveformDigitizationCfg(ConfigFlags))
 
 # Configure verbosity    
-# ConfigFlags.dump()
 if args.verbose:
     acc.foreach_component("*").OutputLevel = VERBOSE
+    ConfigFlags.dump()
 
 else:
     acc.foreach_component("*").OutputLevel = INFO
diff --git a/Control/CalypsoExample/Digitization/scripts/faser_digi.py b/Control/CalypsoExample/Digitization/scripts/faser_digi.py
index 91de9cad754ec1bb739a3f689f27911139e48882..d8cb774850b8b46b94288da30887e3146b7b0c9e 100755
--- a/Control/CalypsoExample/Digitization/scripts/faser_digi.py
+++ b/Control/CalypsoExample/Digitization/scripts/faser_digi.py
@@ -19,8 +19,10 @@ parser.add_argument("file_path",
                     help="Fully qualified path of the raw input file")
 parser.add_argument("run_type", nargs="?", default="",
                     help="Specify run type (if it can't be parsed from path)")
-parser.add_argument("-r", "--reco", default="",
-                    help="Specify reco tag (to append to output filename)")
+parser.add_argument("-t", "--tag", default="",
+                    help="Specify tag (to append to output filename)")
+parser.add_argument("--highCaloGain", action='store_true',
+                    help="Use high gain settings for calo PMTs")
 parser.add_argument("-n", "--nevents", type=int, default=-1,
                     help="Specify number of events to process (default: all)")
 parser.add_argument("-v", "--verbose", action='store_true', 
@@ -95,8 +97,8 @@ else:
 ConfigFlags.Input.Files = [ args.file_path ]
 
 filestem = filepath.stem
-if len(args.reco) > 0:
-    filestem += f"-{args.reco}"
+if len(args.tag) > 0:
+    filestem += f"-{args.tag}"
 
 # ConfigFlags.addFlag("Output.xAODFileName", f"{filestem}-xAOD.root")
 ConfigFlags.Output.RDOFileName = f"{filestem}-RDO.root"
@@ -127,7 +129,11 @@ from FaserSCT_Digitization.FaserSCT_DigitizationConfigNew import FaserSCT_Digiti
 acc.merge(FaserSCT_DigitizationCfg(ConfigFlags))
 
 from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg
-acc.merge(CaloWaveformDigitizationCfg(ConfigFlags))
+if args.highCaloGain:
+    calo_norm = 20.
+else:
+    calo_norm = 4.
+acc.merge(CaloWaveformDigitizationCfg(ConfigFlags, CB_norm=calo_norm))
 
 from ScintDigiAlgs.ScintDigiAlgsConfig import ScintWaveformDigitizationCfg
 acc.merge(ScintWaveformDigitizationCfg(ConfigFlags))
@@ -186,15 +192,9 @@ acc.merge(ScintWaveformDigitizationCfg(ConfigFlags))
 # replicaSvc.UseGeomSQLite = True
 
 # Configure verbosity    
-# ConfigFlags.dump()
 if args.verbose:
     acc.foreach_component("*").OutputLevel = VERBOSE
-
-    #acc.getService("FaserByteStreamInputSvc").DumpFlag = True
-    #acc.getService("FaserEventSelector").OutputLevel = VERBOSE
-    #acc.getService("FaserByteStreamInputSvc").OutputLevel = VERBOSE
-    #acc.getService("FaserByteStreamCnvSvc").OutputLevel = VERBOSE
-    #acc.getService("FaserByteStreamAddressProviderSvc").OutputLevel = VERBOSE
+    ConfigFlags.dump()
 
 else:
     acc.foreach_component("*").OutputLevel = INFO
diff --git a/Control/CalypsoExample/Digitization/scripts/submit_faserMDC_digi.sh b/Control/CalypsoExample/Digitization/scripts/submit_faserMDC_digi.sh
index a8e1b9059aa32d4e91bafe923797d262aa2b84ca..7b1734797cfabcd3a706d389fc8ac396a6ab6248 100755
--- a/Control/CalypsoExample/Digitization/scripts/submit_faserMDC_digi.sh
+++ b/Control/CalypsoExample/Digitization/scripts/submit_faserMDC_digi.sh
@@ -2,7 +2,10 @@
 # Used with a condor file to submit to vanilla universe
 #
 # Usage:
-# submit_faserMDC_digi.sh filepath [release_directory] [working_directory]
+# submit_faserMDC_digi.sh [--highGain] filepath [release_directory] [working_directory]
+#
+# Options:
+#   --highGain - apply high gain settings to the Calorimeter PMTs (for muons)
 # 
 # filepath - full file name (with path)
 # release_directory - optional path to release install directory (default pwd)
@@ -12,13 +15,30 @@
 # (so an unqualified asetup can set up the release properly)
 #
 # Script will use git describe to find the release tag.  
-# If this matches gen/g???? it will be passed to the job
+# If this matches sim/s???? or digi/d???? it will be passed to the job
 #
 #----------------------------------------
 # Keep track of time
 SECONDS=0
 #
 # Parse command-line options
+while [ -n "$1" ]
+do 
+  case "$1" in
+    --highGain) 
+	  echo "Applying high gain settings"
+	  highgain=1 
+	  shift;;  # This 'eats' the argument
+
+    -*) 
+	  echo "Unknown option $1"
+	  shift;;
+
+    *) break;;  # Not an option, don't shift
+  esac
+done
+#
+# Parse command-line options
 file_path=${1}
 release_directory=${2}
 working_directory=${3}
@@ -27,7 +47,7 @@ working_directory=${3}
 if [ -z "$file_path" ]
 then
   echo "No file specified!"
-  echo "Usage: submit_faserMDC_digi.sh file [release dir] [output dir]"
+  echo "Usage: submit_faserMDC_digi.sh [--highGain] file [release dir] [output dir]"
   exit 1
 fi
 #
@@ -124,10 +144,17 @@ fi
 cd "$file_stem"
 #
 # Run job
+#
+if [[ -z "$highgain" ]]; then
+  gainstr=""
+else
+  gainstr="--highCaloGain"
+fi
+# 
 if [[ -z "$tag" ]]; then
-    faserMDC_digi.py "$file_path"
+    faserMDC_digi.py $gainstr "$file_path"
 else
-    faserMDC_digi.py "--tag=$tag" "$file_path"
+    faserMDC_digi.py $gainstr "--tag=$tag" "$file_path"
 fi
 #
 # Print out ending time
diff --git a/Control/CalypsoExample/Digitization/scripts/submit_faserMDC_digi_merge.sh b/Control/CalypsoExample/Digitization/scripts/submit_faserMDC_digi_merge.sh
index 3c500905965fc48e7242838f284226516c03b8f2..7963cda28222ec45e13242e1eaebd96595d62a27 100755
--- a/Control/CalypsoExample/Digitization/scripts/submit_faserMDC_digi_merge.sh
+++ b/Control/CalypsoExample/Digitization/scripts/submit_faserMDC_digi_merge.sh
@@ -2,7 +2,10 @@
 # Used with a condor file to submit to vanilla universe
 #
 # Usage:
-# submit_faserMDC_digi_merge.sh dirpath slice nfiles [release_directory] [working_directory] 
+# submit_faserMDC_digi_merge.sh [--highGain] dirpath slice nfiles [release_directory] [working_directory] 
+#
+# Options:
+#   --highGain - apply high gain settings to the Calorimeter PMTs (for muons)
 # 
 # dirpath - full directory path to HITS files
 # slice - ordinal output file number
@@ -21,6 +24,23 @@
 SECONDS=0
 #
 # Parse command-line options
+while [ -n "$1" ]
+do 
+  case "$1" in
+    --highGain) 
+	  echo "Applying high gain settings"
+	  highgain=1 
+	  shift;;  # This 'eats' the argument
+
+    -*) 
+	  echo "Unknown option $1"
+	  shift;;
+
+    *) break;;  # Not an option, don't shift
+  esac
+done
+#
+# Parse command-line options
 dir_path=${1}
 slice=${2}
 nfiles=${3}
@@ -148,10 +168,16 @@ fi
 cd "$file_stem"
 #
 # Run job
+if [[ -z "$highgain" ]]; then
+  gainstr=""
+else
+  gainstr="--highCaloGain"
+fi
+# 
 if [[ -z "$tag" ]]; then
-    faserMDC_digi_merge.py --slice $slice --files $nfiles $dir_path
+    faserMDC_digi_merge.py $gainstr --slice $slice --files $nfiles $dir_path
 else
-    faserMDC_digi_merge.py --slice $slice --files $nfiles --tag $tag $dir_path
+    faserMDC_digi_merge.py $gainstr --slice $slice --files $nfiles --tag $tag $dir_path
 fi
 #
 # Print out ending time
diff --git a/Control/CalypsoExample/Reconstruction/scripts/faserMDC_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faserMDC_reco.py
index 942340f12c97db7e80850e0ede3be55d44bda54d..19b6564d8fa5f8fb0784918c064049a7c9e4c00f 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faserMDC_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faserMDC_reco.py
@@ -175,6 +175,14 @@ acc.merge(GhostBustersCfg(ConfigFlags))
 from FaserActsKalmanFilter.CKF2Config import CKF2Cfg
 acc.merge(CKF2Cfg(ConfigFlags, noDiagnostics=True))
 
+#
+# Try to write truth for mc
+#if args.isMC:
+#    xAODTruthCnvAlg = CompFactory.xAODMaker.xAODTruthCnvAlg
+#    acc.addEventAlgo(xAODTruthCnvAlg("AOD2xAOD", 
+#                                     AODContainerName="TruthEvent"), 
+#                     primary=True)
+
 #
 # Configure output
 from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
@@ -192,6 +200,13 @@ itemList = [ "xAOD::EventInfo#*"
 if args.isMC:
     # Add truth records here?
     itemList.extend( ["McEventCollection#*", "TrackerSimDataCollection#*"] )
+    # Try adding xAOD versions
+    #itemList.extend( ["xAOD::TruthEventContainer#TruthEvents", 
+    #                  "xAOD::TruthEventAuxContainer#TruthEventsAux.",
+    #                  "xAOD::TruthVertexContainer#TruthVertices", 
+    #                  "xAOD::TruthVertexAuxContainer#TruthVerticesAux.",
+    #                  "xAOD::TruthParticleContainer#TruthParticles", 
+    #                  "xAOD::TruthParticleAuxContainer#TruthParticlesAux."] )
 
 acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList))
 
@@ -216,9 +231,9 @@ if not args.isMC:
     replicaSvc.UseGeomSQLite = True
 
 # Configure verbosity    
-ConfigFlags.dump()
 if args.verbose:
     acc.foreach_component("*").OutputLevel = VERBOSE
+    ConfigFlags.dump()
 else:
     acc.foreach_component("*").OutputLevel = INFO
 
diff --git a/Simulation/G4Extensions/Pythia8Decayer/CMakeLists.txt b/Simulation/G4Extensions/Pythia8Decayer/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..57dcf64bd972fbc0cea03f89edf6a1871694c099
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/CMakeLists.txt
@@ -0,0 +1,31 @@
+################################################################################
+# Package: Pythia8Decayer
+################################################################################
+
+# Declare the package name:
+atlas_subdir( Pythia8Decayer )
+
+#if( NOT GENERATIONBASE )
+
+    # External dependencies:
+    find_package( CLHEP )
+    find_package( Geant4 )
+    find_package( XercesC )
+    find_package( Pythia8 )
+
+    # Component(s) in the package:
+    atlas_add_library( Pythia8DecayerLib
+                       src/*.cxx
+                       NO_PUBLIC_HEADERS
+                       PRIVATE_INCLUDE_DIRS ${GEANT4_INCLUDE_DIRS} ${XERCESC_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${PYTHIA8_INCLUDE_DIRS}
+                       PRIVATE_LINK_LIBRARIES ${GEANT4_LIBRARIES} ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} ${PYTHIA8_LIBRARIES} FaserMCTruth GaudiKernel AthenaBaseComps G4AtlasInterfaces G4AtlasToolsLib Pythia8_iLib )
+
+    atlas_add_component( Pythia8Decayer
+                         src/components/*.cxx
+                         INCLUDE_DIRS ${GEANT4_INCLUDE_DIRS} ${XERCESC_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${PYTHIA8_INCLUDE_DIRS}
+                         LINK_LIBRARIES ${GEANT4_LIBRARIES} ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} ${PYTHIA8_LIBRARIES} Pythia8DecayerLib AthenaBaseComps G4AtlasInterfaces )       
+#endif()
+
+# Install files from the package:
+atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
+#atlas_install_joboptions( share/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} --extend-ignore=F401,F821 )
diff --git a/Simulation/G4Extensions/Pythia8Decayer/python/Pythia8DecayerConfigNew.py b/Simulation/G4Extensions/Pythia8Decayer/python/Pythia8DecayerConfigNew.py
new file mode 100644
index 0000000000000000000000000000000000000000..6085f5cca62f25d1d6fa4e6d97fb6ab30575cccd
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/python/Pythia8DecayerConfigNew.py
@@ -0,0 +1,11 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+
+from AthenaConfiguration.ComponentFactory import CompFactory
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaCommon import Logging
+import os
+
+def Pythia8DecayerPhysicsToolCfg(flags, name='Pythia8DecayerPhysicsTool', **kwargs):
+    result = ComponentAccumulator()
+    result.setPrivateTools( CompFactory.Pythia8DecayerPhysicsTool(name,**kwargs) )
+    return result
diff --git a/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.cxx b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a6f9ea097841a1888761efde8f32303cff44a4ba
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.cxx
@@ -0,0 +1,147 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// Header for my class
+#include "Pythia8Decayer.h"
+
+#include "FaserMCTruth/FaserTrackInformation.h"
+
+// For passing things around
+#include "CLHEP/Vector/LorentzVector.h"
+#include "G4Track.hh"
+#include "G4DynamicParticle.hh"
+#include "G4ParticleTable.hh"
+#include "G4DecayProducts.hh"
+#include "G4VProcess.hh"
+
+Pythia8Decayer::Pythia8Decayer( const std::string s )
+ : G4VExtDecayer(s)
+{
+  // In the constructor, make a decayer instance, so that it's initialized here and not in the event loop
+  std::string docstring = Pythia8_i::xmlpath();
+  m_decayer = std::make_unique<Pythia8::Pythia>(docstring);
+
+  m_decayer->readString("ProcessLevel:all = off"); 
+
+  m_decayer->readString("ProcessLevel:resonanceDecays=on");
+    
+  // shut off Pythia8 (default) verbosty
+  //
+  m_decayer->readString("Init:showAllSettings=false");
+  m_decayer->readString("Init:showChangedSettings=false");
+  m_decayer->readString("Init:showAllParticleData=false");
+  m_decayer->readString("Init:showChangedParticleData=false");
+  //
+  // specify how many Py8 events to print out, at either level
+  // in this particular case print out a maximum of 10 events
+  //
+  m_decayer->readString("Next:numberShowProcess = 0" );
+  m_decayer->readString("Next:numberShowEvent = 10");
+           
+  m_decayer->init();
+   
+  // shut off decays of pi0's as we want Geant4 to handle them
+  // if other immediate decay products should be handled by Geant4,
+  // their respective decay modes should be shut off as well
+  //
+  m_decayer->readString("111:onMode = off");
+}
+
+G4DecayProducts* Pythia8Decayer::ImportDecayProducts(const G4Track& aTrack){
+
+   m_decayer->event.reset();
+   
+   G4DecayProducts* dproducts = nullptr;   
+   
+   G4ParticleDefinition* pd = aTrack.GetDefinition();
+   int    pdgid   = pd->GetPDGEncoding();
+   
+   // check if pdgid is consistent with Pythia8 particle table
+   //   
+   if ( !m_decayer->particleData.findParticle( pdgid ) )
+   {
+      G4cout << " can NOT find pdgid = " << pdgid 
+             << " in Pythia8::ParticleData" << G4endl;
+      return dproducts;
+   }
+   
+   if ( !m_decayer->particleData.canDecay(pdgid) )
+   {
+      G4cout << " Particle of pdgid = " << pdgid 
+             << " can NOT be decayed by Pythia8" << G4endl;
+      return dproducts;
+   }
+   
+   // NOTE: Energy should be in GeV 
+
+   m_decayer->event.append( pdgid, 1, 0, 0, 
+                           aTrack.GetMomentum().x() / CLHEP::GeV, 
+                           aTrack.GetMomentum().y() / CLHEP::GeV,  
+                           aTrack.GetMomentum().z() / CLHEP::GeV,
+                           aTrack.GetDynamicParticle()->GetTotalEnergy() / CLHEP::GeV,
+                           pd->GetPDGMass() / CLHEP::GeV );
+
+   // specify polarization, if any
+   
+   // special logic for primary taus(anti-taus), assumed to have polarization -1(+1), respectively
+   // verified from the polarization info in Genie output
+   double spinup;
+   FaserTrackInformation* info = dynamic_cast<FaserTrackInformation*>(aTrack.GetUserInformation());
+   if (info != nullptr && abs(pdgid) == 15 && (info->GetClassification() == TrackClassification::Primary || info->GetClassification() == TrackClassification::RegeneratedPrimary))
+   {
+       G4cout << "Primary tau decay identified." << G4endl;
+       spinup = (pdgid > 0 ? -1 : +1);
+   }
+   else
+   // NOTE: while in Py8 polarization is a double variable , 
+   //       in reality it's expected to be -1, 0., or 1 in case of "external" tau's, 
+   //       similar to LHA SPINUP; see Particle Decays, Hadron and Tau Decays in docs at
+   //       https://pythia.org/manuals/pythia8305/Welcome.html
+   //       so it's not able to handle anything like 0.99, thus we're rounding off    
+   {
+       spinup = round( std::cos( aTrack.GetPolarization().angle( aTrack.GetMomentumDirection() ) ) );
+   }
+   G4cout << "Using " << aTrack.GetParticleDefinition()->GetParticleName() << " helicity " << spinup << " for Pythia8 decay." << G4endl;
+   m_decayer->event.back().pol( spinup );
+
+   int npart_before_decay = m_decayer->event.size();
+   
+   m_decayer->next();
+   
+   int npart_after_decay = m_decayer->event.size();
+   
+   // create & fill up decay products
+   //
+   dproducts = new G4DecayProducts(*(aTrack.GetDynamicParticle()));
+   
+   // create G4DynamicParticle out of each m_decayer->event entry (except the 1st one)
+   // and push into dproducts
+   
+   for ( int ip=npart_before_decay; ip<npart_after_decay; ++ip )
+   {
+      
+      // only select final state decay products (direct or via subsequent decays);
+      // skip all others
+      //
+      // NOTE: in general, final state decays products will have 
+      //       positive status code between 91 and 99 
+      //       (in case such information could be of interest in the future)
+      //
+      if ( m_decayer->event[ip].status() < 0 ) continue;
+            
+      // should we also skip neutrinos ???
+      // if so, skip products with abs(m_decayer->event[ip].id()) of 12, 14, or 16
+            
+      G4ParticleDefinition* pddec = 
+         G4ParticleTable::GetParticleTable()->FindParticle( m_decayer->event[ip].id() );
+      if ( !pddec ) continue; // maybe we should print out a warning !
+      G4ThreeVector momentum = G4ThreeVector( m_decayer->event[ip].px() * CLHEP::GeV,
+                                              m_decayer->event[ip].py() * CLHEP::GeV,
+                                              m_decayer->event[ip].pz() * CLHEP::GeV ); 
+      dproducts->PushProducts( new G4DynamicParticle( pddec, momentum) ); 
+   }
+   dproducts->DumpInfo();
+   return dproducts;
+
+}
diff --git a/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.h b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.h
new file mode 100644
index 0000000000000000000000000000000000000000..61b52afc41f5b22e42d39052199af5b228d2032e
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.h
@@ -0,0 +1,63 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+//
+/// \file eventgenerator/pythia/pythia8decayer/include/Py8Decayer.hh
+/// \brief Definition of the Py8Decayer class
+///
+/// \author J. Yarba; FNAL
+///
+
+#ifndef Pythia8Decayer_H
+#define Pythia8Decayer_H
+
+#include "G4VExtDecayer.hh"
+#include "globals.hh"
+
+#include "Pythia8_i/Pythia8_i.h"
+
+class G4Track;
+class G4DecayProducts;
+
+class Pythia8Decayer : public G4VExtDecayer
+{
+  
+   public:
+
+      //ctor & dtor
+      Pythia8Decayer( const std::string s );
+
+      virtual G4DecayProducts* ImportDecayProducts(const G4Track&);
+    
+   private:
+   
+      // data members
+      std::unique_ptr<Pythia8::Pythia> m_decayer;
+
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
diff --git a/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.cxx b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..750272b2c5cfa64e85c55acbb991a7b09208683d
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.cxx
@@ -0,0 +1,108 @@
+/*
+  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+*/
+
+// class header
+#include "Pythia8DecayerPhysicsTool.h"
+// package headers
+#include "Pythia8Decayer.h"
+// Geant4 headers
+#include "globals.hh"
+#include "G4ParticleTable.hh"
+#include "G4VProcess.hh"
+#include "G4ProcessManager.hh"
+#include "G4Decay.hh"
+#include "G4DecayTable.hh"
+
+// STL headers
+#include <string>
+
+//=============================================================================
+// Standard constructor, initializes variables
+//=============================================================================
+Pythia8DecayerPhysicsTool::Pythia8DecayerPhysicsTool( const std::string& type,
+                                                      const std::string& nam,
+                                                      const IInterface* parent )
+  : base_class ( type, nam , parent )
+{
+}
+
+//=============================================================================
+// Destructor
+//=============================================================================
+
+Pythia8DecayerPhysicsTool::~Pythia8DecayerPhysicsTool()
+{
+}
+
+//=============================================================================
+// Initialize
+//=============================================================================
+StatusCode Pythia8DecayerPhysicsTool::initialize( )
+{
+  ATH_MSG_DEBUG("Pythia8DecayerPhysicsTool initialize( )");
+  this->SetPhysicsName(name());
+  return StatusCode::SUCCESS;
+}
+
+Pythia8DecayerPhysicsTool* Pythia8DecayerPhysicsTool::GetPhysicsOption()
+{
+  return this;
+}
+
+
+void Pythia8DecayerPhysicsTool::ConstructParticle()
+{
+    // nothing to do here
+}
+void Pythia8DecayerPhysicsTool::ConstructProcess()
+{
+  ATH_MSG_DEBUG("Pythia8DecayerPhysicsTool::ConstructProcess() called ");
+
+   Pythia8Decayer* extDecayer = new Pythia8Decayer(name());
+
+   auto particleIterator=GetParticleIterator();
+   particleIterator->reset();
+   while ((*particleIterator)())
+   {    
+      G4ParticleDefinition* particle = particleIterator->value();
+      G4ProcessManager* pmanager = particle->GetProcessManager();    
+      G4ProcessVector* processVector = pmanager->GetProcessList();
+      for ( size_t i=0; i<processVector->length(); ++i ) 
+      {    
+         G4Decay* decay = dynamic_cast<G4Decay*>((*processVector)[i]);
+         if ( decay ) 
+         {
+            // remove native/existing decay table for
+            // a)tau's 
+            // b) B+/- 
+            // and replace with external decayer
+            if ( std::abs(particle->GetPDGEncoding()) == 15 ||
+                 std::abs(particle->GetPDGEncoding()) == 521 )
+            {
+               if ( particle->GetDecayTable() )
+               {
+                  delete particle->GetDecayTable();
+                  particle->SetDecayTable(nullptr);
+               }
+               decay->SetExtDecayer(extDecayer);
+               G4cout << "Setting ext decayer for: " 
+                <<  particleIterator->value()->GetParticleName()
+                << G4endl;
+            }
+            // now set external decayer to all particles 
+            // that don't yet have a decay table
+            if ( !particle->GetDecayTable() )
+            {
+               decay->SetExtDecayer(extDecayer);
+               G4cout << "Setting ext decayer for: " 
+                <<  particleIterator->value()->GetParticleName()
+                << G4endl;
+            }
+         }
+      }              
+   }
+
+   return;
+   
+}
diff --git a/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.h b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..9acd13bdaee9a2c145e317161bca9f64251c7a9b
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8DecayerPhysicsTool.h
@@ -0,0 +1,43 @@
+/*
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS and FASER collaborations
+*/
+
+#ifndef PYTHIA8DECAYER_PYTHIA8DECAYERPHYSICSTOOL_H
+#define PYTHIA8DECAYER_PYTHIA8DECAYERPHYSICSTOOL_H
+
+// Include files
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "G4AtlasInterfaces/IPhysicsOptionTool.h"
+#include "G4VPhysicsConstructor.hh"
+
+/** @class Pythia8DecayerPhysicsTool Pythia8DecayerPhysicsTool.h "Pythia8Decayer/Pythia8DecayerPhysicsTool.h"
+ *
+ *
+ *
+ *  @author Edoardo Farina (modified for FASER by D. Casper)
+ *  @date  2015-05-14
+ */
+class Pythia8DecayerPhysicsTool :  public G4VPhysicsConstructor, public extends<AthAlgTool, IPhysicsOptionTool>
+{
+public:
+  /// Standard constructor
+  Pythia8DecayerPhysicsTool( const std::string& type , const std::string& name,
+                             const IInterface* parent ) ;
+
+  virtual ~Pythia8DecayerPhysicsTool( ); ///< Destructor
+
+  /// Initialize method
+  virtual StatusCode initialize( ) ;
+  virtual void ConstructParticle();
+  virtual void ConstructProcess();
+
+  /** Implements
+   */
+
+  virtual Pythia8DecayerPhysicsTool* GetPhysicsOption();
+
+};
+
+
+
+#endif
diff --git a/Simulation/G4Extensions/Pythia8Decayer/src/components/Pythia8Decayer_entries.cxx b/Simulation/G4Extensions/Pythia8Decayer/src/components/Pythia8Decayer_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8bbd46608e2fe53980a4338c03dcd06abcf61e59
--- /dev/null
+++ b/Simulation/G4Extensions/Pythia8Decayer/src/components/Pythia8Decayer_entries.cxx
@@ -0,0 +1,4 @@
+#include "../Pythia8DecayerPhysicsTool.h"
+
+DECLARE_COMPONENT( Pythia8DecayerPhysicsTool )
+
diff --git a/Simulation/G4Faser/G4FaserServices/python/G4FaserServicesConfigNew.py b/Simulation/G4Faser/G4FaserServices/python/G4FaserServicesConfigNew.py
index 839a4fda883d0f8e39f04f4efdc7a632918a6176..1b17627a2ccf83cce04907f5784633b597497d50 100644
--- a/Simulation/G4Faser/G4FaserServices/python/G4FaserServicesConfigNew.py
+++ b/Simulation/G4Faser/G4FaserServices/python/G4FaserServicesConfigNew.py
@@ -5,6 +5,8 @@ from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
 
 DetectorGeometrySvc, G4AtlasSvc, G4GeometryNotifierSvc, PhysicsListSvc=CompFactory.getComps("DetectorGeometrySvc","G4AtlasSvc","G4GeometryNotifierSvc","PhysicsListSvc",)
 from G4FaserTools.G4GeometryToolConfig import G4AtlasDetectorConstructionToolCfg
+from G4StepLimitation.G4StepLimitationConfigNew import G4StepLimitationToolCfg
+from Pythia8Decayer.Pythia8DecayerConfigNew import Pythia8DecayerPhysicsToolCfg
 
 def DetectorGeometrySvcCfg(ConfigFlags, name="DetectorGeometrySvc", **kwargs):
     result = ComponentAccumulator()
@@ -34,9 +36,9 @@ def G4GeometryNotifierSvcCfg(ConfigFlags, name="G4GeometryNotifierSvc", **kwargs
 
 def PhysicsListSvcCfg(ConfigFlags, name="PhysicsListSvc", **kwargs):
     result = ComponentAccumulator()
-    G4StepLimitationTool = CompFactory.G4StepLimitationTool
-    PhysOptionList = [G4StepLimitationTool("G4StepLimitationTool")]
+    PhysOptionList = [ result.popToolsAndMerge(G4StepLimitationToolCfg(ConfigFlags)) ]    
     #PhysOptionList += ConfigFlags.Sim.PhysicsOptions # FIXME Missing functionality
+    PhysOptionList += [ result.popToolsAndMerge(Pythia8DecayerPhysicsToolCfg(ConfigFlags)) ]
     PhysDecaysList = []
     kwargs.setdefault("PhysOption", PhysOptionList)
     kwargs.setdefault("PhysicsDecay", PhysDecaysList)