rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

ML script

Browse Source
This commit is contained in:
James Szalkie 2026-08-13 12:41:44 -04:00
parent bdc8949d38
commit 7181b80635
1 changed files with 43 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
Armory/ANASEN_ML.py
View File

@ -15,13 +15,16 @@ from sklearn.metrics import mean_absolute_error
# User settings # User settings
ROOT_FILE = "/Users/jamesszalkie/ANASEN_analysis/Armory/SimAnasen1.root" ROOT_FILE = "/Users/jamesszalkie/ANASEN_analysis/Armory/Ne18Protons.root"
TREE_NAME = "tree1" TREE_NAME = "tree1"
MODEL = "beam_predictor.keras" MODEL = "beam_predictor.keras"
INPUT_SCALER = "input_scaler.pkl" INPUT_SCALER = "input_scaler.pkl"
OUTPUT_SCALER = "output_scaler.pkl" OUTPUT_SCALER = "output_scaler.pkl"
# Candidate excitation energies (MeV) for nearest-state snapping
EX_CANDIDATES = np.array([0, 0.3, 1.7, 2.4, 2.8], dtype=np.float32)
INPUT_BRANCHES = [ INPUT_BRANCHES = [
"Tb", "Tb",
"thetab", "thetab",
@ -30,6 +33,10 @@ INPUT_BRANCHES = [
"MTarget", "MTarget",
"MLight", "MLight",
"MHeavy", "MHeavy",
"ZBeam",
"ZHeavy",
"ZLight",
"ZTarget"
] ]
# Load model # Load model
@ -64,6 +71,10 @@ pred = output_scaler.inverse_transform(pred_scaled)
beam = pred[:,0] beam = pred[:,0]
Ex = pred[:,1] Ex = pred[:,1]
# Snap each predicted excitation to the nearest candidate energy.
nearest_idx = np.argmin(np.abs(Ex[:, None] - EX_CANDIDATES[None, :]), axis=1)
Ex_snapped = EX_CANDIDATES[nearest_idx]
# See whether truth branches exist # See whether truth branches exist
truth_beam = None truth_beam = None
@ -95,8 +106,15 @@ if truth_Ex is not None:
Ex_mae = mean_absolute_error(truth_Ex, Ex) Ex_mae = mean_absolute_error(truth_Ex, Ex)
Ex_rmse = np.sqrt(np.mean((truth_Ex - Ex)**2)) Ex_rmse = np.sqrt(np.mean((truth_Ex - Ex)**2))
Ex_snap_mae = mean_absolute_error(truth_Ex, Ex_snapped)
Ex_snap_rmse = np.sqrt(np.mean((truth_Ex - Ex_snapped)**2))
print(f"Excitation MAE : {Ex_mae:.4f} MeV") print(f"Excitation MAE : {Ex_mae:.4f} MeV")
print(f"Excitation RMSE : {Ex_rmse:.4f} MeV") print(f"Excitation RMSE : {Ex_rmse:.4f} MeV")
print(f"Ex Snapped MAE : {Ex_snap_mae:.4f} MeV")
print(f"Ex Snapped RMSE : {Ex_snap_rmse:.4f} MeV")
print("\nEx candidate states (MeV):", ", ".join(f"{x:.3f}" for x in EX_CANDIDATES))
# Plot Beam Energy # Plot Beam Energy
plt.figure(figsize=(8,6)) plt.figure(figsize=(8,6))
@ -133,7 +151,15 @@ plt.hist(
bins=250, bins=250,
histtype="step", histtype="step",
linewidth=2, linewidth=2,
label="Predicted", label="Predicted (raw)",
)
plt.hist(
Ex_snapped,
bins=250,
histtype="step",
linewidth=2,
label="Predicted (snapped)",
) )
if truth_Ex is not None: if truth_Ex is not None:
@ -175,4 +201,18 @@ plt.plot([mn, mx], [mn, mx], 'k--')
plt.xlabel("True Excitation Energy (MeV)") plt.xlabel("True Excitation Energy (MeV)")
plt.ylabel("Predicted Excitation Energy (MeV)") plt.ylabel("Predicted Excitation Energy (MeV)")
plt.title("Excitation Energy Reconstruction") plt.title("Excitation Energy Reconstruction")
if truth_Ex is not None:
plt.figure(figsize=(6,6))
plt.scatter(truth_Ex, Ex_snapped, s=2)
mn = min(truth_Ex.min(), Ex_snapped.min())
mx = max(truth_Ex.max(), Ex_snapped.max())
plt.plot([mn, mx], [mn, mx], 'k--')
plt.xlabel("True Excitation Energy (MeV)")
plt.ylabel("Snapped Excitation Energy (MeV)")
plt.title("Excitation Energy Reconstruction (Snapped)")