diff --git a/python/experimental/engineParameters.py b/python/experimental/engineParameters.py
new file mode 100644
index 0000000000..aac2d7cb8b
--- /dev/null
+++ b/python/experimental/engineParameters.py
@@ -0,0 +1,354 @@
+#!/usr/bin/env python
+import os
+import sys
+import ROOT
+
+import shipunit as u
+import rootUtils as ut
+import yaml
+
+
+def SetPrimGen(primGen, options, ship_geo, modules):
+    with open(options.yaml_file) as file:
+        config = yaml.safe_load(file)
+
+    #### Cosmics
+    if options.simEngine == "Cosmics":
+        cosmics = AttrDict(config["Cosmics"])
+        primGen.SetTarget(0.0, 0.0)
+        Cosmicsgen = ROOT.CosmicsGenerator()
+        import CMBG_conf
+
+        CMBG_conf.configure(Cosmicsgen, ship_geo)
+        if not Cosmicsgen.Init(cosmics.opt):
+            print("initialization of cosmic background generator failed ", cosmics.opt)
+            sys.exit(0)
+        Cosmicsgen.n_EVENTS = options.nEvents
+        primGen.AddGenerator(Cosmicsgen)
+        print("Process ", options.nEvents, " Cosmic events with option ", cosmics.opt)
+        #
+
+    #####Genie
+    # -----Neutrino Background------------------------
+    if options.simEngine == "Genie":
+        genie = AttrDict(config["Genie"])
+        # Genie
+        ut.checkFileExists(genie.inputFile)
+        primGen.SetTarget(0.0, 0.0)  # do not interfere with GenieGenerator
+        Geniegen = ROOT.GenieGenerator()
+        Geniegen.Init(genie.inputFile, options.firstEvent)
+        Geniegen.SetPositions(
+            ship_geo.target.z0,
+            ship_geo.tauMudet.zMudetC - 5 * u.m,
+            ship_geo.TrackStation2.z,
+        )
+        primGen.AddGenerator(Geniegen)
+        options.nEvents = min(options.nEvents, Geniegen.GetNevents())
+        run.SetPythiaDecayer("DecayConfigNuAge.C")
+        print(
+            "Generate ",
+            options.nEvents,
+            " with Genie input",
+            " first event",
+            options.firstEvent,
+        )
+
+    ##### nuRadiography
+
+    if options.simEngine == "nuRadiography":
+        nuRad = AttrDict(config["NuRadio"])
+        ut.checkFileExists(nuRad.inputFile)
+        primGen.SetTarget(0.0, 0.0)  # do not interfere with GenieGenerator
+        Geniegen = ROOT.GenieGenerator()
+        Geniegen.Init(nuRad.inputFile, options.firstEvent)
+        # Geniegen.SetPositions(ship_geo.target.z0, ship_geo.target.z0, ship_geo.MuonStation3.z)
+        Geniegen.SetPositions(
+            ship_geo.target.z0, ship_geo.tauMudet.zMudetC, ship_geo.MuonStation3.z
+        )
+        Geniegen.NuOnly()
+        primGen.AddGenerator(Geniegen)
+        print(
+            "Generate ",
+            options.nEvents,
+            " for nuRadiography",
+            " first event",
+            options.firstEvent,
+        )
+
+    #####muonDIS
+    if options.simEngine == "muonDIS":
+        muDIS = AttrDict(config["MuonDIS"])
+        ut.checkFileExists(muDIS.inputFile)
+        primGen.SetTarget(0.0, 0.0)
+        DISgen = ROOT.MuDISGenerator()
+        # from nu_tau detector to tracking station 2
+        # mu_start, mu_end =  ship_geo.tauMudet.zMudetC,ship_geo.TrackStation2.z
+        #
+        # in front of UVT up to tracking station 1
+        mu_start, mu_end = (
+            ship_geo.Chamber1.z - ship_geo.chambers.Tub1length - 10.0 * u.cm,
+            ship_geo.TrackStation1.z,
+        )
+        print("MuDIS position info input=", mu_start, mu_end)
+        DISgen.SetPositions(mu_start, mu_end)
+        DISgen.Init(muDIS.inputFile, options.firstEvent)
+        primGen.AddGenerator(DISgen)
+        options.nEvents = min(options.nEvents, DISgen.GetNevents())
+        print(
+            "Generate ",
+            options.nEvents,
+            " with DIS input",
+            " first event",
+            options.firstEvent,
+        )
+
+    ######Nuage
+    # -----neutrino interactions from nuage------------------------
+    if options.simEngine == "Nuage":
+        nuage = AttrDict(config["Nuage"])
+        primGen.SetTarget(0.0, 0.0)
+        Nuagegen = ROOT.NuageGenerator()
+        Nuagegen.EnableExternalDecayer(
+            1
+        )  # with 0 external decayer is disable, 1 is enabled
+
+        # CAMM - This is broken now, need info from dedicated SND geo...
+        print(
+            "Nuage position info input=",
+            ship_geo.EmuMagnet.zC - ship_geo.NuTauTarget.zdim,
+            ship_geo.EmuMagnet.zC + ship_geo.NuTauTarget.zdim,
+        )
+        # --------------------------------
+        # to Generate neutrino interactions in the whole neutrino target
+        # Nuagegen.SetPositions(ship_geo.EmuMagnet.zC, ship_geo.NuTauTarget.zC-ship_geo.NuTauTarget.zdim/2, ship_geo.NuTauTarget.zC+ship_geo.NuTauTarget.zdim/2, -ship_geo.NuTauTarget.xdim/2, ship_geo.NuTauTarget.xdim/2, -ship_geo.NuTauTarget.ydim/2, ship_geo.NuTauTarget.ydim/2)
+        # --------------------------------
+        # to Generate neutrino interactions ONLY in ONE brick
+        ntt = 6
+        nXcells = 7
+        nYcells = 3
+        nZcells = ntt - 1
+        startx = -ship_geo.NuTauTarget.xdim / 2.0 + nXcells * ship_geo.NuTauTarget.BrX
+        endx = (
+            -ship_geo.NuTauTarget.xdim / 2.0 + (nXcells + 1) * ship_geo.NuTauTarget.BrX
+        )
+        starty = -ship_geo.NuTauTarget.ydim / 2.0 + nYcells * ship_geo.NuTauTarget.BrY
+        endy = (
+            -ship_geo.NuTauTarget.ydim / 2.0 + (nYcells + 1) * ship_geo.NuTauTarget.BrY
+        )
+        startz = (
+            ship_geo.EmuMagnet.zC
+            - ship_geo.NuTauTarget.zdim / 2.0
+            + ntt * ship_geo.NuTauTT.TTZ
+            + nZcells * ship_geo.NuTauTarget.CellW
+        )
+        endz = (
+            ship_geo.EmuMagnet.zC
+            - ship_geo.NuTauTarget.zdim / 2.0
+            + ntt * ship_geo.NuTauTT.TTZ
+            + nZcells * ship_geo.NuTauTarget.CellW
+            + ship_geo.NuTauTarget.BrZ
+        )
+        Nuagegen.SetPositions(
+            ship_geo.target.z0, startz, endz, startx, endx, starty, endy
+        )
+        # --------------------------------
+        ut.checkFileExists(nuage.inputFile)
+        Nuagegen.Init(nuage.inputFile, options.firstEvent)
+        primGen.AddGenerator(Nuagegen)
+        options.nEvents = min(options.nEvents, Nuagegen.GetNevents())
+        print(
+            "Generate ",
+            options.nEvents,
+            " with Nuage input",
+            " first event",
+            options.firstEvent,
+        )
+        # -CAMM end broken part
+
+    ######Ntuple
+
+    if options.simEngine == "Ntuple":
+        ntuple = AttrDict(config["Ntuple"])
+        # reading previously processed muon events, [-50m - 50m]
+        ut.checkFileExists(ntuple.inputFile)
+        primGen.SetTarget(ship_geo.target.z0 + 50 * u.m, 0.0)
+        Ntuplegen = ROOT.NtupleGenerator()
+        Ntuplegen.Init(ntuple.inputFile, options.firstEvent)
+        primGen.AddGenerator(Ntuplegen)
+        options.nEvents = min(options.nEvents, Ntuplegen.GetNevents())
+        print("Process ", options.nEvents, " from input file")
+        #
+
+    #### MuonBack
+    if options.simEngine == "MuonBack":
+        # reading muon tracks from previous Pythia8/Geant4 simulation with charm replaced by cascade production
+        muBkg = AttrDict(config["MuonBack"])
+        fileType = ut.checkFileExists(muBkg.inputFile)
+        if fileType == "tree":
+            # 2018 background production
+            primGen.SetTarget(ship_geo.target.z0 + 70.845 * u.m, 0.0)
+        else:
+            primGen.SetTarget(ship_geo.target.z0 + 50 * u.m, 0.0)
+            #
+        MuonBackgen = ROOT.MuonBackGenerator()
+        # MuonBackgen.FollowAllParticles() # will follow all particles after hadron absorber, not only muons
+        MuonBackgen.Init(muBkg.inputFile, options.firstEvent, muBkg.phiRandom)
+        MuonBackgen.SetSmearBeam(5 * u.cm)  # radius of ring, thickness 8mm
+        if muBkg.DownScaleDiMuon:
+            testf = ROOT.TFile.Open(muBkg.inputFile)
+            if not testf.FileHeader.GetTitle().find("diMu100.0") < 0:
+                MuonBackgen.SetDownScaleDiMuon()  # avoid interference with boosted channels
+                print("MuonBackgenerator: set downscale for dimuon on")
+            testf.Close()
+        if muBkg.sameSeed:
+            MuonBackgen.SetSameSeed(muBkg.sameSeed)
+        primGen.AddGenerator(MuonBackgen)
+        options.nEvents = min(options.nEvents, MuonBackgen.GetNevents())
+        print(
+            "Process ",
+            options.nEvents,
+            " from input file, with Phi random=",
+            options.phiRandom,
+            " with MCTracksWithHitsOnly True",
+        )
+        if muBkg.followMuon:
+            muBkg.fastMuon = True
+            modules["Veto"].SetFollowMuon()
+        if muBkg.fastMuon:
+            modules["Veto"].SetFastMuon()
+
+    #### Pythia8
+    if options.simEngine == "Pythia8":
+        p8dict = AttrDict(config["Pythia8"])
+        primGen.SetTarget(ship_geo.target.z0, 0.0)
+        # -----Pythia8--------------------------------------
+        if p8dict.HNL or p8dict.RPVSUSY:
+            P8gen = ROOT.HNLPythia8Generator()
+            import pythia8_conf
+
+            if p8dict.HNL:
+                print("Generating HNL events of mass %.3f GeV" % p8dict.theMass)
+                pythia8_conf.configure(
+                    P8gen,
+                    p8dict.theMass,
+                    p8dict.theProductionCouplings,
+                    p8dict.theDecayCouplings,
+                    p8dict.inclusive,
+                    options.deepCopy,
+                )
+            if p8dict.RPVSUSY:
+                print("Generating RPVSUSY events of mass %.3f GeV" % (theMass * u.GeV))
+                print(
+                    "and with couplings=[%.3f,%.3f]"
+                    % (p8dict.theRPVCouplings[0], p8dict.theRPVCouplings[1])
+                )
+                print("and with stop mass=%.3f GeV\n" % p8dict.theRPVCouplings[2])
+                pythia8_conf.configurerpvsusy(
+                    P8gen,
+                    p8dict.theMass,
+                    [p8dict.theRPVCouplings[0], p8dict.theRPVCouplings[1]],
+                    p8dict.theRPVCouplings[2],
+                    p8dict.RPVSUSYbench,
+                    p8dict.inclusive,
+                    options.deepCopy,
+                )
+                P8gen.SetParameters("ProcessLevel:all = off")
+            if p8dict.inputFile:
+                ut.checkFileExists(p8dict.inputFile)
+                # read from external file
+                P8gen.UseExternalFile(p8dict.inputFile, options.firstEvent)
+        if p8dict.DarkPhoton:
+            P8gen = ROOT.DPPythia8Generator()
+            if inclusive == "qcd":
+                P8gen.SetDPId(4900023)
+            else:
+                P8gen.SetDPId(9900015)
+            import pythia8darkphoton_conf
+
+            passDPconf = pythia8darkphoton_conf.configure(
+                P8gen,
+                p8dict.theMass,
+                p8dict.theDPepsilon,
+                p8dict.inclusive,
+                p8dict.motherMode,
+                options.deepCopy,
+            )
+            if passDPconf != 1:
+                sys.exit()
+        if p8dict.HNL or p8dict.RPVSUSY or p8dict.DarkPhoton:
+            P8gen.SetSmearBeam(1 * u.cm)  # finite beam size
+            P8gen.SetLmin(
+                (ship_geo.Chamber1.z - ship_geo.chambers.Tub1length)
+                - ship_geo.target.z0
+            )
+            P8gen.SetLmax(ship_geo.TrackStation1.z - ship_geo.target.z0)
+        if charmonly:
+            primGen.SetTarget(0.0, 0.0)  # vertex is set in pythia8Generator
+            ut.checkFileExists(p8dict.inputFile)
+            if ship_geo.Box.gausbeam:
+                primGen.SetBeam(
+                    0.0, 0.0, 0.5, 0.5
+                )  # more central beam, for hits in downstream detectors
+                primGen.SmearGausVertexXY(True)  # sigma = x
+            else:
+                primGen.SetBeam(
+                    0.0, 0.0, ship_geo.Box.TX - 1.0, ship_geo.Box.TY - 1.0
+                )  # Uniform distribution in x/y on the target (0.5 cm of margin at both sides)
+                primGen.SmearVertexXY(True)
+            P8gen = ROOT.Pythia8Generator()
+            P8gen.UseExternalFile(p8dict.inputFile, options.firstEvent)
+            P8gen.SetTarget(
+                "volTarget_1", 0.0, 0.0
+            )  # will distribute PV inside target, beam offset x=y=0.
+            # pion on proton 500GeV
+            # P8gen.SetMom(500.*u.GeV)
+            # P8gen.SetId(-211)
+            primGen.AddGenerator(P8gen)
+
+    ####### fixed target
+    if options.simEngine == "FixedTarget":
+        P8gen = ROOT.FixedTargetGenerator()
+        P8gen.SetTarget("volTarget_1", 0.0, 0.0)
+        P8gen.SetMom(400.0 * u.GeV)
+        P8gen.SetEnergyCut(0.0)
+        P8gen.SetHeartBeat(100000)
+        P8gen.SetG4only()
+        primGen.AddGenerator(P8gen)
+
+    ####Pythia6
+    if options.simEngine == "Pythia6":
+        # set muon interaction close to decay volume
+        primGen.SetTarget(ship_geo.target.z0 + ship_geo.muShield.length, 0.0)
+        # -----Pythia6-------------------------
+        test = ROOT.TPythia6()  # don't know any other way of forcing to load lib
+        P6gen = ROOT.tPythia6Generator()
+        P6gen.SetMom(50.0 * u.GeV)
+        P6gen.SetTarget("gamma/mu+", "n0")  # default "gamma/mu-","p+"
+        primGen.AddGenerator(P6gen)
+
+    #### EvtCalc
+
+    # -----EvtCalc--------------------------------------
+    if options.simEngine == "EvtCalc":
+        evtcalc = AttrDict(config["EvtCalc"])
+        primGen.SetTarget(0.0, 0.0)
+        print(f"Opening input file for EvtCalc generator: {evtcalc.inputFile}")
+        ut.checkFileExists(evtcalc.inputFile)
+        EvtCalcGen = ROOT.EvtCalcGenerator()
+        EvtCalcGen.Init(evtcalc.inputFile, options.firstEvent)
+        EvtCalcGen.SetPositions(zTa=ship_geo.target.z, zDV=ship_geo.decayVolume.z)
+        primGen.AddGenerator(EvtCalcGen)
+        options.nEvents = min(options.nEvents, EvtCalcGen.GetNevents())
+        print(
+            f"Generate {options.nEvents} with EvtCalc input. First event: {options.firstEvent}"
+        )
+
+    # -----Particle Gun-----------------------
+    if options.simEngine == "PG":
+        pg = AttrDict(config["ParticleGun"])
+        myPgun = ROOT.FairBoxGenerator(pg.pID, 1)
+        myPgun.SetPRange(pg.Estart, pg.Eend)
+        myPgun.SetPhiRange(0, 360)  # // Azimuth angle range [degree]
+        myPgun.SetXYZ(0.0 * u.cm, 0.0 * u.cm, 0.0 * u.cm)
+        myPgun.SetThetaRange(0, 0)  # // Polar angle in lab system range [degree]
+        primGen.AddGenerator(myPgun)
diff --git a/python/experimental/run_sim_test.py b/python/experimental/run_sim_test.py
new file mode 100644
index 0000000000..30fa384dca
--- /dev/null
+++ b/python/experimental/run_sim_test.py
@@ -0,0 +1,514 @@
+#!/usr/bin/env python
+import os
+import sys
+import ROOT
+
+import shipunit as u
+import shipRoot_conf
+import rootUtils as ut
+from ShipGeoConfig import ConfigRegistry
+from argparse import ArgumentParser
+from array import array
+
+
+mcEngine = "TGeant4"
+
+MCTracksWithHitsOnly = False  # copy particles which produced a hit and their history
+MCTracksWithEnergyCutOnly = True  # copy particles above a certain kin energy cut
+MCTracksWithHitsOrEnergyCut = False  # or of above, factor 2 file size increase compared to MCTracksWithEnergyCutOnly
+
+globalDesigns = {
+    "2023": {"dy": 6.0, "dv": 6, "nud": 4, "caloDesign": 3, "strawDesign": 10}
+}
+default = "2023"
+
+
+parser = ArgumentParser()
+group = parser.add_mutually_exclusive_group()
+
+parser.add_argument(
+    "-E",
+    dest="simEngine",
+    help="Possible options: EvtCalc,Pythia6, Pythia8, PG, Genie, nuRadiography, Ntuple, MuonBack, Nuage, muonDIS, Cosmics",
+    required=True,
+    default="Pythia8",
+)
+parser.add_argument(
+    "-n",
+    "--nEvents",
+    dest="nEvents",
+    help="Number of events to generate",
+    required=False,
+    default=100,
+    type=int,
+)
+parser.add_argument(
+    "-i",
+    "--firstEvent",
+    dest="firstEvent",
+    help="First event of input file to use",
+    required=False,
+    default=0,
+    type=int,
+)
+parser.add_argument(
+    "-F",
+    dest="deepCopy",
+    help="default = False: copy only stable particles to stack, except for HNL events",
+    required=False,
+    action="store_true",
+)
+parser.add_argument(
+    "-s",
+    "--seed",
+    dest="theSeed",
+    help="Seed for random number. Only for experts, see TRrandom::SetSeed documentation",
+    required=False,
+    default=0,
+    type=int,
+)
+parser.add_argument(
+    "-f",
+    dest="yaml_file",
+    help="Yaml configuration file for simEngine parameters",
+    required=True,
+    default="shipgen/shipgen_parameters.yaml",
+)
+parser.add_argument(
+    "-g",
+    dest="geofile",
+    help="geofile for muon shield geometry, for experts only",
+    required=False,
+    default=None,
+)
+parser.add_argument(
+    "-o",
+    "--output",
+    dest="outputDir",
+    help="Output directory",
+    required=False,
+    default=".",
+)
+parser.add_argument(
+    "-Y",
+    dest="dy",
+    help="max height of vacuum tank",
+    required=False,
+    default=globalDesigns[default]["dy"],
+)
+parser.add_argument(
+    "--tankDesign",
+    dest="dv",
+    help="4=TP elliptical tank design, 5 = optimized conical rectangular design, 6=5 without segment-1",
+    required=False,
+    default=globalDesigns[default]["dv"],
+    type=int,
+)
+parser.add_argument(
+    "--nuTauTargetDesign",
+    dest="nud",
+    help="3: emulsion spectrometer and muon filter as in CDS, 4: not magnetized target and muon spectrometer for ECN3",
+    default=globalDesigns[default]["nud"],
+    type=int,
+    choices=[3, 4],
+)
+parser.add_argument(
+    "--caloDesign",
+    help="0=ECAL/HCAL TP 2=splitCal  3=ECAL/ passive HCAL",
+    default=globalDesigns[default]["caloDesign"],
+    type=int,
+    choices=[0, 2, 3],
+)
+parser.add_argument(
+    "--strawDesign",
+    dest="strawDesign",
+    help="simplistic tracker design,  4=sophisticated straw tube design, horizontal wires (default), 10=2cm straw",
+    required=False,
+    default=globalDesigns[default]["strawDesign"],
+    type=int,
+)
+parser.add_argument(
+    "-F",
+    dest="deepCopy",
+    help="default = False: copy only stable particles to stack, except for HNL events",
+    required=False,
+    action="store_true",
+)
+parser.add_argument(
+    "--dry-run",
+    dest="dryrun",
+    help="stop after initialize",
+    required=False,
+    action="store_true",
+)
+parser.add_argument(
+    "-D",
+    "--display",
+    dest="eventDisplay",
+    help="store trajectories",
+    required=False,
+    action="store_true",
+)
+parser.add_argument(
+    "--shieldName",
+    help="The name of the SC shield in the database. SC default: sc_v6, Warm default: warm_opt",
+    default="sc_v6",
+    choices=["sc_v6", "warm_opt"],
+)
+parser.add_argument(
+    "--debug",
+    help="1: print weights and field 2: make overlap check",
+    required=False,
+    default=0,
+    type=int,
+    choices=range(0, 3),
+)
+parser.add_argument("--field_map", default=None, help="Specify spectrometer field map.")
+parser.add_argument(
+    "--helium",
+    dest="decayVolMed",
+    help="Set Decay Volume medium to helium. NOOP, as default is helium",
+    action="store_const",
+    const="helium",
+    default="helium",
+)
+parser.add_argument(
+    "--vacuums",
+    dest="decayVolMed",
+    help="Set Decay Volume medium to vacuum(vessel structure changes)",
+    action="store_const",
+    const="vacuums",
+    default="helium",
+)
+
+parser.add_argument("--SND", dest="SND", help="Activate SND.", action="store_true")
+parser.add_argument(
+    "--noSND",
+    dest="SND",
+    help="Deactivate SND. NOOP, as it currently defaults to off.",
+    action="store_false",
+)
+
+options = parser.parse_args()
+
+
+with open(options.yaml_file) as file:
+    config = yaml.safe_load(file)
+
+
+ROOT.gRandom.SetSeed(
+    options.theSeed
+)  # this should be propagated via ROOT to Pythia8 and Geant4VMC
+shipRoot_conf.configure(0)  # load basic libraries, prepare atexit for python
+# - targetOpt      = 5  # 0=solid   >0 sliced, 5: 5 pieces of tungsten, 4 H20 slits, 17: Mo + W +H2O (default)
+#   nuTauTargetDesign = 3 #3 = 2018 design, 4 = not magnetized target + spectrometer
+ship_geo = ConfigRegistry.loadpy(
+    "$FAIRSHIP/geometry/geometry_config.py",
+    Yheight=options.dy,
+    tankDesign=options.dv,
+    nuTauTargetDesign=options.nud,
+    CaloDesign=options.caloDesign,
+    strawDesign=options.strawDesign,
+    muShieldGeo=options.geofile,
+    shieldName=options.shieldName,
+    DecayVolumeMedium=options.decayVolMed,
+    SND=options.SND,
+)
+
+
+# Output file name, add dy to be able to setup geometry with ambiguities.
+if options.simEngine == "PG":
+    tag = options.simEngine + "_" + str(config["ParticleGun"]["pID"]) + "-" + mcEngine
+else:
+    tag = options.simEngine + "-" + mcEngine
+if config["Pythia8"]["charmonly"] == True:
+    tag = options.simEngine + "CharmOnly-" + mcEngine
+if options.eventDisplay:
+    tag = tag + "_D"
+if options.dv > 4:
+    tag = "conical." + tag
+if not os.path.exists(options.outputDir):
+    os.makedirs(options.outputDir)
+outFile = f"{options.outputDir}/ship.{tag}.root"
+
+# rm older files !!!
+for x in os.listdir(options.outputDir):
+    if not x.find(tag) < 0:
+        os.system(f"rm {options.outputDir}/{x}")
+# Parameter file name
+parFile = f"{options.outputDir}/ship.params.{tag}.root"
+
+
+# In general, the following parts need not be touched
+# ========================================================================
+
+# -----Timer--------------------------------------------------------
+timer = ROOT.TStopwatch()
+timer.Start()
+# ------------------------------------------------------------------------
+# -----Create simulation run----------------------------------------
+run = ROOT.FairRunSim()
+run.SetName(mcEngine)  # Transport engine
+run.SetSink(ROOT.FairRootFileSink(outFile))  # Output file
+run.SetUserConfig("g4Config.C")  # user configuration file default g4Config.C
+rtdb = run.GetRuntimeDb()
+# -----Create geometry----------------------------------------------
+# import shipMuShield_only as shipDet_conf # special use case for an attempt to convert active shielding geometry for use with FLUKA
+# import shipTarget_only as shipDet_conf
+import shipDet_conf
+
+modules = shipDet_conf.configure(run, ship_geo)
+# -----Create PrimaryGenerator--------------------------------------
+primGen = ROOT.FairPrimaryGenerator()
+
+SetPrimGen(primGen, options, ship_geo, modules)
+if options.simEngine == "nuRadiography":
+    #  add tungsten to PDG
+    pdg = ROOT.TDatabasePDG.Instance()
+    pdg.AddParticle("W", "Ion", 1.71350e02, True, 0.0, 74, "XXX", 1000741840)
+    #
+    run.SetPythiaDecayer("DecayConfigPy8.C")
+# this requires writing a C macro, would have been easier to do directly in python!
+
+
+if options.simEngine == "Nuage":
+    run.SetPythiaDecayer("DecayConfigNuAge.C")
+
+if options.simEngine == "MuonBack":
+    MCTracksWithHitsOnly = True  # otherwise, output file becomes too big
+
+
+#
+run.SetGenerator(primGen)
+# ------------------------------------------------------------------------
+
+# ---Store the visualiztion info of the tracks, this make the output file very large!!
+# --- Use it only to display but not for production!
+if options.eventDisplay:
+    run.SetStoreTraj(ROOT.kTRUE)
+else:
+    run.SetStoreTraj(ROOT.kFALSE)
+# -----Initialize simulation run------------------------------------
+run.Init()
+if options.dryrun:  # Early stop after setting up Pythia 8
+    sys.exit(0)
+gMC = ROOT.TVirtualMC.GetMC()
+fStack = gMC.GetStack()
+EnergyCut = 10.0 * u.MeV if options.mudis else 100.0 * u.MeV
+
+if MCTracksWithHitsOnly:
+    fStack.SetMinPoints(1)
+    fStack.SetEnergyCut(-100.0 * u.MeV)
+elif MCTracksWithEnergyCutOnly:
+    fStack.SetMinPoints(-1)
+    fStack.SetEnergyCut(EnergyCut)
+elif MCTracksWithHitsOrEnergyCut:
+    fStack.SetMinPoints(1)
+    fStack.SetEnergyCut(EnergyCut)
+elif options.deepCopy:
+    fStack.SetMinPoints(0)
+    fStack.SetEnergyCut(0.0 * u.MeV)
+
+if options.eventDisplay:
+    # Set cuts for storing the trajectories, can only be done after initialization of run (?!)
+    trajFilter = ROOT.FairTrajFilter.Instance()
+    trajFilter.SetStepSizeCut(1 * u.mm)
+    trajFilter.SetVertexCut(
+        -20 * u.m,
+        -20 * u.m,
+        ship_geo.target.z0 - 1 * u.m,
+        20 * u.m,
+        20 * u.m,
+        200.0 * u.m,
+    )
+    trajFilter.SetMomentumCutP(0.1 * u.GeV)
+    trajFilter.SetEnergyCut(0.0, 400.0 * u.GeV)
+    trajFilter.SetStorePrimaries(ROOT.kTRUE)
+    trajFilter.SetStoreSecondaries(ROOT.kTRUE)
+
+# The VMC sets the fields using the "/mcDet/setIsLocalMagField true" option in "gconfig/g4config.in"
+import geomGeant4
+# geomGeant4.setMagnetField() # replaced by VMC, only has effect if /mcDet/setIsLocalMagField  false
+
+# Define extra VMC B fields not already set by the geometry definitions, e.g. a global field,
+# any field maps, or defining if any volumes feel only the local or local+global field.
+# For now, just keep the fields already defined by the C++ code, i.e comment out the fieldMaker
+if hasattr(ship_geo.Bfield, "fieldMap"):
+    if options.field_map:
+        ship_geo.Bfield.fieldMap = options.field_map
+    fieldMaker = geomGeant4.addVMCFields(ship_geo, verbose=True)
+
+# Print VMC fields and associated geometry objects
+if options.debug == 1:
+    geomGeant4.printVMCFields()
+    geomGeant4.printWeightsandFields(
+        onlyWithField=True,
+        exclude=[
+            "DecayVolume",
+            "Tr1",
+            "Tr2",
+            "Tr3",
+            "Tr4",
+            "Veto",
+            "Ecal",
+            "Hcal",
+            "MuonDetector",
+            "SplitCal",
+        ],
+    )
+# Plot the field example
+# fieldMaker.plotField(1, ROOT.TVector3(-9000.0, 6000.0, 50.0), ROOT.TVector3(-300.0, 300.0, 6.0), 'Bzx.png')
+# fieldMaker.plotField(2, ROOT.TVector3(-9000.0, 6000.0, 50.0), ROOT.TVector3(-400.0, 400.0, 6.0), 'Bzy.png')
+
+# -----Start run----------------------------------------------------
+run.Run(options.nEvents)
+# -----Runtime database---------------------------------------------
+kParameterMerged = ROOT.kTRUE
+parOut = ROOT.FairParRootFileIo(kParameterMerged)
+parOut.open(parFile)
+rtdb.setOutput(parOut)
+rtdb.saveOutput()
+rtdb.printParamContexts()
+getattr(rtdb, "print")()
+# ------------------------------------------------------------------------
+run.CreateGeometryFile(f"{options.outputDir}/geofile_full.{tag}.root")
+# save ShipGeo dictionary in geofile
+import saveBasicParameters
+
+saveBasicParameters.execute(f"{options.outputDir}/geofile_full.{tag}.root", ship_geo)
+
+# checking for overlaps
+if options.debug == 2:
+    fGeo = ROOT.gGeoManager
+    fGeo.SetNmeshPoints(10000)
+    fGeo.CheckOverlaps(0.1)  # 1 micron takes 5minutes
+    fGeo.PrintOverlaps()
+    # check subsystems in more detail
+    for x in fGeo.GetTopNode().GetNodes():
+        x.CheckOverlaps(0.0001)
+        fGeo.PrintOverlaps()
+# -----Finish-------------------------------------------------------
+timer.Stop()
+rtime = timer.RealTime()
+ctime = timer.CpuTime()
+print(" ")
+print("Macro finished successfully.")
+if "P8gen" in globals():
+    if HNL:
+        print("number of retries, events without HNL ", P8gen.nrOfRetries())
+    elif options.DarkPhoton:
+        print("number of retries, events without Dark Photons ", P8gen.nrOfRetries())
+        print(
+            "total number of dark photons (including multiple meson decays per single collision) ",
+            P8gen.nrOfDP(),
+        )
+
+print("Output file is ", outFile)
+print("Parameter file is ", parFile)
+print("Real time ", rtime, " s, CPU time ", ctime, "s")
+
+# remove empty events
+if options.simEngine == "MuonBack":
+    tmpFile = outFile + "tmp"
+    xxx = outFile.split("/")
+    check = xxx[len(xxx) - 1]
+    fin = False
+    for ff in ROOT.gROOT.GetListOfFiles():
+        nm = ff.GetName().split("/")
+        if nm[len(nm) - 1] == check:
+            fin = ff
+    if not fin:
+        fin = ROOT.TFile.Open(outFile)
+    t = fin.Get("cbmsim")
+    fout = ROOT.TFile(tmpFile, "recreate")
+    fSink = ROOT.FairRootFileSink(fout)
+
+    sTree = t.CloneTree(0)
+    nEvents = 0
+    pointContainers = []
+    for x in sTree.GetListOfBranches():
+        name = x.GetName()
+        if not name.find("Point") < 0:
+            pointContainers.append(
+                "sTree." + name + ".GetEntries()"
+            )  # makes use of convention that all sensitive detectors fill XXXPoint containers
+    for n in range(t.GetEntries()):
+        rc = t.GetEvent(n)
+        empty = True
+        for x in pointContainers:
+            if eval(x) > 0:
+                empty = False
+        if not empty:
+            rc = sTree.Fill()
+            nEvents += 1
+
+    branches = ROOT.TList()
+    branches.SetName("BranchList")
+    branches.Add(ROOT.TObjString("MCTrack"))
+    branches.Add(ROOT.TObjString("vetoPoint"))
+    branches.Add(ROOT.TObjString("ShipRpcPoint"))
+    branches.Add(ROOT.TObjString("TargetPoint"))
+    branches.Add(ROOT.TObjString("TTPoint"))
+    branches.Add(ROOT.TObjString("ScoringPoint"))
+    branches.Add(ROOT.TObjString("strawtubesPoint"))
+    branches.Add(ROOT.TObjString("EcalPoint"))
+    branches.Add(ROOT.TObjString("sEcalPointLite"))
+    branches.Add(ROOT.TObjString("smuonPoint"))
+    branches.Add(ROOT.TObjString("TimeDetPoint"))
+    branches.Add(ROOT.TObjString("MCEventHeader"))
+    branches.Add(ROOT.TObjString("UpstreamTaggerPoint"))
+    branches.Add(ROOT.TObjString("sGeoTracks"))
+
+    sTree.AutoSave()
+    fSink.WriteObject(branches, "BranchList", ROOT.TObject.kSingleKey)
+    fSink.SetOutTree(sTree)
+
+    fout.Close()
+    print("removed empty events, left with:", nEvents)
+    rc1 = os.system("rm  " + outFile)
+    rc2 = os.system("mv " + tmpFile + " " + outFile)
+    fin.SetWritable(False)  # bpyass flush error
+
+if options.simEngine == "muonDIS":
+    temp_filename = outFile.replace(".root", "_tmp.root")
+
+    with (
+        ROOT.TFile.Open(outFile, "read") as f_outputfile,
+        ROOT.TFile.Open(inputFile, "read") as f_muonfile,
+        ROOT.TFile.Open(temp_filename, "recreate") as f_temp,
+    ):
+        output_tree = f_outputfile.Get("cbmsim")
+
+        muondis_tree = f_muonfile.Get("DIS")
+
+        new_tree = output_tree.CloneTree(0)
+
+        cross_section = array("f", [0.0])
+        cross_section_leaf = new_tree.Branch(
+            "CrossSection", cross_section, "CrossSection/F"
+        )
+
+        for output_event, muondis_event in zip(output_tree, muondis_tree):
+            mu = muondis_event.InMuon[0]
+            cross_section[0] = mu[10]
+            new_tree.Fill()
+
+        new_tree.Write("", ROOT.TObject.kOverwrite)
+
+    os.replace(temp_filename, outFile)
+    print("Successfully added DISCrossSection to the output file:", outFile)
+
+# ------------------------------------------------------------------------
+import checkMagFields
+
+
+def visualizeMagFields():
+    checkMagFields.run()
+
+
+def checkOverlapsWithGeant4():
+    # after /run/initialize, but prints warning messages, problems with TGeo volume
+    mygMC = ROOT.TGeant4.GetMC()
+    mygMC.ProcessGeantCommand("/geometry/test/recursion_start 0")
+    mygMC.ProcessGeantCommand("/geometry/test/recursion_depth 2")
+    mygMC.ProcessGeantCommand("/geometry/test/run")
diff --git a/shipgen/shipgen_parameters.yaml b/shipgen/shipgen_parameters.yaml
new file mode 100644
index 0000000000..3fa990c272
--- /dev/null
+++ b/shipgen/shipgen_parameters.yaml
@@ -0,0 +1,59 @@
+#### Config parameters for all FairPrimaryGenerators
+
+Cosmics:
+  opt: 1 #Use cosmic generator, argument switch for cosmic generator 0 or 1
+
+Genie:
+  inputFile: "/eos/experiment/ship/data/GenieEvents/genie-nu_mu.root"
+
+NuRadio:
+  inputFile: "/eos/experiment/ship/data/GenieEvents/genie-nu_mu.root"
+
+MuonDIS:
+  inputFile: "/eos/experiment/ship/data/muonDIS/muonDis_1.root"
+
+Nuage:
+  inputFile: 'Numucc.root'
+
+Ntuple:
+  inputFile: "/eos/experiment/ship/data/Mbias/pythia8_Geant4-withCharm_onlyMuons_4magTarget.root"
+
+MuonBack:
+  DownScaleDiMuon: False
+  phiRandom: True
+  sameSeed: 0 #can be set to an integer for the muonBackground simulation with specific seed for each muon, only for experts!"
+  followMuon: True
+  fastMuon: True
+
+Pythia8:
+  inclusive: "c"    # True = all processes if "c" only ccbar -> HNL, if "b" only bbar -> HNL, if "bc" only Bc+/Bc- -> HNL, and for darkphotons: if meson = production through meson decays, pbrem = proton bremstrahlung, qcd = ffbar -> DP.
+  inputFile: "/eos/experiment/ship/data/Charm/Cascade-parp16-MSTP82-1-MSEL4-978Bpot.root"
+
+  #inputFile: "/eos/experiment/ship/data/Beauty/Cascade-run0-19-parp16-MSTP82-1-MSEL5-5338Bpot.root" #for inclusive=b
+  #inputFile = "$FAIRSHIP/files/Cascade-parp16-MSTP82-1-MSEL4-76Mpot_1_5000.root" #test file
+
+  charmonly: False # option to be set with inclusive to enable only charm decays, charm x-sec measurement
+
+
+  HNL: True
+  DarkPhoton: False
+  RPVSUSY: False
+
+  theMass: 1.0 #in GeV, for relevant model (DP, RPV, HNL)
+  #For HNL: Use same values for prod and decay, or set separately.
+  theProductionCouplings: [0.447e-9,7.15e-9,1.88e-9] #"couplings \'U2e,U2mu,U2tau\' or -c \'U2e,U2mu,U2tau\' to set list of HNL couplings. TP default for HNL, ctau=53.3km"
+  theDecayCouplings: [0.447e-9,7.15e-9,1.88e-9] #"couplings \'U2e,U2mu,U2tau\' or -c \'U2e,U2mu,U2tau\' to set list of HNL couplings. TP default for HNL, ctau=53.3km"
+
+  theRPVCouplings: [0,0,0] #CAMM what are good defaults??
+  RPVSUSYbench: 2 #Generate HP Susy
+
+  motherMode: 'pi0' #"Choose DP production meson source: pi0, eta, omega, eta1, eta11"
+  theDPepsilon: 0.00000008
+
+EvtCalc:
+  inputFile:"" #CAMM: need good default?
+
+ParticleGun:
+  pID: 22 #id of particle used by the gun (default=22)
+  Estart: 10 #start of energy range of particle gun (default=10 GeV)
+  Eend: 10 #end of energy range of particle gun (default=10 GeV)