From 5193da6abcb3dcc150bd638768e46eff3bf8b8b0 Mon Sep 17 00:00:00 2001 From: James Szalkie Date: Wed, 30 Sep 2026 14:07:31 -0400 Subject: [PATCH] recoil recon --- Armory/Archives/Claude_recoil_script | 35 ++ Armory/Eloss_and_straggle.py | 37 +- Armory/EnergyLoss.py | 597 +++++++++++++++++++++++++++ Armory/anasenMS.cpp | 72 +++- Armory/testMacro.C | 37 +- 5 files changed, 713 insertions(+), 65 deletions(-) create mode 100644 Armory/Archives/Claude_recoil_script create mode 100644 Armory/EnergyLoss.py diff --git a/Armory/Archives/Claude_recoil_script b/Armory/Archives/Claude_recoil_script new file mode 100644 index 0000000..83df58b --- /dev/null +++ b/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(); + }*/ diff --git a/Armory/Eloss_and_straggle.py b/Armory/Eloss_and_straggle.py index e5d0047..dffb678 100644 --- a/Armory/Eloss_and_straggle.py +++ b/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 - # numerical accuracy of the dedx integration itself. +step_mg_cm2 = 0.001 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 -fine_step_cm = 0.03 # row spacing near the Bragg peak -fine_zone_frac = 0.085 # fraction of the *total* range treated as "near the peak" +coarse_step_cm = 0.2 +fine_step_cm = 0.03 +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, - result.Ein, result.sigma_E, result.sigma_a, result.sigma_r, - result.sigma_x, result.cov, result.tof, result.sp + current_dist_cm, e_total_out, result.sigma_E, result.sigma_x, + result.sigma_a, thickness_g_cm2 * 1000.0, result.Ein, + 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) - energy_mev, dist_cm, sigma_e, tof = data[:, 0], data[:, 2], data[:, 4], data[:, 9] + # Plot in increasing distance order even though the saved table is energy-sorted. + 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") diff --git a/Armory/EnergyLoss.py b/Armory/EnergyLoss.py new file mode 100644 index 0000000..3ba09f5 --- /dev/null +++ b/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 + Ex: >> make table Co60 + Ex: >> make table N17 """ + 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 """ + 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 """ + 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 """ + 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() \ No newline at end of file diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 0eb6ff0..acca090 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -70,7 +70,7 @@ bool IsDeadCathode(int id){ bool IsDeadSX3(int id){ static std::set 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> 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(30, 14); // 30Si* heavy product + transfer.Setb(4, 2); // outgoing proton from the primary transfer + 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 ExAList = {0}; // Beam excited energy - //std::vector 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 ExList = {2.235}; + //std::vector 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 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 - vertexY = gRandom->Gaus(0, vertexRangeY); // mean and standard deviation + vertexX = gRandom->Gaus(0, 0);//vertexRangeX); // 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 aY_sigma = gRandom->Gaus(aY, 10); - double aZ_sigma = gRandom->Gaus(aZ, 10); + double aX_sigma = gRandom->Gaus(aX, 0); + double aY_sigma = gRandom->Gaus(aY, 0); + 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 - - //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; + const TVector3 projectileDirection = TVector3(hitX - aX_sigma, hitY - aY_sigma, hitZ - aZ_sigma); + theta_recon = projectileDirection.Theta() * 180.0 / TMath::Pi(); - Ex_kin = apkin_27Al.getExc(Tb_recon, theta_recon); - /*Checklist: anode smudge, sx3 smudge, beam position off axis, beam angle*/ + if (reactionConfig.recoilLightA == 1) { + //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 { diff --git a/Armory/testMacro.C b/Armory/testMacro.C index b8afd07..b306304 100644 --- a/Armory/testMacro.C +++ b/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");*/ -/* - 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"); -*/ + tree1->Draw("Ex:vZ>>hExVsVZ", "sx3ID >= 0", "box 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:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ)", "same"); + 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) && 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"); - +*/ }