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

range v thickness graph

Browse Source
This commit is contained in:
james szalkie 2026-07-20 10:06:59 -04:00
parent d6e72f1cc7
commit 26be64416c
3 changed files with 10030 additions and 10002 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
.DS_Store vendored
View File

Binary file not shown.

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

File diff suppressed because it is too large Load Diff

32
ELoss/PCEnergyAnalysis.py
View File

@ -503,7 +503,7 @@ def prepare_tree_data(tree, treename, particle, max_events=None, z_max=34.86):
sx3_hit_x,
sx3_hit_y,
sx3_hit_z,
SX3_SI_THICKNESS_CM,) * sx3_sin_theta
SX3_SI_THICKNESS_CM,)
qqq_silicon_path_cm = calculate_qqq_thickness_path(
qqq_vertex_x,
qqq_vertex_y,
@ -1393,6 +1393,35 @@ class MyInteractiveApp(cmd.Cmd):
plt.tight_layout()
plt.savefig(f"{base}/sx3E_vs_theta.png", dpi=300)
plt.show()
# --- Range vs Theta and overlay 0.1/sin(theta) path length ---
try:
mask_range = (~np.isnan(Ei)) & (~np.isnan(thetab)) & sx3_theta_plot_mask & (Ei > 0)
if np.any(mask_range):
theta_deg = thetab[mask_range]
theta_rad = np.deg2rad(theta_deg)
# compute range from initial energy to 0 in He (cm)
ranges_cm = energy_distance(particle, "Si", Ei[mask_range], 0.0)
plt.figure(figsize=(7,6))
plt.hist2d(theta_deg, ranges_cm, bins=200)
plt.xlabel("thetab (deg)")
plt.ylabel("Range (cm)")
plt.title(f"{particle} ({treename}) Range vs Theta (from Ei to 0)")
cb = plt.colorbar(label="Counts")
# overlay 0.1 / sin(theta) path (cm)
sin_th = np.sin(theta_rad)
sin_th = np.where(np.abs(sin_th) < 1e-6, 1e-6, sin_th)
path_cm = 0.1 / sin_th
order = np.argsort(theta_deg)
plt.plot(theta_deg[order], path_cm[order], color='red', linewidth=2, label='0.1/sin(theta)')
plt.legend()
plt.tight_layout()
plt.savefig(f"{base}/Range_vs_theta_with_path.png", dpi=300)
plt.show()
except Exception:
print("Could not compute Range vs Theta plot")
mask1 = (qqqE > 0) & ~np.isnan(thetabqqq) & qqq_theta_plot_mask
plt.figure(figsize=(7,6))
@ -1476,7 +1505,6 @@ class MyInteractiveApp(cmd.Cmd):
plt.legend(loc='upper right')
plt.show()
plt.figure(figsize=(7,5))
plt.hist(dA[sx3_theta_plot_mask], bins=100, label='dA', histtype='step')
plt.hist(dC[sx3_theta_plot_mask], bins=100, label='dC', histtype='step')