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

switching files

Browse Source
This commit is contained in:
James Szalkie 2026-09-29 13:08:35 -04:00
parent 172594a90a
commit dac745fb5c
5 changed files with 634 additions and 100020 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
Armory/Eloss_and_straggle.py
View File

@ -3,7 +3,7 @@ import numpy as np, sys
import matplotlib.pyplot as plt
# --- 1. Constants ---
P_TORR = 250
P_TORR = 250#182
P_CO2 = 3
if(len(sys.argv)==3):
P_TORR=int(sys.argv[1])
@ -43,7 +43,7 @@ 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 = f"cm \tEnergy(MeV) \tmg/cm2 \tEin\tsigE\tsigA\tsigR\tsigX\tcov\ttof\tsp\nStarting Energy: {e_start_mev} MeV"
projectile = catima.Projectile(mass_u, z)
e_u_init = e_start_mev / mass_u
@ -69,7 +69,7 @@ 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, thickness_g_cm2 * 1000.0,
result.Ein, result.sigma_E, result.sigma_a, result.sigma_r,
result.sigma_x, result.cov, result.tof, result.sp
])
@ -107,8 +107,6 @@ def generate_lookup(z, mass_u, e_start_mev, label):
plt.tight_layout()
plt.show()
# --- 3. Run ---
# Format: generate_lookup(Z, mass_u, E_start_MeV, label)
generate_lookup(1, 1.0078, 30, "proton")

618
Armory/EnergyLoss.py Normal file
View File

@ -0,0 +1,618 @@
#!/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"Eloss", 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 energy-loss tables written by Eloss_and_straggle.py or make_E_vs_x.
Returns:
x_array, E_array, sigma_E_array, sigma_x_array (None if column absent)
"""
with open(filename, "r", encoding="utf-8") as table_file:
header = table_file.readline().lstrip("#").split()
data = pd.read_csv(
filename,
sep=r'\s+',
comment="#",
header=None,
skiprows=1
)
# Resolve columns by their headings: the CATIMA table builder writes
# Energy, mg/cm2, cm, Ein, sigE, ..., sigX, so the old positional layout
# interpreted energy as distance and shifted all interpolation inputs.
normalized_header = {
re.sub(r"[^a-z0-9]", "", name.lower()): index
for index, name in enumerate(header)
}
def find_column(names, fallback):
for name in names:
if name in normalized_header:
return normalized_header[name]
return fallback
x_column = find_column(("cm", "distancecm", "xcm"), 0)
energy_column = find_column(("energymev", "energy", "e"), 1)
sigma_E_column = find_column(("sige", "sigmaemev", "sigmae"), 2)
sigma_x_column = find_column(("sigx", "sigmaxcm", "sigmax"), 3)
x = data.iloc[:, x_column].to_numpy()
E = data.iloc[:, energy_column].to_numpy()
sigma_E = data.iloc[:, sigma_E_column].to_numpy() if data.shape[1] > sigma_E_column else None
sigma_x = data.iloc[:, sigma_x_column].to_numpy() if data.shape[1] > sigma_x_column 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()

21
Armory/anasenMS.cpp
View File

@ -40,19 +40,19 @@ void handler(int sig){
// Function to load energy loss table from file
TGraph* LoadELoss(const std::string& filename) {
TGraph* g = new TGraph(filename.c_str(), "%lg %lg");
TGraph* g = new TGraph(filename.c_str(), "%lg %lg"); //column 1 is path length (cm), column 2 is energy (MeV)
return g;
}
// Loads column 2 (Energy_MeV) vs column 4 (Sigma_x_cm) from an E_vs_x_*.dat table
TGraph* LoadSigmaXVsEnergy(const std::string& filename) {
TGraph* g = new TGraph(filename.c_str(), "%*lg %lg %*lg %lg");
TGraph* g = new TGraph(filename.c_str(), "%*lg %lg %*lg %*lg %*lg %*lg %*lg %lg"); // column 2 (Energy_MeV) vs column 4 (Sigma_x_cm)
g->Sort(); // TGraph::Eval requires ascending x (Energy_MeV)
return g;
}
TGraph* LoadSigmaAVsEnergy(const std::string& filename) { //load column 2 vs column 5 (Sigma_a_cm) from an E_vs_a_*.dat table
TGraph* g = new TGraph(filename.c_str(), "%*lg %lg %*lg %*lg %lg");
TGraph* g = new TGraph(filename.c_str(), "%*lg %lg %*lg %*lg %*lg %lg"); // column 2 (Energy_MeV) vs column 6 (Sigma_a_cm)
g->Sort(); // TGraph::Eval requires ascending x (Energy_MeV)
return g;
}
@ -168,8 +168,7 @@ int main(int argc, char **argv){
// 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 = {2.235};
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
const int kMBeam = reactionConfig.beamA; // mass number of beam
const int kMTarget = reactionConfig.targetA; // mass number of target
@ -204,8 +203,8 @@ int main(int argc, char **argv){
double vertexZRange[2] = { -174.3, 174.3}; // -174.3, 174.3 (full length of gas volume, centered at 0)
const double beamEntranceZ = -280 - 174.3; //vertexZRange[0]; // mm, assumed beam entrance into the gas
TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = path length (cm), y = beam energy (MeV)
TGraph* sigmaXBeam = LoadSigmaXVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_x (cm)
TGraph* elossBeam = LoadELoss("Eloss/E_vs_x_aluminum.dat"); // x = path length (cm), y = beam energy (MeV)
TGraph* sigmaXBeam = LoadSigmaXVsEnergy("Eloss/E_vs_x_aluminum.dat"); // x = beam energy (MeV), y = distance straggle sigma_x (cm)
//TGraph* sigmaABeam = LoadSigmaAVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_a (cm)
// Build a temporary inverse (energy -> path) to locate the path at beamE.
TGraph* elossBeamInverseRaw = new TGraph(elossBeam->GetN());
@ -690,15 +689,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,

6
Armory/testMacro.C
View File

@ -38,7 +38,7 @@
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");
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");
@ -66,13 +66,13 @@
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");
*/
}

100001
ELoss/HeLoss/E_vs_x_Al-27.dat
View File

File diff suppressed because it is too large Load Diff