diff --git a/Control/CalypsoExample/BarcodeChecker/python/BarcodeCheckerAlg.py b/Control/CalypsoExample/BarcodeChecker/python/BarcodeCheckerAlg.py
index 0e29c7d2e60b880f22cf383e766ee9ce338f27d3..2aace46b67f79a3d72ea80e19cdf764bc9fa7d61 100644
--- a/Control/CalypsoExample/BarcodeChecker/python/BarcodeCheckerAlg.py
+++ b/Control/CalypsoExample/BarcodeChecker/python/BarcodeCheckerAlg.py
@@ -22,7 +22,7 @@ class BarcodeCheckerAlg(PyAthena.Alg):
     def execute(self):
         evtCollection = self.evtStore[self.MCEventKey] 
         for mcEvt in evtCollection:
-            for mcParticle in mcEvt.particles:
+            for mcParticle in mcEvt.particles():
                 barCode = mcParticle.barcode()
                 if barCode%1000000 <= 200000:
                     self.maxLow = max(self.maxLow, barCode%1000000)
@@ -36,4 +36,4 @@ class BarcodeCheckerAlg(PyAthena.Alg):
         print("Mid part: ", self.maxMid, " out of ", 1000000 - 200000, " (",100*self.maxMid/(1000000-200000),"% of overflow )")
         print("Hi  part: ", self.maxHi, " out of ", (1<<31)//1000000, " (", 100*self.maxHi/((1<<31)//1000000),"% of overflow )")
 
-        return StatusCode.Success
\ No newline at end of file
+        return StatusCode.Success
diff --git a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
index f711d7fdaa1c831ec06106a097d06085bd6baf4d..c2ec26d65d8b55f5cd4f2b2da1c9671dbffba2a9 100755
--- a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
+++ b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py
@@ -255,7 +255,7 @@ if __name__ == '__main__':
 #
     if doShiftLOS:
 
-        import McParticleEvent.Pythonizations
+        #import McParticleEvent.Pythonizations
         from GeneratorUtils.ShiftLOSConfig import ShiftLOSCfg
         cfg.merge(ShiftLOSCfg(configFlags, 
                             xcross = configFlags.Sim.Beam.xangle, ycross = configFlags.Sim.Beam.yangle,
diff --git a/Control/CalypsoExample/TruthEventDumper/python/TruthEventDumperAlg.py b/Control/CalypsoExample/TruthEventDumper/python/TruthEventDumperAlg.py
index 8986c1165ece677bb37a42998bc2e0d05642cd8e..c2b01c8b99bd2b4185e7a0adc78b9e887a3eb6a0 100644
--- a/Control/CalypsoExample/TruthEventDumper/python/TruthEventDumperAlg.py
+++ b/Control/CalypsoExample/TruthEventDumper/python/TruthEventDumperAlg.py
@@ -23,7 +23,7 @@ class TruthEventDumperAlg(PyAthena.Alg):
         evtCollection = self.evtStore[self.MCEventKey] 
         for mcEvt in evtCollection:
             mcEvt.print()
-            # for mcParticle in mcEvt.particles:
+            # for mcParticle in mcEvt.particles():
             #     barCode = mcParticle.barcode()
             #     self.maxLow = max(self.maxLow, barCode%200000)
             #     if barCode%1000000 > 200000:
diff --git a/Control/CalypsoExample/TruthEventDumper/share/TruthEventDumper_jobOptions.py b/Control/CalypsoExample/TruthEventDumper/share/TruthEventDumper_jobOptions.py
index 303f4b8d26db18b0cf96f30e500548ecc1a84208..75ff2533af5ce3b955d0807366a3c71258974c63 100644
--- a/Control/CalypsoExample/TruthEventDumper/share/TruthEventDumper_jobOptions.py
+++ b/Control/CalypsoExample/TruthEventDumper/share/TruthEventDumper_jobOptions.py
@@ -33,6 +33,8 @@ svcMgr.EventSelector.InputCollections = INPUT
 
 from TruthEventDumper.TruthEventDumperAlg import TruthEventDumperAlg
 
+from GeneratorUtils import HepMCFixes    
+
 from AthenaCommon.AlgSequence import AlgSequence
 job = AlgSequence()
 job += TruthEventDumperAlg(MCEventKey=COLLECTION)
diff --git a/Generators/ForeseeGenerator/python/Validate.py b/Generators/ForeseeGenerator/python/Validate.py
index d20a62da5119c8a39085d78512060bf0aa74f15f..b441817f4032c7bfc298741547a774e81f7b41dd 100644
--- a/Generators/ForeseeGenerator/python/Validate.py
+++ b/Generators/ForeseeGenerator/python/Validate.py
@@ -1,22 +1,23 @@
 from AthenaPython.PyAthena import StatusCode, McEventCollection, HepMC, CLHEP
 from GeneratorModules.EvgenAnalysisAlg import EvgenAnalysisAlg
 
+
 import ROOT as R
 import numpy as np
 import os
 from math import sqrt
 
-def fix():
-    "Python Fixes for HepMC"
-    def add(self, other):
-        self.set(self.x() + other.x(), self.y() + other.y(),
-                 self.z() + other.z(), self.t() + other.t())
-        return self
-
-    HepMC.FourVector.__iadd__ = add
-    del add
-
-    return
+# def fix():
+#     "Python Fixes for HepMC"
+#     def add(self, other):
+#         self.set(self.x() + other.x(), self.y() + other.y(),
+#                  self.z() + other.z(), self.t() + other.t())
+#         return self
+# 
+#     HepMC.FourVector.__iadd__ = add
+#     del add
+# 
+#     return
 
 class HistSvc(object):
     "Class to deal with histograms"
@@ -155,16 +156,19 @@ class EvgenValidation(EvgenAnalysisAlg):
     def execute(self):
         evt = self.events()[0]
         self.weight = evt.weights()[0] if evt.weights() else 1
-
+        
         # Loop over all particles in events 
         momenta = []
         vertices = []
         mother = HepMC.FourVector(0,0,0,0)
         llp_vtx = None
+
+        #evt.print()
         
-        for i, p in enumerate(evt.particles):
-            #print("--- ", i)
-            p.print()
+        for i, p in enumerate(evt.particles()):
+            #print(">>> ", i)
+            #p.print()
+            
             self.fillDaughter(p)
 
             if self.mother_stored:
@@ -176,17 +180,12 @@ class EvgenValidation(EvgenAnalysisAlg):
                 momenta.append(p.momentum())    
                 mother += p.momentum()
 
-            if i == 0 and p.production_vertex():
-                #p.production_vertex().print()
-                llp_vtx = p.production_vertex().point3d()
-            elif i == 0 and p.end_vertex():
-                llp_vtx = p.end_vertex().point3d()                            
-
-            if p.production_vertex():
-                vertices.append(p.production_vertex().point3d())
-
+            if i == 0 and p.end_vertex():
+                #p.end_vertex().print()
+                llp_vtx = p.end_vertex().position()#.point3d()                            
+                
             if p.production_vertex():
-                vertices.append(p.production_vertex().point3d())
+                vertices.append(p.production_vertex().position()) #point3d())
 
         # Fill daughter plots
         for i in range(self.ndaughters):
@@ -248,8 +247,9 @@ if __name__ == "__main__":
     from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
     from AthenaConfiguration.ComponentFactory import CompFactory
 
-    import McParticleEvent.Pythonizations
-    fix()
+    from GeneratorUtils import HepMCFixes    
+    #import McParticleEvent.Pythonizations
+    #fix()
     
     acc = ComponentAccumulator()
     valid = EvgenValidation("EvgenValidation", ndaughters =  args.ndaughters, outname = args.output, mother_stored = args.mother_stored)
diff --git a/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py b/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
index 2c6af9baf11cc005b1716ec3a53ace12ed397882..992daaba38d93486fabdcf3651476cdc562b643f 100644
--- a/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
+++ b/Generators/ForeseeGenerator/share/convert_hepmc_to_evnt.py
@@ -27,9 +27,9 @@ if __name__ == "__main__":
 
     from CalypsoConfiguration.AllConfigFlags import initConfigFlags
     configFlags = initConfigFlags()
-    configFlags.Input.RunNumber = [12345]
+    configFlags.Input.RunNumbers = [12345]
     configFlags.Input.OverrideRunNumber = True
-    configFlags.Input.LumiBlockNumber = [1]
+    configFlags.Input.LumiBlockNumbers = [1]
 
     configFlags.Input.isMC = True
     configFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"             # Always needed; must match FaserVersion
@@ -40,6 +40,11 @@ if __name__ == "__main__":
 
     configFlags.Exec.MaxEvents= -1
 
+    configFlags.addFlag("Sim.Beam.xangle", 0)
+    configFlags.addFlag("Sim.Beam.yangle", -160e-6)
+    configFlags.addFlag("Sim.Beam.xshift", 0)
+    configFlags.addFlag("Sim.Beam.yshift", 12)
+    
     import sys
     configFlags.fillFromArgs(sys.argv[1:])
 
@@ -59,14 +64,14 @@ if __name__ == "__main__":
         cfg.merge(HepMCReaderCfg(configFlags))
 
     if doShiftLOS:
-        import McParticleEvent.Pythonizations
+        #import McParticleEvent.Pythonizations
         from GeneratorUtils.ShiftLOSConfig import ShiftLOSCfg
-        cfg.merge(ShiftLOSCfg(configFlags, 
-                              xcross = 0, 
-                              ycross = -150e-6,
-                              xshift = 0,
-                              yshift = 12))
-
+        cfg.merge(ShiftLOSCfg(configFlags,
+                              xcross = configFlags.Sim.Beam.xangle, #0, 
+                              ycross = configFlags.Sim.Beam.yangle, # -160e-6, #+160e-6, #-160e-6,
+                              xshift = configFlags.Sim.Beam.xshift, #0,
+                              yshift = configFlags.Sim.Beam.yshift)) #12))
+                              
     from McEventSelector.McEventSelectorConfig import McEventSelectorCfg
     cfg.merge(McEventSelectorCfg(configFlags))    
 
diff --git a/Generators/ForeseeGenerator/share/generate_forsee_events.py b/Generators/ForeseeGenerator/share/generate_forsee_events.py
index cb8db6c591469f837dd61a83885138b6c9bccbf5..a260f5926ff047bad67b70ad6f3949c7fccc3c32 100644
--- a/Generators/ForeseeGenerator/share/generate_forsee_events.py
+++ b/Generators/ForeseeGenerator/share/generate_forsee_events.py
@@ -4,22 +4,43 @@ import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib
 
+
+PID_TO_MODE = { 
+    (-11, 11) : "e_e",
+    (11, -11) : "e_e",                       
+    (-13, 13) : "mu_mu",
+    (13, -13) : "mu_mu",            
+    (22, 22) : "gamma_gamma",
+    (111, 22) : "pigamma",
+    (22, 111) : "pigamma",            
+    (321, -321) : "KK",
+    (-321, 321) : "KK",
+    (111, 111) : "pi0_pi0",
+    (211, -211) : "pi+pi-",
+    (-211, 211) : "pi+pi-",
+    (None, None) : "Hadrons",
+    None : "Hadrons",    
+    }
+
+MODE_TO_PID = dict((v, k) for k, v in PID_TO_MODE.items())
+
 class ForeseeGenerator(object):
     """
     Generate LLP particles within FASER acceptance from FORESEE
     """
     
-    def __init__(self, modelname, energy, mass, couplings, daughter1_pid, daughter2_pid, outdir = None, path = '.', randomSeed = 12345, t0Shift = 0, zFront = -1.5, notime = False, suffix = ""):
+    def __init__(self, modelname, energy, mass, couplings, daughter1_pid = None, daughter2_pid = None, decay_modes = None,
+                 outdir = None, path = '.', randomSeed = 12345, t0Shift = 0, zFront = -1.5, notime = False, suffix = "", 
+                 tracker_decays = False, yshift = -0.065,
+                 uncertainties = False, allowed_prod_modes = []):
 
         self.modelname = modelname
         self.energy = energy
         self.mass = mass
         self.couplings = [couplings] if isinstance(couplings, (str, int, float)) else couplings
-        self.daughter1_pid = daughter1_pid
-        self.daughter2_pid = daughter2_pid
         self.outdir = outdir
         self.path = path
-        self.version = 3  # Forsee "version"
+        self.version = 4  # Forsee "version" 1.X
         self.seed = randomSeed
         self.zfront = zFront
         self.t0 = t0Shift + zFront
@@ -27,26 +48,43 @@ class ForeseeGenerator(object):
         self.suffix = f"_{suffix}" if suffix else ""
         self.nbinsample = 1
         self.nevents = 0
+        self.tracker_decays = tracker_decays
+        self.uncertainties = uncertainties
+        self.weightnames = []
+        self.prod_modes = {}
+        self.allowed_prod_modes = allowed_prod_modes
+
+        self.distance = 480 # m
+        self.radius = 0.11 # m
+        self.yshift = yshift
+
+        self.wider_decay_volume = False
+
 
         # Set decay mode ...
         
-        self.pid_map = { 
-            (-11, 11) : "e_e",
-            (11, -11) : "e_e",                       
-            (-13, 13) : "mu_mu",
-            (13, -13) : "mu_mu",            
-            (22, 22) : "gamma_gamma",
-            }
+        if decay_modes:
+            if daughter1_pid is not None and daughter2_pid is not None:
+                sys.exit("Cannot suply both decay modes and daughter pids")
+            else:
+                self.decay_modes = decay_modes
+        else:
+            if daughter1_pid is not None and daughter2_pid is not None:
+                self.daughter1_pid = daughter1_pid
+                self.daughter2_pid = daughter2_pid
+                decay_modes = PID_TO_MODE.get((self.daughter1_pid, self.daughter2_pid), None)
+                if decay_modes is None:
+                    sys.exit(f"Undefined decay to {self.daughter1_pid} + {self.daughter2_pid} for {self.modelname}")
+                
+                self.decay_modes = [decay_modes]
 
-        self.mode = self.pid_map.get((self.daughter1_pid, self.daughter2_pid), None)
-        if self.mode is None:
-            sys.exit(f"Undefined decay to {self.daughter1_pid} + {self.daughter2_pid} for {self.modelname}")
+            else:
+                sys.exit("Must suply either decay modes or daughter pids")
 
         # Set detector ...
-        if self.version == 1:
-            self.foresee = Foresee()
-        else:
-            self.foresee = Foresee(path = self.path)
+        self.foresee = Foresee(path = self.path)
+        if self.version >= 4:            
+            self.foresee.rng.seed(self.seed)
 
         # Generate 6.5 cm high to account for translation from ATLAS to FASER coord. system
         # TODO: relax this a bit as daughters may enter even if mother doesn't
@@ -57,34 +95,82 @@ class ForeseeGenerator(object):
         # Generate 6.5 cm high to account for translation from ATLAS to FASER coord. system
         # Relax r to 11 cm to acount for any mismatches with calypso and for duaghter entering if mother doesn't
 
-        if "ALP" in self.modelname:
+        if self.wider_decay_volume:
+            print("*** Generating decays from just before end of magnet to end of calo ***")
+            self.distance = 480 + 1.5 + 2.2 #m
+            self.zfront = 2.2 #m
+            length = 1.3 # m
+            self.radius = 0.15 #m
+            print(f"   d = {self.distance} m -> zfront = {self.zfront} m, l = {length} m, r = {self.radius} m")
+        elif self.tracker_decays or "ALP" in self.modelname or all(m in ["gamma_gamma", "pigamma", "gamma", "pi0_pi0"] for m in self.decay_modes):
             # Since not relying on tracks then extend length to include spectrometer (magents + tracking stations)
-            length = 1.5 + 2 + 0.5 # m
+            print("*** Including decays in tracker volume ***")
+            length = 1.5 + 2 + 0.5 # m (to end of tracker)
+            #print("*** Including decays in tracker volume and calo ***")        
+            #length = 1.5 + 3.5 # m (to end of calo)            
         else:
             length = 1.5 # m
 
-        self.foresee.set_detector(selection="np.sqrt(x.x**2 + (x.y - 0.065)**2)< 0.11",  
-                                  channels=[self.mode], distance=480, length=length ,
-                                  luminosity=1/1000.) # 1 pb-1        
+
+
+#         self.foresee.set_detector(selection="np.sqrt(x.x**2 + (x.y - 0.065)**2)< 0.11",  
+#                                   channels=self.decay_modes, distance=480, length=length ,
+#                                   luminosity=1/1000.) # 1 pb-1        
 
 
         # Set model ...
-        if self.version == 1:
-            self.model = Model(self.modelname)
-        else:
-            self.model = Model(self.modelname, path = f"{self.path}/Models/{self.modelname}/")
+        self.model = Model(self.modelname, path = f"{self.path}/Models/{self.modelname}/")
         
         if self.modelname == "DarkPhoton":
-            self.data = self.darkphoton()
+            if "pigamma" in self.decay_modes:
+                idx = self.decay_modes.index("pigamma")
+                self.decay_modes[idx] = "pi0_gamma" 
+            
+            self.darkphoton()
         elif self.modelname == "ALP-W":
-            self.data = self.alp_W()
+            self.alp_W()
         elif self.modelname == "ALP-photon":
-            self.data = self.alp_photon()
+            if "gamma_gamma" in self.decay_modes:
+                idx = self.decay_modes.index("gamma_gamma")
+                self.decay_modes[idx] = "gamma"             
+            self.alp_photon()
+        elif self.modelname == "ALP-g":
+            if "gamma_gamma" in self.decay_modes:
+                idx = self.decay_modes.index("gamma_gamma")
+                self.decay_modes[idx] = "br_2Gamma" 
+
+            self.alp_gluon()
         elif self.modelname == "U(1)B-L":
-            self.data = self.bminusl()                        
+            self.bminusl()
+        elif self.modelname == "U(1)B":
+            self.decay_modes = [m.replace("_", "") for m in self.decay_modes]
+            self.u1b()
+        elif self.modelname == "UpPhilic":
+            if "gamma_gamma" in self.decay_modes:
+                idx = self.decay_modes.index("gamma_gamma")
+                self.decay_modes[idx] = "gamma" 
+
+            self.upphilic()     
+        elif self.modelname == "DarkHiggs":
+            self.darkhiggs()
         else:
             sys.exit(f"Unknown model {self.modelname}")
 
+
+        # Setup detector along with the decay modes it is sensitive to        
+        self.foresee.set_detector(selection=f"np.sqrt(x.x**2 + (x.y + {self.yshift})**2) < {self.radius}",  
+                                  channels=self.decay_modes, distance=self.distance, length=length ,
+                                  luminosity=1/1000.) # 1 pb-1        
+
+        if self.modelname == "ALP-photon":
+            # Position of TAN
+            self.foresee.distance_prod = 140 # m
+
+            
+
+        # Setup decays
+        self.data = self.decays()
+
         return
 
     def darkphoton(self):
@@ -96,189 +182,389 @@ class ForeseeGenerator(object):
             pid0 =  "111",
             pid1 = "22",
             br = "2.*0.99 * coupling**2 * pow(1.-pow(mass/self.masses('111'),2),3)",
-            generator = "EPOSLHC",
-            energy = self.energy,
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET"] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,            
             nsample = 100)
     
         self.model.add_production_2bodydecay(
             pid0 = "221",
             pid1 = "22",
             br = "2.*0.39 * coupling**2 * pow(1.-pow(mass/self.masses('221'),2),3)",
-            generator = "EPOSLHC",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET"] if self.uncertainties else "EPOSLHC",            
             energy = self.energy,
             nsample = 100)
 
         self.model.add_production_mixing(
             pid = "113",
             mixing = "coupling * 0.3/5. * 0.77545**2/abs(mass**2-0.77545**2+0.77545*0.147*1j)",
-            generator = "EPOSLHC",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET"] if self.uncertainties else "EPOSLHC",            
             energy = self.energy,
             )
 
-        if self.version == 1:
-            self.model.add_production_direct(
-                label = "Brem",
-                energy = self.energy,
-                condition = "p.pt<1",
-                coupling_ref=1,
-                )
-            
-            self.model.add_production_direct(
-                label = "DY",
-                energy = self.energy,
-                coupling_ref=1,
-                massrange=[1.5, 10.]
-                )
-
-            self.model.set_br_1d(
-                modes = [self.mode],
-                filenames=[f"files/models/{self.modelname}/br/{self.mode}.txt"] 
-                )
+        masses_brem = [ 
+            0.01  ,  0.0126,  0.0158,  0.02  ,  0.0251,  0.0316,  0.0398,
+            0.0501,  0.0631,  0.0794,  0.1   ,  0.1122,  0.1259,  0.1413,
+            0.1585,  0.1778,  0.1995,  0.2239,  0.2512,  0.2818,  0.3162,
+            0.3548,  0.3981,  0.4467,  0.5012,  0.5623,  0.6026,  0.631 ,
+            0.6457,  0.6607,  0.6761,  0.6918,  0.7079,  0.7244,  0.7413,
+            0.7586,  0.7762,  0.7943,  0.8128,  0.8318,  0.8511,  0.871 ,
+            0.8913,  0.912 ,  0.9333,  0.955 ,  0.9772,  1.    ,  1.122 ,
+            1.2589,  1.4125,  1.5849,  1.7783,  1.9953,  2.2387,  2.5119,
+            2.8184,  3.1623,  3.9811,  5.0119,  6.3096,  7.9433, 10.    
+            ]
+        
+        self.model.add_production_direct(
+            label = "Brem",
+            energy = self.energy,
+            condition = ["p.pt<1", "p.pt<2", "p.pt<0.5"] if self.uncertainties else "p.pt<1",
+            coupling_ref=1,
+            masses = masses_brem,
+            )
 
+#         # High mass        
+#         masses_dy = [ 
+#             1.5849,  1.7783,  1.9953,  2.2387,  2.5119, 2.8184,  3.1623,  3.9811,  5.0119,  6.3096,  7.9433, 10.    
+#             ]
+#         
+#         self.model.add_production_direct(
+#             label = "DY",
+#             energy = self.energy,
+#             coupling_ref=1,
+#             masses = masses_dy,
+#             )
+
+        # Production modes with order of weights for each (use nominal in case given uncertainty doesn't apply for that mode)
+        if self.uncertainties:        
+            self.prod_modes = {
+                "111"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward", "EPOSLHC", "EPOSLHC"],
+                "221"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward", "EPOSLHC", "EPOSLHC"],
+                "113"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward", "EPOSLHC", "EPOSLHC"],                
+                "Brem" : ["p.pt<1", "p.pt<1", "p.pt<1", "p.pt<1", "p.pt<2", "p.pt<0.5"]
+                }
+            self.weightnames = ["Nominal", "SIBYLL", "QGSJET", "Pythia8-Forward", "BremPtUp", "BremPtDown"]
         else:
-            masses_brem = [ 
-                0.01  ,  0.0126,  0.0158,  0.02  ,  0.0251,  0.0316,  0.0398,
-                0.0501,  0.0631,  0.0794,  0.1   ,  0.1122,  0.1259,  0.1413,
-                0.1585,  0.1778,  0.1995,  0.2239,  0.2512,  0.2818,  0.3162,
-                0.3548,  0.3981,  0.4467,  0.5012,  0.5623,  0.6026,  0.631 ,
-                0.6457,  0.6607,  0.6761,  0.6918,  0.7079,  0.7244,  0.7413,
-                0.7586,  0.7762,  0.7943,  0.8128,  0.8318,  0.8511,  0.871 ,
-                0.8913,  0.912 ,  0.9333,  0.955 ,  0.9772,  1.    ,  1.122 ,
-                1.2589,  1.4125,  1.5849,  1.7783,  1.9953,  2.2387,  2.5119,
-                2.8184,  3.1623,  3.9811,  5.0119,  6.3096,  7.9433, 10.    
-                ]
-            
-            self.model.add_production_direct(
-                label = "Brem",
-                energy = self.energy,
-                condition = "p.pt<1",
-                coupling_ref=1,
-                masses = masses_brem,
-                )
-            
-            masses_dy = [ 
-                1.5849,  1.7783,  1.9953,  2.2387,  2.5119, 2.8184,  3.1623,  3.9811,  5.0119,  6.3096,  7.9433, 10.    
-                ]
-            
-            self.model.add_production_direct(
-                label = "DY",
-                energy = self.energy,
-                coupling_ref=1,
-                masses = masses_dy,
-                )
+            self.prod_modes = {
+                "111"  : ["EPOSLHC"],
+                "221"  : ["EPOSLHC"],
+                "113"  : ["EPOSLHC"],
+                "Brem" : ["p.pt<1"]
+                }                
+            self.weightnames = ["Nominal"]
 
-            self.model.set_br_1d(
-                modes = ["e_e", "mu_mu"],
-                finalstates = [[11, -11], [13, -13]],
-                filenames=["model/br/e_e.txt", "model/br/mu_mu.txt"],
-                #filenames=[f"model/br/all.txt"] 
-                )
 
-        return self.decays()
+        modes = ["e_e", "mu_mu", "pi+_pi-", "pi0_gamma", "pi+_pi-_pi0", "K_K"]    
+        self.model.set_br_1d(
+            modes = modes, # ["e_e", "mu_mu"]
+            finalstates = [[11, -11], [13, -13], [211, -211], [111,22], None, [321, -321]], # [[11, -11], [13, -13]]
+            filenames = [f"model/br/{m}.txt" for m in modes]    # ["model/br/e_e.txt", "model/br/mu_mu.txt"]
+            )
 
+        return #self.decays()
 
-    def alp_W(self):
 
+    def alp_W(self):
+        
+        nsample = 10000
         self.nbinsample = 100 # resample bins to help with smoothness
 
         self.model.add_production_2bodydecay(
-            pid0 = "5",
+            pid0 = "511",
             pid1 = "321",
-            br = "2.2e4 * coupling**2 * np.sqrt((1-(mass+0.495)**2/5.279**2)*(1-(mass-0.495)**2/5.279**2))",
-            generator = "Pythia8",
+            br = "2.3e4 * coupling**2 * (1-mass**2/5.279**2)**2",
+            generator = ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min"] if self.uncertainties else "NLO-P8",
             energy = self.energy,
-            nsample = 500, # Vary over phi and theta -> increase number to improve smoothness for B
+            nsample = nsample,
             )
-            
+
+        self.model.add_production_2bodydecay(
+            pid0 = "-511",
+            pid1 = "321",
+            br = "2.3e4 * coupling**2 * (1-mass**2/5.279**2)**2",
+            generator = ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min"] if self.uncertainties else "NLO-P8",
+            energy = self.energy,
+            nsample = nsample,
+            )
+
+        self.model.add_production_2bodydecay(
+            pid0 = "521",
+            pid1 = "321",
+            br = "2.3e4 * coupling**2 * (1-mass**2/5.279**2)**2",
+            generator = ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min"] if self.uncertainties else "NLO-P8",
+            energy = self.energy,
+            nsample = nsample,
+            )
+
+        self.model.add_production_2bodydecay(
+            pid0 = "-521",
+            pid1 = "333",
+            br = "2.3e4 * coupling**2 * (1-mass**2/5.279**2)**2",
+            generator = ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min"] if self.uncertainties else "NLO-P8",
+            energy = self.energy,
+            nsample = nsample,
+            )
+
+
         self.model.add_production_2bodydecay(
-            pid0 = "-5",
+            pid0 = "531",
+            pid1 = "333",
+            br = "2.3e4 * coupling**2 * (1-mass**2/5.366**2)**2",
+            generator = ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min"] if self.uncertainties else "NLO-P8",
+            energy = self.energy,
+            nsample = nsample,
+            )
+
+        self.model.add_production_2bodydecay(
+            pid0 = "-531",
             pid1 = "321",
-            br = "2.2e4 * coupling**2 * np.sqrt((1-(mass+0.495)**2/5.279**2)*(1-(mass-0.495)**2/5.279**2))",
-            generator = "Pythia8",
+            br = "2.3e4 * coupling**2 * (1-mass**2/5.366**2)**2",
+            generator = ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min"] if self.uncertainties else "NLO-P8",
             energy = self.energy,
-            nsample = 500,
+            nsample = nsample,
             )
         
         self.model.add_production_2bodydecay(
             pid0 = "130",
             pid1 = "111",
             br = "4.5 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.495**2)*(1-(mass-0.135)**2/0.495**2))",
-            generator = "EPOSLHC",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
             energy = self.energy,
-            nsample = 50,
+            nsample = nsample,
             )
 
         self.model.add_production_2bodydecay(
             pid0 = "321",
             pid1 = "211",
             br = "10.5 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.495**2)*(1-(mass-0.135)**2/0.495**2))",
-            generator = "EPOSLHC",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
             energy = self.energy,
-            nsample = 50,
+            nsample = nsample,
             )
 
         self.model.add_production_2bodydecay(
             pid0 = "-321",
             pid1 = "211",
             br = "10.5 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.495**2)*(1-(mass-0.135)**2/0.495**2))",
-            generator = "EPOSLHC",
+            generator =  ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
             energy = self.energy,
-            nsample = 50,
+            nsample = nsample,
             )
 
-
-        if self.version == 1:
-            self.model.set_br_1d(
-                modes = [self.mode],
-                filenames=[f"files/models/{self.modelname}/br/{self.mode}.txt"] 
-                )
+        # Production modes with order of weights for each (use nominal in case given uncertainty doesn't apply for that mode)
+        if self.uncertainties:        
+            self.prod_modes = {
+                "511"  : ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min", "NLO-P8", "NLO-P8", "NLO-P8"],
+                "-511" : ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min", "NLO-P8", "NLO-P8", "NLO-P8"],
+                "521"  : ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min", "NLO-P8", "NLO-P8", "NLO-P8"],
+                "-521" : ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min", "NLO-P8", "NLO-P8", "NLO-P8"],
+                "531"  : ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min", "NLO-P8", "NLO-P8", "NLO-P8"],
+                "-531" : ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min", "NLO-P8", "NLO-P8", "NLO-P8"],
+                "130"  : ["EPOSLHC", "EPOSLHC", "EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"],
+                "321"  : ["EPOSLHC", "EPOSLHC", "EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"],
+                "-321" : ["EPOSLHC", "EPOSLHC", "EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"]
+                }
+            self.weightnames = ["Nominal", "Powheg-Max", "Poweheg-Min", "SIBYLL", "QGSJET", "Pythia8-Forward"]
         else:
-            self.model.set_br_1d(
-                modes = ["gamma_gamma"],
-                finalstates = [[22, 22]],
-                filenames=["model/br/gamma_gamma.txt"]
-                #filenames=[f"model/br/all.txt"] 
-                )   
+            self.prod_modes = {
+                "511"  : ["NLO-P8"],
+                "-511" : ["NLO-P8"],
+                "521"  : ["NLO-P8"],
+                "-521" : ["NLO-P8"],
+                "531"  : ["NLO-P8"],
+                "-531" : ["NLO-P8"],
+                "130"  : ["EPOSLHC"],
+                "321"  : ["EPOSLHC"],
+                "-321" : ["EPOSLHC"]
+                }
+            self.weightnames = ["Nominal"]
+                
+        modes = ["gamma_gamma", "e_e_gamma"]
+        self.model.set_br_1d(
+            modes = modes, # ["gamma_gamma"],
+            finalstates = [[22, 22], [11,-11,22]], # [[22, 22]],
+            filenames = [f"model/br/{m}.txt" for m in modes], #["model/br/gamma_gamma.txt"]
+            )   
 
-        return self.decays()
+        return #self.decays()
         
 
     def alp_photon(self):
 
-        if self.version == 1:
-            sys.exit("ALP-photon model not available in this version of FORESEE")
-
         self.nbinsample = 100 # resample bins to help with smoothness
+        nsample = 10000        
+
+        def conversion_factor(mass, coupling, p0):
+
+            # define quantities
+            alpha, me, ZFe, AFe = 1./137., 0.000511, 26, 56
+            a, d  = 111.*ZFe**(-1./3.) / me,  0.164 * AFe**(-2./3.)
+            SMXSinIGeV2 = 13311.9696379 
+
+            # integration boundaries theta 
+            logthmin, logthmax, nlogth = -12, 0, 20
+            dlogth = (logthmax-logthmin)/float(nlogth)
+            
+            # conversion probability 
+            xs_primakoff = 0
+            k = p0.e
+
+            for ltheta in np.linspace(logthmin+0.5*dlogth,logthmax-0.5*dlogth,nlogth):
+                theta = 10**ltheta
+                t = mass**4 / (2.*k**2) + k**2 * theta**2
+                if t<7.39*me*me: ff = a*a*t/(1.+a*a*t)
+                else: ff = 1./(1.+t/d)
+                xs_primakoff += alpha * ZFe**2 / 4. * ff**2 * k**4 * np.sin(theta)**3 *theta / t**2 * dlogth * np.log(10)
+
+            prob = xs_primakoff / SMXSinIGeV2
+            if np.isnan(prob): prob=0
         
-        masses = [float(x) for x in ['{:0.4e}'.format(m) for m in np.logspace(-2,0,20+1)]]
-        self.model.add_production_direct(
-            label = "Prim",
-            energy = self.energy,
-            coupling_ref = 1,
-            masses=masses,
+            # return
+            return np.sqrt(prob)
+
+        self.model.add_production_mixing(
+            label="Primakoff",
+            pid = "22",
+            mixing = conversion_factor,
+            generator=['EPOSLHC', 'SIBYLL', 'QGSJET', 'Pythia8-Forward'] if self.uncertainties else "EPOSLHC",
+            energy = self.energy
             )
 
+#         # Remove as only relevant above 1.5 TeV or so      
+#         masses = [float(x) for x in ['{:0.4e}'.format(m) for m in np.logspace(-2,0,20+1)]]
+#         self.model.add_production_direct(
+#             label = "Prim",
+#             energy = self.energy,
+#             coupling_ref = 1,
+#             masses=masses,
+#             )
+
+        # Production modes with order of weights for each (use nominal in case given uncertainty doesn't apply for that mode)
+        if self.uncertainties:        
+            self.prod_modes = {
+                "Primakoff"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"]
+                }
+            self.weightnames = ["Nominal", "SIBYLL", "QGSJET", "Pythia8-Forward"]
+        else:
+            self.prod_modes = {
+                "Primakoff"  : ["EPOSLHC"]
+                }
+            self.weightnames = ["Nominal"]
+
+
+        modes = ["gamma", "eegamma"]
         self.model.set_br_1d(
-            modes = ["gamma_gamma"],
-            finalstates = [[22, 22]],
-            filenames=["model/br/gamma_gamma.txt"]
+            modes = modes, # ["gamma_gamma"],
+            finalstates = [[22, 22], [11,-11,22]], # [[22, 22]],
+            filenames = [f"model/br/{m}.txt" for m in modes], #["model/br/gamma_gamma.txt"]
             )   
         
+        return #self.decays()
+
+    def alp_gluon(self):
+
+        self.nbinsample = 100 # resample bins to help with smoothness
+        nsample = 10000
+
+        self.model.add_production_2bodydecay(
+            pid0 = "511",
+            pid1 = "130",
+            br = "76 * coupling**2 * (1-(mass/5.279)**2)**2",
+            generator = ['NLO-P8','NLO-P8-Max', 'NLO-P8-Min'] if self.uncertainties else "NLO-P8",
+            energy = self.energy,
+            nsample = nsample, 
+            ) 
+        
+        self.model.add_production_2bodydecay(
+            pid0 = "-511",
+            pid1 = "130",
+            br = "76 * coupling**2 * (1-(mass/5.279)**2)**2",
+            generator = ['NLO-P8','NLO-P8-Max', 'NLO-P8-Min'] if self.uncertainties else "NLO-P8",
+            energy = self.energy,
+            nsample = nsample, 
+            ) 
+        
+        self.model.add_production_2bodydecay(
+            pid0 = "521",
+            pid1 = "321",
+            br = "76 * coupling**2 * (1-(mass/5.279)**2)**2",
+            generator = ['NLO-P8','NLO-P8-Max', 'NLO-P8-Min'] if self.uncertainties else "NLO-P8",
+            energy = self.energy,
+            nsample = nsample, 
+            ) 
         
-        return self.decays()
+        self.model.add_production_2bodydecay(
+            pid0 = "-521",
+            pid1 = "-321",
+            br = "76 * coupling**2 * (1-(mass/5.279)**2)**2",
+            generator = ['NLO-P8','NLO-P8-Max', 'NLO-P8-Min'] if self.uncertainties else "NLO-P8",
+            energy = self.energy,
+            nsample = nsample, 
+            )
+
+        self.model.add_production_mixing(
+            pid = "111",
+            mixing = "(3.671) * coupling * (-0.180 * mass**2)/(mass**2 - 0.1349**2)",
+            generator = ['EPOSLHC', 'SIBYLL', 'QGSJET', 'Pythia8-Forward'] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            )
+        
+        self.model.add_production_mixing(
+            pid = "221",
+            mixing = "(3.671) * coupling * (-0.3596 * mass**2 + 0.0021)/(mass**2 - 0.5478**2)",
+            generator = ['EPOSLHC', 'SIBYLL', 'QGSJET', 'Pythia8-Forward'] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            )
+        
+        self.model.add_production_mixing(
+            pid = "331",
+            mixing = "(3.671) * coupling * (-0.3362 * mass**2 + 0.0092)/(mass**2 - 0.9574**2)",
+            generator = ['EPOSLHC', 'SIBYLL', 'QGSJET', 'Pythia8-Forward'] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            )
+
+        # Production modes with order of weights for each (use nominal in case given uncertainty doesn't apply for that mode)
+        if self.uncertainties:        
+            self.prod_modes = {
+                "511"  : ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min", "NLO-P8", "NLO-P8", "NLO-P8"],
+                "-511" : ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min", "NLO-P8", "NLO-P8", "NLO-P8"],
+                "521"  : ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min", "NLO-P8", "NLO-P8", "NLO-P8"],
+                "-521" : ["NLO-P8", "NLO-P8-Max", "NLO-P8-Min", "NLO-P8", "NLO-P8", "NLO-P8"],
+                "111"  : ["EPOSLHC", "EPOSLHC", "EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"],                
+                "221"  : ["EPOSLHC", "EPOSLHC", "EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"],
+                "331"  : ["EPOSLHC", "EPOSLHC", "EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"]
+                }
+            self.weightnames = ["Nominal", "Powheg-Max", "Poweheg-Min", "SIBYLL", "QGSJET", "Pythia8-Forward"]
+        else:
+            self.prod_modes = {
+                "511"  : ["NLO-P8"],
+                "-511" : ["NLO-P8"],
+                "521"  : ["NLO-P8"],
+                "-521" : ["NLO-P8"],
+                "111"  : ["EPOSLHC"],
+                "221"  : ["EPOSLHC"],
+                "331"  : ["EPOSLHC"]
+                }
+            self.weightnames = ["Nominal"]
 
-    def bminusl(self):
 
-        if self.version == 1:
-            sys.exit("U(1)B-L model not available in this version of FORESEE")
+        modes =  ["br_2Gamma", "br_3Pi0", "br_Pi+Pi-Gamma", "br_Pi+Pi-Pi0", "br_2PiEta"]     
+        self.model.set_br_1d(
+            modes=modes,
+            finalstates=[[22,22], [111,111,111], [211,-211,22], [211,-211, 111], [211,-211,221]],
+            filenames=[f"model/br/{m}.dat" for m in modes],
+            )
+
+        return #self.decays()
+
 
-        self.nbinsample = 100 # TODO
+    def bminusl(self):
+
+        self.nbinsample = 100 
 
         self.model.add_production_2bodydecay(
             pid0 =  "111",
             pid1 = "22",
             br = "2.*0.99 * (coupling/0.303)**2 * pow(1.-pow(mass/self.masses('111'),2),3)",
-            generator = "EPOSLHC",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET"] if self.uncertainties else "EPOSLHC",
             energy = self.energy,
             nsample = 100)
     
@@ -286,7 +572,7 @@ class ForeseeGenerator(object):
             pid0 = "221",
             pid1 = "22",
             br = "2.*0.25*0.39 * (coupling/0.303)**2 * pow(1.-pow(mass/self.masses('221'),2),3)",
-            generator = "EPOSLHC",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET"] if self.uncertainties else "EPOSLHC",
             energy = self.energy,
             nsample = 100)
 
@@ -306,34 +592,275 @@ class ForeseeGenerator(object):
         self.model.add_production_direct(
             label = "Brem",
             energy = self.energy,
-            condition = "1./0.3**2 * (p.pt<1)",
-            coupling_ref=1,
+            condition = ["p.pt<1", "p.pt<2", "p.pt<0.5"] if self.uncertainties else "p.pt<1",            
+            coupling_ref=0.303, 
             masses = masses_brem,
             )
-            
+
+        if self.uncertainties:        
+            self.prod_modes = {
+                "111"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward", "EPOSLHC", "EPOSLHC"],
+                "221"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward", "EPOSLHC", "EPOSLHC"],
+                "Brem" : ["p.pt<1", "p.pt<1", "p.pt<1", "p.pt<1", "p.pt<2", "p.pt<0.5"]                
+                }
+            self.weightnames = ["Nominal", "SIBYLL", "QGSJET", "Pythia8-Forward", "BremPtUp", "BremPtDown"]
+        else:
+            self.prod_modes = {
+                "111"  : ["EPOSLHC"],
+                "221"  : ["EPOSLHC"],
+                "Brem" : ["p.pt<1"],
+                }
+            self.weightnames = ["Nominal"]
+
+
+        modes = ["e_e", "mu_mu", "nu_nu", "Hadrons"]
         self.model.set_br_1d(
-        modes = ["e_e", "mu_mu", "nu_nu", "Hadrons"],
+            modes = modes, 
             finalstates = [[11, -11], [13, -13], [12, -12], None],
-            filenames=["model/br/e_e.txt", "model/br/mu_mu.txt", "model/br/nu_nu.txt", "model/br/Hadrons.txt"],
-            #filenames=[f"model/br/all.txt"] 
+            filenames=[f"model/br/{m}.txt" for m in modes],
+            )
+        
+        return #self.decays()
+
+
+    def u1b(self):
+
+        self.nbinsample = 100 
+        nsample = 10000
+
+        self.model.add_production_2bodydecay(
+            pid0 = "111",
+            pid1 = "22",
+            br = "2*0.99 * (coupling/0.3)**2 * pow(1.-pow(mass/self.masses('111'),2),3)*np.abs((1-mass**2/0.78265**2 - 0.00849/0.78265*1j)**(-1))**2",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            nsample = nsample
+            )
+        
+        self.model.add_production_2bodydecay(
+            pid0 = "221",
+            pid1 = "22",
+            br = "1./16*0.3941 * (coupling/0.3)**2 * pow(1.-pow(mass/self.masses('221'),2),3)*np.abs((1-mass**2/0.78265**2 - 0.00849/0.78265*1j)**(-1)+(1-mass**2/1.019461**2-0.004249/1.019461*1j)**(-1))**2",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            nsample = nsample
+            )
+
+        self.model.add_production_2bodydecay(
+            pid0 = "331",
+            pid1 = "22",
+            br = "1./49.*0.02307 * (coupling/0.3)**2 * pow(1.-pow(mass/self.masses('331'),2),3) * np.abs((1-mass**2/0.78265**2 - 0.00849/0.78265*1j)**(-1) - 2*(1-mass**2/1.019461**2-0.004249/1.019461*1j)**(-1))**2",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            nsample = nsample
+            )
+
+        self.model.add_production_mixing(
+            pid = "223",
+            mixing = "coupling * 4./17. * 0.78265**2/abs(mass**2-0.78265**2+0.78265*0.00849*1j)",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            )
+        
+        self. model.add_production_mixing(
+            pid = "333",
+            mixing = "coupling /12.88 * 1.019461**2/abs(mass**2-1.019461**2+1.019461*0.004249*1j)",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            )
+        
+        masses_brem = [ 
+            0.01  ,  0.0126,  0.0158,  0.02  ,  0.0251,  0.0316,  0.0398,
+            0.0501,  0.0631,  0.0794,  0.1   ,  0.1122,  0.1259,  0.1413,
+            0.1585,  0.1778,  0.1995,  0.2239,  0.2512,  0.2818,  0.3162,
+            0.3548,  0.3981,  0.4467,  0.5012,  0.5623,  0.6026,  0.631 ,
+            0.6457,  0.6607,  0.6761,  0.6918,  0.7079,  0.7244,  0.7413,
+            0.7586,  0.7762,  0.7943,  0.8128,  0.8318,  0.8511,  0.871 ,
+            0.8913,  0.912 ,  0.9333,  0.955 ,  0.9772,  1.    ,  1.122 ,
+            1.2589,  1.4125,  1.5849,  1.7783,  1.9953,  2.2387,  2.5119,
+            2.8184,  3.1623,  3.9811,  5.0119,  6.3096,  7.9433, 10.    
+            ]
+        
+        self.model.add_production_direct(
+            label = "Brem",
+            energy = self.energy,
+            condition = ["p.pt<1", "p.pt<2", "p.pt<0.5"] if self.uncertainties else "p.pt<1",   
+            coupling_ref = 0.303,          
+            masses = masses_brem,
+            )
+
+        if self.uncertainties:        
+            self.prod_modes = {
+                "111"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward", "EPOSLHC", "EPOSLHC"],
+                "221"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward", "EPOSLHC", "EPOSLHC"],
+                "331"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward", "EPOSLHC", "EPOSLHC"],
+                "223"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward", "EPOSLHC", "EPOSLHC"],
+                "333"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward", "EPOSLHC", "EPOSLHC"],                 
+                "Brem" : ["p.pt<1", "p.pt<1", "p.pt<1", "p.pt<1", "p.pt<2", "p.pt<0.5"]                
+                }
+            self.weightnames = ["Nominal", "SIBYLL", "QGSJET", "Pythia8-Forward", "BremPtUp", "BremPtDown"]
+        else:
+            self.prod_modes = {
+                "111"  : ["EPOSLHC"],
+                "221"  : ["EPOSLHC"],
+                "331"  : ["EPOSLHC"],
+                "323"  : ["EPOSLHC"],
+                "333"  : ["EPOSLHC"],                
+                "Brem" : ["p.pt<1"]
+                }
+            self.weightnames = ["Nominal"]
+        
+        modes = ["pigamma", "3pi", "KK", "ee", "mumu"]
+        self.model.set_br_1d(
+            modes = modes,
+            finalstates = [[111,22], [211,-211,111], [321,-321], [11,-11], [13,-13]],
+            filenames=[f"model/br/{m}.txt" for m in modes],
+            )
+        
+        return #self.decays()
+
+    def upphilic(self):
+
+        self.nbinsample = 100
+        nsample = 10000
+
+        self.model.add_production_2bodydecay(
+            pid0 = "221",
+            pid1 = "111",
+            br = "1.26e5 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.547**2)*(1-(mass-0.135)**2/0.547**2))",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            nsample = nsample,
+            )
+        
+        self.model.add_production_2bodydecay(
+            pid0 = "331",
+            pid1 = "111",
+            br = "273. * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.957**2)*(1-(mass-0.135)**2/0.957**2))",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            nsample = nsample,
+            )
+        
+        self.model.add_production_2bodydecay(
+            pid0 = "321",
+            pid1 = "211",
+            br = "7.42 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.49368**2)*(1-(mass-0.135)**2/0.49368**2))",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            nsample = nsample,
+            )
+        
+        self.model.add_production_2bodydecay(
+            pid0 = "-321",
+            pid1 = "-211",
+            br = "7.42 * coupling**2 * np.sqrt((1-(mass+0.135)**2/0.49368**2)*(1-(mass-0.135)**2/0.49368**2))",
+            generator = ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"] if self.uncertainties else "EPOSLHC",
+            energy = self.energy,
+            nsample = nsample,
             )
+
+        if self.uncertainties:        
+            self.prod_modes = {
+                "221"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"],
+                "331"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"],
+                "321"  : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"],
+                "-321" : ["EPOSLHC", "SIBYLL", "QGSJET", "Pythia8-Forward"],                               
+                }
+            self.weightnames = ["Nominal", "SIBYLL", "QGSJET", "Pythia8-Forward"]
+        else:
+            self.prod_modes = {
+                "221"  : ["EPOSLHC"],
+                "331"  : ["EPOSLHC"],
+                "321"  : ["EPOSLHC"],
+                "-321" : ["EPOSLHC"],                
+                }
+            self.weightnames = ["Nominal"]
+
+
+        modes = ["gamma", "pi+_pi-", "pi0_pi0"]
+        self.model.set_br_1d(
+            modes = modes, 
+            finalstates = [[22,22], [211,-211], [111, 111]],
+            filenames=[f"model/br/{m}.txt" for m in modes],
+            )
+        
+        return #self.decays()
+
+    def darkhiggs(self):
+
+        self.nbinsample = 100 # resample bins to help with asymmetric detector
+
+        self.model.add_production_2bodydecay(
+            pid0 = "511",
+            pid1 = "130",
+            br = "5.6 * coupling**2 * pow(1.-pow(mass/5.279,2),2)",
+            generator = ['NLO-P8','NLO-P8-Max', 'NLO-P8-Min'] if self.uncertainties else ["NLO-P8"],
+            energy = self.energy,
+            nsample = 10000,
+        )
+        
+        self.model.add_production_2bodydecay(
+            pid0 = "-511",
+            pid1 = "130",
+            br = "5.6 * coupling**2 * pow(1.-pow(mass/5.279,2),2)",
+            generator = ['NLO-P8','NLO-P8-Max', 'NLO-P8-Min'] if self.uncertainties else ["NLO-P8"],
+            energy = self.energy,
+            nsample = 10000,
+        )
         
-        return self.decays()
+        self.model.add_production_2bodydecay(
+            pid0 = "521",
+            pid1 = "321",
+            br = "5.6 * coupling**2 * pow(1.-pow(mass/5.279,2),2)",
+            generator = ['NLO-P8','NLO-P8-Max', 'NLO-P8-Min'] if self.uncertainties else ["NLO-P8"],
+            energy = self.energy,
+            nsample = 10000,
+        )
+        
+        self.model.add_production_2bodydecay(
+            pid0 = "-521",
+            pid1 = "321",
+            br = "5.6 * coupling**2 * pow(1.-pow(mass/5.279,2),2)",
+            generator = ['NLO-P8','NLO-P8-Max', 'NLO-P8-Min'] if self.uncertainties else ["NLO-P8"],
+            energy = self.energy,
+            nsample = 10000,
+        )        
+
+        if self.uncertainties:        
+            self.prod_modes = {
+                '511':    ['NLO-P8', 'NLO-P8-max', 'NLO-P8-min'],
+                '-511':   ['NLO-P8', 'NLO-P8-max', 'NLO-P8-min'],
+                '521':    ['NLO-P8', 'NLO-P8-max', 'NLO-P8-min'],
+                '-521':   ['NLO-P8', 'NLO-P8-max', 'NLO-P8-min'],
+                }
+            self.weightnames = ["Nominal", "Powheg-Max", "Poweheg-Min"]
+        else:
+            self.prod_modes = {
+                "511"  : ['NLO-P8'],
+                "-511"  : ['NLO-P8'],
+                "521"  : ['NLO-P8'],
+                "-521" : ['NLO-P8'],                
+                }
+            self.weightnames = ["Nominal"]
+
 
+        modes = ["e_e", "mu_mu", "K_K", "pi_pi"]
+        self.model.set_br_1d(
+            modes=modes,
+            finalstates=[[11,-11], [13,-13], [321,-321], [211,-211]],
+            filenames=["model/br/"+mode+".txt" for mode in modes],
+        )
+
+        return
 
     def decays(self):
         # Set up liftime and BRs
 
-        if self.version == 1:
-            self.model.set_ctau_1d(
-                filename=f"files/models/{self.modelname}/ctau.txt", 
-                coupling_ref=1
-                )        
-        else:
-            self.model.set_ctau_1d(
-                filename=f"model/ctau.txt", 
-                coupling_ref=1
-                )
+        self.model.set_ctau_1d(
+            filename=f"model/ctau.txt", 
+            coupling_ref=1
+            )
             
         # Get LLP spectrum
         self.foresee.set_model(model=self.model)
@@ -349,10 +876,15 @@ class ForeseeGenerator(object):
 
         if self.version >= 3:
             coups, ctaus, nevents, momenta, weights = output
-            fmomenta  = flatten(momenta)
-            fweights  = flatten(weights)            
+            if self.version >= 4:
+                fmomenta = momenta
+                fweights = weights
+            else:
+                fmomenta  = flatten(momenta)
+                fweights  = flatten(weights)            
+
             fenergies = [p.e for p in fmomenta] # Now in MeV ?
-            fthetas   = [p.pt/p.pz for p in fmomenta]            
+            fthetas   = [(p.pt/p.pz if p.pz else 0) for p in fmomenta]            
 
             #self.plot(fthetas, fenergies, fweights)
 
@@ -405,7 +937,7 @@ class ForeseeGenerator(object):
 
     def events(self, lumi, f):
         f.write(f"{self.mass*1000} {self.couplings[0]} {self.nevents[0] * 1000 * lumi}\n")
-        print(f"{self.mass*1000} {self.couplings[0]} {self.nevents[0] * 1000 * lumi}")        
+        print(f"--> {self.mass*1000} {self.couplings[0]} {self.nevents[0] * 1000 * lumi}")        
         return
         
 
@@ -415,19 +947,23 @@ class ForeseeGenerator(object):
         energies, thetas, weights = self.data
 
         if self.outdir is None:
-            if self.version == 1:
-                self.outdir = f"files/models/{self.modelname}/events"
-            else:
-                self.outdir = f"{self.foresee.dirpath}/Models/{self.modelname}/model/events"
+            self.outdir = f"{self.foresee.dirpath}/Models/{self.modelname}/model/events"
 
         if not os.path.exists(self.outdir):
             os.mkdir(self.outdir)
 
-        if len(self.couplings) == 1:
-            filename = f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]}to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
-        else:
-            filename = f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
-
+        try:
+            if len(self.couplings) == 1:
+                filename = f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]}to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
+            else:
+                filename = f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_to_{self.daughter1_pid}_{self.daughter2_pid}.npy"
+        except AttributeError:
+            if len(self.couplings) == 1:
+                filename = f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]}to_{'+'.join(self.decay_modes)}.npy"
+            else:
+                filename = f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_to_{'+'.join(self.decay_modes)}.npy"
+            
+                
         print(f"Generated {len(thetas)} events")
         print(f"save data to file: {filename}")
         np.save(filename,[energies,thetas, weights])
@@ -446,10 +982,28 @@ class ForeseeGenerator(object):
         elif not os.path.exists(self.outdir):
             os.mkdir(self.outdir)
 
-        filename =  f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]:.1e}to_{self.daughter1_pid}_{self.daughter2_pid}{self.suffix}.hepmc"
+        try:
+            filename =  f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]:.1e}to_{self.daughter1_pid}_{self.daughter2_pid}{self.suffix}.hepmc"
+        except AttributeError:
+            filename =  f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]:.1e}to_{'+'.join(self.decay_modes)}{self.suffix}.hepmc"            
+
+        if self.version >= 4:
+            _, weights, _ = self.foresee.write_events(self.mass, self.couplings[0], self.energy, filename, nevents,
+                                                      zfront = self.zfront, # seed = self.seed, # seed now set differently
+                                                      modes = ({k: v for k,v in self.prod_modes.items()
+                                                                if (not self.allowed_prod_modes or k in self.allowed_prod_modes)}),
+                                                      weightnames = self.weightnames,
+                                                      notime = self.notime, t0 = self.t0, nsample = self.nbinsample,
+                                                      return_data = True)
+        else:
+            _, weights, _ = self.foresee.write_events(self.mass, self.couplings[0], self.energy, filename, nevents,
+                                                      zfront = self.zfront, seed = self.seed,
+                                                      modes = ({k: v for k,v in self.prod_modes.items()
+                                                                if (not self.allowed_prod_modes or k in self.allowed_prod_modes)}),
+                                                      weightnames = self.weightnames,
+                                                      notime = self.notime, t0 = self.t0, nsample = self.nbinsample,
+                                                      return_data = True)
 
-        _, weights, _ = self.foresee.write_events(self.mass, self.couplings[0], self.energy, filename, nevents, zfront = self.zfront, seed = self.seed, decaychannels = [self.mode],
-                                  notime = self.notime, t0 = self.t0, nsample = self.nbinsample, return_data = True)
 
         cfgname = f"{self.foresee.dirpath}/Models/{self.modelname}/" + filename.replace(".hepmc", ".cfg")
         print(f"save config to file: {cfgname}")
@@ -458,7 +1012,7 @@ class ForeseeGenerator(object):
 
         filename = f"{self.foresee.dirpath}/Models/{self.modelname}/{filename}"
 
-        return filename, weights
+        return filename, weights.T[0]
 
 
 def setup_foresee(path):
@@ -537,9 +1091,18 @@ def production(args):
         couplings = parse_couplings([coup], args.hepmc)
 
         for m in masses:
-            print(f">>> {nrun}: {data['model']} ({m} GeV) -> {data['pid1']} + {data['pid2']} with couplings = {couplings[0]:.2e} @ Ecom = {args.Ecom} TeV") 
+            if "modes" in data:
+                if isinstance(data['modes'], str):
+                    data['modes'] = [data['modes']]
+                print(f">>> {nrun}: {data['model']} ({m} GeV) -> {','.join(data['modes'])} with couplings = {couplings[0]:.2e} @ Ecom = {args.Ecom} TeV")
+            else:
+                print(f">>> {nrun}: {data['model']} ({m} GeV) -> {data['pid1']} + {data['pid2']} with couplings = {couplings[0]:.2e} @ Ecom = {args.Ecom} TeV") 
 
-            f = ForeseeGenerator(data["model"], args.Ecom, m, couplings, data["pid1"], data["pid2"], path = args.path, randomSeed = args.randomSeed, t0Shift = args.t0, zFront = args.zfront, notime = args.notime, suffix = data["name"])
+            f = ForeseeGenerator(data["model"], args.Ecom, m, couplings, data.get("pid1", None), data.get("pid2", None),
+                                 data.get("modes", None), path = args.path, randomSeed = args.randomSeed, t0Shift = args.t0, 
+                                 zFront = args.zfront, yshift = data.get("yshift", args.yshift),
+                                 notime = args.notime, suffix = data["name"], tracker_decays = data.get("tracker_decays", False),
+                                 uncertainties = data.get("uncertainties", False), allowed_prod_modes=args.production)
 
             with open("events.txt", "a") as fout:
                 f.events(40, fout)
@@ -556,15 +1119,20 @@ def production(args):
                     proc.run(["rm", filename])
 
             xsect[nrun] = {
-                "Description" : f"{data['model']} ({m} GeV) -> {data['pid1']} + {data['pid2']} with couplings = {couplings[0]:.2e} @ Ecom = {args.Ecom} TeV",
                 "Mass (GeV)" : m,
                 "Coupling" : f"{couplings[0]:e}",
-                "Cross-section (pb)" : sum(weights[0]), # / data["nevents"],
+                "Cross-section (pb)" : sum(weights), # / data["nevents"],
                 "ECom (TeV)" : float(args.Ecom),
-                "NEvents" : float(data["nevents"])                
+                "NEvents" : float(data["nevents"]),
+                "WtNames" : f.weightnames,
                 }
 
-            print(xsect)
+            if "modes" in data:
+                xsect[nrun]["Description"] = f"{data['model']} ({m} GeV) -> {','.join(data['modes'])} with couplings = {couplings[0]:.2e} @ Ecom = {args.Ecom} TeV",                            
+            else:
+                xsect[nrun]["Description"] = f"{data['model']} ({m} GeV) -> {data['pid1']} + {data['pid2']} with couplings = {couplings[0]:.2e} @ Ecom = {args.Ecom} TeV",            
+
+            print("X-sect", xsect)
             with open (outdir + "/cross-sections.json", "w") as f:
                 json.dump(xsect, f, indent = 4)
             
@@ -585,6 +1153,7 @@ if __name__ == "__main__":
     parser.add_argument("--couplings", "-c", required = False, nargs = "+", help = "Couplings of mother (either single/mulitple values or tuple to pass to np.logspace)")
     parser.add_argument("--pid1", default = None, required = False, type = int, help = "PID of daughter 1")
     parser.add_argument("--pid2", default = None, type = int, help = "PID of daughter 2 (if not set then will be -PID1)")
+    parser.add_argument("--modes", default = None, required = False, nargs = "*", help = "Space-separated list of foresee decay modes")   
     parser.add_argument("--Ecom", default = "13.6", help = "Center of mass energy [TeV]")
     parser.add_argument("--outdir", "-o", default = None, help = "Output path")    
     parser.add_argument("--path", default = ".", help = "Path to foresee installation")
@@ -595,6 +1164,12 @@ if __name__ == "__main__":
     parser.add_argument("--zfront", "-z", default = -1.5, help = "Minimum FASER z coordinate for decays (in meters)")
     parser.add_argument("--notime", action = "store_true", help = "Set all vertex times to 0 rather than calculating from start of decay volume")
     parser.add_argument("--suffix", default = "", help = "Filename suffix")
+    parser.add_argument("--tracker", default = False, action = "store_true", help = "Include decays in the tracker")    
+    parser.add_argument("--uncertainties", "-u", default = False, action = "store_true", help = "Store weights for uncertainties")    
+    parser.add_argument("--yshift", "-y", default = -0.065, 
+                        help = "Shift due to crossing angle (minus any detector offset) in m with negative being down")
+    parser.add_argument("--production", default = None, required = False, nargs = "*", help = "Space-separated list of foresee production modes")   
+    
     args = parser.parse_args()
 
     add_to_python_path(f"{args.path}/src")
@@ -616,7 +1191,10 @@ if __name__ == "__main__":
         print(f"   decay = {args.pid1} {args.pid2}")
         print(f"   couplings = {couplings}")    
         
-        f = ForeseeGenerator(args.model, args.Ecom, args.mass, couplings, args.pid1, args.pid2, outdir = args.outdir, path = args.path, randomSeed = args.randomSeed, t0Shift = args.t0, zFront = args.zfront, notime = args.notime, suffix = args.suffix)
+        f = ForeseeGenerator(args.model, args.Ecom, args.mass, couplings, args.pid1, args.pid2, args.modes,
+                             outdir = args.outdir, path = args.path, randomSeed = args.randomSeed, t0Shift = args.t0, zFront = args.zfront,
+                             notime = args.notime, suffix = args.suffix, tracker_decays = args.tracker, yshift = args.yshift,
+                             uncertainties = args.uncertainties, allowed_prod_modes=args.production)
         
         if args.hepmc:
             f.write_hepmc(args.nevents)
diff --git a/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid.json b/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid.json
new file mode 100644
index 0000000000000000000000000000000000000000..d717f91e17352dd50d249f1554f7d04d011d4d60
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid.json
@@ -0,0 +1,27 @@
+{
+"name" : "ALPWGrid13p6_65mm_110mm_3p5m_fullgrid",
+"model" : "ALP-W",
+"pid1" : 22,
+"pid2" : 22,
+"nevents" : 50000,
+"nrun" : 110210,
+"samples" : {
+	"3e-03"   : [0.06],
+	"2e-03"   : [0.06],
+	"1.1e-03" : [      0.08],
+	"1e-03"   : [0.06, 0.08, 0.1],
+	"6e-04"   : [0.06,       0.1, 0.12],
+	"3e-04"   : [            0.1, 0.12, 0.14],  
+	"2e-04"   : [      0.08,            0.14, 0.2],
+	"1.1e-04" : [0.06,       0.1],
+	"1e-04"   : [                 0.12, 0.14, 0.2], 
+	"6e-05"   : [0.06,       0.1,             0.2, 0.3],
+	"4e-05"   : [      0.08,            0.14, 0.2, 0.3],
+	"3e-05"   : [0.06, 0.08, 0.1, 0.12,                 0.4],
+	"2e-05"   : [            0.1, 0.12, 0.14, 0.2, 0.3, 0.4, 0.43],
+	"1.1e-05" : [                       0.14,                0.43, 0.5],
+	"1e-05"   : [                             0.2, 0.3, 0.4,       0.5],
+	"7e-06"   : [                                  0.3, 0.4, 0.43]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid_nlo.json b/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid_nlo.json
new file mode 100644
index 0000000000000000000000000000000000000000..aa01c89bd04cbf939bf32030a018e6225da3f253
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid_nlo.json
@@ -0,0 +1,27 @@
+{
+"name" : "ALPWGrid13p6_65mm_110mm_3p5m_fullgrid_nlo",
+"model" : "ALP-W",
+"pid1" : 22,
+"pid2" : 22,
+"nevents" : 50000,
+"nrun" : 110210,
+"samples" : {
+	"3e-03"   : [0.06],
+	"2e-03"   : [0.06],
+	"1.1e-03" : [      0.08],
+	"1e-03"   : [0.06, 0.08, 0.1],
+	"6e-04"   : [0.06,       0.1, 0.12],
+	"3e-04"   : [            0.1, 0.12, 0.14],  
+	"2e-04"   : [      0.08,            0.14, 0.2],
+	"1.1e-04" : [0.06,       0.1],
+	"1e-04"   : [                 0.12, 0.14, 0.2], 
+	"6e-05"   : [0.06,       0.1,             0.2, 0.3],
+	"4e-05"   : [      0.08,            0.14, 0.2, 0.3],
+	"3e-05"   : [0.06, 0.08, 0.1, 0.12,                 0.4],
+	"2e-05"   : [            0.1, 0.12, 0.14, 0.2, 0.3, 0.4, 0.43],
+	"1.1e-05" : [                       0.14,                0.43, 0.5],
+	"1e-05"   : [                             0.2, 0.3, 0.4,       0.5],
+	"7e-06"   : [                                  0.3, 0.4, 0.43]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid_nlo_uncertainties.json b/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid_nlo_uncertainties.json
new file mode 100644
index 0000000000000000000000000000000000000000..b3dbafac46e60f939b7e1e36eaa75281dc8fcf2e
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid_nlo_uncertainties.json
@@ -0,0 +1,28 @@
+{
+"name" : "ALPWGrid13p6_65mm_110mm_3p5m_fullgrid_nlo_uncertainties",
+"model" : "ALP-W",
+"pid1" : 22,
+"pid2" : 22,
+"nevents" : 50000,
+"nrun" : 110210,
+"uncertainties" : true,
+"samples" : {
+	"3e-03"   : [0.06],
+	"2e-03"   : [0.06],
+	"1.1e-03" : [      0.08],
+	"1e-03"   : [0.06, 0.08, 0.1],
+	"6e-04"   : [0.06,       0.1, 0.12],
+	"3e-04"   : [            0.1, 0.12, 0.14],  
+	"2e-04"   : [      0.08,            0.14, 0.2],
+	"1.1e-04" : [0.06,       0.1],
+	"1e-04"   : [                 0.12, 0.14, 0.2], 
+	"6e-05"   : [0.06,       0.1,             0.2, 0.3],
+	"4e-05"   : [      0.08,            0.14, 0.2, 0.3],
+	"3e-05"   : [0.06, 0.08, 0.1, 0.12,                 0.4],
+	"2e-05"   : [            0.1, 0.12, 0.14, 0.2, 0.3, 0.4, 0.43],
+	"1.1e-05" : [                       0.14,                0.43, 0.5],
+	"1e-05"   : [                             0.2, 0.3, 0.4,       0.5],
+	"7e-06"   : [                                  0.3, 0.4, 0.43]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid_nlo_uncertainties_extra.json b/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid_nlo_uncertainties_extra.json
new file mode 100644
index 0000000000000000000000000000000000000000..b3d438eb40f39c15dad4a225d731804af4d4e23f
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_ALP_W_fullgrid_nlo_uncertainties_extra.json
@@ -0,0 +1,17 @@
+{
+"name" : "ALPWGrid13p6_65mm_110mm_3p5m_fullgrid_nlo_uncertainties_extra",
+"model" : "ALP-W",
+"pid1" : 22,
+"pid2" : 22,
+"nevents" : 50000,
+"nrun" : 110265,
+"uncertainties" : true,
+"samples" : {
+	"1e-04"   : [0.26],
+	"6e-05"   : [0.26, 0.32],
+	"4e-05"   : [0.26, 0.32, 0.36],
+	"2e-05"   : [0.26, 0.32, 0.36],
+	"1e-05"   : [0.26, 0.32, 0.36]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_ALP_W_larger_decay_volume.json b/Generators/ForeseeGenerator/share/points_ALP_W_larger_decay_volume.json
new file mode 100644
index 0000000000000000000000000000000000000000..d77938b77b9ed093b032c213381795e16b1b9eb8
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_ALP_W_larger_decay_volume.json
@@ -0,0 +1,23 @@
+{
+"name" : "ALPWGrid13p6_65mm_150mm_483p7m_1p3m",
+"model" : "ALP-W",
+"pid1" : 22,
+"pid2" : 22,
+"nevents" : 50000,
+"nrun" : 110265,
+"uncertainties" : true,
+"samples" : {
+	"6e-04"   : [0.12],
+	"3e-04"   : [0.12, 0.14],  
+	"2e-04"   : [      0.14, 0.2],
+	"1e-04"   : [0.12, 0.14, 0.2], 
+	"6e-05"   : [            0.2, 0.3],
+	"4e-05"   : [      0.14, 0.2, 0.3],
+	"3e-05"   : [                      0.4],
+	"2e-05"   : [            0.2, 0.3, 0.4, 0.43],
+	"1.1e-05" : [                           0.43, 0.5],
+	"1e-05"   : [                 0.3, 0.4,       0.5],
+	"7e-06"   : [                 0.3, 0.4, 0.43]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_ALP_gluon_uncertainties.json b/Generators/ForeseeGenerator/share/points_ALP_gluon_uncertainties.json
new file mode 100644
index 0000000000000000000000000000000000000000..ad7de32ceef4dd3d62df99b3f80c87871962c474
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_ALP_gluon_uncertainties.json
@@ -0,0 +1,60 @@
+{
+"name" : "ALPGluonGrid3p6_65mm_110mm_3p5m_nlo_uncertainties",
+"model" : "ALP-g",
+"modes" : ["gamma_gamma", "br_3Pi0", "br_Pi+Pi-Gamma", "br_Pi+Pi-Pi0", "br_2PiEta"],
+"nevents" : 50000,
+"nrun" : 111000,
+"tracker_decays" : true,
+"uncertainties" : true,
+"samples" : {
+	"0.0001": [0.05],
+	"0.00015": [0.3, 0.4],
+	"0.0002": [0.125, 0.14, 0.15, 0.16],
+	"0.00025": [0.05],
+	"0.0003": [0.25, 0.3],
+	"0.0006": [0.125, 0.2, 0.25],
+	"0.001": [0.2],
+	"0.0015": [0.1, 0.125],
+	"0.002": [0.14, 0.15, 0.16],
+	"0.003": [0.075, 0.1],
+	"0.004": [0.14, 0.15, 0.16],
+	"0.006": [0.05, 0.075],
+	"0.01": [0.15, 0.16],
+	"0.015": [0.05],
+	"0.04": [0.15],
+	"1.5e-05": [0.1,
+            0.125,
+            0.135,
+            0.14,
+            0.15,
+            0.16,
+            0.3,
+            0.5,
+            0.6,
+            0.65,
+            0.7,
+            0.75,
+            0.8,
+            0.85],
+	"1.5e-06": [0.85, 0.9],
+	"1e-05": [0.2, 0.525, 0.575, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9],
+	"1e-06": [0.5, 0.525, 0.547, 0.575, 0.6],
+	"2.5e-05": [0.5, 0.65, 0.7, 0.75, 0.8],
+	"2.5e-06": [0.5],
+	"2e-05": [0.2, 0.25],
+	"2e-06": [0.125, 0.135, 0.14, 0.525, 0.575, 0.6, 0.65, 0.7, 0.75, 0.8],
+	"2e-07": [0.525, 0.547, 0.575],
+	"3e-05": [0.125, 0.135, 0.14, 0.15, 0.16],
+	"4e-05": [0.075],
+	"4e-06": [0.4, 0.85, 0.9],
+	"5e-05": [0.1],
+	"5e-06": [0.525, 0.575, 0.6, 0.65, 0.7, 0.75, 0.8],
+	"5e-07": [0.525, 0.547, 0.575],
+	"6e-06": [0.125, 0.135, 0.14, 0.5],
+	"7e-05": [0.4],
+	"7e-06": [0.3],
+	"8e-05": [0.075],
+	"8e-06": [0.25, 0.4]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_ALP_photon_fullgrid_nlo_uncertainties.json b/Generators/ForeseeGenerator/share/points_ALP_photon_fullgrid_nlo_uncertainties.json
new file mode 100644
index 0000000000000000000000000000000000000000..53757d49e18453dfb476e2d9fa46a8eee38031b6
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_ALP_photon_fullgrid_nlo_uncertainties.json
@@ -0,0 +1,29 @@
+{
+"name" : "ALPPhotonGrid13p6_65mm_110mm_3p5m_fullgrid_nlo_uncertainties_fixdist",
+"model" : "ALP-photon",
+"pid1" : 22,
+"pid2" : 22,
+"nevents" : 50000,
+"nrun" : 110900,
+"uncertainties" : true,
+"samples" : {
+	"1e-02"   : [0.01, 0.013],
+	"6e-03"   : [      0.013, 0.02],
+	"4e-03"   : [0.01,        0.02],
+	"2e-03"   : [0.01, 0.013, 0.02, 0.03],
+	"1.3e-03" : [                   0.03, 0.04],
+	"1e-03"   : [                         0.04],
+	"8e-04"   : [0.01, 0.013, 0.02, 0.03, 0.04],
+	"6e-04"   : [                         0.04, 0.05],
+	"4e-04"   : [0.01, 0.013, 0.02, 0.03, 0.04, 0.05, 0.06],
+	"3e-04"   : [                                     0.06, 0.07],
+	"2e-04"   : [0.01, 0.013, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07],
+	"1.3e-04" : [0.01, 0.013, 0.02,                   0.06, 0.07, 0.08],
+	"9e-05"   : [             0.02, 0.03,             0.06, 0.07, 0.08, 0.09],
+	"7e-05"   : [                   0.03, 0.04],
+	"6e-05"   : [                               0.05, 0.06, 0.07, 0.08, 0.09],
+	"5e-05"   : [                         0.04],
+	"4e-05"   : [                               0.05, 0.06, 0.07, 0.08, 0.09]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_BminusL_extra.json b/Generators/ForeseeGenerator/share/points_BminusL_extra.json
new file mode 100644
index 0000000000000000000000000000000000000000..8438c4f0406923c03ed6a4b4f54abc6d43d74675
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_BminusL_extra.json
@@ -0,0 +1,25 @@
+{
+"name" : "U1B-LGrid13p6_65mm_110mm_extra",
+"model" : "U(1)B-L",
+"pid1" : 11,
+"pid2" : -11,
+"nevents" : 50000,
+"nrun" : 110329,
+"samples" : {
+	"5e-5"   : [0.01],
+	"4e-5"   : [      0.011],
+	"3e-5"   : [      0.011, 0.014, 0.015],
+	"2e-5"   : [      0.011,        0.015, 0.02],
+	"1.4e-5" : [                    0.015,                    0.034],
+	"1e-5"   : [                                       0.03],
+	"8e-6"   : [                                0.024,        0.034],	
+	"7e-6"   : [                                0.024,        0.034],	
+	"6e-6"   : [                                0.024,        0.034, 0.04],	
+	"5e-6"   : [                                0.024,               0.04, 0.044],	
+	"4e-6"   : [                                0.024,                     0.044], 
+	"3e-6"   : [                                                           0.044],	
+	"2.4e-6" : [                                                     0.04, 0.044],	
+	"2e-6"   : [0.01,        0.014,       0.02, 0.024, 0.03,  0.034]	
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_DP22_fullgrid.json b/Generators/ForeseeGenerator/share/points_DP22_fullgrid.json
new file mode 100644
index 0000000000000000000000000000000000000000..6bcafaad67b94c47b5a039d6c6f3baa8f0e4046e
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_DP22_fullgrid.json
@@ -0,0 +1,21 @@
+{
+    "name" : "DarkPhotonGrid13p6_2022_110mm_tracker_fullgrid",
+    "model" : "DarkPhoton",
+    "modes" : ["e_e", "mu_mu", "pi+_pi-"],
+    "nevents" : 20000,
+    "nrun" : 113000,
+    "yshift" : -0.065, 
+    "tracker_decays" : true,
+    "uncertainties" : true,
+    "samples" : {
+        "3e-04"  : [0.01, 0.013, 0.016, 0.025],
+        "1.5e-04": [0.01, 0.013, 0.016, 0.025, 0.032, 0.036, 0.05],
+        "5e-05"  : [0.01,        0.016, 0.025, 0.032, 0.036, 0.05, 0.063, 0.079],
+        "3e-05"  : [0.01,               0.025, 0.032, 0.036, 0.05, 0.063, 0.079, 0.1, 0.11], 
+	"2e-05"  : [                                         0.05, 0.063, 0.079, 0.1, 0.11, 0.2, 0.22, 0.3],
+        "1e-05"  : [0.01, 0.013, 0.016, 0.025, 0.032, 0.036, 0.05, 0.063, 0.079, 0.1, 0.11, 0.2, 0.22, 0.3],
+	"5e-06"  : [0.01, 0.013, 0.016, 0.025, 0.032, 0.036, 0.05, 0.063, 0.079, 0.1, 0.11, 0.2, 0.22, 0.3],
+	"3e-06"  : [0.01, 0.013, 0.016, 0.025, 0.032, 0.036, 0.05, 0.063, 0.079, 0.1, 0.11, 0.2, 0.22, 0.3],
+	"1e-06"  : [0.01, 0.013, 0.016, 0.025, 0.032, 0.036, 0.05, 0.063, 0.079, 0.1, 0.11, 0.2]
+    } 
+}
diff --git a/Generators/ForeseeGenerator/share/points_DP24_fullgrid.json b/Generators/ForeseeGenerator/share/points_DP24_fullgrid.json
new file mode 100644
index 0000000000000000000000000000000000000000..61c6c60dd03d6fd192446018b9a22a01b6dc83d7
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_DP24_fullgrid.json
@@ -0,0 +1,21 @@
+{
+    "name" : "DarkPhotonGrid13p6_2024_110mm_tracker_fullgrid",
+    "model" : "DarkPhoton",
+    "modes" : ["e_e", "mu_mu", "pi+_pi-"],
+    "nevents" : 30000,
+    "nrun" : 114000,
+    "yshift" : 0.089, 
+    "tracker_decays" : true,
+    "uncertainties" : true,
+    "samples" : {
+        "3e-04"  : [0.01, 0.013, 0.016, 0.025],
+        "1.5e-04": [0.01, 0.013, 0.016, 0.025, 0.032, 0.036, 0.05],
+        "5e-05"  : [0.01,        0.016, 0.025, 0.032, 0.036, 0.05, 0.063, 0.079],
+        "3e-05"  : [0.01,               0.025, 0.032, 0.036, 0.05, 0.063, 0.079, 0.1, 0.11], 
+	"2e-05"  : [                                         0.05, 0.063, 0.079, 0.1, 0.11, 0.2, 0.22, 0.3],
+        "1e-05"  : [0.01, 0.013, 0.016, 0.025, 0.032, 0.036, 0.05, 0.063, 0.079, 0.1, 0.11, 0.2, 0.22, 0.3],
+	"5e-06"  : [0.01, 0.013, 0.016, 0.025, 0.032, 0.036, 0.05, 0.063, 0.079, 0.1, 0.11, 0.2, 0.22, 0.3],
+	"3e-06"  : [0.01, 0.013, 0.016, 0.025, 0.032, 0.036, 0.05, 0.063, 0.079, 0.1, 0.11, 0.2, 0.22, 0.3],
+	"1e-06"  : [0.01, 0.013, 0.016, 0.025, 0.032, 0.036, 0.05, 0.063, 0.079, 0.1, 0.11, 0.2]
+    } 
+}
diff --git a/Generators/ForeseeGenerator/share/points_DP_fullgrid.json b/Generators/ForeseeGenerator/share/points_DP_fullgrid.json
new file mode 100644
index 0000000000000000000000000000000000000000..b3060e95d1c48e971fd7bd8cd4ce8882867a976a
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_DP_fullgrid.json
@@ -0,0 +1,22 @@
+{
+    "name" : "DarkPhotonGrid13p6_110mm_5m_fullgrid",
+    "model" : "DarkPhoton",
+    "pid1" : 11,
+    "pid2" : -11,
+    "nevents" : 50000,
+    "nrun" : 112000,
+    "tracker_decays" : true,
+    "uncertainties" : true,
+    "samples" : {
+        "3e-04"  : [0.01, 0.0251, 0.0126, 0.0158],
+        "1.5e-04": [0.01, 0.0126, 0.0158, 0.0251],
+        "1e-04"  : [0.01, 0.0251, 0.0501, 0.0126, 0.0158],    
+        "5e-05"  : [0.01, 0.0251, 0.0501, 0.0631, 0.0158, 0.0316, 0.0358],
+        "3e-05"  : [0.01, 0.0251, 0.0501, 0.0631, 0.0794, 0.1, 0.1122, 0.0316, 0.0358],
+        "2e-05"  : [0.0501, 0.0631, 0.0794],
+        "1e-05"  : [0.01, 0.0251, 0.0501, 0.0631, 0.0794, 0.1, 0.1122, 0.0126, 0.0158, 0.0316, 0.0358],
+        "5e-06"  : [0.01, 0.0251, 0.0501, 0.0631, 0.0794, 0.1, 0.1122, 0.0126, 0.0158, 0.0316, 0.0358],
+        "3e-06"  : [0.01, 0.0126, 0.0158, 0.0251, 0.0316, 0.0358, 0.0501, 0.0631, 0.0794],
+        "1e-06"  : [0.01, 0.0251, 0.0501, 0.0631, 0.0794]
+    } 
+}
diff --git a/Generators/ForeseeGenerator/share/points_DP_incTracker.json b/Generators/ForeseeGenerator/share/points_DP_incTracker.json
new file mode 100644
index 0000000000000000000000000000000000000000..595f790e032592fa7fe2c2666b8a534317b7ea14
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_DP_incTracker.json
@@ -0,0 +1,14 @@
+{
+    "name" : "DarkPhotonGrid13p6_110mm_5m",
+    "model" : "DarkPhoton",
+    "pid1" : 11,
+    "pid2" : -11,
+    "nevents" : 50000,
+    "nrun" : 112000,
+    "tracker_decays" : true,
+    "uncertainties" : false,
+    "samples" : {
+        "1e-04"  : [0.01],
+        "1e-05"  : [0.01, 0.1]
+    } 
+}
diff --git a/Generators/ForeseeGenerator/share/points_DarkHiggs22.json b/Generators/ForeseeGenerator/share/points_DarkHiggs22.json
new file mode 100644
index 0000000000000000000000000000000000000000..965a631929121e045ab4f6dd8bd57fdccd9b8b63
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_DarkHiggs22.json
@@ -0,0 +1,17 @@
+{
+    "name" : "DarkHiggsGrid13p6_2022_110mm_tracker",
+    "model" : "DarkHiggs",
+    "modes" : ["mu_mu", "pi_pi"],
+    "nevents" : 20000,
+    "nrun" : 113900,
+    "yshift" : -0.065, 
+    "tracker_decays" : true,
+    "uncertainties" : false,
+    "samples" : {
+        "7e-04"  : [0.23, 0.26, 0.3, 0.33, 0.4],
+        "6e-04"  : [0.23, 0.26, 0.3, 0.33, 0.4],
+        "5e-04"  : [      0.26, 0.3, 0.33, 0.4],
+        "4e-04"  : [            0.3, 0.33, 0.4],
+        "3e-04"  : [                 0.33     ]
+    } 
+}
diff --git a/Generators/ForeseeGenerator/share/points_DarkHiggs24.json b/Generators/ForeseeGenerator/share/points_DarkHiggs24.json
new file mode 100644
index 0000000000000000000000000000000000000000..b54bdc5530482bf60b72927bd11e24a6201e5a02
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_DarkHiggs24.json
@@ -0,0 +1,17 @@
+{
+    "name" : "DarkHiggsGrid13p6_2024_110mm_tracker",
+    "model" : "DarkHiggs",
+    "modes" : ["mu_mu", "pi_pi"],
+    "nevents" : 30000,
+    "nrun" : 114900,
+    "yshift" : 0.089, 
+    "tracker_decays" : true,
+    "uncertainties" : false,
+    "samples" : {
+        "7e-04"  : [0.23, 0.26, 0.3, 0.33, 0.4],
+        "6e-04"  : [0.23, 0.26, 0.3, 0.33, 0.4],
+        "5e-04"  : [      0.26, 0.3, 0.33, 0.4],
+        "4e-04"  : [            0.3, 0.33, 0.4],
+        "3e-04"  : [                 0.33     ]
+    } 
+}
diff --git a/Generators/ForeseeGenerator/share/points_U1B.json b/Generators/ForeseeGenerator/share/points_U1B.json
new file mode 100644
index 0000000000000000000000000000000000000000..5c5324abb25ad2005c800f553fbb4e16b0166add
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_U1B.json
@@ -0,0 +1,14 @@
+{
+"name" : "U1BGrid3p6_65mm_110mm_3p5m",
+"model" : "U(1)B",
+"pid1" : 111,
+"pid2" : 22,
+"nevents" : 50000,
+"nrun" : 110400,
+"samples" : {
+	"1e-04" : [0.3],
+	"5e-05" : [0.4],
+	"3e-05" : [0.5]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_U1B_extra.json b/Generators/ForeseeGenerator/share/points_U1B_extra.json
new file mode 100644
index 0000000000000000000000000000000000000000..48afc408a49a30fa40e6065f0fb468d1bb18f8f2
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_U1B_extra.json
@@ -0,0 +1,13 @@
+{
+"name" : "U1BGrid3p6_65mm_110mm_3p5m",
+"model" : "U(1)B",
+"pid1" : 111,
+"pid2" : 22,
+"nevents" : 50000,
+"nrun" : 110403,
+"samples" : {
+	"3e-04" : [0.2],
+	"3e-05" : [0.2, 0.3]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_U1B_fullgrid_nlo.json b/Generators/ForeseeGenerator/share/points_U1B_fullgrid_nlo.json
new file mode 100644
index 0000000000000000000000000000000000000000..8e846bcd6c6f257efb94843c6e038b1b93a7b7c9
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_U1B_fullgrid_nlo.json
@@ -0,0 +1,24 @@
+{
+"name" : "U1BGrid3p6_65mm_110mm_3p5m_fullgrid_nlo",
+"model" : "U(1)B",
+"modes" : ["ee", "pigamma"],
+"nevents" : 50000,
+"nrun" : 110410,
+"tracker_decays" : true,
+"samples" : {
+	"1e-02" : [0.15],
+	"3e-03" : [0.15, 0.2],
+	"2e-03" : [      0.2],
+	"1e-03" : [           0.23],
+	"5e-04" : [0.15, 0.2, 0.23, 0.3],
+	"3e-04" : [                 0.3, 0.33],
+	"2e-04" : [                      0.33, 0.4],
+	"1e-04" : [0.15, 0.2, 0.23, 0.3, 0.33, 0.4, 0.43],
+	"6e-05" : [                                       0.5],
+	"4e-05" : [0.15,                       0.4, 0.43, 0.5, 0.53],
+	"2e-05" : [0.15, 0.2, 0.23, 0.3,       0.4,       0.5, 0.53, 0.6],
+	"1e-05" : [      0.2, 0.23, 0.3, 0.33, 0.4, 0.43, 0.5, 0.53, 0.6],
+	"6e-06" : [                 0.3, 0.33, 0.4, 0.43, 0.5, 0.53, 0.6]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_U1B_fullgrid_nlo_uncertainties.json b/Generators/ForeseeGenerator/share/points_U1B_fullgrid_nlo_uncertainties.json
new file mode 100644
index 0000000000000000000000000000000000000000..0f0e9f1488889831adcc1e74271709e32b8e9f57
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_U1B_fullgrid_nlo_uncertainties.json
@@ -0,0 +1,25 @@
+{
+"name" : "U1BGrid3p6_65mm_110mm_3p5m_fullgrid_nlo_uncertainties",
+"model" : "U(1)B",
+"modes" : ["ee", "pigamma", "3pi"],
+"nevents" : 50000,
+"nrun" : 110410,
+"uncertainties" : true,
+"tracker_decays" : true,
+"samples" : {
+	"1e-02" : [0.15],
+	"3e-03" : [0.15, 0.2],
+	"2e-03" : [      0.2],
+	"1e-03" : [           0.23],
+	"5e-04" : [0.15, 0.2, 0.23, 0.3],
+	"3e-04" : [                 0.3, 0.33],
+	"2e-04" : [                      0.33, 0.4],
+	"1e-04" : [0.15, 0.2, 0.23, 0.3, 0.33, 0.4, 0.43],
+	"6e-05" : [                                       0.5],
+	"4e-05" : [0.15,                       0.4, 0.43, 0.5, 0.53],
+	"2e-05" : [0.15, 0.2, 0.23, 0.3,       0.4,       0.5, 0.53, 0.6],
+	"1e-05" : [      0.2, 0.23, 0.3, 0.33, 0.4, 0.43, 0.5, 0.53, 0.6],
+	"6e-06" : [                 0.3, 0.33, 0.4, 0.43, 0.5, 0.53, 0.6]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_UpPhilic.json b/Generators/ForeseeGenerator/share/points_UpPhilic.json
new file mode 100644
index 0000000000000000000000000000000000000000..51ba8e94610c575e8f251c22f36829316df565e8
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_UpPhilic.json
@@ -0,0 +1,12 @@
+{
+"name" : "UpPhilicGrid3p6_65mm_110mm_3p5m",
+"model" : "UpPhilic",
+"pid1" : 111,
+"pid2" : 111,
+"nevents" : 50000,
+"nrun" : 110500,
+"samples" : {
+	"5e-08" : [0.28]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_UpPhilicCharged2022_fullgrid.json b/Generators/ForeseeGenerator/share/points_UpPhilicCharged2022_fullgrid.json
new file mode 100644
index 0000000000000000000000000000000000000000..b04275eea474351e7e539ff62d04b56540fd9c79
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_UpPhilicCharged2022_fullgrid.json
@@ -0,0 +1,18 @@
+{
+"name" : "UpPhilicChargedGrid3p6_2022_110mm_3p5m_tracker_fullgrid",
+"model" : "UpPhilic",
+"modes" : ["pi+_pi-"],
+"nevents" : 20000,
+"nrun" : 113500,
+"yshift" : -0.065, 
+"tracker_decays" : true,
+"uncertainties" : true,
+"samples" : {
+	"1e-07" : [0.28, 0.3],  
+	"5e-08" : [0.28, 0.3, 0.303],  
+	"4e-08" : [                  0.31],  
+	"3e-08" : [0.28, 0.3, 0.303, 0.31],  
+	"2e-08" : [0.28, 0.3, 0.303, 0.31]  
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_UpPhilicCharged2024_fullgrid.json b/Generators/ForeseeGenerator/share/points_UpPhilicCharged2024_fullgrid.json
new file mode 100644
index 0000000000000000000000000000000000000000..54a3f4cf429f2019bd8b4cbe48df88341d1d4b5f
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_UpPhilicCharged2024_fullgrid.json
@@ -0,0 +1,18 @@
+{
+"name" : "UpPhilicChargedGrid3p6_2024_110mm_3p5m_tracker_fullgrid",
+"model" : "UpPhilic",
+"modes" : ["pi+_pi-"],
+"nevents" : 30000,
+"nrun" : 114500,
+"yshift" : 0.089, 
+"tracker_decays" : true,
+"uncertainties" : true,
+"samples" : {
+	"1e-07" : [0.28, 0.3],  
+	"5e-08" : [0.28, 0.3, 0.303],  
+	"4e-08" : [                  0.31],  
+	"3e-08" : [0.28, 0.3, 0.303, 0.31],  
+	"2e-08" : [0.28, 0.3, 0.303, 0.31]  
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_UpPhilic_extra.json b/Generators/ForeseeGenerator/share/points_UpPhilic_extra.json
new file mode 100644
index 0000000000000000000000000000000000000000..ed3b93212da6875fc1b1e8c68d20137d4a35471b
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_UpPhilic_extra.json
@@ -0,0 +1,17 @@
+{
+"name" : "UpPhilicGrid3p6_65mm_110mm_3p5m_extra_uncertainties",
+"model" : "UpPhilic",
+"modes" : ["gamma_gamma", "pi0_pi0"],
+"nevents" : 50000,
+"nrun" : 110527,
+"tracker_decays" : true,
+"uncertainties" : true,
+"samples" : {
+    "4e-07" : [0.275 ],
+    "2e-07" : [0.275, 0.28],
+    "1e-07" : [0.275, 0.285, 0.29, 0.295],
+    "5e-08" : [0.275, 0.2875, 0.295, 0.303, 0.31],
+    "3e-08" : [0.275, 0.2875, 0.295]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_UpPhilic_fullgrid_uncertainties.json b/Generators/ForeseeGenerator/share/points_UpPhilic_fullgrid_uncertainties.json
new file mode 100644
index 0000000000000000000000000000000000000000..ad5f6d32ed29a24913c90fbefb01b128635ab141
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_UpPhilic_fullgrid_uncertainties.json
@@ -0,0 +1,20 @@
+{
+"name" : "UpPhilicGrid3p6_65mm_110mm_3p5m_fullgrid_uncertainties",
+"model" : "UpPhilic",
+"modes" : ["gamma_gamma", "pi0_pi0"],
+"nevents" : 50000,
+"nrun" : 110502,
+"tracker_decays" : true,
+"uncertainties" : true,
+"samples" : {
+	"4e-07" : [0.25, 0.26],  
+	"3e-07" : [            0.27],  
+	"2e-07" : [0.25, 0.26, 0.27],  
+	"1e-07" : [            0.27, 0.28, 0.3],  
+	"5e-08" : [0.25, 0.26,       0.28, 0.3, 0.303],  
+	"4e-08" : [                                    0.31],  
+	"3e-08" : [                  0.28, 0.3, 0.303, 0.31],  
+	"2e-08" : [0.25, 0.26,       0.28, 0.3, 0.303, 0.31]  
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_UpPhilic_incTracker.json b/Generators/ForeseeGenerator/share/points_UpPhilic_incTracker.json
new file mode 100644
index 0000000000000000000000000000000000000000..456cb2b008aba69d0e39800bf7ba58d7566084fa
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_UpPhilic_incTracker.json
@@ -0,0 +1,12 @@
+{
+"name" : "UpPhilicGrid3p6_65mm_110mm_3p5m_inTracker",
+"model" : "UpPhilic",
+"pid1" : 111,
+"pid2" : 111,
+"nevents" : 50001,
+"nrun" : 110501,
+"samples" : {
+	"5e-08" : [0.28]
+   }
+}
+
diff --git a/Generators/GeneratorUtils/CMakeLists.txt b/Generators/GeneratorUtils/CMakeLists.txt
index d471fbe183ae3da80c2a0ce1a49785ca0b907387..ee21848694bce3abefd48ae794a4299452b17379 100644
--- a/Generators/GeneratorUtils/CMakeLists.txt
+++ b/Generators/GeneratorUtils/CMakeLists.txt
@@ -8,3 +8,4 @@ atlas_subdir( GeneratorUtils )
 # Install files from the package:
 atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
 
+atlas_install_joboptions( share/*.py )
diff --git a/Generators/GeneratorUtils/python/HepMCFixes.py b/Generators/GeneratorUtils/python/HepMCFixes.py
new file mode 100644
index 0000000000000000000000000000000000000000..49ee98b7fc0e89843e34aa2725563d7bf8656c39
--- /dev/null
+++ b/Generators/GeneratorUtils/python/HepMCFixes.py
@@ -0,0 +1,34 @@
+from AthenaPython.PyAthena import McEventCollection, HepMC, HepMC3  # MCEventCollection is needed for HepMC3 -> HepMC alias
+import ROOT as R
+
+def add(self, other):
+    self.set(self.x() + other.x(), self.y() + other.y(),
+             self.z() + other.z(), self.t() + other.t())
+    return self
+
+HepMC.FourVector.__iadd__ = add
+HepMC3.FourVector.__iadd__ = add
+del add
+
+def print_particle(self):
+    HepMC.Print.line(R.std.shared_ptr["HepMC3::GenParticle"](self))
+    return
+    
+HepMC.GenParticle.print = print_particle
+HepMC3.GenParticle.print = print_particle
+del print_particle
+
+def print_vertex(self):
+    HepMC.Print.line(R.std.shared_ptr["HepMC3::GenVertex"](self))
+    return
+
+HepMC.GenVertex.print = print_vertex
+HepMC3.GenVertex.print = print_vertex
+del print_vertex
+
+def print_event(self):
+    HepMC.Print.content(self)
+    
+HepMC.GenEvent.print = print_event
+HepMC3.GenEvent.print = print_event
+del print_event
diff --git a/Generators/GeneratorUtils/python/ShiftLOS.py b/Generators/GeneratorUtils/python/ShiftLOS.py
index 3c882e269ff9d6369953fa3ebb8a6f06ef0b5f86..922e75ecf4353b399b24819a9ac7dfeb589de7c8 100644
--- a/Generators/GeneratorUtils/python/ShiftLOS.py
+++ b/Generators/GeneratorUtils/python/ShiftLOS.py
@@ -29,7 +29,7 @@ class ShiftLOS(PyAthena.Alg):
             return evt
 
         # Loop over all vertices
-        for v in evt.vertices:
+        for v in evt.vertices():
             # Get position
             pos = v.position()
             x = pos.x()
@@ -66,7 +66,7 @@ class ShiftLOS(PyAthena.Alg):
         pxsum_orig, pysum_orig = 0,0
 
         # Loop over all particles
-        for p in evt.particles:
+        for p in evt.particles():
             # Get momentum
             mom = p.momentum()
 
@@ -75,7 +75,7 @@ class ShiftLOS(PyAthena.Alg):
 
             # Boost in x or y using CLHEP
             boost = CLHEP.Hep3Vector(self.xcross, self.ycross, 0.0)
-            tmp = CLHEP.HepLorentzVector(mom.px(), mom.py(), mom.pz(), mom.e())
+            tmp = CLHEP.HepLorentzVector(mom.px(), mom.py(), mom.pz(), mom.e())  # this assumes the angle is small so theta = tan(theta)
             tmp.boost(boost)
             
             pxsum += tmp.x() - mom.x()
@@ -88,6 +88,42 @@ class ShiftLOS(PyAthena.Alg):
 
         return evt
 
+    def rotate_particles(self):
+        if self.xcross == self.ycross == 0:
+            return evt
+
+        pxsum, pysum = 0,0
+        pxsum_orig, pysum_orig = 0,0
+
+        # Loop over all particles
+        for p in evt.particles():
+            # Get momentum
+            mom = p.momentum()
+
+            pxsum_orig += mom.x()
+            pysum_orig += mom.y()            
+
+            # Rotate in x or y using CLHEP
+
+            tmp = CLHEP.HepLorentzVector(mom.px(), mom.py(), mom.pz(), mom.e())
+            
+            rotate = CLHEP.HepRotation()
+            rotate.rotateX(self.xcross)
+            rotate.rotateY(self.ycross)
+            
+            lorentzRoation = CLHEP.HepLorentzRotation()
+            tmp = mom * lotentzRotation
+            
+            pxsum += tmp.x() - mom.x()
+            pysum += tmp.y() - mom.y()            
+
+            # Convert back to HepMC
+            p.set_momentum(HepMC.FourVector(tmp.px(), tmp.py(), tmp.pz(), tmp.e()))
+            
+            self.msg.debug(f"Change in total px = {pxsum:.1f} MeV ({pxsum/pxsum_orig * 100: .3f} %), change in total py = {pysum:.1f} MeV ({pysum/pysum_orig * 100: .3f} %)")
+
+
+        
     def execute(self):
         self.msg.debug(f"Exectuing {self.getName()}")
 
diff --git a/Generators/GeneratorUtils/share/convert_hepmc_2to3.py b/Generators/GeneratorUtils/share/convert_hepmc_2to3.py
new file mode 100644
index 0000000000000000000000000000000000000000..5efb96d922db7a3df1a946193aca3b148e641fdb
--- /dev/null
+++ b/Generators/GeneratorUtils/share/convert_hepmc_2to3.py
@@ -0,0 +1,27 @@
+from pyHepMC3 import HepMC3
+
+def convert(infile, outfile):
+    reader = HepMC3.ReaderAsciiHepMC2(infile)
+    writer = HepMC3.WriterAscii(outfile)
+
+    event = HepMC3.GenEvent()
+
+    while not reader.failed():
+        reader.read_event(event)
+        if reader.failed():
+            break
+
+        writer.write_event(event)
+
+    reader.close()
+    writer.close()
+
+if __name__ == "__main__":
+    
+    import argparse, sys
+    parser = argparse.ArgumentParser(description="Conver HEPMC2 ASCII files to HEPMC3 ASCII files")
+    parser.add_argument("file", help = "full path to input file")
+    parser.add_argument("--output", "-o", default = "output_HepMC3.hepmc")
+    args = parser.parse_args()
+                        
+    convert(args.file, args.output)