recoil recon

Armory/Archives/Claude_recoil_script

@@ -0,0 +1,35 @@
+          /*
+          // Missing-mass reconstruction for 27Al + alpha -> alpha + 27Al*.
+          // The beam travels along +z; the target alpha is at rest.
+          if (std::isfinite(EBeam_track) && std::isfinite(Tb_recon) &&
+              EBeam_track >= 0.0 && Tb_recon >= 0.0 &&
+              projectileDirection.Mag2() > 0.0) {
+            const double beamMass = al27Iso.Mass;
+            const double alphaMass = alphaIso.Mass;
+            const double beamMomentum = std::sqrt(EBeam_track *
+                                                  (EBeam_track + 2.0 * beamMass));
+            const double alphaMomentum = std::sqrt(Tb_recon *
+                                                   (Tb_recon + 2.0 * alphaMass));
+
+            TLorentzVector beam4;
+            beam4.SetPxPyPzE(0.0, 0.0, beamMomentum, beamMass + EBeam_track);
+
+            TLorentzVector target4;
+            target4.SetPxPyPzE(0.0, 0.0, 0.0, alphaMass);
+
+            TLorentzVector outgoingAlpha4;
+            outgoingAlpha4.SetVectM(projectileDirection.Unit() * alphaMomentum,
+                                    alphaMass);
+
+            const TLorentzVector missingAl27 = beam4 + target4 - outgoingAlpha4;
+            const double missingMass2 = missingAl27.M2();
+            Ex_kin = missingMass2 >= 0.0
+                         ? std::sqrt(missingMass2) - beamMass
+                         : TMath::QuietNaN();
+            EBeam_kin = EBeam_track;
+            vX_kin = vX_track;
+            vY_kin = vY_track;
+            vZ_kin = vZ_track;
+          } else {
+            Ex_kin = TMath::QuietNaN();
+          }*/

Armory/Eloss_and_straggle.py

@@ -4,7 +4,7 @@ import matplotlib.pyplot as plt
 
 # --- 1. Constants ---
 P_TORR = 250
-P_CO2 = 3
+P_CO2 = 0
 if(len(sys.argv)==3):
     P_TORR=int(sys.argv[1])
     P_CO2 = int(sys.argv[2])
@@ -27,14 +27,13 @@ gas_mix = catima.Material(material_def)
 gas_mix.density(rho_g_cm3)
 
 # Thickness step settings
-step_mg_cm2 = 0.001           # 1 ug/cm2 steps as per your example -- kept fine for
+step_mg_cm2 = 0.001
-                               # numerical accuracy of the dedx integration itself.
 step_g_cm2 = step_mg_cm2 / 1000.0
 max_steps = 1000000000             # Adjust based on how far you want to track
 
-coarse_step_cm = 0.2   # row spacing over most of the track
+coarse_step_cm = 0.2  
-fine_step_cm = 0.03     # row spacing near the Bragg peak
+fine_step_cm = 0.03     
-fine_zone_frac = 0.085   # fraction of the *total* range treated as "near the peak"
+fine_zone_frac = 0.085  
 
 # Set relative integration tolerance (lower = higher precision, slower calculation)
 catima.Config.epsrel = 1e-6
@@ -43,7 +42,10 @@ catima.Config.epsabs = 1e-9
 
 def generate_lookup(z, mass_u, e_start_mev, label):
     filename = f"Eloss/E_vs_x_{label}.dat"
-    header = f"Energy(MeV) \tmg/cm2 \tcm\tEin\tsigE\tsigA\tsigR\tsigX\tcov\ttof\tsp\nStarting Energy: {e_start_mev} MeV"
+    header = (
+        "Distance_cm\tEnergy_MeV\tSigma_E_MeV\tSigma_x_cm\tSigma_a_cm\t"
+        "Thickness_mg_cm2\tEin_MeV\tSigma_r\tCov\tToF\tStoppingPower"
+    )
 
     projectile = catima.Projectile(mass_u, z)
     e_u_init = e_start_mev / mass_u
@@ -69,9 +71,9 @@ def generate_lookup(z, mass_u, e_start_mev, label):
         # Save checkpoint row
         e_total_out = result.Eout * mass_u
         output.append([
-            e_total_out, thickness_g_cm2 * 1000.0, current_dist_cm,
+            current_dist_cm, e_total_out, result.sigma_E, result.sigma_x,
-            result.Ein, result.sigma_E, result.sigma_a, result.sigma_r,
+            result.sigma_a, thickness_g_cm2 * 1000.0, result.Ein,
-            result.sigma_x, result.cov, result.tof, result.sp
+            result.sigma_r, result.cov, result.tof, result.sp
         ])
 
         if result.Eout == 0:
@@ -81,11 +83,18 @@ def generate_lookup(z, mass_u, e_start_mev, label):
         step = fine_step_cm if current_dist_cm >= fine_zone_start_cm else coarse_step_cm
         current_dist_cm += step
 
-    np.savetxt(filename, output, fmt='%.6f', delimiter='\t', header=header)
+    data = np.asarray(output)
+    data = data[np.argsort(data[:, 1])]
+    np.savetxt(
+        filename, data, fmt='%.6f', delimiter='\t', header=header, comments=""
+    )
     print(f"Lookup table created: {filename} ({len(output)} rows)")
 
-    data = np.array(output)
+    # Plot in increasing distance order even though the saved table is energy-sorted.
-    energy_mev, dist_cm, sigma_e, tof = data[:, 0], data[:, 2], data[:, 4], data[:, 9]
+    plot_data = data[::-1]
+    dist_cm, energy_mev, sigma_e, tof = (
+        plot_data[:, 0], plot_data[:, 1], plot_data[:, 2], plot_data[:, 9]
+    )
 
     fig, (ax_e, ax_sig, ax_tof) = plt.subplots(3, 1, figsize=(8, 12), sharex=True)
 
@@ -114,7 +123,7 @@ def generate_lookup(z, mass_u, e_start_mev, label):
 generate_lookup(1, 1.0078, 30, "proton")
 generate_lookup(1, 2.01355, 30, "deuteron")
 generate_lookup(2, 4.0026, 50, "alpha")
-generate_lookup(13,26.9815, 80, "aluminum")
+generate_lookup(13,26.9815, 56.1, "aluminum")
 #generate_lookup(9,17.0021, 70, "fluorine")
 #generate_lookup(8,15.9949, 70, "oxygen")
 

Armory/EnergyLoss.py

@@ -0,0 +1,597 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Aug 20 12:46:47 2026
+
+@author: jamesszalkie
+"""
+
+
+import gc
+import numpy as np
+import pandas as pd
+from scipy.interpolate import interp1d
+import uproot
+import pycatima as catima
+from scipy.integrate import cumulative_trapezoid
+import matplotlib.pyplot as plt
+import cmd
+import shlex
+import textwrap
+import os
+import periodictable as pt
+import re
+from matplotlib.colors import LinearSegmentedColormap
+from mpl_toolkits.mplot3d import Axes3D
+import nbformat as nbf
+import shutil
+
+alpha_data = [2, 4.0015, 40, "alpha"]
+proton_data = [1, 1.0073, 20, "proton"]
+deuteron_data = [1, 2.014102, 30, "deuteron"]
+
+interp_cache = {}
+
+particles = {
+    "alpha": alpha_data,
+    "proton": proton_data,
+    "deuteron": deuteron_data
+}
+
+
+def get_loss_table_path(medium, particle_label):
+    script_dir = os.path.dirname(os.path.abspath(__file__))
+    return os.path.join(script_dir, f"{medium}Loss", f"E_vs_x_{particle_label}.dat")
+
+
+def normalize_particle_label(particle):
+    """Return the canonical table label for a particle or isotope name."""
+    if isinstance(particle, str):
+        try:
+            _, _, _, label = resolve_particle(particle)
+            return label
+        except Exception:
+            return particle.strip()
+    return str(particle)
+
+
+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,
+    emax_mev,
+    medium,
+    label,
+    npoints,
+    P_TORR,
+    TEMP_K
+):
+    """
+    Builds energy vs depth table using Catima stopping powers.
+    Output:
+        x [cm], E [MeV]
+    """
+    R = 8.3144  # J/mol/K
+
+    p_pa = P_TORR * 133.322
+    molar_density = p_pa / (R * TEMP_K)  # mol/m^3
+
+    # Medium definition
+    if medium == "He":
+
+        m_he = 4.0026
+        m_c = 12.0000
+        m_o = 15.9949
+        """
+        material_def = [
+            (m_he, 2, 0.96),
+            (m_c, 6, 0.04),
+            (m_o, 8, 0.08)
+        ] """
+        material_def = [(m_he, 2, 1.0)]
+
+        m_mix_avg = 0.96 * m_he + 0.04 * (m_c + 2 * m_o)
+        #m_mix_avg = 1.0 * m_he
+        rho_g_cm3 = (molar_density * m_mix_avg) / 1e6
+        
+        gas = catima.Material(material_def)
+        gas.density(rho_g_cm3)
+
+    elif medium == "Si":
+
+        m_si = 28.084
+        rho_g_cm3 = 2.33
+
+        gas = catima.Material([(m_si, 14, 1.0)])
+        gas.density(rho_g_cm3)
+
+    elif medium == "kapton":
+        # Kapton (C22H10N2O5)
+        # Density: 1.42 g/cm3
+        m_h = 1.0078
+        m_n = 14.0067
+        m_c = 12.0000
+        m_o = 15.9949
+        kapton_molar_mass = 22*m_c + 10*m_h + 2*m_n + 5*m_o
+        
+        material_def = [
+            (m_c, 6, 22/kapton_molar_mass * kapton_molar_mass / m_c),
+            (m_h, 1, 10/kapton_molar_mass * kapton_molar_mass / m_h),
+            (m_n, 7, 2/kapton_molar_mass * kapton_molar_mass / m_n),
+            (m_o, 8, 5/kapton_molar_mass * kapton_molar_mass / m_o)
+        ]
+        
+        rho_g_cm3 = 1.42
+        gas = catima.Material(material_def)
+        gas.density(rho_g_cm3)
+
+    elif medium == "mylar":
+        # Mylar (C10H8O4, polyethylene terephthalate)
+        # Density: 1.39 g/cm3
+        m_h = 1.0078
+        m_c = 12.0000
+        m_o = 15.9949
+        mylar_molar_mass = 10*m_c + 8*m_h + 4*m_o
+        
+        material_def = [
+            (m_c, 6, 10/mylar_molar_mass * mylar_molar_mass / m_c),
+            (m_h, 1, 8/mylar_molar_mass * mylar_molar_mass / m_h),
+            (m_o, 8, 4/mylar_molar_mass * mylar_molar_mass / m_o)
+        ]
+        
+        rho_g_cm3 = 1.39
+        gas = catima.Material(material_def)
+        gas.density(rho_g_cm3)
+
+    else:
+        raise ValueError("Unsupported medium")
+
+    print(f"[INFO] density = {rho_g_cm3:.3e} g/cm^3")
+
+    projectile = catima.Projectile(mass_u, z)
+    E = np.linspace(0.1, emax_mev, npoints)
+    S_mass = np.zeros_like(E)
+    for i, energy in enumerate(E):
+        projectile.T(energy / mass_u)
+        S_mass[i] = catima.dedx(projectile, gas)  # MeV/(g/cm^2)
+
+    S_linear = S_mass * rho_g_cm3
+    invS = 1.0 / np.clip(S_linear, 1e-30, None)
+    x = cumulative_trapezoid(invS[::-1], E[::-1], initial=0)
+    x = x[::-1]
+    x = -x
+
+    # Energy and range straggling at each cumulative depth, starting from emax_mev
+    e_u_init = emax_mev / mass_u
+    sigma_E = np.zeros_like(E)
+    sigma_x = np.zeros_like(E)
+    sigma_a = np.zeros_like(E)
+    for i, depth_cm in enumerate(x):
+        gas.density(rho_g_cm3).thickness(depth_cm * rho_g_cm3)
+        projectile.T(e_u_init)
+        result = catima.calculate(projectile, gas)
+        sigma_E[i] = result.sigma_E
+        sigma_x[i] = result.sigma_x
+        sigma_a[i] = result.sigma_a
+
+    df = pd.DataFrame({
+        "Distance_cm": x,
+        "Energy_MeV": E,
+        "Sigma_E_MeV": sigma_E,
+        "Sigma_x_cm": sigma_x,
+        "Sigma_a_cm": sigma_a
+    })
+
+    outfile = get_loss_table_path(medium, label)
+    os.makedirs(os.path.dirname(outfile), exist_ok=True)
+    df.to_csv(outfile, sep="\t", index=False)
+
+    print(f"[INFO] saved: {outfile}")
+
+    interp_cache.pop((label.lower(), medium.lower()), None)
+    gc.collect()
+
+    return x, E
+
+   #Generate energy loss tables from file
+def load_table(filename):
+    """
+    Load table with columns:
+        x(cm)   E(MeV)   Sigma_E(MeV) [optional]   Sigma_x(cm) [optional]
+
+    Returns:
+        x_array, E_array, sigma_E_array, sigma_x_array (None if column absent)
+    """
+
+    data = pd.read_csv(
+        filename,
+        sep=r'\s+',
+        comment="#",
+        header=None,
+        skiprows=1
+    )
+
+    x = data.iloc[:, 0].to_numpy()
+    E = data.iloc[:, 1].to_numpy()
+    sigma_E = data.iloc[:, 2].to_numpy() if data.shape[1] > 2 else None
+    sigma_x = data.iloc[:, 3].to_numpy() if data.shape[1] > 3 else None
+
+    return x, E, sigma_E, sigma_x
+
+def get_interpolators(particle, medium):
+
+    canonical_particle = normalize_particle_label(particle)
+    cache_key = (canonical_particle.lower(), medium.lower())
+
+    if cache_key in interp_cache:
+        return interp_cache[cache_key]
+
+    candidate_labels = [canonical_particle]
+    if isinstance(particle, str):
+        raw_particle = particle.strip()
+        candidate_labels.extend([
+            raw_particle,
+            canonical_particle.replace("-", "")
+        ])
+
+    filename = None
+    seen_candidates = set()
+    for label in candidate_labels:
+        normalized_label = label.strip()
+        if not normalized_label:
+            continue
+        key = normalized_label.lower()
+        if key in seen_candidates:
+            continue
+        seen_candidates.add(key)
+
+        candidate_filename = get_loss_table_path(medium, normalized_label)
+        if os.path.exists(candidate_filename):
+            filename = candidate_filename
+            break
+
+    if filename is None:
+        filename = get_loss_table_path(medium, canonical_particle)
+
+    x, E, sigma_E, sigma_x = load_table(filename)
+
+    E_of_x = interp1d(
+        x,
+        E,
+        bounds_error=False,
+        fill_value="extrapolate"
+    )
+
+    x_of_E = interp1d(
+        E[::-1],
+        x[::-1],
+        bounds_error=False,
+        fill_value="extrapolate"
+    )
+
+    if sigma_E is not None:
+        sigma_of_x = interp1d(
+            x,
+            sigma_E,
+            bounds_error=False,
+            fill_value="extrapolate"
+        )
+    else:
+        sigma_of_x = None
+
+    if sigma_x is not None:
+        sigma_x_of_x = interp1d(
+            x,
+            sigma_x,
+            bounds_error=False,
+            fill_value="extrapolate"
+        )
+    else:
+        sigma_x_of_x = None
+
+    interp_cache[cache_key] = (E_of_x, x_of_E, sigma_of_x, sigma_x_of_x)
+
+    return E_of_x, x_of_E, sigma_of_x, sigma_x_of_x
+
+def energy_loss(particle, medium, Ei, dl):
+
+    E_of_x, x_of_E, _, _ = get_interpolators(particle, medium)
+
+    xi = x_of_E(Ei)
+
+    xf = xi + dl
+
+    xmax = E_of_x.x.max()  # maximum tabulated range
+
+    Ef = E_of_x(xf)
+
+    Ef = np.where(xf >= xmax, 0.0, Ef)
+
+    return np.maximum(Ef, 0.0)
+    
+def energy_reconstruction(particle, medium, Ef, dl):
+
+    E_of_x, x_of_E, _, _ = get_interpolators(particle, medium)
+
+    xf = x_of_E(Ef)
+
+    xi = xf - dl
+
+    Ei = E_of_x(xi)
+
+    return np.maximum(Ei, 0.0)
+
+def energy_distance(particle, medium, Ei, Ef):
+
+    _, x_of_E, _, _ = get_interpolators(particle, medium)
+
+    xi = x_of_E(Ei)
+
+    xf = x_of_E(Ef)
+
+    return np.abs(xf - xi)
+
+def resolve_particle(name):
+    name = name.lower().strip().rstrip("s")
+    if name in particles:
+        return particles[name]
+    match = re.match(r"([a-zA-Z]+)[-\s]?(\d+)$", name)
+    if match:
+        element_symbol = match.group(1).capitalize()
+        A = int(match.group(2))
+
+        try:
+            element = pt.elements.symbol(element_symbol)
+            isotope = element[A]
+
+            return (
+                isotope.number,   # Z
+                isotope.mass,     # mass in u
+                30.0,
+                f"{element_symbol}-{A}"
+            )
+        except Exception:
+            raise ValueError(f"Unknown isotope: {name}")
+
+    match = re.match(r"(\d+)[-\s]?([a-zA-Z]+)$", name)
+    if match:
+        A = int(match.group(1))
+        element_symbol = match.group(2).capitalize()
+
+        try:
+            element = pt.elements.symbol(element_symbol)
+            isotope = element[A]
+
+            return (
+                isotope.number,   # Z
+                isotope.mass,     # mass in u
+                30.0,
+                f"{element_symbol}-{A}"
+            )
+        except Exception:
+            raise ValueError(f"Unknown isotope: {name}")
+    try:
+        elem = pt.elements.symbol(name.capitalize())
+        return elem.number, elem.mass, 30.0, name
+    except Exception:
+        raise ValueError(f"Unknown particle/isotope: {name}")
+        
+        
+class MyInteractiveApp(cmd.Cmd):
+    def __init__(self):
+        super().__init__()
+        # Initial value set when the script starts
+        self.T = 293.15
+        self.P = 250#379
+        self.temp_particle = [0, 0.0, 0.0, ""]
+        print("-" * 30)
+        print("INTERACTIVE SHELL STARTED")
+        self.print_params()
+        print("Type 'help' for commands.")
+        print("Type 'exit' to end program")
+        print("-" * 30)
+        
+    def print_params(self):
+        """Helper method to display current state"""
+        print(f"Current Parameters: T={self.T} K, P={self.P} Torr")
+        
+    #intro = "Interactive Shell Started. Type 'help' to see commands."
+    prompt = ">> "
+    
+    def default(self, line):
+        # Check if the command starts with our multi-word phrase
+        if line.startswith("make table "):
+            # Extract everything after "make table "
+            args = line[len("make table "):].strip()
+            self.do_make_table(args)
+        elif line.startswith("set t") or line.startswith("Set T") or line.startswith("set T") or line.startswith("Set t"):
+            args = line[len("set t "):].strip()
+            self.do_set_T(args)
+        elif line.startswith("set p") or line.startswith("Set P") or line.startswith("set P") or line.startswith("Set p"):
+            args = line[len("set p "):].strip()
+            self.do_set_P(args)
+        elif line.startswith("energy loss") or line.startswith("Energy Loss") or line.startswith("Energy loss"):
+            args = line[len("energy loss "):].strip()
+            self.do_energy_loss(args)
+        elif line.startswith("energy reconstruction") or line.startswith("Energy Reconstruction") or line.startswith("Energy reconstruction"):
+            args = line[len("energy reconstruction "):].strip()
+            self.do_energy_reconstruction(args)
+        elif line.startswith("energy distance") or line.startswith("Energy Distance") or line.startswith("Energy distance"):
+            args = line[len("energy distance "):].strip()
+            self.do_energy_distance(args)
+        else:
+            print(f"*** Unknown syntax: {line}")
+            
+    def do_exit(self, arg):
+        """Exits the application."""
+        print("Closing application...")
+        return True  # Returning True stops the cmdloop()
+            
+    def do_T(self, arg):
+        """Print value of T"""
+        print(self.T)
+        
+    def do_P(self, arg):
+        """Print value of P (pressure)"""
+        print(self.P)
+        
+    def do_set_T(self, arg):
+        """Changes the value of T. Usage: set_t 300"""
+        try:
+            self.T = float(arg)
+            print(f"T has been updated to {self.T}")
+        except ValueError:
+            print("Please enter a valid number for T.")
+            
+    def do_set_P(self, arg):
+        """Changes the value of P in Torr. Usage: set_ 400"""
+        try:
+            self.P = float(arg)
+            print(f"P has been updated to {self.P}")
+        except ValueError:
+            print("Please enter a valid number for P.")
+            
+    
+    def do_make_table(self, arg):
+        """Create E vs X tables for particle, or isotopes
+        Ex: >> make table proton <max energy (optional) >
+        Ex: >> make table Co60 <max energy (optional) >
+        Ex: >> make table N17 <max energy (optional) >"""
+        try:
+            args = shlex.split(arg)
+    
+            if not args:
+                print("Please enter desired reaction particle")
+                return
+    
+            name = args[0]
+    
+            if len(args) > 1:
+                emax_mev = float(args[1])
+            else:
+                emax_mev = None
+            
+            if len(args) > 2:
+                medium = args[2]
+            else:
+                medium = "He"
+    
+            z, mass_u, default_emax, label = resolve_particle(name)
+    
+            if emax_mev is None:
+                emax_mev = default_emax
+    
+            x, E = make_E_vs_x(
+                z,
+                mass_u,
+                emax_mev,
+                medium,
+                label,
+                100000,
+                self.P,
+                self.T
+            )
+    
+            plt.figure(figsize=(8,6))
+            plt.plot(x, E)
+            
+            plt.xlabel("Distance (cm)")
+            plt.ylabel("Energy (MeV)")
+            plt.title(f"Energy Loss Curve {label.capitalize()} {medium}")
+            plt.grid(True)
+            textstr = f"T = {self.T:.2f} K\nP = {self.P:.2f} Torr"
+            
+            plt.gca().text(
+                0.02, 0.02,
+                textstr,
+                transform=plt.gca().transAxes,
+                fontsize=10,
+                verticalalignment='bottom',
+                bbox=dict(boxstyle="round", facecolor="white", alpha=0.7)
+            )
+            
+            plt.tight_layout()
+            filename = f"{medium}Loss/Energy_Loss_Curve_{label}.png"
+            plt.savefig(filename, dpi=300, bbox_inches="tight")
+            plt.show()
+            
+            print(f"Saved plot: {filename}")
+    
+        except Exception as e:
+            print(f"Error in make_table: {e}")
+            
+    def do_energy_loss(self, arg):
+        """Find a final energy given an initial energy and distance travelled
+        Ex: >> energy loss <particle> <medium> <initial energy MeV> <distance travelled cm>"""
+        args = shlex.split(arg)
+        try:
+            particle = args[0]
+            medium = args[1]
+            Ei = float(args[2])
+            dl = float(args[3])
+            #Ei_offset = Ei * 1.1
+            #table_specs = f"{particle} {Ei_offset}"
+            try:
+                #self.do_make_table(table_specs)
+                Ef = energy_loss(particle, medium, Ei, dl)
+                print(f"\nFinal energy: {Ef:.6f} MeV\n")
+            except:
+                return
+        except IndexError:
+            print("Please input particle, initial energy, and distance travelled")
+            
+    def do_energy_reconstruction(self, arg):
+        """Find a vertex energy given an final energy and distance travelled
+        Ex: >> energy reconstruction <particle> <medium> <final energy MeV> <distance travelled cm>"""
+        args = shlex.split(arg)
+        try:
+            particle = args[0]
+            medium = args[1]
+            Ef = float(args[2])
+            dl = float(args[3])
+            try:
+                Ei = energy_reconstruction(particle, medium, Ef, dl)
+                if Ei > 0:
+                    print(f"\nInitial energy: {Ei:.6f} MeV\n")
+                else:
+                    print("Error: remake table with larger value, fallen off map")
+            except:
+                print("Particle energy table not made yet, please do so using 'make table'")
+        except IndexError:
+            print("Please input particle, final energy from detector, and distance travelled")
+            
+            
+    def do_energy_distance(self, arg):
+        """Find a distance travelled given an initial and final energy
+        Ex: >> energy distance <particle> <medium> <initial energy> <final energy MeV>"""
+        args = shlex.split(arg)
+        try:
+            particle = args[0]
+            medium = args[1]
+            Ei = float(args[2])
+            Ef = float(args[3])
+            dE = Ei - Ef
+            try:
+                dl = energy_distance(particle, medium, Ei, Ef)
+                if dl > 0:
+                    print(f"\nChange in energy: {dE:.6f} MeV")
+                    print(f"Distance travelled: {dl:.6f} cm\n")
+                    print(f"{(dl * 10):.6f} mm")
+                else:
+                    print("Error: remake table with larger value, fallen off map")
+            except:
+                print("Particle energy table not made yet, please do so using 'make table'")
+        except IndexError:
+            print("Please input particle, final energy from detector, and distance travelled")
+            
+            
+
+
+if __name__ == "__main__":
+    MyInteractiveApp().cmdloop()

Armory/anasenMS.cpp

@@ -70,7 +70,7 @@ bool IsDeadCathode(int id){
 
 bool IsDeadSX3(int id){ 
     static std::set<int> dead = {
-      //0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23
+      0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23
     }; // add dead SX3 IDs here, 0-23 1,7,9,3
     return dead.count(id);
 }
@@ -91,7 +91,7 @@ bool IsDeadSX3FrontDnChannel(int sx3ID, int chDn){
 
 bool IsDeadSX3BackChannel(int sx3ID, int chBk){
   static std::set<std::pair<int, int>> dead = {
-    //{1, 10}
+    {1, 10}
     // {sx3ID, back-channel} 
   };
   return dead.count({sx3ID, chBk});
@@ -160,22 +160,26 @@ int main(int argc, char **argv){
   //To set beam energy loss, use energy loss app, and create table with target isotope, set Initial beam energy as max energy 
   transfer.SetA(27, 13, 0);            // 22Mg projectile
   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);               // 30Si* heavy product
+  transfer.SetB(27, 13);               // 30Si* heavy product
   double beamE = 56.1; //56.1
   const ReactionConfig reactionConfig = transfer.GetRectionConfig();
   const double beamA = reactionConfig.beamA;                   // mass number of 14N beam
 
   // Excited state lists (projectile and heavy-product excitation states)
   std::vector<float> ExAList = {0};              // Beam excited energy
-  //std::vector<float> ExList = {0, 2.235, 3.498, 3.774, 4.831, 4.809, 5.279, 5.614, 5.951, 6.550};             // Heavy product excited energy
+  //std::vector<float> ExList = {0, 2.235, 3.498, 3.774, 4.831, 4.809, 5.279, 5.614, 5.951, 6.550};  // a/p aluminum
-  std::vector<float> ExList = {2.235};  
+  std::vector<float> ExList = {0.84, 1.01, 2.21, 2.73, 2.98, 3}; 
 
   const int kMBeam = reactionConfig.beamA; // mass number of beam
   const int kMTarget = reactionConfig.targetA; // mass number of target
   const int kMLight = reactionConfig.recoilLightA; // mass number of light ejectile
   const int kMHeavy = reactionConfig.recoilHeavyA; // mass number of heavy product
 
+  // Nuclear masses are in MeV/c^2; cache them rather than rereading mass20.txt per event.
+  const Isotope al27Iso(27, 13);
+  const Isotope alphaIso(4, 2);
+
   //const int kZBeam = reactionConfig.beamZ; // atomic number of beam
   //const int kZTarget = reactionConfig.targetZ; // atomic number of target
   //const int kZLight = reactionConfig.recoilLightZ; // atomic number of light ejectile
@@ -450,6 +454,9 @@ int main(int argc, char **argv){
   tree1->Branch("EBeam_kin", &EBeam_kin, "EBeam_kin/D");
   tree1->Branch("Ex_kin", &Ex_kin, "Ex_kin/D");
 
+  double recoilEnergy;
+  tree1->Branch("recoilEnergy", &recoilEnergy, "recoilEnergy/D");
+
   // reconstructed vertex Z from PW fit
   double z0;
   tree1->Branch("z0",     &z0, "reconstucted_Z/D");
@@ -538,8 +545,8 @@ int main(int argc, char **argv){
     double vertexRangeX = std::sqrt((vertexXRange[1] * vertexXRange[1]) + (sigmaX_mm * sigmaX_mm));
     double vertexRangeY = std::sqrt((vertexYRange[1] * vertexYRange[1]) + (sigmaX_mm * sigmaX_mm));
     //std::cout << "vertexRangeX: " << vertexRangeX << ", vertexRangeY: " << vertexRangeY << std::endl;
-    vertexX = gRandom->Gaus(0, vertexRangeX); // mean and standard deviation
+    vertexX = gRandom->Gaus(0, 0);//vertexRangeX); // mean and standard deviation
-    vertexY = gRandom->Gaus(0, vertexRangeY); // mean and standard deviation
+    vertexY = gRandom->Gaus(0, 0);//vertexRangeY); // mean and standard deviation
 
     //vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() +  vertexZRange[0];
 
@@ -690,15 +697,15 @@ int main(int argc, char **argv){
       //TVector3 hitSigma = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
 
       //smudging
-      TVector3 hitSigma = sx3->GetHitPosWithSigma(10, 10); // args are sigmaSX3_W, sigmaSX3_L, currently set to 0 for no smearing
+      TVector3 hitSigma = sx3->GetHitPosWithSigma(0, 0); // args are sigmaSX3_W, sigmaSX3_L, currently set to 0 for no smearing
 
       double hitX = hitSigma.X();
       double hitY = hitSigma.Y();
       double hitZ = hitSigma.Z();
 
-      double aX_sigma = gRandom->Gaus(aX, 10);
+      double aX_sigma = gRandom->Gaus(aX, 0);
-      double aY_sigma = gRandom->Gaus(aY, 10);
+      double aY_sigma = gRandom->Gaus(aY, 0);
-      double aZ_sigma = gRandom->Gaus(aZ, 10);
+      double aZ_sigma = gRandom->Gaus(aZ, 0);
                                                                                              
       double distance_A;
       Eanode = CalculateEnergyLoss(vertexX, vertexY, vertexZ,
@@ -749,15 +756,40 @@ int main(int argc, char **argv){
                                             b, "He", Esx3,
                                             distance_sx3);
 
-        Kinematics apkin_27Al(26.981538408,4.00260325413,1.00782503224,29.973770136,EBeam_track/26.981538408); //m3 is proton
+        const TVector3 projectileDirection = TVector3(hitX - aX_sigma, hitY - aY_sigma, hitZ - aZ_sigma);
-        
+        theta_recon = projectileDirection.Theta() * 180.0 / TMath::Pi();
-        //kinematic reconstruction
-        theta_recon = TVector3(hitX - aX_sigma, hitY - aY_sigma, hitZ - aZ_sigma).Theta() * 180.0 / TMath::Pi();
-        EBeam_kin = apkin_27Al.getEbeam_givenQ(Tb_recon, Ex, theta_recon); //args are (reaction energy, Ex, thetab) still tracking
-        vZ_kin = elossBeamInverse->Eval(EBeam_kin) * 10 - 454.3;
 
-        Ex_kin = apkin_27Al.getExc(Tb_recon, theta_recon);
+        if (reactionConfig.recoilLightA == 1) {
-        /*Checklist: anode smudge, sx3 smudge, beam position off axis, beam angle*/
+          //Kinematics apkin_27Al(26.981538408,4.00260325413,1.00782503224,29.973770136,EBeam_track/26.981538408); //m3 is proton
+          Kinematics apkin_27Al(26.981538408,4.00260325413,4.00260325413,26.981538408,EBeam_track/26.981538408); //args are (m1, m2, m3, m4, Ebeam/m1)
+          
+          //kinematic reconstruction
+          EBeam_kin = apkin_27Al.getEbeam_givenQ(Tb_recon, Ex, theta_recon); //args are (reaction energy, Ex, thetab) still tracking
+          vZ_kin = elossBeamInverse->Eval(EBeam_kin) * 10 - 454.3;
+
+          Ex_kin = apkin_27Al.getExc(Tb_recon, theta_recon);
+          /*Checklist: anode smudge, sx3 smudge, beam position off axis, beam angle*/
+
+        } else if (reactionConfig.recoilLightA == 4) {
+            double beamMomentum = std::sqrt((EBeam_track + al27Iso.Mass) * (EBeam_track + al27Iso.Mass) - al27Iso.Mass * al27Iso.Mass);
+            TLorentzVector beam4(0.0, 0.0, beamMomentum, EBeam_track + al27Iso.Mass);
+
+            TLorentzVector helium4(0, 0, 0, alphaIso.Mass);
+
+            double alphaEnergy = Tb_recon + alphaIso.Mass;
+            double alphaMomentum = std::sqrt(alphaEnergy * alphaEnergy - alphaIso.Mass * alphaIso.Mass);
+            TLorentzVector alpha4(alphaMomentum * std::sin(theta_recon * TMath::Pi() / 180.0), 
+                                    0, 
+                                    alphaMomentum * std::cos(theta_recon * TMath::Pi() / 180.0), 
+                                    alphaEnergy);
+            TLorentzVector recoil4 = beam4 + helium4 - alpha4;
+
+            double recoilMass = recoil4.M();
+
+            Ex_kin = recoilMass - al27Iso.Mass;
+            vZ_kin = vZ_track;
+
+        }
 
 
       } else {

Armory/testMacro.C

@@ -4,44 +4,18 @@
 	TTree *tree1 = (TTree*)(file0->Get("tree1"));
 	new TBrowser();
 
-	/*
-	new TCanvas();
-	TH2F *h3 = new TH2F(
-		"hExReconVsVZ",
-		"Ex_recon vs vZ;vZ (mm);Excitation energy (MeV)",
-		100, 0, 0,
-		100, -1, 0
-	);
-	tree1->Draw("Ex_recon:vZ>>hExReconVsVZ", "sx3ID >= 0", "colz");
-
 	TH2F *h4 = new TH2F(
 		"hExVsVZ",
 		"",
 		100, 0, 0,
 		100, -1, 0
 	);
-	tree1->Draw("Ex:vZ>>hExVsVZ", "sx3ID >= 0", "box same");*/
+	tree1->Draw("Ex:vZ>>hExVsVZ", "sx3ID >= 0", "box same");
-/*
-	new TCanvas();
-	tree1->Draw("originalEnergy:Tb", "sx3ID >=0", "");
-	TF1 *f_diagonal = new TF1("f_diagonal", "x", 0, 200); // assuming the range for Tb is 0 to 200
-	f_diagonal->SetLineColor(kRed);
-	f_diagonal->SetLineStyle(2);
-	f_diagonal->Draw("same");*/
-/*
-	TH2F *h5 = new TH2F("Excitation Energy vs vZ", "Excitation Energy vs vZ;vZ (mm);Excitation energy (MeV)", 300, 0, 0, 300, 0, 0); //arguments are (name, title, nbinsX, xlow, xup, nbinsY, ylow, yup)
-	tree1->Draw("Ex_recon:vZ >> hExReconVsVZ", "sx3ID >= 0", "colz");
-	tree1->Draw("Ex:vZ >> h5", "sx3ID >= 0 && EPC >= 0.03", "same");
-*/
 
 	new TCanvas("c1", "Canvas", 900, 600);
-	tree1->Draw("Ex_kin:vZ_kin>>h2(800,-500,300,110,-1,10)", "sx3ID >=0 && !TMath::IsNaN(vZ_kin)", "colz");
+	tree1->Draw("Ex_kin:vZ_kin>>h2(800,-500,300,110,-1,10)", "sx3ID >=0 && !TMath::IsNaN(vZ_kin) && EPC > .02", "colz");
-	tree1->Draw("Ex:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ)", "same");
+	tree1->Draw("Ex:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ) && EPC > .02", "same");
 	c1->SaveAs("Ex_kin_vz.png");
-	/*
-	new TCanvas("c2", "Canvas", 900, 600);
-	tree1->Draw("theta_recon:thetab", "sx3ID >=0 && Esx3 > 0", "");
-	c2->SaveAs("theta_recon_vs_thetab.png");
 
 	new TCanvas("c3", "Canvas", 900, 600);
 	tree1->Draw("vZ_track:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ_track)", "");
@@ -66,13 +40,14 @@
 	TF1 *line_kin_kin = new TF1("line_beam_kin", "x", gPad->GetUxmin(), gPad->GetUxmax());
 	line_kin_kin->SetLineColor(kRed); line_kin_kin->Draw("SAME");
 	c6->SaveAs("Ebeam_kin_vs_beamEnergy.png");
-*/
+
+/*
 	new TCanvas("c7", "Canvas", 900, 600);
 	//c7->SetLogy();
 	tree1->Draw("Ex_kin", "sx3ID >= 0");
 	tree1->Draw("Ex", "sx3ID >= 0", "same");
 	tree1->GetHistogram()->SetLineColor(kRed);
 	c7->SaveAs("Ex_kin.png");
-
+*/
 
 }