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/constants.py

@@ -55,3 +55,29 @@
 # are not a factor of two apart in protons per spill.
 ONEBB1500W_POT_PER_CYCLE = 1.556e12 # protons per 1.33 s cycle at 1500 W
 ONEBB1500W_PROTONS_PER_SPILL = 6.14e6 # protons per 1695 ns spill at 1500 W
+
+# --- RMC (Radiative Muon Capture) Physics Constants ---
+# Branching ratio for muon capture on nucleus
+RMC_BR_MUON_CAPTURE = 0.609
+# RMC rate above 57 MeV, relative to OMC
+RMC_RATE_GT_57 = 1.41e-5
+# Branching ratio fraction for 0-nucleon knockout above 57 MeV
+RMC_BR_0N_FRAC_GT_57 = 0.099  # BR(0 knockout | E > 57) / BR(RMC | E > 57)
+# Branching ratio fraction for 1-nucleon knockout above 57 MeV
+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
+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)
+RMC_INTERNAL_EXTERNAL_RATIO = 0.0069  # rho = BR(internal) / BR(external)

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