more plotting troubleshooting

Armory/anasenMS.cpp

@@ -100,13 +100,13 @@ int main(int argc, char **argv){
   // then sequential decay of 21Na* -> p + 20Ne
   TransferReaction transfer;
 
-  transfer.SetA(27, 13, 0);            // 18Ne projectile
+  transfer.SetA(18, 10, 0);            // 18Ne projectile
-  transfer.SetIncidentEnergyAngle((72/27.0), 0, 0); // KEA in MeV/u, theta and phi in degree
+  transfer.SetIncidentEnergyAngle((4.4), 0, 0); // KEA in MeV/u, theta and phi in degree
   transfer.Seta(4, 2);                 // 4He target
-  transfer.Setb(4, 2);                 // outgoing proton from the primary transfer
+  transfer.Setb(1, 1);                 // outgoing proton from the primary transfer
-  transfer.SetB(27, 13);               // 21Na* heavy product
+  transfer.SetB(21, 11);               // 21Na* heavy product
 
-  bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable
+  bool enableSequentialDecay = true; // turning to false to disable sequential decay for now, can be set to true to enable
   const int decayDaughterA = 20;
   const int decayDaughterZ = 10;
   const int decayEjectA = 1;
@@ -239,7 +239,7 @@ int main(int argc, char **argv){
 
   // reconstructed SX3 hit position
   double sx3X, sx3Y, sx3Z;
-  tree1->Branch("sx3X",     &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing)
+  tree1->Branch("sx3X",     &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing) in mm
   tree1->Branch("sx3Y",     &sx3Y, "sx3Y/D"); // reconstructed Y position from SX3 (with optional smearing)
   tree1->Branch("sx3Z",     &sx3Z, "sx3Z/D"); // reconstructed Z position from SX3 (with optional smearing)
 

ELoss/PCEnergyAnalysis.py

@@ -11,30 +11,54 @@ from scipy.interpolate import interp1d
 import uproot
 import pycatima as catima
 from scipy.integrate import cumulative_trapezoid
-import matplotlib
+#matplotlib.use("Agg")
-matplotlib.use("Agg")
 import matplotlib.pyplot as plt
 import cmd
 import shlex
 import textwrap
-import readline
-import atexit
 import os
 import periodictable as pt
 import re
-from matplotlib.colors import PowerNorm
+from matplotlib.colors import LinearSegmentedColormap
+import matplotlib as mpl
 
-#Run program from terminal or IDE, and prompts will provide user steps
+# ROOT-like styling
+plt.rcParams.update({
 
-#histfile = os.path.expanduser("~/.uproot_shell_history")
+    # Figure
+    "figure.facecolor": "white",
+    "axes.facecolor": "white",
 
-#try:
+    # ROOT-like axes
-#    readline.read_history_file(histfile)
+    "axes.linewidth": 1.2,
-#except FileNotFoundError:
-#    pass
 
-#atexit.register(readline.write_history_file, histfile)
+    # Ticks on all sides
+    "xtick.top": True,
+    "ytick.right": True,
 
+    # Tick direction inward
+    "xtick.direction": "in",
+    "ytick.direction": "in",
+
+    # Major/minor ticks
+    "xtick.major.size": 10,
+    "ytick.major.size": 10,
+    "xtick.minor.size": 5,
+    "ytick.minor.size": 5,
+
+    # Smaller ROOT-like fonts
+    "font.size": 12,
+})
+
+base_cmap = plt.cm.jet
+
+# Stop before the red region
+colors = base_cmap(np.linspace(0.2, 0.65, 256))
+
+yellow_jet = LinearSegmentedColormap.from_list(
+    'yellow_jet',
+    colors
+)
 #data = [z, mass_u, maximum MeV, name]
 alpha_data = [2, 4.0026, 40, "alpha"]
 proton_data = [1, 1.0078, 20, "proton"]
@@ -172,7 +196,7 @@ def get_interpolators(particle):
         x,
         E,
         bounds_error=False,
-        fill_value=0.0
+        fill_value="extrapolate"
     )
 
     x_of_E = interp1d(
@@ -258,43 +282,67 @@ def update_plot_data(name, values):
             return
     plot_data.append((name, values))
     
-def process_file(filename):
+def calculate_distance_tree1(vx, vy, vz, sx3x, sx3y, sx3z):
+    """Calculate 3D distance from vertex to SX3 hit position for tree1"""
+    dx = (sx3x - vx) / 10
+    dy = (sx3y - vy) / 10
+    dz = (sx3z - vz) / 10
+    return np.sqrt(dx*dx + dy*dy + dz*dz)
 
-    tree = uproot.open(filename)["tree"]
+def calculate_distance_tree2(vz, theta, z_max=34.86):
+    """Calculate distance along trajectory from vertex to z_max for tree2
+    theta is polar angle from z-axis in radians"""
+    cos_theta = np.cos(theta)
+    cos_theta = np.where(np.abs(cos_theta) < 1e-10, 1e-10, cos_theta)
+    return np.abs((z_max - vz) / cos_theta)
 
-    branches = ["Tb", "thetab", "sx3Z"]
 
-    data = tree.arrays(branches, library="np")
+def load_tree_arrays(tree, treename, max_events=None):
+    branches = ["Tb", "thetab", "sx3Z", "vX", "vY", "vZ"]
+    if treename == 'tree1':
+        branches.extend(["sx3X", "sx3Y"])
+    return tree.arrays(branches, library="np", entry_stop=max_events)
+
+
+def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
+    data = load_tree_arrays(tree, treename, max_events)
 
     Ei = data["Tb"]
     theta = np.radians(data["thetab"])
-    sx3Z = data["sx3Z"]
+    vX = data["vX"]
+    vY = data["vY"]
+    vZ = data["vZ"]
 
-    mask = np.sin(theta) != 0
+    if treename == 'tree1':
+        sx3X = data["sx3X"]
+        sx3Y = data["sx3Y"]
+        sx3Z = data["sx3Z"]
+        mask = ~np.isnan(sx3X) & ~np.isnan(sx3Y) & ~np.isnan(sx3Z)
+    else:
+        sx3Z = np.full_like(vZ, z_max) # Assume sx3Z is at z_max for tree2
+        mask = ~np.isnan(Ei) & ~np.isnan(theta) & ~np.isnan(vZ)
 
     Ei = Ei[mask]
     theta = theta[mask]
-    sx3Z = sx3Z[mask]
+    vX = vX[mask]
+    vY = vY[mask]
+    vZ = vZ[mask]
+
+    if treename == 'tree1':
+        sx3X = sx3X[mask]
+        sx3Y = sx3Y[mask]
+        sx3Z = sx3Z[mask]
+        dsx3 = calculate_distance_tree1(vX, vY, vZ, sx3X, sx3Y, sx3Z)
+    else:
+        dsx3 = calculate_distance_tree2(vZ, theta, z_max=z_max)
 
     sin_theta = np.sin(theta)
-
+    sin_theta = np.where(sin_theta != 0, sin_theta, 1e-10)
-    radii = np.array([3.2, 4.2, 6.6])
+    radii = np.array([3.7, 4.2])
-
     dA = radii[0] / sin_theta
     dC = radii[1] / sin_theta
-    dsx3 = radii[2] / sin_theta
 
-    # Determine particle from filename
+    print(f"Computing energies for {particle} ({treename})...")
-    lower = filename.lower()
-
-    if "proton" in lower:
-        particle = "proton"
-    elif "alpha" in lower:
-        particle = "alpha"
-    else:
-        particle = "proton"
-
-    print(f"Computing energies for {particle}...")
 
     EA = energy_loss(particle, Ei, dA)
     EC = energy_loss(particle, Ei, dC)
@@ -314,6 +362,21 @@ def process_file(filename):
         "Eprop": Eprop
     }
 
+
+def process_file(filename, treename, particle=None, max_events=None):
+    tree = uproot.open(filename)[f"{treename}"]
+
+    if particle is None:
+        lower = filename.lower()
+        if "proton" in lower:
+            particle = "proton"
+        elif "alpha" in lower:
+            particle = "alpha"
+        else:
+            particle = "proton"
+    print(f"File {filename} particle {particle}, tree {treename}")
+    return prepare_tree_data(tree, treename, particle, max_events=max_events)
+
 class MyInteractiveApp(cmd.Cmd):
     def __init__(self):
         super().__init__()
@@ -484,7 +547,7 @@ class MyInteractiveApp(cmd.Cmd):
             plt.tight_layout()
             filename = f"Energy_Loss_Curve_{label}.png"
             plt.savefig(filename, dpi=300, bbox_inches="tight")
-            plt.close()
+            plt.show()
             
             print(f"Saved plot: {filename}")
     
@@ -601,7 +664,7 @@ class MyInteractiveApp(cmd.Cmd):
     
         print("Available trees:", keys)
     
-        treeName = arg if len(arg) > 0 else "tree"
+        treeName = arg if len(arg) > 0 else "tree1"
     
         try:
             # uproot automatically resolves ";1"
@@ -696,18 +759,20 @@ class MyInteractiveApp(cmd.Cmd):
             self.do_set_tree("")
             print(f"Using TTree: {self.tree}")
     
-        branches = [
+        treename = self.tree.name if hasattr(self.tree, 'name') else "tree1"
-            "Tb",
+        
-            "thetab",
+
-            "sx3Z"
+        branches = ["Tb", "thetab", "sx3Z", "vX", "vY", "vZ", "sx3X", "sx3Y"]
-        ]
+        Tb_key = "Tb"
+        theta_key = "thetab"
+
     
         if max_events:
             n_events = max_events
         else:
             n_events = self.tree.num_entries
     
-        print(f"Loading {n_events} events...")
+        print(f"Loading {n_events} events from {treename}...")
     
         data = self.tree.arrays(
             branches,
@@ -715,25 +780,42 @@ class MyInteractiveApp(cmd.Cmd):
             entry_stop=max_events
         )
         global Ei
-        Ei = data["Tb"]
+        Ei = data[Tb_key]
         global sx3Z
         sx3Z = data["sx3Z"]
     
-        theta = np.radians(data["thetab"])
+        theta = np.radians(data[theta_key])
-    
+        vX = data["vX"]
-        # Remove theta = 0 events
+        vY = data["vY"]
-        mask = np.sin(theta) != 0
+        vZ = data["vZ"]
+        
+        if treename == 'tree1':
+            sx3X = data["sx3X"]
+            sx3Y = data["sx3Y"]
+            mask = ~np.isnan(sx3X) & ~np.isnan(sx3Y) & ~np.isnan(sx3Z)
+        else:
+            mask = ~np.isnan(Ei) & ~np.isnan(theta)
     
         Ei = Ei[mask]
         theta = theta[mask]
-    
+        sx3Z = sx3Z[mask]
+        vZ_masked = vZ[mask]
+        
         sin_theta = np.sin(theta)
+        sin_theta = np.where(sin_theta != 0, sin_theta, 1e-10)
     
-        radii = np.array([3.2, 4.2, 6.6])
+        if treename == 'tree1':
-    
+            vX_masked = vX[mask]
+            vY_masked = vY[mask]
+            sx3X_masked = sx3X[mask]
+            sx3Y_masked = sx3Y[mask]
+            dsx3 = calculate_distance_tree1(vX_masked, vY_masked, vZ_masked, sx3X_masked, sx3Y_masked, sx3Z)
+        else:
+            dsx3 = calculate_distance_tree2(vZ_masked, theta, z_max=34.86)
+        
+        radii = np.array([3.7, 4.2])
         dA = radii[0] / sin_theta
         dC = radii[1] / sin_theta
-        dsx3 = radii[2] / sin_theta
     
         print("Calculating energy losses...")
         global EA
@@ -830,162 +912,265 @@ class MyInteractiveApp(cmd.Cmd):
         print("Finished writing augmented ROOT file")   
         
     def do_make_plots(self, arg):
-    
+
         import os
-        
+
         global plot_data
-        
         plot_data = []
-    
+
         args = shlex.split(arg)
-    
         particle = args[0] if len(args) > 0 else "proton"
         max_events = int(args[1]) if len(args) > 1 else None
-    
+
         if self.tree is None:
             self.do_set_tree("")
             print(f"Using TTree: {self.tree}")
-    
+
-        branches = ["Tb", "thetab", "sx3Z"]
+        treename = self.tree.name if hasattr(self.tree, 'name') else "tree1"
-    
+
-        data = self.tree.arrays(
+        data = prepare_tree_data(self.tree, treename, particle, max_events=max_events)
-            branches,
+
-            library="np",
+        Ei = data["Ei"]
-            entry_stop=max_events
-        )
-    
-        Ei = data["Tb"]
-        theta = np.radians(data["thetab"])
         sx3Z = data["sx3Z"]
-    
+        EA = data["EA"]
-        mask = np.sin(theta) != 0
+        EC = data["EC"]
-    
+        Esx3 = data["Esx3"]
-        Ei = Ei[mask]
+        Elost = data["Elost"]
-        theta = theta[mask]
+        Eprop = data["Eprop"]
-        sx3Z = sx3Z[mask]
+
-        
+        update_plot_data(f"{particle}_{treename}_Ei", Ei)
-        update_plot_data(f"{particle}_Ei", Ei)
+        update_plot_data(f"{particle}_{treename}_sx3Z", sx3Z)
-        update_plot_data(f"{particle}_sx3Z", sx3Z)
+        update_plot_data(f"{particle}_{treename}_EA", EA)
-    
+        update_plot_data(f"{particle}_{treename}_EC", EC)
-        sin_theta = np.sin(theta)
+        update_plot_data(f"{particle}_{treename}_Esx3", Esx3)
-    
+        update_plot_data(f"{particle}_{treename}_Elost", Elost)
-        radii = np.array([3.2, 4.2, 6.6])
+        update_plot_data(f"{particle}_{treename}_Eprop", Eprop)
-    
+
-        dA = radii[0] / sin_theta
+        base = f"{particle}_{treename}_plots"
-        dC = radii[1] / sin_theta
-        dsx3 = radii[2] / sin_theta
-    
-        print("Computing energies...")
-    
-        EA = energy_loss(particle, Ei, dA)
-        EC = energy_loss(particle, Ei, dC)
-        Esx3 = energy_loss(particle, Ei, dsx3)
-    
-        Elost = Ei - Esx3
-        Eprop = EA - EC
-        
-        update_plot_data(f"{particle}_EA", EA)
-        update_plot_data(f"{particle}_EC", EC)
-        update_plot_data(f"{particle}_Esx3", Esx3)
-        update_plot_data(f"{particle}_Elost", Elost)
-        update_plot_data(f"{particle}_Eprop", Eprop)
-    
-        base = f"{particle}_plots"
         os.makedirs(base, exist_ok=True)
-    
+
-        print(f"Saving plots to folder: {base}")
+        print(f"Saving plots to folder: {base} ({treename})")
-    
+
-        # 1. Histogram: energy loss
+        branch_names = []
+        for key in self.tree.keys():
+            if isinstance(key, bytes):
+                branch_names.append(key.decode("ascii"))
+            elif hasattr(key, "name"):
+                branch_names.append(key.name)
+            else:
+                branch_names.append(str(key))
+        branch_names = list(dict.fromkeys(branch_names))
+
+        if branch_names:
+            print(f"Creating histograms for {len(branch_names)} branches...")
+            all_branches = self.tree.arrays(branch_names, library="np", entry_stop=max_events)
+
+            for branch in branch_names:
+                values = all_branches[branch]
+                try:
+                    values = np.asarray(values, dtype=float)
+                except Exception:
+                    print(f"Skipping non-numeric branch: {branch}")
+                    continue
+                if values.ndim != 1:
+                    print(f"Skipping non-scalar branch: {branch}")
+                    continue
+                values = values[~np.isnan(values)]
+                if values.size == 0:
+                    print(f"Skipping empty branch: {branch}")
+                    continue
+
+                plt.figure(figsize=(7,5))
+                plt.hist(values, bins=100)
+                plt.xlabel(branch)
+                plt.ylabel("Counts")
+                plt.title(f"{particle} ({treename}) {branch} distribution")
+                plt.grid(True)
+                plt.tight_layout()
+                safe_name = re.sub(r"[^0-9A-Za-z_-]", "_", branch)
+                plt.savefig(f"{base}/{safe_name}_hist.png", dpi=300)
+                plt.close()
+        else:
+            print("No branches found to histogram.")
+
         plt.figure(figsize=(7,5))
         plt.hist(Elost, bins=200)
         plt.xlabel("Energy Loss (MeV)")
         plt.ylabel("Counts")
-        plt.title(f"{particle} Total Energy Loss Distribution")
+        plt.title(f"{particle} ({treename}) Total Energy Loss Distribution")
         plt.grid(True)
         plt.tight_layout()
         plt.savefig(f"{base}/Elost_hist.png", dpi=300)
-        plt.close()
+        plt.show()
-    
+
-        #Histogram: sx3Z
         plt.figure(figsize=(7,5))
-        plt.hist(sx3Z, bins=100)
+        plt.hist(Ei, bins=100, histtype='step', label='Vertex Energy')
-        plt.xlabel("SX3 Z")
+        plt.hist(Esx3[Esx3 != 0], bins=100, histtype='step', label='SX3 Energy')
+        plt.xlabel("Energies (MeV)")
         plt.ylabel("Counts")
-        plt.title(f"{particle} SX3 Position Distribution")
+        plt.title(f"{particle} ({treename}) Vertex and SX3 Energy")
         plt.grid(True)
         plt.tight_layout()
-        plt.savefig(f"{base}/sx3Z_hist.png", dpi=300)
+        plt.savefig(f"{base}/Vertex_vs_SX3_Energy.png", dpi=300)
-        plt.close()
+        plt.legend(loc='upper right')
+        plt.show()
+
+        try:
+            plt.figure(figsize=(7,5))
+            plt.hist(sx3Z, bins=100)
+            plt.xlabel("SX3 Z")
+            plt.ylabel("Counts")
+            plt.title(f"{particle} ({treename}) SX3 Position Distribution")
+            plt.grid(True)
+            plt.tight_layout()
+            plt.savefig(f"{base}/sx3Z_hist.png", dpi=300)
+            plt.show()
+        except ValueError:
+            print("Value error")
+
+        try:
+            plt.figure(figsize=(7,6))
+            plt.hist2d(sx3Z, Elost, bins=200)
+            plt.ylabel("Energy Loss (MeV)")
+            plt.xlabel("SX3 Z")
+            plt.title(f"{particle} ({treename}) Energy Loss vs SX3 Position")
+            plt.colorbar(label="Counts")
+            plt.tight_layout()
+            plt.savefig(f"{base}/Elost_vs_sx3Z.png", dpi=300)
+            plt.show()
+        except ValueError:
+            print("Value error")
+
+        try:
+            plt.figure(figsize=(7,6))
+            plt.hist2d(EA, Esx3, bins=200)
+            plt.xlabel("EA (MeV)")
+            plt.ylabel("Esx3 (MeV)")
+            plt.title(f"{particle} ({treename}) Anode vs SX3 Energy")
+            plt.colorbar(label="Counts")
+            plt.tight_layout()
+            plt.savefig(f"{base}/EA_vs_Esx3.png", dpi=300)
+            plt.show()
+        except ValueError:
+            print("Value error")
+
+        try:
+            plt.figure(figsize=(7,6))
+            plt.hist2d(sx3Z, Eprop, bins=200)
+            plt.ylabel("EA - EC (MeV)")
+            plt.xlabel("SX3 Z")
+            plt.title(f"{particle} ({treename}) Energy Propagation Difference vs Position")
+            plt.colorbar(label="Counts")
+            plt.tight_layout()
+            plt.savefig(f"{base}/Eprop_vs_sx3Z.png", dpi=300)
+            plt.show()
+        except ValueError:
+            print("Value error")
+
+        try:
+            mask1 = ~np.isnan(Esx3) & (Esx3 > 0)
+            plt.figure(figsize=(7,6))
+            plt.hist2d(Esx3[mask1], Eprop[mask1], bins=200)
+            plt.ylabel("PCEnergy")
+            plt.xlabel("SX3 Energy (MeV)")
+            plt.title(f"{particle} ({treename}) Energy Propagation Difference vs SX3 Energy")
+            plt.colorbar(label="Counts")
+            plt.xlim(0,30)
+            plt.ylim(0,0.45)
+            plt.tight_layout()
+            plt.savefig(f"{base}/Eprop_vs_Esx3.png", dpi=300)
+            plt.show()
+        except ValueError:
+            print("Value error")
     
-        #2D: Elost vs sx3Z
+        print(f"Plotting complete ({treename}).") 
-        plt.figure(figsize=(7,6))
-        plt.hist2d(sx3Z, Elost, bins=200)
-        plt.ylabel("Energy Loss (MeV)")
-        plt.xlabel("SX3 Z")
-        plt.title(f"{particle} Energy Loss vs SX3 Position")
-        plt.colorbar(label="Counts")
-        plt.tight_layout()
-        plt.savefig(f"{base}/Elost_vs_sx3Z.png", dpi=300)
-        plt.close()
-    
-        #Anode energy vs sx3 energy
-        plt.figure(figsize=(7,6))
-        plt.hist2d(EA, Esx3, bins=200)
-        plt.xlabel("EA (MeV)")
-        plt.ylabel("Esx3 (MeV)")
-        plt.title(f"{particle} Anode vs SX3 Energy")
-        plt.colorbar(label="Counts")
-        plt.tight_layout()
-        plt.savefig(f"{base}/EA_vs_Esx3.png", dpi=300)
-        plt.close()
-        
-        #Prop counter energy loss
-        plt.figure(figsize=(7,6))
-        plt.hist2d(Eprop, sx3Z, bins=200)
-        plt.xlabel("EA - EC (MeV)")
-        plt.ylabel("SX3 Z")
-        plt.title(f"{particle} Energy Propagation Difference vs Position")
-        plt.colorbar(label="Counts")
-        plt.tight_layout()
-        plt.savefig(f"{base}/Eprop_vs_sx3Z.png", dpi=300)
-        plt.close()
-        
-        plt.figure(figsize=(7,6))
-        plt.hist2d(Esx3, Eprop, bins=200)
-        plt.ylabel("PCEnergy")
-        plt.xlabel("SX3 Energy (MeV)")
-        plt.title(f"{particle} Energy Propagation Difference vs sx3 Energy")
-        plt.colorbar(label="Counts")
-        plt.xlim(0,30)
-        plt.ylim(0,.45)
-        plt.tight_layout()
-        plt.savefig(f"{base}/Eprop_vs_Esx3.png", dpi=300)
-        plt.close()
-    
-        print("Plotting complete.")
                 
+    def do_hist_comp(self, arg):
+        particle = arg
+        self.do_set_tree("tree1")
+        data1 = prepare_tree_data(self.tree, "tree1", particle, max_events=None)
+        self.do_set_tree("tree2")
+        data2 = prepare_tree_data(self.tree, "tree2", particle, max_events=None)
+        
+        base = f"{particle}_comp_plots"
+        os.makedirs(base, exist_ok=True)
+        
+        branch_names = []
+        for key in self.tree.keys():
+            if isinstance(key, bytes):
+                branch_names.append(key.decode("ascii"))
+            elif hasattr(key, "name"):
+                branch_names.append(key.name)
+            else:
+                branch_names.append(str(key))
+        branch_names = list(dict.fromkeys(branch_names))
+        
+        print(branch_names)
+        
+        if branch_names:
+            print(f"Creating histograms for {len(branch_names)} branches...")
+            self.do_set_tree("tree1")
+            all_branches1 = self.tree.arrays(branch_names, library="np", entry_stop=None)
+            self.do_set_tree("tree2")
+            all_branches2 = self.tree.arrays(branch_names, library="np", entry_stop=None)
+            for branch in branch_names:
+                values1 = all_branches1[branch]
+                values2 = all_branches2[branch]
+                try:
+                    values1 = np.asarray(values1, dtype=float)
+                    values2 = np.asarray(values2, dtype=float)
+                except Exception:
+                    print(f"Skipping non-numeric branch: {branch}")
+                    continue
+                if values1.ndim != 1 or values2.ndim != 1:
+                    print(f"Skipping non-scalar branch: {branch}")
+                    continue
+                values1 = values1[~np.isnan(values1)]
+                values2 = values2[~np.isnan(values2)]
+                if values1.size == 0 or values2.size == 0:
+                    print(f"Skipping empty branch: {branch}")
+                    continue
+
+                plt.figure(figsize=(7,5))
+                plt.hist(values1, bins=100, histtype='step')
+                plt.hist(values2, bins=100, histtype='step')
+                plt.xlabel(branch)
+                plt.ylabel("Counts")
+                plt.title(f"{particle} {branch} distribution")
+                plt.grid(True)
+                plt.tight_layout()
+                safe_name = re.sub(r"[^0-9A-Za-z_-]", "_", branch)
+                plt.savefig(f"{base}/{safe_name}_hist.png", dpi=300)
+                plt.show()
+        else:
+            print("No branches found to histogram.")
+        
+        
     def do_dual_plotter(self, arg):
         args = shlex.split(arg)
     
         if len(args) < 2:
-            print("Usage: make dual plots proton_data.root alpha_data.root")
+            try:
-            return
+                args = ['SimAnasenProton.root', 'SimAnasenAlpha.root']
+            except:
+                print("Usage: make dual plots proton_data.root alpha_data.root")
+                return
     
         file1 = args[0]
         file2 = args[1]
-    
+        #tree1_name = args[2] if len(args) > 2 else 'tree2'
-        #Helper function
+        #tree2_name = args[3] if len(args) > 3 else 'tree1'
-
-    
-        #Process both files
-        data1 = process_file(f"../Armory/{file1}")
-        data2 = process_file(f"../Armory/{file2}")
-    
         outdir = "dual_plots"
+        
+        data1 = process_file(os.path.join("..", "Armory", file1), "tree1")
+        data2 = process_file(os.path.join("..", "Armory", file2), "tree1")
+
+        #print(f"File one {file1} ({tree1_name}) length {len(data1['Ei'])} \nFile two {file2} ({tree2_name}) length {len(data2['Ei'])}")
         os.makedirs(outdir, exist_ok=True)
     
         print(f"Saving plots to: {outdir}")
     
+    
         #Overlay histogram: Elost
         plt.figure(figsize=(8,6))
     
@@ -1015,7 +1200,7 @@ class MyInteractiveApp(cmd.Cmd):
         plt.tight_layout()
     
         plt.savefig(f"{outdir}/Elost_overlay.png", dpi=300)
-        plt.close()
+        plt.show()
     
         #Overlay histogram: sx3Z
         plt.figure(figsize=(8,6))
@@ -1046,7 +1231,7 @@ class MyInteractiveApp(cmd.Cmd):
         plt.tight_layout()
     
         plt.savefig(f"{outdir}/sx3Z_overlay.png", dpi=300)
-        plt.close()
+        plt.show()
     
         #Side-by-side 2D plots
         fig, axes = plt.subplots(1, 2, figsize=(14,6))
@@ -1077,7 +1262,7 @@ class MyInteractiveApp(cmd.Cmd):
         plt.tight_layout()
     
         plt.savefig(f"{outdir}/Elost_vs_sx3_comparison.png", dpi=300)
-        plt.close()
+        plt.show()
         
         #EA vs Esx3 overlay scatter
         plt.figure(figsize=(8,6))
@@ -1106,48 +1291,148 @@ class MyInteractiveApp(cmd.Cmd):
         plt.tight_layout()
     
         plt.savefig(f"{outdir}/EA_vs_Esx3_overlay.png", dpi=300)
-        plt.close()
+        plt.show()
         
         #PCE vs SiE
         
+        mask1 = data1["Esx3"] > 0
+        mask2 = data2["Esx3"] > 0
         combined_Esx3 = np.concatenate([
-            data1["Esx3"],
+            (data1["Esx3"][mask1]),
-            data2["Esx3"]
+            data2["Esx3"][mask2]
         ])
-        
         combined_Eprop = np.concatenate([
-            data1["Eprop"],
+            data1["Eprop"][mask1],
-            data2["Eprop"]
+            data2["Eprop"][mask2]
         ])
         
-        plt.figure(figsize=(7,6))
+        plt.figure(figsize=(8,6))
-        
         plt.hist2d(
             combined_Esx3,
             combined_Eprop,
-            bins=100,
+            bins=200
-            range=[[0,30],[0,0.45]],
-            cmap='viridis'   # same default matplotlib style
         )
-        
-        plt.ylabel("PCEnergy")
         plt.xlabel("SX3 Energy (MeV)")
-        
+        plt.ylabel("PCEnergy")
-        plt.title(
+        plt.xlim(0,35)
-            f'{data1["particle"]} + {data2["particle"]} '
+        plt.ylim(0,0.6)
-            'Energy Propagation Difference vs SX3 Energy'
-        )
-        
         plt.colorbar(label="Counts")
-        
-        plt.xlim(0,30)
-        plt.ylim(0,0.45)
-        
         plt.tight_layout()
+        plt.savefig(f"{outdir}/Eprop_vs_Esx3.png", dpi=300)
+        plt.show()
+        
+        
+        
+        plt.figure(figsize=(12,6), facecolor='white')
+        plt.hist2d(
+            combined_Esx3,
+            combined_Eprop,
+            bins=[500,500],
+            #range=[[0,35],[0,0.6]],
+            cmap=yellow_jet,
+            cmin=1
+        )
+        plt.xlabel("SX3 Energy (MeV)")
+        plt.ylabel("PCEnergy")
+        plt.xlim(0,35)
+        plt.ylim(0,0.6)
+        cbar = plt.colorbar()
+        cbar.set_label("Counts")
+        plt.minorticks_on()
+        plt.grid(False)
+        plt.tight_layout()
+        plt.savefig(
+            f"{outdir}/ROOT_style_plot.png",
+            dpi=300,
+            facecolor='white'
+        )
+        plt.show()
+        if False:
+            plt.figure(figsize=(7,6))
+            plt.hist2d(
+                combined_Esx3,
+                combined_Eprop,
+                bins=100,
+                range=[[0,30],[0,0.45]],
+                cmap='viridis'
+            )
+            
+            def power_fit_and_plot(x, y, label, color=None):
+            
+                x = np.array(x)
+                y = np.array(y)
+            
+                mask = (
+                    np.isfinite(x) &
+                    np.isfinite(y) &
+                    (x > 0) &
+                    (y > 0)
+                )
+            
+                x = x[mask]
+                y = y[mask]
+            
+                logx = np.log(x)
+                logy = np.log(y)
+            
+                b, loga = np.polyfit(logx, logy, 1)
+            
+                a = np.exp(loga)
+            
+                y_pred = a * x**b
+            
+                ss_res = np.sum((y - y_pred)**2)
+                ss_tot = np.sum((y - np.mean(y))**2)
+            
+                r2 = 1 - (ss_res / ss_tot)
+            
+                x_fit = np.linspace(np.min(x), np.max(x), 1000)
+                y_fit = a * x_fit**b
+            
+                plt.plot(
+                    x_fit,
+                    y_fit,
+                    linewidth=3,
+                    color=color,
+                    label=f'{label} fit ($R^2$={r2:.4f})'
+                )
+            
+                print(f"\n{label} power-law fit:")
+                print(f"y = {a:.6f} * x^{b:.6f}")
+                print(f"R^2 = {r2:.6f}")
+            
+                return a, b, r2
+        
+            
+            power_fit_and_plot(
+                data1["Esx3"]["Esx3" > 0],
+                data1["Eprop"]["Esx3" > 0],
+                label=data1["particle"],
+                color='red'
+            )
+            
+            power_fit_and_plot(
+                data2["Esx3"],
+                data2["Eprop"],
+                label=data2["particle"],
+                color='cyan'
+            )
+            
+            plt.ylabel("PCEnergy")
+            plt.xlabel("SX3 Energy (MeV)")
+            plt.title(
+                f'{data1["particle"]} + {data2["particle"]} '
+                'Energy Propagation Difference vs SX3 Energy'
+            )
+            plt.colorbar(label="Counts")
+            #plt.xlim(0,30)
+            #plt.ylim(0,0.45)
+            plt.legend()
+            plt.tight_layout()
+            plt.savefig(f"{outdir}/Combined_Eprop_vs_Esx3_fit.png",dpi=300)
+            plt.show()
         
-        plt.savefig(f"{outdir}/Combined_Eprop_vs_Esx3.png", dpi=300)
         
-        plt.close()
     
         print("Dual plotting complete.")
         

ELoss/anasenMS(ElossVersion).cpp

@@ -0,0 +1,518 @@
+#include "TRandom.h"      // ROOT random number generators, gRandom
+#include "TFile.h"        // ROOT file I/O
+#include "TTree.h"        // ROOT tree storage
+#include "TH1.h"          // 1D histograms
+#include "TH2.h"          // 2D histograms
+#include "TStyle.h"       // ROOT plotting style controls
+#include "TCanvas.h"      // ROOT canvas drawing
+#include "TBenchmark.h"   // timing measurement
+#include "TGraph.h"       // for energy loss interpolation
+#include <cstring>
+#include "TApplication.h" // ROOT app loop
+#include "ClassTransfer.h" // Reaction kinematics and MC event generation
+#include "ClassAnasen.h"   // ANASEN detector model classes (SX3, PW, etc.)
+#include <stdio.h>
+#include <stdlib.h>
+#include <set>
+#include "TLegend.h"
+#include "TH1D.h"
+#include "TObjArray.h"
+#include "TBranch.h"
+#include <iostream>
+#include <fstream>
+
+//======== Generate light particle based on reaction
+// calculate real and reconstructed tracks and Q-value uncertainty
+
+// Function to load energy loss table from file
+TGraph* LoadELoss(const char* filename) {
+  TGraph* g = new TGraph(filename, "%lg %lg");
+  return g;
+}
+
+bool IsDeadAnode(int id){
+    static std::set<int> dead = {}; // add dead anode IDs here, 0-23
+    return dead.count(id);
+}
+
+bool IsDeadCathode(int id){
+    static std::set<int> dead = {}; // add dead cathode IDs here, 0-23
+    return dead.count(id);
+}
+
+bool IsDeadSX3(int id){ 
+    static std::set<int> dead = {}; // add dead SX3 IDs here, 0-23 1,7,9,3
+    return dead.count(id);
+}
+
+static std::set<pair<int,int>> ReactionProductb = { {1,1} }; // add reaction product b (light particle) A,Z pairs here, e.g. {1,1} for proton, {4,2} for alpha 
+
+int main(int argc, char **argv){
+
+  printf("=========================================\n");
+  printf("===       ANASEN Monte Carlo          ===\n");
+  printf("=========================================\n");
+
+  // number of events can be overridden from command line
+  int numEvent = 1000000;
+  if( argc >= 2 ) numEvent = atoi(argv[1]);
+
+  // load energy loss tables (assume units: E in MeV, dE/dx in MeV/(mg/cm^2), density in mg/cm^3)
+  TGraph* elossAlpha = LoadELoss("../ELoss/E_vs_x_alpha.dat"); // for light particle (alpha)
+  TGraph* elossProton = LoadELoss("../ELoss/E_vs_x_proton.dat");     // for heavy particle (proton)
+  TGraph *invgAlpha = new TGraph(elossAlpha->GetN(), elossAlpha->GetY(), elossAlpha->GetX());
+  TGraph *invgProton = new TGraph(elossProton->GetN(), elossProton->GetY(), elossProton->GetX());
+
+  
+  //Plot energy loss tables (sanity check), vis will not work if this is ran without X11 display (e.g. on cluster), so comment out if running in headless mode
+  auto c1 = new TCanvas("c1", "Graph Example", 800, 600);
+  auto g = elossAlpha;
+  g->SetTitle("Energy Loss Table (Alpha);cm;Kinetic Energy (MeV)");
+  g->Draw("ALP");
+  g->SetLineColor(kRed);
+  //c1->SetLogy();
+  //c1->SetLogx();
+  c1->Print("eloss_alpha.png");
+
+  auto c2 = new TCanvas("c2", "Graph Example", 800, 600);
+  auto g2 = elossProton;
+  g2->SetTitle("Energy Loss Table (Proton);cm;Kinetic Energy (MeV)");
+  g2->Draw("ALP");
+  g2->SetLineColor(kBlue);
+  c2->Print("eloss_proton.png");
+
+  // Reaction setup: projectile + target configuration, energy, and product IDs
+  TransferReaction transfer;
+
+  transfer.SetA(18, 10, 0);           // e.g., 24Mg (Z=12) with 0 excitation
+  transfer.SetIncidentEnergyAngle((4.4), 0, 0); // arguments are KEA in MeV/u, theta and phi in degree
+  transfer.Seta( 4, 2);              // identify reaction product a in internal indexing e.g., 4He (alpha)
+  transfer.Setb(ReactionProductb.begin()->first, ReactionProductb.begin()->second);              // identify reaction product b e.g., 1H (proton)
+  transfer.SetB(21, 11);              // identify reaction product B e.g., 23Na (Z=11)
+
+  // TODO add alpha source or alternative reaction channel selection
+
+  // Excited state lists (target and projectile/excited products)
+  std::vector<float> ExAList = {0}; // projectile excitation states in MeV
+  std::vector<float> ExList = {0}; // target excitation states in MeV
+
+  // define vertex position uniform distribution ranges (mm)
+  double vertexXRange[2] = {   -5,    5}; // mm
+  double vertexYRange[2] = {   -5,    5};
+  double vertexZRange[2] = { -100,  100};
+
+  // detector resolution / uncertainty parameters
+  double sigmaSX3_W = -1; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension)
+  double sigmaSX3_L = 3;  // mm, vertical spread for SX3
+  double sigmaPW_A  = 0;  // normalized anode uncertainty term (0-1)
+  double sigmaPW_C  = 0;  // normalized cathode uncertainty term (0-1)
+
+  // status printout
+  printf("------------ Vertex :\n");
+  printf("X : %7.2f - %7.2f mm\n", vertexXRange[0], vertexXRange[1]);
+  printf("Y : %7.2f - %7.2f mm\n", vertexYRange[0], vertexYRange[1]);
+  printf("Z : %7.2f - %7.2f mm\n", vertexZRange[0], vertexZRange[1]);
+  printf("------------ Uncertainty :\n");
+  printf(" SX3 horizontal : %.1f\n", sigmaSX3_W);
+  printf(" SX3   vertical : %.1f\n", sigmaSX3_L);
+  printf("          Anode : %.1f mm\n", sigmaPW_A);
+  printf("        Cathode : %.1f mm\n", sigmaPW_C);
+  printf("        num_eve : %d \n",numEvent);
+
+  // calculates energy/momentum/kinematics constants for transfer reaction
+  transfer.CalReactionConstant();
+
+  int nExA = ExAList.size();
+  int nEx  = ExList.size();
+
+  // optional visualization control: pass "vis" as 3rd arg
+  bool enableVis = (argc >= 3 && strcmp(argv[2], "vis") == 0);
+  TApplication *app = nullptr;
+  if(enableVis){
+    app = new TApplication("anasenVis", &argc, argv);
+  }
+
+  // storage for tracks during simulation (for visualization)
+  std::vector<TVector3> visTrackVertex, visTrackDir, visTrackHitPos;
+  std::vector<std::pair<int,int>> visTrackWires;  // {anodeID, cathodeID}
+
+  // create detector representation in memory
+  ANASEN * anasen = new ANASEN();       // top-level detector object
+  SX3 * sx3 = anasen->GetSX3();         // silicon array part
+  PW * pw = anasen->GetPW();            // proportional wire chamber part
+
+  // output file + tree
+  TString saveFileName = "SimAnasen1.root";
+  printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
+  TFile * saveFile = new TFile(saveFileName, "recreate");
+  TTree * tree = new TTree("tree", "tree");
+   
+
+  // beam and CM variables saved in tree
+  double KEA;
+  tree->Branch("beamKEA",     &KEA, "beamKEA/D");
+
+  double thetaCM, phiCM;
+  tree->Branch("thetaCM", &thetaCM, "thetaCM/D");
+  tree->Branch("phiCM",     &phiCM, "phiCM/D");
+
+  // outgoing particles in lab frame (light/heavy)
+  double thetab, phib, Tb;
+  double thetaB, phiB, TB;
+  double dEb;
+  double dEB;
+  std::array<double, 2> T;
+  tree->Branch("thetab", &thetab, "thetab/D"); // polar angle of light particle in lab frame
+  tree->Branch("phib",     &phib, "phib/D"); // azimuthal angle of light particle in lab frame
+  tree->Branch("Tb",         &Tb, "Tb/D"); // kinetic energy of light particle at vertex (before energy loss)
+  tree->Branch("thetaB", &thetaB, "thetaB/D");
+  tree->Branch("phiB",     &phiB, "phiB/D");
+  tree->Branch("TB",         &TB, "TB/D"); // kinetic energy of heavy particle at vertex (before energy loss)
+  tree->Branch("dEb",        &dEb, "dEb/D");
+  tree->Branch("dEB",        &dEB, "dEB/D"); // placeholder for heavy particle energy loss, currently set equal to light particle loss for simplicity
+  tree->Branch("T",          &T, "T/D");   // placeholder for true Q-value, currently set to 0 for simplicity
+
+  // excitation state identifiers
+  int ExAID;
+  double ExA;
+  tree->Branch("ExAID", &ExAID, "ExAID/I"); // projectile excitation state ID
+  tree->Branch("ExA",     &ExA, "ExA/D"); // projectile excitation energy in MeV
+
+  int ExID;
+  double Ex;
+  tree->Branch("ExID", &ExID, "ExID/I"); // target excitation state ID
+  tree->Branch("Ex",     &Ex, "Ex/D"); // target excitation energy in MeV
+
+  // true vertex position in target volume
+  double vertexX, vertexY, vertexZ;
+  tree->Branch("vX",     &vertexX, "VertexX/D"); // true vertex X position in mm
+  tree->Branch("vY",     &vertexY, "VertexY/D"); // true vertex Y position in mm
+  tree->Branch("vZ",     &vertexZ, "VertexZ/D"); // true vertex Z position in mm
+
+  // reconstructed SX3 hit position
+  double sx3X, sx3Y, sx3Z;
+  tree->Branch("sx3X",     &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing)
+  tree->Branch("sx3Y",     &sx3Y, "sx3Y/D"); // reconstructed Y position from SX3 (with optional smearing)
+  tree->Branch("sx3Z",     &sx3Z, "sx3Z/D"); // reconstructed Z position from SX3 (with optional smearing)
+
+  // PW nearest and next nearest wires
+  int anodeID[2], cathodeID[2];
+  tree->Branch("aID", anodeID, "anodeID/I"); // anodeID[0] is nearest anode wire, anodeID[1] is next nearest anode wire
+  tree->Branch("cID", cathodeID, "cathodeID/I"); // cathodeID[0] is nearest cathode wire, cathodeID[1] is next nearest cathode wire
+
+  // distances to nearest wires
+  double anodeDist[2], cathodeDist[2];
+  tree->Branch("aDist",   anodeDist, "anodeDist/D");
+  tree->Branch("cDist",  cathodeDist, "cathodeDist/D");
+
+  // SX3 channel assignment and Z fraction (depth) information
+  int sx3ID, sx3Up, sx3Dn, sx3Bk;
+  double sx3ZFrac;
+  tree->Branch("sx3ID",   &sx3ID, "sx3ID/I");
+  tree->Branch("sx3Up",   &sx3Up, "sx3Up/I");
+  tree->Branch("sx3Dn",   &sx3Dn, "sx3Dn/I");
+  tree->Branch("sx3Bk",   &sx3Bk, "sx3Bk/I");
+  tree->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D");
+
+  // reconstructed angles from PW track fit, method 1 and 2
+  double reTheta, rePhi;
+  tree->Branch("reTheta", &reTheta, "reconstucted_theta/D");
+  tree->Branch("rePhi",     &rePhi, "reconstucted_phi/D");
+
+  double reTheta1, rePhi1;
+  tree->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D");
+  tree->Branch("rePhi1",     &rePhi1, "reconstucted_phi1/D");
+
+  // reconstructed vertex Z from PW fit
+  double z0;
+  tree->Branch("z0",     &z0, "reconstucted_Z/D");
+
+  TTree* tree2 = tree->CloneTree(0);
+  tree2->SetName("tree2");
+
+  //========timer
+  TBenchmark clock;
+  bool shown ;
+  clock.Reset();
+  clock.Start("timer");
+  shown = false;
+  int ELossTotal = 0;
+
+  //================================= Calculate event loop
+  for( int i = 0; i < numEvent ; i++){
+
+    // randomly sample target/projectile excitations
+    ExAID = gRandom->Integer(nExA);
+    ExA = ExAList[ExAID];
+    transfer.SetExA(ExA);
+
+    ExID = gRandom->Integer(nEx);
+    Ex = ExList[ExID];
+    transfer.SetExB(Ex);
+
+    // recalc kinematic constants for chosen states
+    transfer.CalReactionConstant();
+
+    // isotropic CM direction
+    thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ;
+    phiCM   = (gRandom->Rndm() - 0.5) * TMath::TwoPi();
+
+    //==== Calculate reaction kinematics in lab frame
+    TLorentzVector * output = transfer.Event(thetaCM, phiCM); // returns array of outputs
+    TLorentzVector Pb = output[2]; // light particle or product A
+    TLorentzVector PB = output[3]; // heavy particle or product B
+
+    thetab = Pb.Theta() * TMath::RadToDeg();
+    thetaB = PB.Theta() * TMath::RadToDeg();
+
+    Tb = (Pb.E() - Pb.M()); // kinetic energy of light particle at vertex (before energy loss) units of MeV
+    TB = (PB.E() - PB.M());
+    T[0] = Tb;
+    T[1] = TB;
+    //if (Tb < 1.5) {
+    //  //skip event if light particle energy after loss is below detection threshold of 1.5 MeV
+    //  continue;
+    //}
+
+    phib = Pb.Phi() * TMath::RadToDeg();
+    phiB = PB.Phi() * TMath::RadToDeg();
+
+    // vertex position in target volume
+    vertexX = (vertexXRange[1]- vertexXRange[0])*gRandom->Rndm() +  vertexXRange[0];
+    vertexY = (vertexYRange[1]- vertexYRange[0])*gRandom->Rndm() +  vertexYRange[0];
+    vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() +  vertexZRange[0];
+
+    TVector3 vertex(vertexX, vertexY, vertexZ);
+
+    // set direction vector from lab angle
+    TVector3 dir(1, 0, 0);
+    dir.SetTheta(thetab * TMath::DegToRad());
+    dir.SetPhi(phib * TMath::DegToRad());
+
+    // run detector response models for PW and SX3
+    pw->FindWireID(vertex, dir, false);
+    sx3->FindSX3Pos(vertex, dir, false);
+
+    PWHitInfo hitInfo = pw->GetHitInfo();
+
+    anodeID[0]   = hitInfo.nearestWire.first; // nearest anode wire ID
+    cathodeID[0] = hitInfo.nearestWire.second; // nearest cathode wire ID
+    anodeID[1]   = hitInfo.nextNearestWire.first; // next nearest anode wire ID
+    cathodeID[1] = hitInfo.nextNearestWire.second; // next nearest cathode wire ID
+
+    anodeDist[1]   = hitInfo.nextNearestDist.first; // distance to next nearest anode wire
+    cathodeDist[1] = hitInfo.nextNearestDist.second; // distance to next nearest cathode wire
+
+    if(IsDeadAnode(anodeID[0])) continue;
+    if(IsDeadCathode(cathodeID[0])) continue;
+
+    // SX3 hit channel info and depth fraction
+    sx3ID = sx3->GetID();
+
+    if(IsDeadSX3(sx3ID)) continue;
+
+    anodeDist[0]   = hitInfo.nearestDist.first; // distance to nearest anode wire
+    cathodeDist[0] = hitInfo.nearestDist.second; // distance to nearest cathode wire
+
+    if( sx3ID >= 0 ){
+      sx3Up    = sx3->GetChUp();
+      sx3Dn    = sx3->GetChDn();
+      sx3Bk    = sx3->GetChBk();
+      sx3ZFrac = sx3->GetZFrac();
+      
+      // apply intrinsic detector resolution to true SX3 hit position
+      // for no smearing comment out and use GetHitPos();
+      TVector3 hitPos = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
+
+      sx3X = hitPos.X();
+      sx3Y = hitPos.Y();
+      sx3Z = hitPos.Z();
+
+      // store track data for visualization if enabled
+      if(enableVis){
+        visTrackVertex.push_back(vertex);
+        visTrackDir.push_back(dir);
+        visTrackHitPos.push_back(hitPos);
+        visTrackWires.push_back({anodeID[0], cathodeID[0]});
+      }
+      // reconstruct track from PW readings + SX3 hit
+      pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
+      reTheta = pw->GetTrackTheta() * TMath::RadToDeg();
+      rePhi   = pw->GetTrackPhi() * TMath::RadToDeg();
+
+      // alternative track algorithm with uncertainty parameters
+      pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false);
+      reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg();
+      rePhi1   = pw->GetTrackPhi() * TMath::RadToDeg();
+
+      z0 = pw->GetZ0();
+      dEb = 0;
+      dEB = 0;
+      tree->Fill();
+
+      //Energy loss
+      double dl = (hitPos - vertex).Mag(); // path length in units of cm
+      if (numEvent <= 100){
+        //printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl);
+        //printf("Total T before loss: %f MeV\n", T);
+      }
+      double tb_temp = Tb;
+
+      dEb = tb_temp - Tb; // total energy loss
+      if (ReactionProductb.count({4, 2})){ // if light particle is alpha, use alpha energy loss table
+        double x0b = invgAlpha->Eval(Tb);
+        x0b = x0b + dl;
+        Tb = elossAlpha->Eval(x0b);
+      } else if (ReactionProductb.count({1, 1})){ // if light particle is proton, use proton energy loss table
+        double x0b = invgProton->Eval(Tb);
+        x0b = x0b + dl;
+        Tb = elossProton->Eval(x0b);
+      } else {
+        // for other particle types, can add additional energy loss tables or use a generic approximation
+        // for now, we will just apply a simple linear energy loss as a placeholder
+        double dE_dx = 5; // MeV/cm, placeholder value for energy loss per unit length
+        Tb = Tb - dE_dx * dl;
+      }
+
+      //if (Tb < 0) {
+      //  Tb = TMath::QuietNaN();
+      //}
+
+      dEb = tb_temp - Tb; // total energy loss
+      
+      // fill tree2 with energy loss adjusted data
+      //Fill T so it can make a histogram of both Tb and TB in root script
+      T[0] = Tb;
+      T[1] = 0; 
+      //to plot both as one histogram in root, can use tree2->Draw("T(0)"); for light particle and tree2->Draw("T(1)") for heavy particle
+
+      tree2->Fill();
+      
+      if (numEvent <= 10){
+        //printf("Event %d: Tb after energy loss = %f MeV, energy loss = %f MeV\n", i, Tb, tb_temp - Tb);
+      } //to give in scientific notation, use %e instead of %f in the printf format string. For example: printf("Event %d: Tb after energy loss = %e MeV, energy loss = %e MeV\n", i, Tb, tb_temp - Tb);
+      ELossTotal += (tb_temp - Tb);
+
+    }else{
+      // no valid SX3 hit: mark clearly invalid
+      sx3Up = -1;
+      sx3Dn = -1;
+      sx3Bk = -1;
+      sx3ZFrac = TMath::QuietNaN();
+
+      sx3X = TMath::QuietNaN();
+      sx3Y = TMath::QuietNaN();
+      sx3Z = TMath::QuietNaN();
+
+      reTheta  = TMath::QuietNaN();
+      rePhi    = TMath::QuietNaN();
+      reTheta1 = TMath::QuietNaN();
+      rePhi1   = TMath::QuietNaN();
+      z0       = TMath::QuietNaN();
+      dEb      = TMath::QuietNaN();
+      dEB      = TMath::QuietNaN();
+      //Tb = -12354567; // mark kinetic energy as invalid for no hit case
+      // fill tree with original data (no energy loss for these events)
+      //comment out tree fill for no hit case  
+      //tree->Fill();
+      //tree2->Fill();
+    }
+
+    //#################################################################### Timer
+    // measure elapsed real time and print progress roughly every 10 sec
+    clock.Stop("timer");
+    Double_t time = clock.GetRealTime("timer");
+    clock.Start("timer");
+
+    if ( !shown ) {
+      if (fmod(time, 10) < 1 ){
+        printf( "%10d[%2d%%]| %8.2f sec | expect: %5.1f min \n", i, TMath::Nint((i+1)*100./numEvent), time , numEvent*time/(i+1)/60);
+        shown = 1;
+      }
+    } else {
+      if (fmod(time, 10) > 9 ){
+        shown = 0;
+      }
+    }
+
+  }
+
+  // write results to ROOT file and close
+  //tree->Write();
+  //tree2->Write();
+  tree->Write("", TObject::kOverwrite);
+  tree2->Write("", TObject::kOverwrite);
+  int count = tree->GetEntries();
+  int count2 = tree2->GetEntries();
+  saveFile->Close();
+
+  printf("=============== done. saved as %s. tree entries: %d, tree2 entries: %d\n", saveFileName.Data(), count, count2);
+  printf("Total energy loss across all events: %f MeV\n", (double)ELossTotal);
+  printf("Average energy loss across events: %f MeV\n", (double)ELossTotal / count);
+
+  if(enableVis){ // to enable visualization, run with 3rd argument "vis", e.g. "./anasenMC 1000 vis"
+    printf("Displaying geometry with %zu tracks from simulation\n", visTrackVertex.size());
+    
+    // Build full geometry with all wires
+    anasen->DrawAnasen(0, 23, 0, 23, -1, true);
+    
+    // Add all stored tracks to the geometry
+    TGeoManager *geom = anasen->GetGeoManager();
+    TGeoVolume *worldBox = anasen->GetWorldBox();
+    
+    if(geom && worldBox && visTrackVertex.size() > 0){
+      int trackNodeID = 500;  // start node IDs for tracks
+      
+      for(size_t iTrack = 0; iTrack < visTrackVertex.size(); ++iTrack){
+        TVector3 vertex = visTrackVertex[iTrack];
+        TVector3 dir = visTrackDir[iTrack];
+        TVector3 hitPos = visTrackHitPos[iTrack];
+        
+        double theta = dir.Theta() * TMath::RadToDeg();
+        double phi = dir.Phi() * TMath::RadToDeg();
+        
+        // Add a line marker at the vertex
+        TGeoVolume *startMarker = geom->MakeSphere("startMarker", 0, 0, 2.0);
+        startMarker->SetLineColor(kBlack);
+        worldBox->AddNode(startMarker, trackNodeID, 
+                         new TGeoCombiTrans(vertex.X(), vertex.Y(), vertex.Z(), 
+                                          new TGeoRotation("rot", 0, 0, 0)));
+        trackNodeID++;
+        
+        // Add track line from vertex toward hit position
+        TGeoVolume *trackLine = geom->MakeTube("trackLine", 0, 0, 0.08, 150.0);
+        trackLine->SetLineColor(kBlue);
+        worldBox->AddNode(trackLine, trackNodeID,
+                         new TGeoCombiTrans(vertex.X(), vertex.Y(), vertex.Z(),
+                                          new TGeoRotation("rotTrack", phi + 90, theta, 0)));
+        trackNodeID++;
+        
+        // Add hit position marker
+        TGeoVolume *hitMarker = geom->MakeSphere("hitMarker", 0, 0, 2.0);
+        hitMarker->SetLineColor(kRed);
+        worldBox->AddNode(hitMarker, trackNodeID,
+                         new TGeoCombiTrans(hitPos.X(), hitPos.Y(), hitPos.Z(),
+                                          new TGeoRotation("rotHit", 0, 0, 0)));
+        trackNodeID++;
+      }
+      
+      // Redraw geometry with all tracks
+      geom->CloseGeometry();
+      geom->SetVisLevel(4);
+      worldBox->Draw("ogle");
+    }
+    
+    if(app){
+      printf("Entering ROOT event loop\n");
+      app->Run();
+    }
+  }
+
+  delete anasen;
+  delete elossAlpha;
+  delete elossProton;
+
+
+  return 0;
+
+}