diff --git a/JobConfig/ensemble/fcl/prolog.fcl b/JobConfig/ensemble/fcl/prolog.fcl
index 2b75cde8..37be233a 100644
--- a/JobConfig/ensemble/fcl/prolog.fcl
+++ b/JobConfig/ensemble/fcl/prolog.fcl
@@ -5,7 +5,8 @@ Ensemble : {
     InputCommands : [ "keep *", 
           "keep mu2e::CosmicLivetime_*_*_*",
           "drop *_genCounter_*_*",
-          "drop *_protonBunchIntensity_*_*"]
+          "drop *_protonBunchIntensity_*_*",
+          "drop mu2e::SpectrumConfig_*_*_*"]
     OutputCommandsMC : [   "keep *_*_*_*"]
     OutputCommandsData : [   "drop *_*_*_*",
           "keep mu2e::KalSeeds_*_*_*",
diff --git a/JobConfig/ensemble/python/calculateEvents.py b/JobConfig/ensemble/python/calculateEvents.py
index d20a0091..f13827ad 100755
--- a/JobConfig/ensemble/python/calculateEvents.py
+++ b/JobConfig/ensemble/python/calculateEvents.py
@@ -1,4 +1,5 @@
 #! /usr/bin/env python
+import argparse
 from normalizations import *
 
 def main(args):
@@ -39,10 +40,18 @@ def main(args):
       print("ExternalRPC_yield=",Yield)
     if(args.prc == "RMC" and int(args.internal) == 1):
       Yield = rmc_normalization(float(args.livetime),  str(args.internal), float(args.rmcemin))
-      print("InternalRMC_yield=",Yield)
-    if(args.prc == "RMC" and int(args.internal) == 0):
-      Yield = rmc_normalization(float(args.livetime),  str(args.internal), float(args.rmcemin))
-      print("ExternalRMC_yield=",Yield)
+    if(args.prc == "RMCPhaseSpace0NExternal"):
+      Yield = rmc_0n_normalization(float(args.livetime), float(args.rmcn0emin), internal=0, run_mode=str(args.BB))
+      print("ExternalRMCPhaseSpace0N_yield=",Yield)
+    if(args.prc == "RMCPhaseSpace0NInternal"):
+      Yield = rmc_0n_normalization(float(args.livetime), float(args.rmcn0emin), internal=1, run_mode=str(args.BB))
+      print("InternalRMCPhaseSpace0N_yield=",Yield)
+    if(args.prc == "RMCPhaseSpace1NExternal"):
+      Yield = rmc_1n_normalization(float(args.livetime), float(args.rmcn1emin), internal=0, run_mode=str(args.BB))
+      print("ExternalRMCPhaseSpace1N_yield=",Yield)
+    if(args.prc == "RMCPhaseSpace1NInternal"):
+      Yield = rmc_1n_normalization(float(args.livetime), float(args.rmcn1emin), internal=1, run_mode=str(args.BB))
+      print("InternalRMCPhaseSpace1N_yield=",Yield)
     if(args.prc == "IPAMichel"):
       Yield = ipaMichel_normalization(float(args.livetime), float(args.ipaemin), str(args.BB))
       print("IPAMichel_yield=",Yield)
@@ -58,6 +67,8 @@ def main(args):
     parser.add_argument("--ipaemin", help="min energy cut dio ipa")
     parser.add_argument("--rpcemin", help="rpcemin", default=0)
     parser.add_argument("--rmcemin", help="min energy cut rmc")
+    parser.add_argument("--rmcn0emin", help="min energy cut rmc 0N")
+    parser.add_argument("--rmcn1emin", help="min energy cut rmc 1N")
     parser.add_argument("--prc", help="process")
     parser.add_argument("--printpot", help="print pot", default="no")
     parser.add_argument("--tmin", help="tmin", default=0)
diff --git a/JobConfig/ensemble/python/make_template_fcl.py b/JobConfig/ensemble/python/make_template_fcl.py
index d829af83..c2ec6406 100755
--- a/JobConfig/ensemble/python/make_template_fcl.py
+++ b/JobConfig/ensemble/python/make_template_fcl.py
@@ -32,17 +32,48 @@ def main(args):
 
   ROOT.gRandom.SetSeed(0)
 
-  # extract normalization of each background/signal process:
-  norms = {
-          "CRYCosmic": cry_onspill_normalization(livetime, args.BB),
-          "CORSIKACosmic": corsika_onspill_normalization(livetime, args.BB),
-          "DIO": dio_normalization(livetime, dioemin, args.BB),
-          "RPCInternal": rpc_normalization(livetime, args.tmin, 1, args.rpcemin, args.BB),
-          "RPCExternal": rpc_normalization(livetime, args.tmin, 0, args.rpcemin, args.BB),
-          "RMCInternal": rmc_normalization(livetime, 1, args.rmcemin, args.rmckmax, args.BB),
-          "RMCExternal": rmc_normalization(livetime, 0, args.rmcemin, args.rmckmax, args.BB),
-          "IPAMichel": ipaMichel_normalization(livetime, args.ipaemin, args.BB)
-          }
+  # Convert args.prc into a set for fast lookup
+  requested_processes = set(args.prc)
+  
+  # Initialize an empty dictionary
+  norms = {}
+
+  # Only call functions and define keys if they are in the requested processes
+  if "CRYCosmic" in requested_processes:
+      norms["CRYCosmic"] = cry_onspill_normalization(livetime, args.BB)
+      
+  if "CORSIKACosmic" in requested_processes:
+      norms["CORSIKACosmic"] = corsika_onspill_normalization(livetime, args.BB)
+      
+  if "DIO" in requested_processes:
+      norms["DIO"] = dio_normalization(livetime, dioemin, args.BB)
+      
+  if "RPCInternal" in requested_processes:
+      norms["RPCInternal"] = rpc_normalization(livetime, tmin, 1, args.rpcemin, args.BB)
+      
+  if "RPCExternal" in requested_processes:
+      norms["RPCExternal"] = rpc_normalization(livetime, tmin, 0, args.rpcemin, args.BB)
+      
+  if "RMCInternal" in requested_processes:
+      norms["RMCInternal"] = rmc_normalization(livetime, 1, args.rmcemin, args.rmckmax, args.BB)
+      
+  if "RMCExternal" in requested_processes:
+      norms["RMCExternal"] = rmc_normalization(livetime, 0, args.rmcemin, args.rmckmax, args.BB)
+      
+  if "RMCN0External" in requested_processes or "RMCPhaseSpace0NExternal" in requested_processes:
+      norms["RMCN0External"] = norms["RMCPhaseSpace0NExternal"] = rmc_0n_normalization(livetime, args.rmcn0emin, internal=0, run_mode=args.BB)
+      
+  if "RMCN0Internal" in requested_processes or "RMCPhaseSpace0NInternal" in requested_processes:
+      norms["RMCN0Internal"] = norms["RMCPhaseSpace0NInternal"] = rmc_0n_normalization(livetime, args.rmcn0emin, internal=1, run_mode=args.BB)
+      
+  if "RMCN1External" in requested_processes or "RMCPhaseSpace1NExternal" in requested_processes:
+      norms["RMCN1External"] = norms["RMCPhaseSpace1NExternal"] = rmc_1n_normalization(livetime, args.rmcn1emin, internal=0, run_mode=args.BB)
+      
+  if "RMCN1Internal" in requested_processes or "RMCPhaseSpace1NInternal" in requested_processes:
+      norms["RMCN1Internal"] = norms["RMCPhaseSpace1NInternal"] = rmc_1n_normalization(livetime, args.rmcn1emin, internal=1, run_mode=args.BB)
+      
+  if "IPAMichel" in requested_processes:
+      norms["IPAMichel"] = ipaMichel_normalization(livetime, args.ipaemin, args.BB)
 
   starting_event_num = {}
   max_possible_events = {}
@@ -210,6 +241,8 @@ def main(args):
     parser.add_argument("--rpcemin", help="min energy cut rpc")
     parser.add_argument("--ipaemin", help="min energy cut ipa")
     parser.add_argument("--rmcemin", help="min energy cut rmc")
+    parser.add_argument("--rmcn0emin", help="min energy cut rmc 0N")
+    parser.add_argument("--rmcn1emin", help="min energy cut rmc 1N")
     parser.add_argument("--rmckmax", help="kmax theory value")
     parser.add_argument("--run", help="run number")
     parser.add_argument("--samplingseed", help="samplingseed")
diff --git a/JobConfig/ensemble/python/normalizations.py b/JobConfig/ensemble/python/normalizations.py
index d768fd67..e016a4b2 100755
--- a/JobConfig/ensemble/python/normalizations.py
+++ b/JobConfig/ensemble/python/normalizations.py
@@ -58,6 +58,32 @@ def set_verbose(verbose=True):
 dutyfactor = 1.0
 total_pot = 0.
 
+# RMC 0N and 1N Physics Constants
+RMC_BR_MUON_CAPTURE = 0.609
+RMC_RATE_GT_57 = 1.41e-5  # RMC rate above 57 MeV, relative to OMC
+RMC_BR_0N_FRAC_GT_57 = 0.099  # BR(0 knockout | E > 57) / BR(RMC | E > 57)
+RMC_BR_1N_FRAC_GT_57 = 0.901  # BR(1 knockout | E > 57) / BR(RMC | E > 57)
+
+# RMC K_max values: Energy endpoints for each knockout mode on Al-27
+RMC_KMAX_0N = 101.8667  # MeV, 0-nucleon knockout endpoint on Al-27
+RMC_KMAX_1N = 95.4489   # MeV, 1-nucleon knockout endpoint on Al-27
+
+# RMC Spectrum Fractions - from experimental/theoretical physics literature
+# These represent: R(* knockout | E > threshold) / R(* knockout) 
+# i.e., the fraction of the full spectrum above the given energy threshold
+# Used to scale branching ratios as a function of energy cut:
+#   BR(* knockout | E > E_min) = BR(muon capture) * RMC_RATE_GT_57 * BR(* frac | E > 57) * (R_*_E_min / R_*_57)
+# Note: These are the RAW spectrum fractions (not pre-multiplied by RMC_BR_MUON_CAPTURE)
+RMC_SPECTRUM_FRAC_0N_57 = 0.22887  # Fraction of 0-knockout spectrum above 57 MeV
+RMC_SPECTRUM_FRAC_1N_57 = 0.061620  # Fraction of 1-knockout spectrum above 57 MeV
+RMC_SPECTRUM_FRAC_0N_80 = 0.03319  # Fraction of 0-knockout spectrum above 80 MeV
+RMC_SPECTRUM_FRAC_1N_80 = 0.0013175  # Fraction of 1-knockout spectrum above 80 MeV
+
+# Internal/external conversion ratio for RMC
+# rho = BR(internal) / BR(external)
+# Note: Should ideally use Plestid-Hill or Kroll-Wada-Joseph integrals for precision
+RMC_INTERNAL_EXTERNAL_RATIO = 0.0069  # rho = BR(internal) / BR(external)
+
 #-------------------------------------------------------------------------------------#
 
 # --- Database Interaction ---
@@ -118,6 +144,120 @@ def set_verbose(verbose=True):
         ipa_stopping_rate = ipa_stopping_rate * float(words[3])
         ipa_stopped_mu_per_POT = ipa_stopping_rate
 #print("IPAStopMuonRate=", ipa_stopped_mu_per_POT)
+
+#-------------------------------------------------------------------------------------#
+# Plestid spectrum functions for RMC 0N and 1N
+
+def plestid_integral(K_1, K_2, KMax, knockout):
+    """
+    Calculates the integral of the Plestid phase-space approximation spectrum
+    between two energy points K_1 and K_2.
+    
+    This implements the Plestid spectrum shape for RMC with different knockout modes.
+    
+    Args:
+        K_1 (float): Lower energy bound for integration (MeV).
+        K_2 (float): Upper energy bound for integration (MeV).
+        KMax (float): Maximum possible RMC energy (MeV).
+        knockout (int): Knockout mode (0 for 0N knockout, 1 for 1N knockout).
+    
+    Returns:
+        float: The integral of the spectrum between K_1 and K_2.
+    """
+    if KMax <= 0.0:
+        return 0.0
+    if knockout < 0:
+        return 0.0
+    
+    K_1 = max(0.0, min(KMax, K_1))
+    K_2 = max(0.0, min(KMax, K_2))
+    
+    if K_1 >= K_2:
+        return 0.0
+    
+    power = 2.0 + 1.5 * knockout
+    x_1 = K_1 / KMax
+    x_2 = K_2 / KMax
+    
+    val_1 = (x_1 - 1.0) * pow(1.0 - x_1, power) * (power * x_1 + x_1 + 1.0)
+    val_2 = (x_2 - 1.0) * pow(1.0 - x_2, power) * (power * x_2 + x_2 + 1.0)
+    
+    integral = val_2 - val_1
+    return integral
+
+
+def plestid_spectrum(energy, kmax, knockout):
+    """
+    Calculates the Plestid phase-space approximation spectrum value at a given energy.
+    
+    Args:
+        energy (float): Energy point to evaluate the spectrum (MeV).
+        kmax (float): Maximum possible RMC energy (MeV).
+        knockout (int): Knockout mode (0 for 0N knockout, 1 for 1N knockout).
+    
+    Returns:
+        float: The spectrum value at the given energy.
+    """
+    if energy <= 0.0 or energy >= kmax:
+        return 0.0
+    
+    power = 2.0 + 1.5 * knockout
+    norm = (power + 1.0) * (power + 2.0) / kmax
+    x = energy / kmax
+    p = norm * x * pow(1.0 - x, power)
+    
+    return p
+
+
+def compute_rmc_spectrum_fractions():
+    """
+    Compute RMC spectrum fractions dynamically using plestid_integral.
+    
+    Uses the correct K_max values for each knockout mode:
+    - kmax_0n = 101.8667 MeV (0-nucleon knockout endpoint on Al-27)
+    - kmax_1n = 95.4489 MeV (1-nucleon knockout endpoint on Al-27)
+    
+    These fractions represent: R(* knockout | E > threshold) / R(* knockout)
+    """
+    global RMC_SPECTRUM_FRAC_0N_57, RMC_SPECTRUM_FRAC_1N_57
+    global RMC_SPECTRUM_FRAC_0N_80, RMC_SPECTRUM_FRAC_1N_80
+    
+    # Compute integrals for each knockout mode using its own K_max
+    frac_0_0n  = plestid_integral(0.0, RMC_KMAX_0N, RMC_KMAX_0N, 0)
+    frac_0_1n  = plestid_integral(0.0, RMC_KMAX_1N, RMC_KMAX_1N, 1)
+    frac_57_0n = plestid_integral(57.0, RMC_KMAX_0N, RMC_KMAX_0N, 0)
+    frac_57_1n = plestid_integral(57.0, RMC_KMAX_1N, RMC_KMAX_1N, 1)
+    frac_80_0n = plestid_integral(80.0, RMC_KMAX_0N, RMC_KMAX_0N, 0)
+    frac_80_1n = plestid_integral(80.0, RMC_KMAX_1N, RMC_KMAX_1N, 1)
+    
+    # Compute ratios: integral above threshold / integral from 0 to K_max
+    RMC_SPECTRUM_FRAC_0N_57 = frac_57_0n / frac_0_0n if frac_0_0n != 0 else 0.0
+    RMC_SPECTRUM_FRAC_1N_57 = frac_57_1n / frac_0_1n if frac_0_1n != 0 else 0.0
+    RMC_SPECTRUM_FRAC_0N_80 = frac_80_0n / frac_0_0n if frac_0_0n != 0 else 0.0
+    RMC_SPECTRUM_FRAC_1N_80 = frac_80_1n / frac_0_1n if frac_0_1n != 0 else 0.0
+    
+    if VERBOSE:
+        print(f"RMC Spectrum Fractions (computed via Plestid integral):")
+        print(f"  RMC_SPECTRUM_FRAC_0N_57 = {RMC_SPECTRUM_FRAC_0N_57:.6f}")
+        print(f"  RMC_SPECTRUM_FRAC_1N_57 = {RMC_SPECTRUM_FRAC_1N_57:.6f}")
+        print(f"  RMC_SPECTRUM_FRAC_0N_80 = {RMC_SPECTRUM_FRAC_0N_80:.6f}")
+        print(f"  RMC_SPECTRUM_FRAC_1N_80 = {RMC_SPECTRUM_FRAC_1N_80:.6f}")
+    
+    return (RMC_SPECTRUM_FRAC_0N_57, RMC_SPECTRUM_FRAC_1N_57, 
+            RMC_SPECTRUM_FRAC_0N_80, RMC_SPECTRUM_FRAC_1N_80)
+
+
+def rmc_spectrum_fraction(e_min, k_max, knockout):
+    """Return the RMC spectrum fraction above ``e_min`` relative to E > 57 MeV."""
+    reference_fraction = plestid_integral(57.0, k_max, k_max, knockout)
+    requested_fraction = plestid_integral(float(e_min), k_max, k_max, knockout)
+
+    if reference_fraction == 0.0:
+        return 0.0
+    return requested_fraction / reference_fraction
+
+# Compute RMC spectrum fractions from Plestid integral using correct K_max values
+compute_rmc_spectrum_fractions()
     
 #-------------------------------------------------------------------------------------#    
 def get_duty_factor(run_mode='1BB'):
@@ -398,7 +538,6 @@ def rpc_normalization(on_spill_time, t_min, internal, e_min, run_mode='1BB'):
     base_physics_events = (
         total_pot *
         target_stopped_pions_per_pot *
-        filter_efficiency *
         survival_probability_weight *
         RPC_PER_STOPPED_PION *
         rpc_e_sample_frac
@@ -500,7 +639,112 @@ def rmc_normalization(on_spill_time, internal, e_min, k_max=90.1, run_mode='1BB'
         
     return base_physics_events
 
-# get IPA Michel normalization:
+
+def rmc_0n_normalization(on_spill_time, e_min=80.0, internal=1, run_mode='1BB'):
+    """
+    Calculates the expected number of RMC 0-nucleon knockout (0N) events
+    above a given energy threshold.
+    
+    Uses the Plestid phase-space approximation spectrum shape.
+    
+    Args:
+        on_spill_time (float): Time the beam was on spill (seconds).
+        e_min (float): Minimum energy threshold for the spectrum cut (MeV).
+        internal (int/bool): Flag (1 or 0) to include internal conversion scaling.
+                            Defaults to 1.
+        run_mode (str): The operational mode ('1BB' or '2BB'). Defaults to '1BB'.
+    
+    Returns:
+        float: The expected number of RMC 0N physics events passing the cuts.
+    """
+    # 1. Calculate total Protons on Target (POT)
+    total_pot = get_pot(on_spill_time, run_mode)
+    
+    # 2. Determine the spectrum fraction above the requested threshold.
+    R_spectrum = rmc_spectrum_fraction(e_min, RMC_KMAX_0N, knockout=0)
+    
+    # 3. Calculate the branching ratio for 0N events above the energy threshold
+    br_0n_above_emin = (
+        RMC_BR_MUON_CAPTURE *
+        RMC_RATE_GT_57 *
+        RMC_BR_0N_FRAC_GT_57 *
+        R_spectrum
+    )
+    
+    # 4. Calculate base physics events
+    # Note: br_0n_above_emin already includes RMC_BR_MUON_CAPTURE, so we do NOT multiply by CAPTURES_PER_STOPPED_MUON
+    base_physics_events = (
+        total_pot *
+        target_stopped_muons_per_pot *
+        br_0n_above_emin
+    )
+    
+    # 5. Apply internal conversion scaling if requested
+    is_internal_conversion = bool(int(internal))
+    
+    if is_internal_conversion:
+        if VERBOSE:
+            print("RMC_0N_emin=", e_min)
+            print("RMC_0N_spectrum_frac=", R_spectrum)
+            print("RMC_0N_BR=", br_0n_above_emin)
+        
+        base_physics_events *= RMC_INTERNAL_EXTERNAL_RATIO
+    
+    return base_physics_events
+
+
+def rmc_1n_normalization(on_spill_time, e_min=80.0, internal=1, run_mode='1BB'):
+    """
+    Calculates the expected number of RMC 1-nucleon knockout (1N) events
+    above a given energy threshold.
+    
+    Uses the Plestid phase-space approximation spectrum shape.
+    
+    Args:
+        on_spill_time (float): Time the beam was on spill (seconds).
+        e_min (float): Minimum energy threshold for the spectrum cut (MeV).
+        internal (int/bool): Flag (1 or 0) to include internal conversion scaling.
+                            Defaults to 1.
+        run_mode (str): The operational mode ('1BB' or '2BB'). Defaults to '1BB'.
+    
+    Returns:
+        float: The expected number of RMC 1N physics events passing the cuts.
+    """
+    # 1. Calculate total Protons on Target (POT)
+    total_pot = get_pot(on_spill_time, run_mode)
+    
+    # 2. Determine the spectrum fraction above the requested threshold.
+    R_spectrum = rmc_spectrum_fraction(e_min, RMC_KMAX_1N, knockout=1)
+    
+    # 3. Calculate the branching ratio for 1N events above the energy threshold
+    br_1n_above_emin = (
+        RMC_BR_MUON_CAPTURE *
+        RMC_RATE_GT_57 *
+        RMC_BR_1N_FRAC_GT_57 *
+        R_spectrum
+    )
+    
+    # 4. Calculate base physics events
+    # Note: br_1n_above_emin already includes RMC_BR_MUON_CAPTURE, so we do NOT multiply by CAPTURES_PER_STOPPED_MUON
+    base_physics_events = (
+        total_pot *
+        target_stopped_muons_per_pot *
+        br_1n_above_emin
+    )
+    
+    # 5. Apply internal conversion scaling if requested
+    is_internal_conversion = bool(int(internal))
+    
+    if is_internal_conversion:
+        if VERBOSE:
+            print("RMC_1N_emin=", e_min)
+            print("RMC_1N_spectrum_frac=", R_spectrum)
+            print("RMC_1N_BR=", br_1n_above_emin)
+        
+        base_physics_events *= RMC_INTERNAL_EXTERNAL_RATIO
+    
+    return base_physics_events
+
 def ipaMichel_normalization(on_spill_time, ipa_de_min, run_mode='1BB'):
     """
     Calculates the expected number of IPA (Incoming Particle Decay After Stopping)
@@ -589,4 +833,4 @@ def corsika_onspill_normalization(livetime, run_mode = '1BB'):
   tst_1BB = get_pot(9.52e6)
   tst_2BB = get_pot(1.58e6)
   tst_rpc = rpc_normalization(3.77e19,350,1,1)
-  print("SU2020", tst_1BB, tst_2BB)
+  print("SU2020", tst_1BB, tst_2BB)
\ No newline at end of file
diff --git a/JobConfig/ensemble/python/test_normalizations.py b/JobConfig/ensemble/python/test_normalizations.py
index df54c929..2f147252 100644
--- a/JobConfig/ensemble/python/test_normalizations.py
+++ b/JobConfig/ensemble/python/test_normalizations.py
@@ -25,8 +25,8 @@ def setUpClass(cls):
             os.environ['MUSE_WORK_DIR'] = '/exp/mu2e/app/users/sophie/newOffline'
 
         # Test livetimes from MDS3 (cosmic livetime)
-        cls.on_spill_time_1BB = 4.4e6  # seconds
-        cls.on_spill_time_2BB = 4.4e6  # seconds
+        cls.on_spill_time_1BB = 4.4e6 * (88/496) # seconds
+        cls.on_spill_time_2BB = 4.4e6 * (88/496)  # seconds
         
         cls.tst_1BB_cycle = normalizations.get_pot(cls.on_spill_time_1BB, run_mode='1BB',printout=False,method='cycle')
         cls.tst_2BB_cycle = normalizations.get_pot(cls.on_spill_time_2BB, run_mode='2BB',printout=False,method='cycle')
@@ -201,6 +201,78 @@ def test_rmc_internal_yields(self):
         self.assertGreater(rmc_yield_1bb, 0)
         self.assertGreater(rmc_yield_2bb, 0)
 
+    def test_rmc_0n_external_yields(self):
+        """Calculate RMC 0-nucleon knockout (external) yields for both 1BB and 2BB POT values."""
+        e_min = 80  # MeV
+        
+        rmc_0n_yield_1bb = normalizations.rmc_0n_normalization(
+            self.on_spill_time_1BB, e_min, internal=0, run_mode='1BB'
+        )
+        rmc_0n_yield_2bb = normalizations.rmc_0n_normalization(
+            self.on_spill_time_2BB, e_min, internal=0, run_mode='2BB'
+        )
+        
+        print(f"\nRMC 0N External Yields (above {e_min} MeV):")
+        print(f"  1BB: {rmc_0n_yield_1bb:.2e}")
+        print(f"  2BB: {rmc_0n_yield_2bb:.2e}")
+        
+        self.assertGreater(rmc_0n_yield_1bb, 0)
+        self.assertGreater(rmc_0n_yield_2bb, 0)
+
+    def test_rmc_0n_internal_yields(self):
+        """Calculate RMC 0-nucleon knockout (internal conversion) yields for both 1BB and 2BB POT values."""
+        e_min = 80  # MeV
+        
+        rmc_0n_yield_1bb = normalizations.rmc_0n_normalization(
+            self.on_spill_time_1BB, e_min, internal=1, run_mode='1BB'
+        )
+        rmc_0n_yield_2bb = normalizations.rmc_0n_normalization(
+            self.on_spill_time_2BB, e_min, internal=1, run_mode='2BB'
+        )
+        
+        print(f"\nRMC 0N Internal Yields (above {e_min} MeV):")
+        print(f"  1BB: {rmc_0n_yield_1bb:.2e}")
+        print(f"  2BB: {rmc_0n_yield_2bb:.2e}")
+        
+        self.assertGreater(rmc_0n_yield_1bb, 0)
+        self.assertGreater(rmc_0n_yield_2bb, 0)
+
+    def test_rmc_1n_external_yields(self):
+        """Calculate RMC 1-nucleon knockout (external) yields for both 1BB and 2BB POT values."""
+        e_min = 80  # MeV
+        
+        rmc_1n_yield_1bb = normalizations.rmc_1n_normalization(
+            self.on_spill_time_1BB, e_min, internal=0, run_mode='1BB'
+        )
+        rmc_1n_yield_2bb = normalizations.rmc_1n_normalization(
+            self.on_spill_time_2BB, e_min, internal=0, run_mode='2BB'
+        )
+        
+        print(f"\nRMC 1N External Yields (above {e_min} MeV):")
+        print(f"  1BB: {rmc_1n_yield_1bb:.2e}")
+        print(f"  2BB: {rmc_1n_yield_2bb:.2e}")
+        
+        self.assertGreater(rmc_1n_yield_1bb, 0)
+        self.assertGreater(rmc_1n_yield_2bb, 0)
+
+    def test_rmc_1n_internal_yields(self):
+        """Calculate RMC 1-nucleon knockout (internal conversion) yields for both 1BB and 2BB POT values."""
+        e_min = 80  # MeV
+        
+        rmc_1n_yield_1bb = normalizations.rmc_1n_normalization(
+            self.on_spill_time_1BB, e_min, internal=1, run_mode='1BB'
+        )
+        rmc_1n_yield_2bb = normalizations.rmc_1n_normalization(
+            self.on_spill_time_2BB, e_min, internal=1, run_mode='2BB'
+        )
+        
+        print(f"\nRMC 1N Internal Yields (above {e_min} MeV):")
+        print(f"  1BB: {rmc_1n_yield_1bb:.2e}")
+        print(f"  2BB: {rmc_1n_yield_2bb:.2e}")
+        
+        self.assertGreater(rmc_1n_yield_1bb, 0)
+        self.assertGreater(rmc_1n_yield_2bb, 0)
+
     def test_ipa_michel_normalization_yields(self):
         """Calculate IPA Michel yields for both 1BB and 2BB POT values."""
         ipa_de_min = 50.0  # MeV
@@ -236,16 +308,24 @@ def test_all_yields_summary(self):
         rpc_int_1bb = normalizations.rpc_normalization(self.on_spill_time_1BB, 350, 1, 50.0, run_mode='1BB')
         rmc_ext_1bb = normalizations.rmc_normalization(self.on_spill_time_1BB, 0, 85, run_mode='1BB')
         rmc_int_1bb = normalizations.rmc_normalization(self.on_spill_time_1BB, 1, 85, run_mode='1BB')
+        rmc_0n_ext_1bb = normalizations.rmc_0n_normalization(self.on_spill_time_1BB, 80, internal=0, run_mode='1BB')
+        rmc_0n_int_1bb = normalizations.rmc_0n_normalization(self.on_spill_time_1BB, 80, internal=1, run_mode='1BB')
+        rmc_1n_ext_1bb = normalizations.rmc_1n_normalization(self.on_spill_time_1BB, 80, internal=0, run_mode='1BB')
+        rmc_1n_int_1bb = normalizations.rmc_1n_normalization(self.on_spill_time_1BB, 80, internal=1, run_mode='1BB')
         ipa_1bb = normalizations.ipaMichel_normalization(self.on_spill_time_1BB, 50.0, run_mode='1BB')
         
         print(f"\n  Process Yields:")
-        print(f"    CE (RUE=1e-13):            {ce_1bb:.3e}")
-        print(f"    DIO (E>95 MeV):            {dio_1bb:.3e}")
-        print(f"    RPC External (E>50 MeV):   {rpc_ext_1bb:.3e}")
-        print(f"    RPC Internal (E>50 MeV):   {rpc_int_1bb:.3e}")
-        print(f"    RMC External (E>85 MeV):   {rmc_ext_1bb:.3e}")
-        print(f"    RMC Internal (E>85 MeV):   {rmc_int_1bb:.3e}")
-        print(f"    IPA Michel (E>50 MeV):     {ipa_1bb:.3e}")
+        print(f"    CE (RUE=1e-13):                  {ce_1bb:.3e}")
+        print(f"    DIO (E>95 MeV):                  {dio_1bb:.3e}")
+        print(f"    RPC External (E>50 MeV):        {rpc_ext_1bb:.3e}")
+        print(f"    RPC Internal (E>50 MeV):        {rpc_int_1bb:.3e}")
+        print(f"    RMC External (E>85 MeV):        {rmc_ext_1bb:.3e}")
+        print(f"    RMC Internal (E>85 MeV):        {rmc_int_1bb:.3e}")
+        print(f"    RMC 0N External (E>80 MeV):     {rmc_0n_ext_1bb:.3e}")
+        print(f"    RMC 0N Internal (E>80 MeV):     {rmc_0n_int_1bb:.3e}")
+        print(f"    RMC 1N External (E>80 MeV):     {rmc_1n_ext_1bb:.3e}")
+        print(f"    RMC 1N Internal (E>80 MeV):     {rmc_1n_int_1bb:.3e}")
+        print(f"    IPA Michel (E>50 MeV):          {ipa_1bb:.3e}")
         
         # 2BB mode
         print(f"\n2BB Mode (POT: {self.tst_2BB_spill:.3e}):")
@@ -258,16 +338,24 @@ def test_all_yields_summary(self):
         rpc_int_2bb = normalizations.rpc_normalization(self.on_spill_time_2BB, 350, 1, 50.0, run_mode='2BB')
         rmc_ext_2bb = normalizations.rmc_normalization(self.on_spill_time_2BB, 0, 85, run_mode='2BB')
         rmc_int_2bb = normalizations.rmc_normalization(self.on_spill_time_2BB, 1, 85, run_mode='2BB')
+        rmc_0n_ext_2bb = normalizations.rmc_0n_normalization(self.on_spill_time_2BB, 80, internal=0, run_mode='2BB')
+        rmc_0n_int_2bb = normalizations.rmc_0n_normalization(self.on_spill_time_2BB, 80, internal=1, run_mode='2BB')
+        rmc_1n_ext_2bb = normalizations.rmc_1n_normalization(self.on_spill_time_2BB, 80, internal=0, run_mode='2BB')
+        rmc_1n_int_2bb = normalizations.rmc_1n_normalization(self.on_spill_time_2BB, 80, internal=1, run_mode='2BB')
         ipa_2bb = normalizations.ipaMichel_normalization(self.on_spill_time_2BB, 50.0, run_mode='2BB')
         
         print(f"\n  Process Yields:")
-        print(f"    CE (RUE=1e-13):            {ce_2bb:.3e}")
-        print(f"    DIO (E>95 MeV):            {dio_2bb:.3e}")
-        print(f"    RPC External (E>50 MeV):   {rpc_ext_2bb:.3e}")
-        print(f"    RPC Internal (E>50 MeV):   {rpc_int_2bb:.3e}")
-        print(f"    RMC External (E>85 MeV):   {rmc_ext_2bb:.3e}")
-        print(f"    RMC Internal (E>85 MeV):   {rmc_int_2bb:.3e}")
-        print(f"    IPA Michel (E>50 MeV):     {ipa_2bb:.3e}")
+        print(f"    CE (RUE=1e-13):                  {ce_2bb:.3e}")
+        print(f"    DIO (E>95 MeV):                  {dio_2bb:.3e}")
+        print(f"    RPC External (E>50 MeV):        {rpc_ext_2bb:.3e}")
+        print(f"    RPC Internal (E>50 MeV):        {rpc_int_2bb:.3e}")
+        print(f"    RMC External (E>85 MeV):        {rmc_ext_2bb:.3e}")
+        print(f"    RMC Internal (E>85 MeV):        {rmc_int_2bb:.3e}")
+        print(f"    RMC 0N External (E>80 MeV):     {rmc_0n_ext_2bb:.3e}")
+        print(f"    RMC 0N Internal (E>80 MeV):     {rmc_0n_int_2bb:.3e}")
+        print(f"    RMC 1N External (E>80 MeV):     {rmc_1n_ext_2bb:.3e}")
+        print(f"    RMC 1N Internal (E>80 MeV):     {rmc_1n_int_2bb:.3e}")
+        print(f"    IPA Michel (E>50 MeV):          {ipa_2bb:.3e}")
         
         print(f"\n{'='*70}\n")
 
diff --git a/JobConfig/ensemble/python/test_rmc_fractions.py b/JobConfig/ensemble/python/test_rmc_fractions.py
new file mode 100644
index 00000000..954c7839
--- /dev/null
+++ b/JobConfig/ensemble/python/test_rmc_fractions.py
@@ -0,0 +1,76 @@
+#!/usr/bin/env python3
+"""
+Test RMC spectrum fractions computation using correct K_max values.
+"""
+
+import math
+
+RMC_KMAX_0N = 101.8667  # MeV, 0-nucleon knockout endpoint on Al-27
+RMC_KMAX_1N = 95.4489   # MeV, 1-nucleon knockout endpoint on Al-27
+
+def plestid_integral(K_1, K_2, KMax, knockout):
+    """Calculates the integral of the Plestid phase-space approximation spectrum."""
+    if KMax <= 0.0:
+        return 0.0
+    if knockout < 0:
+        return 0.0
+    
+    K_1 = max(0.0, min(KMax, K_1))
+    K_2 = max(0.0, min(KMax, K_2))
+    
+    if K_1 >= K_2:
+        return 0.0
+    
+    power = 2.0 + 1.5 * knockout
+    x_1 = K_1 / KMax
+    x_2 = K_2 / KMax
+    
+    val_1 = (x_1 - 1.0) * pow(1.0 - x_1, power) * (power * x_1 + x_1 + 1.0)
+    val_2 = (x_2 - 1.0) * pow(1.0 - x_2, power) * (power * x_2 + x_2 + 1.0)
+    
+    integral = val_2 - val_1
+    return integral
+
+# Compute integrals for each knockout mode using its own K_max
+frac_0_0n  = plestid_integral(0.0, RMC_KMAX_0N, RMC_KMAX_0N, 0)
+frac_0_1n  = plestid_integral(0.0, RMC_KMAX_1N, RMC_KMAX_1N, 1)
+frac_57_0n = plestid_integral(57.0, RMC_KMAX_0N, RMC_KMAX_0N, 0)
+frac_57_1n = plestid_integral(57.0, RMC_KMAX_1N, RMC_KMAX_1N, 1)
+frac_80_0n = plestid_integral(80.0, RMC_KMAX_0N, RMC_KMAX_0N, 0)
+frac_80_1n = plestid_integral(80.0, RMC_KMAX_1N, RMC_KMAX_1N, 1)
+
+# Compute ratios: integral above threshold / integral from 0 to K_max
+RMC_SPECTRUM_FRAC_0N_57 = frac_57_0n / frac_0_0n if frac_0_0n != 0 else 0.0
+RMC_SPECTRUM_FRAC_1N_57 = frac_57_1n / frac_0_1n if frac_0_1n != 0 else 0.0
+RMC_SPECTRUM_FRAC_0N_80 = frac_80_0n / frac_0_0n if frac_0_0n != 0 else 0.0
+RMC_SPECTRUM_FRAC_1N_80 = frac_80_1n / frac_0_1n if frac_0_1n != 0 else 0.0
+
+# Expected values from the hardcoded constants
+expected = {
+    'RMC_SPECTRUM_FRAC_0N_57': 0.22887,
+    'RMC_SPECTRUM_FRAC_1N_57': 0.061620,
+    'RMC_SPECTRUM_FRAC_0N_80': 0.03319,
+    'RMC_SPECTRUM_FRAC_1N_80': 0.0013175
+}
+
+print("=" * 80)
+print("RMC Spectrum Fractions: Computed vs. Expected")
+print("=" * 80)
+print(f"{'Constant':<25} {'Computed':<15} {'Expected':<15} {'Difference':<15}")
+print("-" * 80)
+
+computed_values = {
+    'RMC_SPECTRUM_FRAC_0N_57': RMC_SPECTRUM_FRAC_0N_57,
+    'RMC_SPECTRUM_FRAC_1N_57': RMC_SPECTRUM_FRAC_1N_57,
+    'RMC_SPECTRUM_FRAC_0N_80': RMC_SPECTRUM_FRAC_0N_80,
+    'RMC_SPECTRUM_FRAC_1N_80': RMC_SPECTRUM_FRAC_1N_80
+}
+
+for key in computed_values:
+    comp = computed_values[key]
+    exp = expected[key]
+    diff = abs(comp - exp)
+    rel_diff = diff / exp * 100 if exp != 0 else 0
+    print(f"{key:<25} {comp:<15.8f} {exp:<15.8f} {diff:.2e} ({rel_diff:.2f}%)")
+
+print("=" * 80)
diff --git a/JobConfig/ensemble/scripts/Stage1_initiate_ensemble.sh b/JobConfig/ensemble/scripts/Stage1_initiate_ensemble.sh
index c3430fde..d82e74db 100755
--- a/JobConfig/ensemble/scripts/Stage1_initiate_ensemble.sh
+++ b/JobConfig/ensemble/scripts/Stage1_initiate_ensemble.sh
@@ -21,6 +21,8 @@ STOPS="MDC2025ac"
 RELEASE="MDC2025"
 VERSION="ac"
 GEN="Signal" #cosmic generator name CRY or CORSIKA only Cat = "Signal"
+INCLUDE_RMCN0=1 # Include RMC 0N processes (default: no)
+INCLUDE_RMCN1=1 # Include RMC 1N processes (default: no)
 # Loop: Get the next option;
 while getopts ":-:" options; do
   case "${options}" in
@@ -59,6 +61,12 @@ while getopts ":-:" options; do
         gen)
           GEN=${!OPTIND} OPTIND=$(( $OPTIND + 1 ))
           ;;
+        rmcn0)
+          INCLUDE_RMCN0=${!OPTIND} OPTIND=$(( $OPTIND + 1 ))
+          ;;
+        rmcn1)
+          INCLUDE_RMCN1=${!OPTIND} OPTIND=$(( $OPTIND + 1 ))
+          ;;
         *)
           echo "Unknown option " ${OPTARG}
           exit_abnormal
@@ -114,7 +122,8 @@ BEAM_NMOT=$(echo "${BEAM_INFO}" | grep "^NMOT=" | awk '{print $NF}')
 echo "      • POT: ${BEAM_POT}"
 # Energy cut parameters (hardcoded for now)
 RPC_EMIN=50
-RMC_EMIN=85
+RMC_N0_EMIN=80
+RMC_N1_EMIN=80
 RMC_kmax=90.1
 IPA_EMIN=70
 # Extract just the numeric values from event yields (remove labels and spaces)
@@ -126,10 +135,18 @@ echo "      • Calculating RPC Internal events (emin=${RPC_EMIN})..."
 RPC_INTERNAL_EVENTS=$(calculateEvents.py --livetime ${LIVETIME} --prc "RPC" --tmin ${TMIN} --internal 1 --rpcemin ${RPC_EMIN} --BB ${BB} --printpot "no" --verbose false 2>/dev/null | tail -1 | awk '{print $NF}')
 echo "      • Calculating RPC External events (emin=${RPC_EMIN})..."
 RPC_EXTERNAL_EVENTS=$(calculateEvents.py --livetime ${LIVETIME} --prc "RPC" --tmin ${TMIN} --internal 0 --rpcemin ${RPC_EMIN} --BB ${BB} --printpot "no" --verbose false 2>/dev/null | tail -1 | awk '{print $NF}')
-echo "      • Calculating RMC Internal events (emin=${RMC_EMIN})..."
-RMC_INTERNAL_EVENTS=$(calculateEvents.py --livetime ${LIVETIME} --prc "RMC" --tmin ${TMIN} --internal 1 --rmcemin ${RMC_EMIN} --BB ${BB} --printpot "no" --verbose false 2>/dev/null | tail -1 | awk '{print $NF}')
-echo "      • Calculating RMC External events (emin=${RMC_EMIN})..."
-RMC_EXTERNAL_EVENTS=$(calculateEvents.py --livetime ${LIVETIME} --prc "RMC" --tmin ${TMIN} --internal 0 --rmcemin ${RMC_EMIN} --BB ${BB} --printpot "no" --verbose false 2>/dev/null | tail -1 | awk '{print $NF}')
+if [[ ${INCLUDE_RMCN0} -eq 1 ]]; then
+  echo "      • Calculating RMC 0N External events (emin=${RMC_N0_EMIN})..."
+  RMC_N0_EXTERNAL_EVENTS=$(calculateEvents.py --livetime ${LIVETIME} --prc "RMCPhaseSpace0NExternal" --internal 0 --rmcn0emin ${RMC_N0_EMIN} --BB ${BB} --printpot "no" --verbose false 2>/dev/null | tail -1 | awk '{print $NF}')
+  echo "      • Calculating RMC 0N Internal events (emin=${RMC_N0_EMIN})..."
+  RMC_N0_INTERNAL_EVENTS=$(calculateEvents.py --livetime ${LIVETIME} --prc "RMCPhaseSpace0NInternal" --internal 1 --rmcn0emin ${RMC_N0_EMIN} --BB ${BB} --printpot "no" --verbose false 2>/dev/null | tail -1 | awk '{print $NF}')
+fi
+if [[ ${INCLUDE_RMCN1} -eq 1 ]]; then
+  echo "      • Calculating RMC 1N External events (emin=${RMC_N1_EMIN})..."
+  RMC_N1_EXTERNAL_EVENTS=$(calculateEvents.py --livetime ${LIVETIME} --prc "RMCPhaseSpace1NExternal" --internal 0 --rmcn1emin ${RMC_N1_EMIN} --BB ${BB} --printpot "no" --verbose false 2>/dev/null | tail -1 | awk '{print $NF}')
+  echo "      • Calculating RMC 1N Internal events (emin=${RMC_N1_EMIN})..."
+  RMC_N1_INTERNAL_EVENTS=$(calculateEvents.py --livetime ${LIVETIME} --prc "RMCPhaseSpace1NInternal" --internal 1 --rmcn1emin ${RMC_N1_EMIN} --BB ${BB} --printpot "no" --verbose false 2>/dev/null | tail -1 | awk '{print $NF}')
+fi
 echo "   ✓ All event yields calculated"
 echo ""
 echo "💾 [4/4] Writing configuration file..."
@@ -149,7 +166,12 @@ echo "   Output: ${TAG}.txt"
   echo "DEM_emin=\"${DEM_EMIN}\""
   echo "RPC_TMIN=\"${TMIN}\""
   echo "RPC_emin=\"${RPC_EMIN}\""
-  echo "RMC_emin=\"${RMC_EMIN}\""
+  if [[ ${INCLUDE_RMCN0} -eq 1 ]]; then
+    echo "RMC_N0_emin=\"${RMC_N0_EMIN}\""
+  fi
+  if [[ ${INCLUDE_RMCN1} -eq 1 ]]; then
+    echo "RMC_N1_emin=\"${RMC_N1_EMIN}\""
+  fi
   echo "RMC_kmax=\"${RMC_kmax}\""
   echo "IPA_emin=\"${IPA_EMIN}\"" 
   echo ""
@@ -165,8 +187,14 @@ echo "   Output: ${TAG}.txt"
   echo "ipa_events=\"${IPA_EVENTS}\""
   echo "rpc_internal_events=\"${RPC_INTERNAL_EVENTS}\""
   echo "rpc_external_events=\"${RPC_EXTERNAL_EVENTS}\""
-  echo "rmc_internal_events=\"${RMC_INTERNAL_EVENTS}\""
-  echo "rmc_external_events=\"${RMC_EXTERNAL_EVENTS}\""
+  if [[ ${INCLUDE_RMCN0} -eq 1 ]]; then
+    echo "rmc_n0_internal_events=\"${RMC_N0_INTERNAL_EVENTS}\""
+    echo "rmc_n0_external_events=\"${RMC_N0_EXTERNAL_EVENTS}\""
+  fi
+  if [[ ${INCLUDE_RMCN1} -eq 1 ]]; then
+    echo "rmc_n1_internal_events=\"${RMC_N1_INTERNAL_EVENTS}\""
+    echo "rmc_n1_external_events=\"${RMC_N1_EXTERNAL_EVENTS}\""
+  fi
 } > ${TAG}.txt
 
 echo "   ✓ Configuration file written successfully"
diff --git a/JobConfig/ensemble/scripts/Stage2_build_sampler.sh b/JobConfig/ensemble/scripts/Stage2_build_sampler.sh
index 2ba40068..ee9c50de 100755
--- a/JobConfig/ensemble/scripts/Stage2_build_sampler.sh
+++ b/JobConfig/ensemble/scripts/Stage2_build_sampler.sh
@@ -11,12 +11,13 @@ exit_abnormal() {
 }
 OWNER="mu2e"
 RELEASE=MDC2025
-CURRENT="af"
+CURRENT="au"
 TAG=""
 VERBOSE=1
 
 DIOVERSION=af
-RMCVERSION=af
+RMCVERSIONINT=au
+RMCVERSIONEXT=at
 RPCVERSION=af
 IPAVERSION=af
 
@@ -64,8 +65,11 @@ while getopts ":-:" options; do
         dioversion)
           DIOVERSION=${!OPTIND} OPTIND=$(( $OPTIND + 1 ))
           ;;
-        rmcversion)
-          RMCVERSION=${!OPTIND} OPTIND=$(( $OPTIND + 1 ))
+        rmcversionext)
+          RMCVERSIONEXT=${!OPTIND} OPTIND=$(( $OPTIND + 1 ))
+          ;;
+        rmcversionint)
+          RMCVERSIONINT=${!OPTIND} OPTIND=$(( $OPTIND + 1 ))
           ;;
         rpcversion)
           RPCVERSION=${!OPTIND} OPTIND=$(( $OPTIND + 1 ))
@@ -109,8 +113,11 @@ fi
 if [[ ! -z ${RPC_emin} ]]; then
   RPC_EMIN=${RPC_emin}
 fi
-if [[ ! -z ${RMC_emin} ]]; then
-  RMC_EMIN=${RMC_emin}
+if [[ ! -z ${RMC_N0_emin} ]]; then
+  RMC_N0_EMIN=${RMC_N0_emin}
+fi
+if [[ ! -z ${RMC_N1_emin} ]]; then
+  RMC_N1_EMIN=${RMC_N1_emin}
 fi
 if [[ ! -z ${RMC_kmax} ]]; then
   kmax=${RMC_kmax}
@@ -128,15 +135,18 @@ echo "     • CosmicGen: ${CosmicGen}"
 echo "     • CosmicJob: ${CosmicJob}"
 echo "     • DIO_EMIN: ${DEM_emin}"
 echo "     • RPC_EMIN: ${RPC_emin}"
-echo "     • RMC_EMIN: ${RMC_emin}"
+echo "     • RMC_N0_EMIN: ${RMC_N0_emin}"
+echo "     • RMC_N1_EMIN: ${RMC_N1_emin}"
 echo "     • IPA_EMIN: ${IPA_emin}"
 echo ""
 echo "   Event Yields from Config:"
 echo "     • DIO: ${dio_events:-N/A}"
 echo "     • RPC Internal: ${rpc_internal_events:-N/A}"
 echo "     • RPC External: ${rpc_external_events:-N/A}"
-echo "     • RMC Internal: ${rmc_internal_events:-N/A}"
-echo "     • RMC External: ${rmc_external_events:-N/A}"
+echo "     • RMC 0N Internal: ${rmc_n0_internal_events:-N/A}"
+echo "     • RMC 0N External: ${rmc_n0_external_events:-N/A}"
+echo "     • RMC 1N Internal: ${rmc_n1_internal_events:-N/A}"
+echo "     • RMC 1N External: ${rmc_n1_external_events:-N/A}"
 echo "     • IPA Michel: ${ipa_events:-N/A}"
 echo ""
 
@@ -162,13 +172,21 @@ VALIDATION_FAILED=0
 # Define datasets to check with their corresponding yield variables
 declare -a DATASETS=(
   "dts.mu2e.DIOtail${DIO_EMIN}.${RELEASE}${DIOVERSION}.art:DIO:dio_events"
-  "dts.mu2e.RMCInternal.${RELEASE}${RMCVERSION}.art:RMCInternal:rmc_internal_events"
-  "dts.mu2e.RMCExternal.${RELEASE}${RMCVERSION}.art:RMCExternal:rmc_external_events"
   "dts.mu2e.RPCInternalPhysical.${RELEASE}${RPCVERSION}.art:RPCInternal:rpc_internal_events"
   "dts.mu2e.RPCExternalPhysical.${RELEASE}${RPCVERSION}.art:RPCExternal:rpc_external_events"
   "dts.mu2e.IPAMuminusMichel.${RELEASE}${IPAVERSION}.art:IPAMichel:ipa_events"
 )
 
+# Add conditional RMCPhaseSpace checks based on config
+if [[ ! -z ${RMC_N0_emin} ]]; then
+  DATASETS+=("dts.mu2e.RMCPhaseSpace0NInternal.${RELEASE}${RMCVERSIONINT}.art:RMCPhaseSpace0NInternal:rmc_n0_internal_events")
+  DATASETS+=("dts.mu2e.RMCPhaseSpace0NExternal.${RELEASE}${RMCVERSIONEXT}.art:RMCPhaseSpace0NExternal:rmc_n0_external_events")
+fi
+if [[ ! -z ${RMC_N1_emin} ]]; then
+  DATASETS+=("dts.mu2e.RMCPhaseSpace1NInternal.${RELEASE}${RMCVERSIONINT}.art:RMCPhaseSpace1NInternal:rmc_n1_internal_events")
+  DATASETS+=("dts.mu2e.RMCPhaseSpace1NExternal.${RELEASE}${RMCVERSIONEXT}.art:RMCPhaseSpace1NExternal:rmc_n1_external_events")
+fi
+
 # Also check CRYCosmic for files only (no yield check needed)
 echo "   Checking CRYCosmic files (${NJOBS} files needed)..."
 cry_file_count=$(mu2eDatasetFileList "dts.mu2e.CosmicSignal.${COSMICTAG}.art" 2>/dev/null | wc -l)
@@ -193,7 +211,7 @@ for dataset_pair in "${DATASETS[@]}"; do
   IFS=':' read -r dataset_name dataset_label yield_var <<< "$dataset_pair"
   
   # Check file count
-  file_count=$(mu2eDatasetFileList "$dataset_name" 2>/dev/null | wc -l)
+  file_count=$(mu2eDatasetFileList "$dataset_name" --disk 2>/dev/null | wc -l)
   file_count=$((file_count + 0))  # Ensure it's a number
   if [[ -z "$file_count" ]] || [[ $file_count -lt $NJOBS ]]; then
     echo "   ❌ ${dataset_label}: Only ${file_count} files available (need ${NJOBS})"
@@ -243,19 +261,60 @@ fi
 echo "   ✓ All datasets validated"
 echo ""
 
+# Function: Check if filename lists are empty
+check_file_lists() {
+  local filelist_check_failed=0
+  for file in "$@"; do
+    if [[ ! -f "$file" ]] || [[ ! -s "$file" ]]; then
+      echo "   ⚠️  WARNING: File list is empty or missing: $file"
+      filelist_check_failed=1
+    fi
+  done
+  return $filelist_check_failed
+}
+
 echo "🔨 [3/6] Building file lists (${NJOBS} files per process)..."
 mu2eDatasetFileList "dts.mu2e.CosmicSignal.${COSMICTAG}.art" | head -${NJOBS} > filenames_CRYCosmic
 mu2eDatasetFileList "dts.mu2e.DIOtail${DIO_EMIN}.${RELEASE}${DIOVERSION}.art"| head -${NJOBS} > filenames_DIO
-mu2eDatasetFileList "dts.mu2e.RMCInternal.${RELEASE}${RMCVERSION}.art" | head -${NJOBS} > filenames_RMCInternal
-mu2eDatasetFileList "dts.mu2e.RMCExternal.${RELEASE}${RMCVERSION}.art" | head -${NJOBS} > filenames_RMCExternal
-mu2eDatasetFileList "dts.mu2e.RPCInternalPhysical.${RELEASE}${RPCVERSION}.art" | head -${NJOBS} > filenames_RPCInternal
+if [[ ! -z ${RMC_N0_emin} ]]; then
+  mu2eDatasetFileList "dts.mu2e.RMCPhaseSpace0NInternal.${RELEASE}${RMCVERSIONINT}.art" | head -${NJOBS} > filenames_RMCPhaseSpace0NInternal
+  mu2eDatasetFileList "dts.mu2e.RMCPhaseSpace0NExternal.${RELEASE}${RMCVERSIONEXT}.art" | head -${NJOBS} > filenames_RMCPhaseSpace0NExternal
+fi
+if [[ ! -z ${RMC_N1_emin} ]]; then
+  mu2eDatasetFileList "dts.mu2e.RMCPhaseSpace1NInternal.${RELEASE}${RMCVERSIONINT}.art" | head -${NJOBS} > filenames_RMCPhaseSpace1NInternal
+  mu2eDatasetFileList "dts.mu2e.RMCPhaseSpace1NExternal.${RELEASE}${RMCVERSIONEXT}.art" | head -${NJOBS} > filenames_RMCPhaseSpace1NExternal
+fi
+mu2eDatasetFileList "dts.mu2e.RPCInternalPhysical.${RELEASE}${RPCVERSION}.art" --tape | head -${NJOBS} > filenames_RPCInternal
 mu2eDatasetFileList "dts.mu2e.RPCExternalPhysical.${RELEASE}${RPCVERSION}.art" | head -${NJOBS} > filenames_RPCExternal
 mu2eDatasetFileList "dts.mu2e.IPAMuminusMichel.${RELEASE}${IPAVERSION}.art" | head -${NJOBS} > filenames_IPAMichel
 echo "   ✓ File lists created"
 echo ""
 
+echo "🔍 Validating file list contents..."
+FILES_TO_CHECK="filenames_CRYCosmic filenames_DIO filenames_RPCInternal filenames_RPCExternal filenames_IPAMichel"
+if [[ ! -z ${RMC_N0_emin} ]]; then
+  FILES_TO_CHECK="${FILES_TO_CHECK} filenames_RMCPhaseSpace0NInternal filenames_RMCPhaseSpace0NExternal"
+fi
+if [[ ! -z ${RMC_N1_emin} ]]; then
+  FILES_TO_CHECK="${FILES_TO_CHECK} filenames_RMCPhaseSpace1NInternal filenames_RMCPhaseSpace1NExternal"
+fi
+if check_file_lists ${FILES_TO_CHECK}; then
+  echo "   ✓ All file lists populated"
+else
+  echo "   ⚠️  WARNING: Some file lists may be empty or incomplete"
+fi
+echo ""
+
 echo "📝 [4/6] Generating template FCL files..."
-make_template_fcl.py --BB=${BB} --release=${RELEASE}${CURRENT}  --tag=${TAG} --verbose=${VERBOSE} --livetime=${LIVETIME} --run=${RUN} --dioemin=${DIO_EMIN} --rpcemin=${RPC_EMIN} --rmcemin=${RMC_EMIN} --rmckmax=${RMC_kmax} --ipaemin=${IPA_EMIN} --tmin=${TMIN} --samplingseed=${SAMPLINGSEED} --prc "DIO" "CRYCosmic" "RPCInternal" "RPCExternal" "RMCInternal" "RMCExternal" "IPAMichel"
+# Build process list dynamically
+PRC_LIST="DIO CRYCosmic RPCInternal RPCExternal IPAMichel"
+if [[ ! -z ${RMC_N0_emin} ]]; then
+  PRC_LIST="${PRC_LIST} RMCPhaseSpace0NInternal RMCPhaseSpace0NExternal"
+fi
+if [[ ! -z ${RMC_N1_emin} ]]; then
+  PRC_LIST="${PRC_LIST} RMCPhaseSpace1NInternal RMCPhaseSpace1NExternal"
+fi
+make_template_fcl.py --BB=${BB} --release=${RELEASE}${CURRENT}  --tag=${TAG} --verbose=${VERBOSE} --livetime=${LIVETIME} --run=${RUN} --dioemin=${DIO_EMIN} --rpcemin=${RPC_EMIN} --rmcn0emin=${RMC_N0_EMIN} --rmcn1emin=${RMC_N1_EMIN} --rmckmax=${RMC_kmax} --ipaemin=${IPA_EMIN} --tmin=${TMIN} --samplingseed=${SAMPLINGSEED} --prc ${PRC_LIST}
 echo "   ✓ Template FCL files generated"
 echo ""
 
@@ -266,26 +325,62 @@ rm -f filenames_CRYCosmic_${NJOBS}.txt
 rm -f filenames_DIO_${NJOBS}.txt
 rm -f filenames_RPCInternal_${NJOBS}.txt
 rm -f filenames_RPCExternal_${NJOBS}.txt
-rm -f filenames_RMCInternal_${NJOBS}.txt
-rm -f filenames_RMCExternal_${NJOBS}.txt
+if [[ ! -z ${RMC_N0_emin} ]]; then
+  rm -f filenames_RMCPhaseSpace0NInternal_${NJOBS}.txt
+  rm -f filenames_RMCPhaseSpace0NExternal_${NJOBS}.txt
+fi
+if [[ ! -z ${RMC_N1_emin} ]]; then
+  rm -f filenames_RMCPhaseSpace1NInternal_${NJOBS}.txt
+  rm -f filenames_RMCPhaseSpace1NExternal_${NJOBS}.txt
+fi
 rm -f filenames_IPAMichel_${NJOBS}.txt
 
 echo "   Creating SAM file lists..."
 samweb list-files "dh.dataset=dts.mu2e.CosmicSignal.${COSMICTAG}.art" | head -${NJOBS} > filenames_CRYCosmic_${NJOBS}.txt
 samweb list-files "dh.dataset=dts.mu2e.DIOtail${DIO_EMIN}.${RELEASE}${DIOVERSION}.art"  | head -${NJOBS} > filenames_DIO_${NJOBS}.txt
-samweb list-files "dh.dataset=dts.mu2e.RMCInternal.${RELEASE}${RMCVERSION}.art  and availability:anylocation"  | head -${NJOBS}  >  filenames_RMCInternal_${NJOBS}.txt
-samweb list-files "dh.dataset=dts.mu2e.RMCExternal.${RELEASE}${RMCVERSION}.art  and availability:anylocation"  | head -${NJOBS}  >  filenames_RMCExternal_${NJOBS}.txt
+if [[ ! -z ${RMC_N0_emin} ]]; then
+  samweb list-files "dh.dataset=dts.mu2e.RMCPhaseSpace0NInternal.${RELEASE}${RMCVERSIONINT}.art  and availability:anylocation"  | head -${NJOBS}  >  filenames_RMCPhaseSpace0NInternal_${NJOBS}.txt
+  samweb list-files "dh.dataset=dts.mu2e.RMCPhaseSpace0NExternal.${RELEASE}${RMCVERSIONEXT}.art  and availability:anylocation"  | head -${NJOBS}  >  filenames_RMCPhaseSpace0NExternal_${NJOBS}.txt
+fi
+if [[ ! -z ${RMC_N1_emin} ]]; then
+  samweb list-files "dh.dataset=dts.mu2e.RMCPhaseSpace1NInternal.${RELEASE}${RMCVERSIONINT}.art  and availability:anylocation"  | head -${NJOBS}  >  filenames_RMCPhaseSpace1NInternal_${NJOBS}.txt
+  samweb list-files "dh.dataset=dts.mu2e.RMCPhaseSpace1NExternal.${RELEASE}${RMCVERSIONEXT}.art  and availability:anylocation"  | head -${NJOBS}  >  filenames_RMCPhaseSpace1NExternal_${NJOBS}.txt
+fi
 samweb list-files "dh.dataset=dts.mu2e.RPCInternalPhysical.${RELEASE}${RPCVERSION}.art  and availability:anylocation"  | head -${NJOBS}  >  filenames_RPCInternal_${NJOBS}.txt
 samweb list-files "dh.dataset=dts.mu2e.RPCExternalPhysical.${RELEASE}${RPCVERSION}.art  and availability:anylocation"  | head -${NJOBS}  >  filenames_RPCExternal_${NJOBS}.txt
 samweb list-files "dh.dataset=dts.mu2e.IPAMuminusMichel.${RELEASE}${IPAVERSION}.art  and availability:anylocation"  | head -${NJOBS}  >  filenames_IPAMichel_${NJOBS}.txt
 echo "   ✓ SAM file lists ready"
 echo ""
 
+echo "🔍 Validating SAM file list contents..."
+SAM_FILES_TO_CHECK="filenames_CRYCosmic_${NJOBS}.txt filenames_DIO_${NJOBS}.txt filenames_RPCInternal_${NJOBS}.txt filenames_RPCExternal_${NJOBS}.txt filenames_IPAMichel_${NJOBS}.txt"
+if [[ ! -z ${RMC_N0_emin} ]]; then
+  SAM_FILES_TO_CHECK="${SAM_FILES_TO_CHECK} filenames_RMCPhaseSpace0NInternal_${NJOBS}.txt filenames_RMCPhaseSpace0NExternal_${NJOBS}.txt"
+fi
+if [[ ! -z ${RMC_N1_emin} ]]; then
+  SAM_FILES_TO_CHECK="${SAM_FILES_TO_CHECK} filenames_RMCPhaseSpace1NInternal_${NJOBS}.txt filenames_RMCPhaseSpace1NExternal_${NJOBS}.txt"
+fi
+if check_file_lists ${SAM_FILES_TO_CHECK}; then
+  echo "   ✓ All SAM file lists populated"
+else
+  echo "   ⚠️  WARNING: Some SAM file lists may be empty or incomplete"
+fi
+echo ""
+
 echo "🚀 [6/6] Submitting ensemble jobs..."
 DSCONF=${RELEASE}${CURRENT}
 # note change setup to code to use a custom tarball
 echo "   Running mu2ejobdef..."
-cmd="mu2ejobdef --desc=ensemble${TAG} --dsconf=${DSCONF} --run=${RUN} --setup ${SETUP} --sampling=1:DIO:filenames_DIO_${NJOBS}.txt --sampling=1:CRYCosmic:filenames_CRYCosmic_${NJOBS}.txt --sampling=1:RPCInternal:filenames_RPCInternal_${NJOBS}.txt  --embed SamplingInput_sr0.fcl  --sampling=1:RPCExternal:filenames_RPCExternal_${NJOBS}.txt --sampling=1:RMCInternal:filenames_RMCInternal_${NJOBS}.txt --sampling=1:RMCExternal:filenames_RMCExternal_${NJOBS}.txt --sampling=1:IPAMichel:filenames_IPAMichel_${NJOBS}.txt --verb "
+# Build sampling options dynamically
+SAMPLING_OPTS="--sampling=1:DIO:filenames_DIO_${NJOBS}.txt --sampling=1:CRYCosmic:filenames_CRYCosmic_${NJOBS}.txt --sampling=1:RPCInternal:filenames_RPCInternal_${NJOBS}.txt --sampling=1:RPCExternal:filenames_RPCExternal_${NJOBS}.txt"
+if [[ ! -z ${RMC_N0_emin} ]]; then
+  SAMPLING_OPTS="${SAMPLING_OPTS} --sampling=1:RMCPhaseSpace0NInternal:filenames_RMCPhaseSpace0NInternal_${NJOBS}.txt --sampling=1:RMCPhaseSpace0NExternal:filenames_RMCPhaseSpace0NExternal_${NJOBS}.txt"
+fi
+if [[ ! -z ${RMC_N1_emin} ]]; then
+  SAMPLING_OPTS="${SAMPLING_OPTS} --sampling=1:RMCPhaseSpace1NInternal:filenames_RMCPhaseSpace1NInternal_${NJOBS}.txt --sampling=1:RMCPhaseSpace1NExternal:filenames_RMCPhaseSpace1NExternal_${NJOBS}.txt"
+fi
+SAMPLING_OPTS="${SAMPLING_OPTS} --sampling=1:IPAMichel:filenames_IPAMichel_${NJOBS}.txt"
+cmd="mu2ejobdef --desc=ensemble${TAG} --dsconf=${DSCONF} --run=${RUN} --setup ${SETUP} ${SAMPLING_OPTS} --embed SamplingInput_sr0.fcl --verb "
 
 $cmd
 parfile=$(ls cnf.*.tar)
diff --git a/JobConfig/ensemble/scripts/Stage3_configure_ensemble_campaign.sh b/JobConfig/ensemble/scripts/Stage3_configure_ensemble_campaign.sh
index d918f439..b9925fd0 100755
--- a/JobConfig/ensemble/scripts/Stage3_configure_ensemble_campaign.sh
+++ b/JobConfig/ensemble/scripts/Stage3_configure_ensemble_campaign.sh
@@ -16,7 +16,7 @@ exit_abnormal() {
 # Default values
 TAG=""
 RELEASE="MDC2025"
-VERSION="ai"
+VERSION="au"
 OWNER="mu2e"
 APPEND=0
 
@@ -124,7 +124,7 @@ if [[ ${APPEND} -eq 1 ]]; then
     },
     "inloc": "tape",
     "outloc": {"*.root": "disk"},
-    "simjob_setup": "/cvmfs/mu2e.opensciencegrid.org/Musings/AnalysisMDC2025/v01_01_03/setup.sh"
+    "simjob_setup": "/cvmfs/mu2e.opensciencegrid.org/Musings/AnalysisMDC2025/v02_00_00/setup.sh"
   }
 ]
 CAMPAIGN_EOF
@@ -228,7 +228,7 @@ echo "   Tape location: ${TAPE_PNFS}"
 echo ""
 
 echo "   Adding file location to SAM..."
-samweb add-file-location cnf.${OWNER}.ensemble${TAG}.${RELEASE}${VERSION}.tar ${TAPE_PATH}
+samweb add-file-location cnf.${OWNER}.ensemble${TAG}.${RELEASE}${VERSION}.0.tar ${TAPE_PATH}
 if [[ $? -ne 0 ]]; then
   echo "   ✗ Error: Failed to add file location"
   exit 1
@@ -236,6 +236,71 @@ fi
 echo "   ✓ File location registered"
 echo ""
 
+echo "📋 [3b/4] Processing tag information file..."
+TAG_INFO_FILE="cnf.mu2e.${TAG}-info.${RELEASE}${VERSION}.0.txt"
+TAG_INFO_JSON_FILE="${TAG_INFO_FILE}.json"
+
+echo "   Renaming config file to tag-info..."
+mv ${CONFIG_FILE} ${TAG_INFO_FILE}
+if [[ ! -f ${TAG_INFO_FILE} ]]; then
+  echo "   ✗ Error: Failed to rename config file to ${TAG_INFO_FILE}"
+  exit 1
+fi
+echo "   ✓ Renamed to ${TAG_INFO_FILE}"
+echo ""
+
+echo "   Generating tag-info metadata JSON..."
+printJson --no-parents ${TAG_INFO_FILE} > ${TAG_INFO_JSON_FILE}
+
+if [[ ! -f "${TAG_INFO_JSON_FILE}" ]] || [[ ! -s "${TAG_INFO_JSON_FILE}" ]]; then
+  echo "   ✗ Error: Failed to generate tag-info JSON file"
+  exit 1
+fi
+echo "   ✓ Generated ${TAG_INFO_JSON_FILE}"
+echo ""
+
+echo "   Declaring tag-info JSON file to SAM..."
+ls ${TAG_INFO_JSON_FILE} | mu2eFileDeclare
+if [[ $? -ne 0 ]]; then
+  echo "   ✗ Error: Failed to declare tag-info file"
+  exit 1
+fi
+echo "   ✓ Tag-info file declared"
+echo ""
+
+echo "   Uploading tag-info file to tape..."
+UPLOAD_TAG_OUTPUT=$(ls ${TAG_INFO_FILE} | mu2eFileUpload --disk 2>&1)
+if [[ ${PIPESTATUS[1]} -ne 0 ]]; then
+  echo "   ✗ Error: Failed to upload tag-info file"
+  exit 1
+fi
+
+# Extract tape path from upload output
+TAG_TAPE_FILE=$(echo "${UPLOAD_TAG_OUTPUT}" | sed -n 's/.*to \(\/pnfs[^ ]*\).*/\1/p')
+
+if [[ -z ${TAG_TAPE_FILE} ]]; then
+  echo "   ✗ Error: Could not extract tape path from upload output"
+  echo "   Output was: ${UPLOAD_TAG_OUTPUT}"
+  exit 1
+fi
+
+# Get directory path (remove filename)
+TAG_TAPE_PNFS=$(dirname "${TAG_TAPE_FILE}")
+TAG_TAPE_PATH="enstore:${TAG_TAPE_PNFS}"
+
+echo "   ✓ Tag-info file uploaded"
+echo "   Tape location: ${TAG_TAPE_PNFS}"
+echo ""
+
+echo "   Adding tag-info file location to SAM..."
+samweb add-file-location ${TAG_INFO_FILE} ${TAG_TAPE_PATH}
+if [[ $? -ne 0 ]]; then
+  echo "   ✗ Error: Failed to add tag-info file location"
+  exit 1
+fi
+echo "   ✓ Tag-info file location registered"
+echo ""
+
 echo "📋 [4/4] Generating campaign JSON (multipart stages)..."
 # Create JSON for the campaign with all stages
 cat > ${CAMPAIGN_JSON_FILE} << EOF
@@ -279,7 +344,12 @@ echo "════════════════════════
 echo "✅ Stage 3 Complete!"
 echo "   TAR file: cnf.${OWNER}.ensemble${TAG}.${RELEASE}${VERSION}.tar"
 echo "   TAR metadata: cnf.${OWNER}.ensemble${TAG}.${RELEASE}${VERSION}.tar.json (declared and uploaded)"
-echo "   Tape path: ${TAPE_PATH}"
+echo "   Tag-info file: cnf.mu2e.tag-info.0.txt (declared and uploaded)"
+echo "   Tag-info metadata: cnf.mu2e.tag-info.0.txt.json (declared and uploaded)"
+echo "   Tape path (TAR): ${TAPE_PATH}"
+echo "   Tape path (Tag-info): ${TAG_TAPE_PATH}"
 echo "   Campaign: ${CAMPAIGN_JSON_FILE} (pending declaration)"
+echo "   For ensemble generation enter mu2epro and launch: e.g. mkidxdef --jobdefs /exp/mu2e/app/users/mu2epro/production_manager/poms_map/MDC2025-MDS3b.json --prod"
+echo "   For digi/mix/reco/ntuple enter mu2epro and launch: e.g. json2jobdef --json digi.json --dsconf MDC2025af_best_v1_3 --desc ensembleMDS3aOnSpill --jobdefs /exp/mu2e/app/users/mu2epro/production_manager/poms_map/MDC2025-002.json --prod"
 echo "═══════════════════════════════════════════════════════════════"
 echo ""
diff --git a/JobConfig/ensemble/scripts/Stage5_signal.sh b/JobConfig/ensemble/scripts/Stage5_signal.sh
index 081f6354..e9cd641e 100755
--- a/JobConfig/ensemble/scripts/Stage5_signal.sh
+++ b/JobConfig/ensemble/scripts/Stage5_signal.sh
@@ -138,7 +138,10 @@ if [[ ${#PARTS[@]} -ge 3 ]]; then
   RELEASE="${PARTS[0]}"
   DBPURPOSE="${PARTS[1]}"
   # Everything after PURPOSE is VERSION (in case it has multiple underscores like v1_3)
-  DBVERSION=$(IFS='_'; echo "${PARTS[@]:2}")
+  DBVERSION="${PARTS[2]}"
+  for ((j=3; j<${#PARTS[@]}; j++)); do
+    DBVERSION+="_${PARTS[$j]}"
+  done
 else
   echo "ERROR: Could not parse RELEASE, DBPURPOSE, and DBVERSION from dataset name"
   echo "Expected format: mcs.mu2e.NAME.RELEASE_PURPOSE_VERSION.art"