fixed alpha punch-through plots

Armory/ClassQQQ.h

@@ -100,13 +100,13 @@ inline void QQQ::Clear(){
 }
 
 inline void QQQ::ConstructGeo(){
-  TGeoVolume *qqq = geom->MakeTubs("qqq", Al, qqqR1, qqqR2, 0.5, 5, 85);
+  TGeoVolume *qqq = geom->MakeTubs("qqq", Al, qqqR1, qqqR2, 0.5, 5, 85); // thickness 0.5 mm, phi from 5 to 90 deg in each quadrant
   qqq->SetLineColor(7);
   for( int i = 0; i < 4; i++){
-    worldBox->AddNode(qqq, i+1, new TGeoCombiTrans( 0, 
+    worldBox->AddNode(qqq, i+1, new TGeoCombiTrans( 0,
                                                     0, 
                                                     qqqZPos, 
-                                                    new TGeoRotation("rot1", 360/4 * (i), 0., 0.)));
+                                                    new TGeoRotation("rot1", 360/4 * (i), 0., 0.))); //arguments are (name, material, inner radius, outer radius, half length in z, start phi, delta phi
   }
 }
 
@@ -120,10 +120,6 @@ inline void QQQ::FindQQQPos(TVector3 pos,
   chDn = -1;
   chBk = -1;
 
-  //--------------------------------------------
-  // Intersect trajectory with QQQ plane
-  //--------------------------------------------
-
   if( TMath::Abs(direction.Z()) < 1e-10 ) return;
 
   double t = (qqqZPos - pos.Z()) / direction.Z();
@@ -166,10 +162,6 @@ inline void QQQ::FindQQQPos(TVector3 pos,
 
   if( id < 0 ) return;
 
-  //--------------------------------------------
-  // Ring number (32 strips)
-  //--------------------------------------------
-
   const double ringWidth =
       (qqqR2 - qqqR1)/32.0;
 

Armory/anasenMS.cpp

@@ -101,8 +101,8 @@ int main(int argc, char **argv){
   transfer.SetA(27, 13, 0);            // 18Ne projectile
   TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Al-27.dat"); 
   transfer.Seta(4, 2);                 // 4He target
-  transfer.Setb(1, 1);                 // outgoing proton from the primary transfer
+  transfer.Setb(4, 2);                 // outgoing proton from the primary transfer
-  transfer.SetB(30, 14);               // 21Na* heavy product
+  transfer.SetB(27, 13);               // 21Na* heavy product
   const double beamA = 27;                   // mass number of 27Al beam
 
   bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable

ELoss/PCEnergyAnalysis.py

@@ -5,6 +5,7 @@ Created on Wed May 20 13:32:14 2026
 
 @author: jamesszalkie
 """
+import gc
 import numpy as np
 import pandas as pd
 from scipy.interpolate import interp1d
@@ -74,6 +75,13 @@ particles = {
 interp_cache = {}
 
 
+def clear_interpolator_cache():
+    """Drop cached energy-loss interpolators and force Python to release memory."""
+    global interp_cache
+    interp_cache.clear()
+    gc.collect()
+
+
 def make_E_vs_x(
     z,
     mass_u,
@@ -121,7 +129,7 @@ def make_E_vs_x(
         if z == 2:
             rho_g_cm3 = 0.00017
             #rho_g_cm3 = 8.2928e-5
-        if z == 1:
+        else:
             rho_g_cm3 = 8.2928e-5
         
         gas = catima.Material(material_def)
@@ -164,6 +172,7 @@ def make_E_vs_x(
     print(f"[INFO] saved: {outfile}")
 
     interp_cache.pop((label.lower(), medium.lower()), None)
+    gc.collect()
 
     return x, E
 
@@ -365,7 +374,7 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
     radii = np.array([3.7, 4.3])
     dA = (radii[0] - np.sqrt((vXsx3/10)**2 + (vYsx3/10)**2))/ sin_theta
     dC = (radii[1] - np.sqrt((vXsx3/10)**2 + (vYsx3/10)**2))/ sin_theta
-    lsx3 = (0.1 / sin_theta) * 10
+    lsx3 = (.1 / sin_theta)
 
     # Filter out unphysical distances (negative or unreasonably small)
     # These typically occur at large angles where trajectory doesn't properly intersect proportional counters
@@ -386,13 +395,16 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
     print(f"Computing energies for {particle} ({treename})...")
     #print(f"  Retained {np.sum(distance_mask)} / {len(distance_mask)} events after distance filter")
 
+    clear_interpolator_cache()
     EA = energy_loss(particle, "He", Eisx3, dA)
     EC = energy_loss(particle, "He", Eisx3, dC)
     Esx3 = energy_loss(particle, "He", Eisx3, dsx3)
     Eqqq = energy_loss(particle, "He", Eiqqq, dqqq)
+
+    clear_interpolator_cache()
     Efinal = energy_loss(particle, "Si", Esx3, lsx3)
-    
+    Edet = Esx3 - Efinal
-    Edet = np.where(Efinal <= 0, Esx3, Esx3 - Efinal)
+    #Edet = np.where(Efinal <= 0, Esx3, Esx3 - Efinal)
 
     Elost = Eisx3 - Esx3
     Elostqqq = Eiqqq - Eqqq
@@ -421,19 +433,24 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
         "Edet": Edet
     }
 
+def infer_particle_from_filename(filename):
+    lower = os.path.basename(filename).lower()
+    if "proton" in lower:
+        return "proton"
+    if "alpha" in lower:
+        return "alpha"
+    if "deuteron" in lower:
+        return "deuteron"
+    return None
+
+
 def process_file(filename, treename, particle=None, max_events=None):
     tree = uproot.open(filename)[f"{treename}"]
 
     if particle is None:
-        lower = filename.lower()
+        particle = infer_particle_from_filename(filename)
-        if "proton" in lower:
+    if particle is None:
-            particle = "proton"
+        particle = "proton"
-        elif "alpha" in lower:
-            particle = "alpha"
-        elif "deuteron" in lower:
-            particle = "deuteron"
-        else:
-            particle = "proton"
     print(f"File {filename} particle {particle}, tree {treename}")
     return prepare_tree_data(tree, treename, particle, max_events=max_events)
 
@@ -1281,13 +1298,7 @@ class MyInteractiveApp(cmd.Cmd):
             plt.show()
         except ValueError:
             print("Value error")
-            
+    
-        #try:
-
-        #except ValueError:
-        #    print("Value error")
-            
-            
         branch_names = []
         for key in self.tree.keys():
             if isinstance(key, bytes):
@@ -1301,9 +1312,6 @@ class MyInteractiveApp(cmd.Cmd):
         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)
-        
-            # Keep only events with thetab >= 45
-            #mask = (all_branches["Tb"] > 0) & (all_branches["qqqTb"] > 0)
 
             for branch in branch_names:
                 values = all_branches[branch]
@@ -1341,8 +1349,6 @@ class MyInteractiveApp(cmd.Cmd):
         else:
             print("No branches found to histogram.")
 
-
-    
         print(f"Plotting complete ({treename}).") 
                 
     def do_hist_comp(self, arg):
@@ -1360,7 +1366,6 @@ class MyInteractiveApp(cmd.Cmd):
             else:
                 branch_names.append(str(key))
         branch_names = list(dict.fromkeys(branch_names))
-        
         print(branch_names)
         
         if branch_names:
@@ -1401,7 +1406,6 @@ class MyInteractiveApp(cmd.Cmd):
         else:
             print("No branches found to histogram.")
         
-        
     def do_dual_plotter(self, arg):
     
         args = shlex.split(arg)
@@ -1425,21 +1429,25 @@ class MyInteractiveApp(cmd.Cmd):
         for file in files:
     
             try:
-                # If you want to combine tree1 + tree2:
+                resolved_path = os.path.join("..", "Armory", file)
+                resolved_particle = infer_particle_from_filename(resolved_path)
+                
                 tree1 = process_file(
-                    os.path.join("..", "Armory", file),
+                    resolved_path,
-                    "tree1"
+                    "tree1",
+                    particle=resolved_particle
                 )
     
                 try:
                     tree2 = process_file(
-                        os.path.join("..", "Armory", file),
+                        resolved_path,
-                        "tree2"
+                        "tree2",
+                        particle=resolved_particle
                     )
     
                     data = {
                         "particle":
-                        f"{tree1['particle']}_combined"
+                        f"{resolved_particle or tree1['particle']}_combined"
                     }
     
                     for key in tree1:
@@ -1474,8 +1482,9 @@ class MyInteractiveApp(cmd.Cmd):
         plt.figure(figsize=(8, 6))
     
         for data in datasets:
+            dep_values = data.get("Edet", data.get("Elost", []))
             plt.hist(
-                data["Elost"],
+                dep_values,
                 bins=200,
                 histtype="step",
                 linewidth=2,
@@ -1483,15 +1492,15 @@ class MyInteractiveApp(cmd.Cmd):
                 label=data["particle"]
             )
     
-        plt.xlabel("Energy Loss (MeV)")
+        plt.xlabel("Energy Deposition (MeV)")
         plt.ylabel("Normalized Counts")
-        plt.title("Energy Loss Comparison")
+        plt.title("Energy Deposition Comparison")
         plt.legend()
         plt.grid(True)
         plt.tight_layout()
     
         plt.savefig(
-            f"{outdir}/Elost_overlay.png",
+            f"{outdir}/Edet_overlay.png",
             dpi=300
         )
     
@@ -1536,18 +1545,19 @@ class MyInteractiveApp(cmd.Cmd):
     
         for ax, data in zip(axes, datasets):
     
+            dep_values = data.get("Edet", data.get("Elost", []))
             h = ax.hist2d(
                 data["sx3Z"],
-                data["Elost"],
+                dep_values,
                 bins=200
             )
     
             ax.set_title(
-                f'{data["particle"]}\nElost vs SX3'
+                f'{data["particle"]}\nEnergy Deposition vs SX3'
             )
     
             ax.set_xlabel("SX3 Z")
-            ax.set_ylabel("Energy Loss (MeV)")
+            ax.set_ylabel("Energy Deposition (MeV)")
     
             fig.colorbar(
                 h[3],
@@ -1558,11 +1568,40 @@ class MyInteractiveApp(cmd.Cmd):
         plt.tight_layout()
     
         plt.savefig(
-            f"{outdir}/Elost_vs_sx3_comparison.png",
+            f"{outdir}/Edet_vs_sx3_comparison.png",
             dpi=300
         )
     
         plt.show()
+
+        for data in datasets:
+            mask = (
+                ~np.isnan(data.get("Edet", [])) &
+                (data.get("Edet", []) > 0) &
+                ~np.isnan(data.get("Eprop", [])) &
+                ~np.isnan(data.get("thetab", []))
+            )
+            if not np.any(mask):
+                print(f"No valid events for {data['particle']} when plotting Eprop vs Edet.")
+                continue
+
+            x = np.asarray(data["Edet"])[mask]
+            y = np.asarray(data["Eprop"])[mask] * np.asarray(data["thetab"])[mask]
+
+            plt.figure(figsize=(7, 6))
+            plt.hist2d(x, y, bins=200)
+            plt.ylabel("PCEnergy")
+            plt.xlabel("SX3 Energy Detected (MeV)")
+            plt.title(f"{data['particle']} Energy Propagation Difference vs SX3 Energy Detected")
+            plt.colorbar(label="Counts")
+            plt.tight_layout()
+
+            safe_name = re.sub(r"[^0-9A-Za-z_-]", "_", data["particle"])
+            plt.savefig(
+                f"{outdir}/{safe_name}_Eprop_vs_Edet.png",
+                dpi=300
+            )
+            plt.show()
     
         plt.figure(figsize=(8, 6))
     
@@ -1592,29 +1631,26 @@ class MyInteractiveApp(cmd.Cmd):
     
         all_Esx3 = []
         all_Eprop = []
-        all_Edet =  []
+        all_Edet = []
     
         for data in datasets:
-    
+            mask = (
-            mask = data["Esx3"] > 1
+                ~np.isnan(data["Edet"])
-    
+                & (data["Edet"] > 0)
-            thetab = np.deg2rad(
+                & ~np.isnan(data["Eprop"])
-                data["thetab"][mask]
+                & ~np.isnan(data["thetab"])
             )
-    
+            thetab = data["thetab"][mask]
             all_Esx3.append(
                 data["Esx3"][mask]
             )
-    
             all_Eprop.append(
                 data["Eprop"][mask]
-                * np.sin(thetab)
+                * thetab
-                * 3
             )
-            
             all_Edet.append(
                 data["Edet"][mask]
-                * np.sin(thetab))
+            )
     
         combined_Esx3 = np.concatenate(all_Esx3)
         combined_Eprop = np.concatenate(all_Eprop)
@@ -1660,6 +1696,7 @@ class MyInteractiveApp(cmd.Cmd):
         )
         plt.xlabel("SX3 Energy Deposition (MeV)")
         plt.ylabel("PCEnergy x Sin(theta)")
+        plt.title("PCEnergy vs Energy Deposition in sx3")
         plt.colorbar(label="Counts")
         plt.tight_layout()
         plt.savefig(
@@ -1705,8 +1742,6 @@ class MyInteractiveApp(cmd.Cmd):
             f"Completed plotting "
             f"{len(datasets)} datasets."
         )
-            
-#exec(open("PCEnergyAnalysis.py").read())
 
 if __name__ == "__main__":
     MyInteractiveApp().cmdloop()