new file: anasen_labels.csv
new file: detector_geometry.dat new file: eloss_calculations/Eloss.py new file: eloss_calculations/alpha_lookup_20MeV.dat new file: eloss_calculations/alpha_lookup_6.0MeV.dat new file: eloss_calculations/make_eloss_table.C new file: eloss_calculations/proton_lookup_20MeV.dat new file: eloss_calculations/proton_lookup_6.0MeV.dat new file: grid_generate.py new file: shadowplay.py
This commit is contained in:
parent
c8923ca870
commit
9d581fa72d
5
.gitignore
vendored
5
.gitignore
vendored
|
|
@ -10,6 +10,7 @@ EventBuilder*
|
||||||
*.png
|
*.png
|
||||||
*.gif
|
*.gif
|
||||||
*.msh
|
*.msh
|
||||||
|
*.vtk
|
||||||
Mapper
|
Mapper
|
||||||
AnasenMS
|
AnasenMS
|
||||||
|
|
||||||
|
|
@ -26,12 +27,8 @@ MakePlotsQQQ.C
|
||||||
MakePlotsQQQ.h
|
MakePlotsQQQ.h
|
||||||
MakePlotsSX3.C
|
MakePlotsSX3.C
|
||||||
MakePlotsSX3.h
|
MakePlotsSX3.h
|
||||||
qqq_gains_det3.dat
|
|
||||||
qqq_relative_gains.dat
|
|
||||||
Armory/CorrelateQQQ.h
|
Armory/CorrelateQQQ.h
|
||||||
QQQStage2.C
|
QQQStage2.C
|
||||||
anasen_fem/scalars.dat.names
|
anasen_fem/scalars.dat.names
|
||||||
myenv/
|
|
||||||
eloss_calculations/
|
|
||||||
anasen_fem/He96_CO2_4_260Torr.gas
|
anasen_fem/He96_CO2_4_260Torr.gas
|
||||||
anasen_fem/heavy_ion_track.csv
|
anasen_fem/heavy_ion_track.csv
|
||||||
|
|
|
||||||
97
anasen_labels.csv
Normal file
97
anasen_labels.csv
Normal file
|
|
@ -0,0 +1,97 @@
|
||||||
|
X,Y,Z,Label
|
||||||
|
-28.779245994292484,28.779245994292488,189.3,A0
|
||||||
|
-6.173888892008441,46.89534194298123,204.3,C0
|
||||||
|
-35.24723393402666,20.349999999999998,189.3,A1
|
||||||
|
6.1738888920084465,46.89534194298123,204.3,C1
|
||||||
|
-39.31318112996508,10.533935135672605,189.3,A2
|
||||||
|
18.10092635086875,43.69950188778387,204.3,C2
|
||||||
|
-40.7,4.984312472529728e-15,189.3,A3
|
||||||
|
28.794415592112486,37.525612995775425,204.3,C3
|
||||||
|
-39.31318112996508,-10.533935135672598,189.3,A4
|
||||||
|
37.525612995775425,28.794415592112486,204.3,C4
|
||||||
|
-35.24723393402665,-20.350000000000005,189.3,A5
|
||||||
|
43.69950188778387,18.10092635086875,204.3,C5
|
||||||
|
-28.77924599429249,-28.779245994292484,189.3,A6
|
||||||
|
46.89534194298123,6.17388889200844,204.3,C6
|
||||||
|
-20.350000000000023,-35.247233934026646,189.3,A7
|
||||||
|
46.89534194298123,-6.17388889200844,204.3,C7
|
||||||
|
-10.53393513567259,-39.31318112996508,189.3,A8
|
||||||
|
43.69950188778387,-18.10092635086875,204.3,C8
|
||||||
|
-7.476468708794592e-15,-40.7,189.3,A9
|
||||||
|
37.525612995775425,-28.794415592112486,204.3,C9
|
||||||
|
10.533935135672577,-39.31318112996509,189.3,A10
|
||||||
|
28.794415592112486,-37.525612995775425,204.3,C10
|
||||||
|
20.350000000000005,-35.24723393402665,189.3,A11
|
||||||
|
18.10092635086875,-43.69950188778387,204.3,C11
|
||||||
|
28.77924599429248,-28.77924599429249,189.3,A12
|
||||||
|
6.1738888920084465,-46.89534194298123,204.3,C12
|
||||||
|
35.247233934026646,-20.350000000000023,189.3,A13
|
||||||
|
-6.173888892008441,-46.89534194298123,204.3,C13
|
||||||
|
39.31318112996508,-10.533935135672593,189.3,A14
|
||||||
|
-18.100926350868747,-43.69950188778387,204.3,C14
|
||||||
|
40.7,-9.968624945059456e-15,189.3,A15
|
||||||
|
-28.794415592112486,-37.525612995775425,204.3,C15
|
||||||
|
39.31318112996509,10.533935135672575,189.3,A16
|
||||||
|
-37.52561299577542,-28.7944155921125,204.3,C16
|
||||||
|
35.24723393402665,20.35,189.3,A17
|
||||||
|
-43.69950188778387,-18.100926350868754,204.3,C17
|
||||||
|
28.77924599429249,28.77924599429248,189.3,A18
|
||||||
|
-46.89534194298123,-6.17388889200844,204.3,C18
|
||||||
|
20.34999999999999,35.24723393402666,189.3,A19
|
||||||
|
-46.89534194298124,6.173888892008428,204.3,C19
|
||||||
|
10.533935135672595,39.31318112996508,189.3,A20
|
||||||
|
-43.699501887783875,18.100926350868743,204.3,C20
|
||||||
|
1.246078118132432e-14,40.7,189.3,A21
|
||||||
|
-37.525612995775425,28.794415592112486,204.3,C21
|
||||||
|
-10.533935135672573,39.31318112996509,189.3,A22
|
||||||
|
-28.7944155921125,37.52561299577542,204.3,C22
|
||||||
|
-20.34999999999997,35.247233934026674,189.3,A23
|
||||||
|
-18.100926350868736,43.699501887783875,204.3,C23
|
||||||
|
8.677040182954032,-100.82732255526474,0,D0_S0
|
||||||
|
20.363633298595417,-99.13002793745362,0,D0_S1
|
||||||
|
31.929590218792463,-96.03093912099365,0,D0_S2
|
||||||
|
42.89912404953582,-91.657542819904,0,D0_S3
|
||||||
|
57.928198505728936,-82.98050263694998,0,D1_S0
|
||||||
|
67.20043771866115,-75.66730582239823,0,D1_S1
|
||||||
|
75.66730582239825,-67.20043771866115,0,D1_S2
|
||||||
|
82.98050263694998,-57.928198505728936,0,D1_S3
|
||||||
|
91.657542819904,-42.89912404953581,0,D2_S0
|
||||||
|
96.03093912099365,-31.929590218792463,0,D2_S1
|
||||||
|
99.13002793745362,-20.363633298595417,0,D2_S2
|
||||||
|
100.82732255526474,-8.677040182954027,0,D2_S3
|
||||||
|
100.82732255526474,8.677040182954027,0,D3_S0
|
||||||
|
99.13002793745362,20.363633298595417,0,D3_S1
|
||||||
|
96.03093912099365,31.929590218792466,0,D3_S2
|
||||||
|
91.657542819904,42.89912404953581,0,D3_S3
|
||||||
|
82.98050263694998,57.928198505728936,0,D4_S0
|
||||||
|
75.66730582239825,67.20043771866115,0,D4_S1
|
||||||
|
67.20043771866115,75.66730582239823,0,D4_S2
|
||||||
|
57.928198505728936,82.98050263694998,0,D4_S3
|
||||||
|
42.89912404953582,91.657542819904,0,D5_S0
|
||||||
|
31.929590218792463,96.03093912099365,0,D5_S1
|
||||||
|
20.363633298595417,99.13002793745362,0,D5_S2
|
||||||
|
8.677040182954055,100.82732255526474,0,D5_S3
|
||||||
|
-8.677040182954043,100.82732255526474,0,D6_S0
|
||||||
|
-20.36363329859545,99.1300279374536,0,D6_S1
|
||||||
|
-31.929590218792452,96.03093912099366,0,D6_S2
|
||||||
|
-42.89912404953579,91.65754281990401,0,D6_S3
|
||||||
|
-57.92819850572895,82.98050263694996,0,D7_S0
|
||||||
|
-67.20043771866116,75.6673058223982,0,D7_S1
|
||||||
|
-75.66730582239823,67.20043771866116,0,D7_S2
|
||||||
|
-82.98050263694995,57.92819850572897,0,D7_S3
|
||||||
|
-91.657542819904,42.89912404953581,0,D8_S0
|
||||||
|
-96.03093912099366,31.929590218792452,0,D8_S1
|
||||||
|
-99.1300279374536,20.363633298595445,0,D8_S2
|
||||||
|
-100.82732255526474,8.67704018295404,0,D8_S3
|
||||||
|
-100.82732255526474,-8.677040182954059,0,D9_S0
|
||||||
|
-99.13002793745362,-20.36363329859542,0,D9_S1
|
||||||
|
-96.03093912099366,-31.929590218792423,0,D9_S2
|
||||||
|
-91.65754281990401,-42.899124049535786,0,D9_S3
|
||||||
|
-82.98050263694996,-57.92819850572895,0,D10_S0
|
||||||
|
-75.66730582239822,-67.20043771866116,0,D10_S1
|
||||||
|
-67.20043771866116,-75.66730582239823,0,D10_S2
|
||||||
|
-57.928198505728936,-82.98050263694998,0,D10_S3
|
||||||
|
-42.899124049535814,-91.657542819904,0,D11_S0
|
||||||
|
-31.929590218792413,-96.03093912099368,0,D11_S1
|
||||||
|
-20.363633298595452,-99.1300279374536,0,D11_S2
|
||||||
|
-8.67704018295409,-100.82732255526474,0,D11_S3
|
||||||
|
48
detector_geometry.dat
Normal file
48
detector_geometry.dat
Normal file
|
|
@ -0,0 +1,48 @@
|
||||||
|
=========================================================
|
||||||
|
SX3 BARREL AZIMUTHAL ANGLES (Degrees)
|
||||||
|
=========================================================
|
||||||
|
Det ID | Strip 0 | Strip 1 | Strip 2 | Strip 3 | Det Center
|
||||||
|
---------------------------------------------------------
|
||||||
|
0 | 85.08 | 78.39 | 71.61 | 64.92 | 75.00
|
||||||
|
1 | 55.08 | 48.39 | 41.61 | 34.92 | 45.00
|
||||||
|
2 | 25.08 | 18.39 | 11.61 | 4.92 | 15.00
|
||||||
|
3 | -4.92 | -11.61 | -18.39 | -25.08 | -15.00
|
||||||
|
4 | -34.92 | -41.61 | -48.39 | -55.08 | -45.00
|
||||||
|
5 | -64.92 | -71.61 | -78.39 | -85.08 | -75.00
|
||||||
|
6 | -94.92 | -101.61 | -108.39 | -115.08 | -105.00
|
||||||
|
7 | -124.92 | -131.61 | -138.39 | -145.08 | -135.00
|
||||||
|
8 | -154.92 | -161.61 | -168.39 | -175.08 | -165.00
|
||||||
|
9 | 175.08 | 168.39 | 161.61 | 154.92 | 165.00
|
||||||
|
10 | 145.08 | 138.39 | 131.61 | 124.92 | 135.00
|
||||||
|
11 | 115.08 | 108.39 | 101.61 | 94.92 | 105.00
|
||||||
|
|
||||||
|
|
||||||
|
=========================================================
|
||||||
|
PROPORTIONAL COUNTER WIRE ANGLES (Degrees)
|
||||||
|
=========================================================
|
||||||
|
Wire ID | Anode (-z) | Anode (+z) | Cathode (-z) | Cathode (+z)
|
||||||
|
----------------------------------------------------------------
|
||||||
|
0 | -135.00 | -90.00 | -97.50 | -142.50
|
||||||
|
1 | -150.00 | -105.00 | -82.50 | -127.50
|
||||||
|
2 | -165.00 | -120.00 | -67.50 | -112.50
|
||||||
|
3 | -180.00 | -135.00 | -52.50 | -97.50
|
||||||
|
4 | 165.00 | -150.00 | -37.50 | -82.50
|
||||||
|
5 | 150.00 | -165.00 | -22.50 | -67.50
|
||||||
|
6 | 135.00 | -180.00 | -7.50 | -52.50
|
||||||
|
7 | 120.00 | 165.00 | 7.50 | -37.50
|
||||||
|
8 | 105.00 | 150.00 | 22.50 | -22.50
|
||||||
|
9 | 90.00 | 135.00 | 37.50 | -7.50
|
||||||
|
10 | 75.00 | 120.00 | 52.50 | 7.50
|
||||||
|
11 | 60.00 | 105.00 | 67.50 | 22.50
|
||||||
|
12 | 45.00 | 90.00 | 82.50 | 37.50
|
||||||
|
13 | 30.00 | 75.00 | 97.50 | 52.50
|
||||||
|
14 | 15.00 | 60.00 | 112.50 | 67.50
|
||||||
|
15 | 0.00 | 45.00 | 127.50 | 82.50
|
||||||
|
16 | -15.00 | 30.00 | 142.50 | 97.50
|
||||||
|
17 | -30.00 | 15.00 | 157.50 | 112.50
|
||||||
|
18 | -45.00 | 0.00 | 172.50 | 127.50
|
||||||
|
19 | -60.00 | -15.00 | -172.50 | 142.50
|
||||||
|
20 | -75.00 | -30.00 | -157.50 | 157.50
|
||||||
|
21 | -90.00 | -45.00 | -142.50 | 172.50
|
||||||
|
22 | -105.00 | -60.00 | -127.50 | -172.50
|
||||||
|
23 | -120.00 | -75.00 | -112.50 | -157.50
|
||||||
62
eloss_calculations/Eloss.py
Normal file
62
eloss_calculations/Eloss.py
Normal file
|
|
@ -0,0 +1,62 @@
|
||||||
|
import pycatima as catima
|
||||||
|
import numpy as np
|
||||||
|
|
||||||
|
# --- 1. Constants ---
|
||||||
|
P_TORR = 250
|
||||||
|
TEMP_K = 293.15
|
||||||
|
R = 8.3144
|
||||||
|
MEV2U = 1.0 / 931.494
|
||||||
|
|
||||||
|
# Gas Density Calculations
|
||||||
|
p_pa = P_TORR * 133.322
|
||||||
|
molar_density = p_pa / (R * TEMP_K)
|
||||||
|
m_he, m_c, m_o = 4.0026, 12.0000, 15.9949
|
||||||
|
m_mix_avg = (0.96 * m_he) + (0.04 * (m_c + 2*m_o))
|
||||||
|
rho_g_cm3 = (molar_density * m_mix_avg) / 1e6
|
||||||
|
print(f"Gas density at {P_TORR} Torr: {rho_g_cm3:.6e} g/cm^3")
|
||||||
|
|
||||||
|
# --- 2. Material & Step Setup ---
|
||||||
|
material_def = [(m_he, 2, 0.96), (m_c, 6, 0.04), (m_o, 8, 0.08)]
|
||||||
|
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
|
||||||
|
step_g_cm2 = step_mg_cm2 / 1000.0
|
||||||
|
max_steps = 100000 # Adjust based on how far you want to track
|
||||||
|
|
||||||
|
def generate_lookup(z, mass_u, e_start_mev, label):
|
||||||
|
filename = f"{label}_lookup_{e_start_mev}MeV.dat"
|
||||||
|
projectile = catima.Projectile(mass_u, z)
|
||||||
|
|
||||||
|
current_e_total = e_start_mev
|
||||||
|
current_thickness_g_cm2 = 0.0
|
||||||
|
|
||||||
|
output = []
|
||||||
|
header = f"Energy(MeV) \tmg/cm2 \tcm\nStarting Energy: {e_start_mev} MeV"
|
||||||
|
|
||||||
|
for i in range(max_steps):
|
||||||
|
# 1. Record current state
|
||||||
|
dist_cm = current_thickness_g_cm2 / rho_g_cm3
|
||||||
|
output.append([current_e_total, current_thickness_g_cm2 * 1000.0, dist_cm])
|
||||||
|
|
||||||
|
# 2. Calculate energy loss for the NEXT step
|
||||||
|
e_u = current_e_total / mass_u
|
||||||
|
if e_u < 0.01: # Stop at ATIMA limit
|
||||||
|
break
|
||||||
|
|
||||||
|
projectile.T(e_u)
|
||||||
|
# dedx returns MeV / (g/cm2)
|
||||||
|
loss_mev = catima.dedx(projectile, gas_mix) * step_g_cm2
|
||||||
|
|
||||||
|
# 3. Update values
|
||||||
|
current_e_total -= loss_mev
|
||||||
|
current_thickness_g_cm2 += step_g_cm2
|
||||||
|
|
||||||
|
np.savetxt(filename, output, fmt='%.6f', delimiter='\t', header=header)
|
||||||
|
print(f"Lookup table created: {filename}")
|
||||||
|
|
||||||
|
# --- 3. Run ---
|
||||||
|
# Format: generate_lookup(Z, mass_u, E_start_MeV, label)
|
||||||
|
generate_lookup(1, 1.0078, 20, "proton")
|
||||||
|
generate_lookup(2, 4.0026, 20, "alpha")
|
||||||
31981
eloss_calculations/alpha_lookup_20MeV.dat
Normal file
31981
eloss_calculations/alpha_lookup_20MeV.dat
Normal file
File diff suppressed because it is too large
Load Diff
4263
eloss_calculations/alpha_lookup_6.0MeV.dat
Normal file
4263
eloss_calculations/alpha_lookup_6.0MeV.dat
Normal file
File diff suppressed because it is too large
Load Diff
29
eloss_calculations/make_eloss_table.C
Normal file
29
eloss_calculations/make_eloss_table.C
Normal file
|
|
@ -0,0 +1,29 @@
|
||||||
|
#include "/home/sud/Desktop/Software2/propagator/elastcaller.h"
|
||||||
|
void make_eloss_table_protons() {
|
||||||
|
double einput = 20.0, estepnow;
|
||||||
|
double target_thickness_unit = 4e-2; //mg/cm2.
|
||||||
|
double density = 0.0711;//mg/cm3
|
||||||
|
long i=0;
|
||||||
|
while(einput > 0.001) {
|
||||||
|
std::cout << "After " << i << " steps, 1H is at " << einput << " MeV after penetrating " << i*target_thickness_unit << " mg/cm2 " << i*target_thickness_unit/density << " cm of HeCO2" << std::endl;
|
||||||
|
estepnow = slowmedown("1H",einput,"3(12C)6(16O)97(4He)",target_thickness_unit);
|
||||||
|
|
||||||
|
einput = estepnow;
|
||||||
|
i+=1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void make_eloss_table() {
|
||||||
|
double einput = 20.0, estepnow;
|
||||||
|
double target_thickness_unit = 1e-3; //mg/cm2.
|
||||||
|
double density = 0.0711;//mg/cm3
|
||||||
|
long i=0;
|
||||||
|
while(einput > 0.001) {
|
||||||
|
std::cout << "After " << i << " steps, 4He is at " << einput << " MeV after penetrating " << i*target_thickness_unit << " mg/cm2 " << i*target_thickness_unit/density << " cm of HeCO2" << std::endl;
|
||||||
|
estepnow = slowmedown("4He",einput,"3(12C)6(16O)97(4He)",target_thickness_unit);
|
||||||
|
|
||||||
|
|
||||||
|
einput = estepnow;
|
||||||
|
i+=1;
|
||||||
|
}
|
||||||
|
}
|
||||||
397468
eloss_calculations/proton_lookup_20MeV.dat
Normal file
397468
eloss_calculations/proton_lookup_20MeV.dat
Normal file
File diff suppressed because it is too large
Load Diff
44120
eloss_calculations/proton_lookup_6.0MeV.dat
Normal file
44120
eloss_calculations/proton_lookup_6.0MeV.dat
Normal file
File diff suppressed because it is too large
Load Diff
56
grid_generate.py
Normal file
56
grid_generate.py
Normal file
|
|
@ -0,0 +1,56 @@
|
||||||
|
import math
|
||||||
|
|
||||||
|
def wrap_phi(angle):
|
||||||
|
"""Wraps an angle to be between -180 and +180 degrees"""
|
||||||
|
return (angle + 180) % 360 - 180
|
||||||
|
|
||||||
|
# The name of the file you want to generate
|
||||||
|
filename = "detector_geometry.dat"
|
||||||
|
|
||||||
|
# Open the file in 'write' mode
|
||||||
|
with open(filename, "w") as f:
|
||||||
|
|
||||||
|
# --- 1. SX3 Geometry ---
|
||||||
|
f.write("=========================================================\n")
|
||||||
|
f.write(" SX3 BARREL AZIMUTHAL ANGLES (Degrees) \n")
|
||||||
|
f.write("=========================================================\n")
|
||||||
|
f.write(" Det ID | Strip 0 | Strip 1 | Strip 2 | Strip 3 | Det Center\n")
|
||||||
|
f.write("---------------------------------------------------------\n")
|
||||||
|
|
||||||
|
for det_id in range(12): # Assuming 12 barrel detectors
|
||||||
|
strips = []
|
||||||
|
for stripF in range(4):
|
||||||
|
stripF_rev = 3 - stripF
|
||||||
|
num = (2 * stripF_rev - 3) * 40.30
|
||||||
|
den = 8.0 * 88.0 * math.cos(math.radians(15.0))
|
||||||
|
beta_n = 15.0 + math.degrees(math.atan2(num, den))
|
||||||
|
phi_n = ((-det_id + 0.5) * 30.0 + beta_n) + 45.0
|
||||||
|
strips.append(wrap_phi(phi_n))
|
||||||
|
|
||||||
|
det_center = wrap_phi(((-det_id + 0.5) * 30.0 + 15.0) + 45.0)
|
||||||
|
f.write(f" {det_id:2d} | {strips[0]:6.2f} | {strips[1]:6.2f} | {strips[2]:6.2f} | {strips[3]:6.2f} | {det_center:6.2f}\n")
|
||||||
|
|
||||||
|
|
||||||
|
# --- 2. PC Wire Geometry ---
|
||||||
|
f.write("\n\n=========================================================\n")
|
||||||
|
f.write(" PROPORTIONAL COUNTER WIRE ANGLES (Degrees) \n")
|
||||||
|
f.write("=========================================================\n")
|
||||||
|
f.write(" Wire ID | Anode (-z) | Anode (+z) | Cathode (-z) | Cathode (+z)\n")
|
||||||
|
f.write("----------------------------------------------------------------\n")
|
||||||
|
|
||||||
|
k = 360.0 / 24.0
|
||||||
|
offset_a1 = -6*k - 3*k
|
||||||
|
offset_c1 = -4*k - 2*k - (360.0/48.0)
|
||||||
|
wireShift = 3
|
||||||
|
offset_a2 = offset_a1 + wireShift*k
|
||||||
|
offset_c2 = offset_c1 - wireShift*k
|
||||||
|
|
||||||
|
for i in range(24):
|
||||||
|
phi_a1 = wrap_phi(-k * i + offset_a1)
|
||||||
|
phi_a2 = wrap_phi(-k * i + offset_a2)
|
||||||
|
phi_c1 = wrap_phi(k * i + offset_c1)
|
||||||
|
phi_c2 = wrap_phi(k * i + offset_c2)
|
||||||
|
|
||||||
|
f.write(f" {i:2d} | {phi_a1:7.2f} | {phi_a2:7.2f} | {phi_c1:7.2f} | {phi_c2:7.2f}\n")
|
||||||
|
|
||||||
|
print(f"Success! The geometry lookup table has been saved to '{filename}'.")
|
||||||
152
shadowplay.py
Normal file
152
shadowplay.py
Normal file
|
|
@ -0,0 +1,152 @@
|
||||||
|
import math
|
||||||
|
import csv
|
||||||
|
|
||||||
|
# ==========================================
|
||||||
|
# CONFIGURATION
|
||||||
|
# ==========================================
|
||||||
|
VTK_FILE = "anasen_geometry.vtk"
|
||||||
|
CSV_FILE = "anasen_labels.csv"
|
||||||
|
|
||||||
|
Z_SOURCE = 53.0 # Source position
|
||||||
|
R_A = 37.0 # Anode radius (mm)
|
||||||
|
R_C = 43.0 # Cathode radius (mm)
|
||||||
|
R_SX3 = 88.0 # SX3 Barrel radius (mm)
|
||||||
|
Z_LEN = 348.6 # Total Z length (174.3 * 2 from ClassPW.h)
|
||||||
|
Z_HALF = Z_LEN / 2.0
|
||||||
|
|
||||||
|
points = []
|
||||||
|
lines = []
|
||||||
|
polygons = []
|
||||||
|
cell_colors = []
|
||||||
|
labels = []
|
||||||
|
|
||||||
|
def add_point(x_code, y_code, z_code):
|
||||||
|
# TRANSFORM: Map the C++ Frame (+Y Down, +Z Into Page)
|
||||||
|
# to the Visual ParaView Frame (+Y Up, +Z Out of Page)
|
||||||
|
points.append((x_code, -y_code, -z_code))
|
||||||
|
return len(points) - 1
|
||||||
|
|
||||||
|
def project_shadow(r_wire, phi_minusZ_deg, phi_plusZ_deg, color_code):
|
||||||
|
shadow_pts = []
|
||||||
|
steps = 50
|
||||||
|
for i in range(steps + 1):
|
||||||
|
t = i / float(steps)
|
||||||
|
# 1. Coordinate on the wire in CODE frame
|
||||||
|
z_w = -Z_HALF + (t * Z_LEN)
|
||||||
|
phi_w_deg = phi_minusZ_deg * (1 - t) + phi_plusZ_deg * t
|
||||||
|
phi_w = math.radians(phi_w_deg)
|
||||||
|
|
||||||
|
# 2. Ray from Source through the wire
|
||||||
|
dx = r_wire * math.cos(phi_w)
|
||||||
|
dy = r_wire * math.sin(phi_w)
|
||||||
|
dz = z_w - Z_SOURCE
|
||||||
|
|
||||||
|
# 3. Scale ray until it hits SX3 (R = 88)
|
||||||
|
r_current = math.sqrt(dx**2 + dy**2)
|
||||||
|
alpha = R_SX3 / r_current
|
||||||
|
|
||||||
|
x_proj = alpha * dx
|
||||||
|
y_proj = alpha * dy
|
||||||
|
z_proj = Z_SOURCE + alpha * dz
|
||||||
|
|
||||||
|
shadow_pts.append(add_point(x_proj, y_proj, z_proj))
|
||||||
|
|
||||||
|
lines.append(shadow_pts)
|
||||||
|
cell_colors.append(color_code)
|
||||||
|
|
||||||
|
# ==========================================
|
||||||
|
# 1. GENERATE PC WIRES & SHADOWS (LINES)
|
||||||
|
# ==========================================
|
||||||
|
k = 360.0 / 24.0
|
||||||
|
|
||||||
|
for i in range(24):
|
||||||
|
# --- ANODES ---
|
||||||
|
# From ClassPW.h: offset_a1 = -135 deg, offset_a2 = -90 deg
|
||||||
|
phi_a_minusZ = -k * i - 135.0
|
||||||
|
phi_a_plusZ = -k * i - 90.0
|
||||||
|
|
||||||
|
x_a1, y_a1 = R_A * math.cos(math.radians(phi_a_minusZ)), R_A * math.sin(math.radians(phi_a_minusZ))
|
||||||
|
x_a2, y_a2 = R_A * math.cos(math.radians(phi_a_plusZ)), R_A * math.sin(math.radians(phi_a_plusZ))
|
||||||
|
|
||||||
|
pa1 = add_point(x_a1, y_a1, -Z_HALF)
|
||||||
|
pa2 = add_point(x_a2, y_a2, Z_HALF)
|
||||||
|
lines.append([pa1, pa2])
|
||||||
|
cell_colors.append(1) # Color 1 = Anodes
|
||||||
|
project_shadow(R_A, phi_a_minusZ, phi_a_plusZ, 4)
|
||||||
|
|
||||||
|
# Place label at +Z_visual (which maps from -Z_code in the transform)
|
||||||
|
labels.append((x_a1 * 1.1, -y_a1 * 1.1, Z_HALF + 15, f"A{i}"))
|
||||||
|
|
||||||
|
# --- CATHODES ---
|
||||||
|
# From ClassPW.h: offset_c1 = -97.5 deg, offset_c2 = -142.5 deg
|
||||||
|
phi_c_minusZ = k * i - 97.5
|
||||||
|
phi_c_plusZ = k * i - 142.5
|
||||||
|
|
||||||
|
x_c1, y_c1 = R_C * math.cos(math.radians(phi_c_minusZ)), R_C * math.sin(math.radians(phi_c_minusZ))
|
||||||
|
x_c2, y_c2 = R_C * math.cos(math.radians(phi_c_plusZ)), R_C * math.sin(math.radians(phi_c_plusZ))
|
||||||
|
|
||||||
|
pc1 = add_point(x_c1, y_c1, -Z_HALF)
|
||||||
|
pc2 = add_point(x_c2, y_c2, Z_HALF)
|
||||||
|
lines.append([pc1, pc2])
|
||||||
|
cell_colors.append(2) # Color 2 = Cathodes
|
||||||
|
project_shadow(R_C, phi_c_minusZ, phi_c_plusZ, 5)
|
||||||
|
|
||||||
|
labels.append((x_c1 * 1.1, -y_c1 * 1.1, Z_HALF + 30, f"C{i}"))
|
||||||
|
|
||||||
|
# ==========================================
|
||||||
|
# 2. GENERATE SX3 BARREL (POLYGONS)
|
||||||
|
# ==========================================
|
||||||
|
# Calculate exact strip positions based on MakeVertex.C formulas
|
||||||
|
for det_id in range(12):
|
||||||
|
for stripF in range(4):
|
||||||
|
stripF_rev = 3 - stripF
|
||||||
|
num = (2 * stripF_rev - 3) * 40.30
|
||||||
|
den = 8.0 * R_SX3 * math.cos(math.radians(15.0))
|
||||||
|
beta_n = 15.0 + math.degrees(math.atan2(num, den))
|
||||||
|
|
||||||
|
# phi_n perfectly computes the C++ angle
|
||||||
|
phi_code_deg = ((-det_id + 0.5) * 30.0 + beta_n) + 45.0
|
||||||
|
phi_code = math.radians(phi_code_deg)
|
||||||
|
|
||||||
|
d_phi = 10.0 / R_SX3 # ~6.5 degrees width per strip
|
||||||
|
phi_L = phi_code - (d_phi / 2.0)
|
||||||
|
phi_R = phi_code + (d_phi / 2.0)
|
||||||
|
|
||||||
|
p1 = add_point(R_SX3 * math.cos(phi_L), R_SX3 * math.sin(phi_L), -Z_HALF)
|
||||||
|
p2 = add_point(R_SX3 * math.cos(phi_R), R_SX3 * math.sin(phi_R), -Z_HALF)
|
||||||
|
p3 = add_point(R_SX3 * math.cos(phi_R), R_SX3 * math.sin(phi_R), Z_HALF)
|
||||||
|
p4 = add_point(R_SX3 * math.cos(phi_L), R_SX3 * math.sin(phi_L), Z_HALF)
|
||||||
|
|
||||||
|
polygons.append([p1, p2, p3, p4])
|
||||||
|
cell_colors.append(3) # SX3 Color
|
||||||
|
|
||||||
|
# Place Labels correctly in visual space
|
||||||
|
x_lbl, y_lbl = R_SX3 * 1.15 * math.cos(phi_code), R_SX3 * 1.15 * math.sin(phi_code)
|
||||||
|
labels.append((x_lbl, -y_lbl, 0, f"D{det_id}_S{stripF}"))
|
||||||
|
|
||||||
|
# ==========================================
|
||||||
|
# 3. WRITE VTK AND CSV
|
||||||
|
# ==========================================
|
||||||
|
with open(VTK_FILE, "w") as f:
|
||||||
|
f.write("# vtk DataFile Version 2.0\nANASEN Geometry\nASCII\nDATASET POLYDATA\n")
|
||||||
|
f.write(f"POINTS {len(points)} float\n")
|
||||||
|
for pt in points: f.write(f"{pt[0]:.4f} {pt[1]:.4f} {pt[2]:.4f}\n")
|
||||||
|
|
||||||
|
line_data_size = sum(len(l) + 1 for l in lines)
|
||||||
|
f.write(f"\nLINES {len(lines)} {line_data_size}\n")
|
||||||
|
for l in lines: f.write(f"{len(l)} " + " ".join(map(str, l)) + "\n")
|
||||||
|
|
||||||
|
poly_data_size = sum(len(p) + 1 for p in polygons)
|
||||||
|
f.write(f"\nPOLYGONS {len(polygons)} {poly_data_size}\n")
|
||||||
|
for p in polygons: f.write(f"{len(p)} " + " ".join(map(str, p)) + "\n")
|
||||||
|
|
||||||
|
f.write(f"\nCELL_DATA {len(lines) + len(polygons)}\n")
|
||||||
|
f.write("SCALARS EntityType int 1\nLOOKUP_TABLE default\n")
|
||||||
|
for c in cell_colors: f.write(f"{c}\n")
|
||||||
|
|
||||||
|
with open(CSV_FILE, "w", newline='') as f:
|
||||||
|
writer = csv.writer(f)
|
||||||
|
writer.writerow(["X", "Y", "Z", "Label"])
|
||||||
|
writer.writerows(labels)
|
||||||
|
|
||||||
|
print(f"Generated {VTK_FILE} and {CSV_FILE}.")
|
||||||
Loading…
Reference in New Issue
Block a user