Updates for ensembling

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_*_*_*",

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)

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")

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,110 @@ 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)
+
+# Compute RMC spectrum fractions from Plestid integral using correct K_max values
+compute_rmc_spectrum_fractions()
 
 #-------------------------------------------------------------------------------------#    
 def get_duty_factor(run_mode='1BB'):
@@ -500,7 +630,126 @@ 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, 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 to use based on energy threshold
+    # Default: fraction from E > 57 MeV
+    e_threshold = 57.0
+    R_spectrum = RMC_SPECTRUM_FRAC_0N_57
+
+    # If threshold is higher, use interpolated fraction
+    # For E > 80 MeV: use the E > 80 spectrum fraction
+    if float(e_min) >80.0:
+        R_spectrum = RMC_SPECTRUM_FRAC_0N_80 / RMC_SPECTRUM_FRAC_0N_57
+
+    # 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, 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 to use based on energy threshold
+    # Default: fraction from E > 57 MeV
+    e_threshold = 57.0
+    R_spectrum = RMC_SPECTRUM_FRAC_1N_57
+
+    # If threshold is higher, use interpolated fraction
+    # For E > 80 MeV: use the E > 80 spectrum fraction
+    if float(e_min) >80.0:
+        R_spectrum = RMC_SPECTRUM_FRAC_1N_80 / RMC_SPECTRUM_FRAC_1N_57
+
+    # 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 +838,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)