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

modified: Armory/ClassPW.h changes made to explain the andoe and cathode shifts in a more intuitive geometric manner

Browse Source 
modified:   MakeVertex.C changed code to make nA plots instead of nA0C to add more stats to n anode analysis
	modified:   anasen_fem/paraview_plotter.py changed the analysis to plot every 32nd field line, dure to indreaed field density from the new wire shift
	modified:   anasen_fem/run.py going from 0 to 174.6 in 17.43 increments to account for the remeasured lenth of the pc
	modified:   anasen_fem/scalars.dat.names
	modified:   anasen_fem/wires_gmsh2d_bc.py length of pc changed, 4 wire shift incorporated, set mesh order to 2 to decrease the meashing density
    pccal folder seems to make 17F analysis work, not sure why still trying to figure it out
	new file:   pccal/anode_gainmatch.C
	new file:   pccal/anode_gm_coeffs.dat
	new file:   pccal/cathode_gainmatch.C
	new file:   pccal/cathode_gm_coeffs.dat
	new file:   pccal/pc_gm_coeffs.dat
	new file:   pccal/slope_intercept_26Al.dat
    scratch folder analysis for pcz vs pczfix analysis from sudarsan's branch
	new file:   scratch/sx3z_vs_phiz/scan_offset.C
	new file:   scratch/sx3z_vs_phiz/scan_offset_fix.C
This commit is contained in:
Vignesh Sitaraman 2026-05-19 17:57:54 -04:00
parent d8593b2f55
commit bb05baf89d
14 changed files with 587 additions and 126 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Armory/ClassPW.h
View File

@ -154,8 +154,8 @@ inline void PW::ConstructGeo()
std::pair<TVector3, TVector3> q1; // cathode std::pair<TVector3, TVector3> q1; // cathode
double k = TMath::TwoPi() / 24.; // 48 solder thru holes, wires in every other one double k = TMath::TwoPi() / 24.; // 48 solder thru holes, wires in every other one
double offset_a1 = -5 * k - 5 * k; double offset_a1 = -6 * k - 4 * k; // -6 to go from 0,0 to 90degree up, and 4 for the wire offset //old version -5 * k - 5 * k;
double offset_c1 = -5 * k - 2 * k - TMath::TwoPi() / 48; // correct for a half-turn double offset_c1 = -6 * k + k / 2.0; // correct for a half-turn
// std::cerr << "Here!" << std::endl; // std::cerr << "Here!" << std::endl;
// #include "../scratch/testing.h" // #include "../scratch/testing.h"
double offset_a2 = offset_a1 + wireShift * k; double offset_a2 = offset_a1 + wireShift * k;

66
MakeVertex.C
View File

@ -3,8 +3,8 @@
// Comment out any line below to disable that diagnostic section // Comment out any line below to disable that diagnostic section
#define DIAG_1WIRE // raw per-wire dE vs Si (no PC required) #define DIAG_1WIRE // raw per-wire dE vs Si (no PC required)
#define DIAG_PC_SX3 // PC-SX3 coincidence analysis #define DIAG_PC_SX3 // PC-SX3 coincidence analysis
// #define DIAG_NA0C_SX3 // nA0C (n>=1) pseudo-wire vertex with SX3 #define DIAG_NA_SX3 // nA (n>=1) pseudo-wire vertex with SX3
// #define DIAG_NA0C_QQQ // nA0C (n>=1) pseudo-wire vertex with QQQ #define DIAG_NA_QQQ // nA (n>=1) pseudo-wire vertex with QQQ
#define DIAG_PC_QQQ // PC-QQQ coincidence analysis #define DIAG_PC_QQQ // PC-QQQ coincidence analysis
Int_t colors[40] = { Int_t colors[40] = {
@ -1269,12 +1269,12 @@ Bool_t MakeVertex::Process(Long64_t entry)
} }
}*/ }*/
///////////////////nA0C analysis using pseudo-wire (GetPseudoWire + getClosestWirePosAtWirePhi)/////////////////// ///////////////////nA analysis using pseudo-wire (GetPseudoWire + getClosestWirePosAtWirePhi)///////////////////
#if defined(DIAG_NA0C_SX3) || defined(DIAG_NA0C_QQQ) #if defined(DIAG_NA_SX3) || defined(DIAG_NA_QQQ)
if (cClusters.size() == 0 && aClusters.size() >= 1) if (aClusters.size() >= 1)
{ {
std::string nA0C_label = std::to_string(aClusters.size()) + "A0C"; std::string nA_label = std::to_string(aClusters.size()) + "A";
// Flatten all anode clusters into a combined hit list for the pseudo-wire // Flatten all anode clusters into a combined hit list for the pseudo-wire
std::vector<std::tuple<int, double, double>> allAnodeHits; std::vector<std::tuple<int, double, double>> allAnodeHits;
@ -1284,7 +1284,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(allAnodeHits, "ANODE"); auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(allAnodeHits, "ANODE");
#ifdef DIAG_NA0C_SX3 #ifdef DIAG_NA_SX3
for (auto sx3event : SX3_Events) for (auto sx3event : SX3_Events)
{ {
if (sx3event.Time1 - apTSMaxE < -150) if (sx3event.Time1 - apTSMaxE < -150)
@ -1328,26 +1328,26 @@ Bool_t MakeVertex::Process(Long64_t entry)
double Ex_from_proton = apkin_p.getExc(sx3Efix, theta_recon * 180. / M_PI); double Ex_from_proton = apkin_p.getExc(sx3Efix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(sx3Efixalpha, theta_recon * 180. / M_PI); double Ex_from_alpha = apkin_a.getExc(sx3Efixalpha, theta_recon * 180. / M_PI);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_SX3", 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA0C_label); plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_SX3", 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA_label);
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label); plotter->Fill1D(nA_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA_label);
plotter->Fill1D(nA0C_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label); plotter->Fill1D(nA_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA_label);
plotter->Fill2D(nA0C_label + "_sx3_E_vs_theta_raw_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3event.Energy1, nA0C_label); plotter->Fill2D(nA_label + "_sx3_E_vs_theta_raw_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3event.Energy1, nA_label);
plotter->Fill2D(nA0C_label + "_sx3_E_vs_theta_corr_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3Efix, nA0C_label); plotter->Fill2D(nA_label + "_sx3_E_vs_theta_corr_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3Efix, nA_label);
if (vtx_gate != "") if (vtx_gate != "")
{ {
plotter->Fill1D(nA0C_label + "_twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA0C_label); plotter->Fill1D(nA_label + "_twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA_label);
plotter->Fill1D(nA0C_label + "_twisted_vertex_recon_SX3" + vtx_gate, 600, -300, 300, vertex_recon, nA0C_label); plotter->Fill1D(nA_label + "_twisted_vertex_recon_SX3" + vtx_gate, 600, -300, 300, vertex_recon, nA_label);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_SX3" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA0C_label); plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_SX3" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA_label);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_SX3_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efixalpha, apSumE * sinTheta, nA0C_label); plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_SX3_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efixalpha, apSumE * sinTheta, nA_label);
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label); plotter->Fill1D(nA_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA_label);
plotter->Fill1D(nA0C_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label); plotter->Fill1D(nA_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA_label);
} }
} }
} }
#endif // DIAG_NA0C_SX3 #endif // DIAG_nA_SX3
#ifdef DIAG_NA0C_QQQ #ifdef DIAG_NA_QQQ
for (auto qqqevent : QQQ_Events) for (auto qqqevent : QQQ_Events)
{ {
if (qqqevent.Time1 - apTSMaxE < -150) if (qqqevent.Time1 - apTSMaxE < -150)
@ -1391,26 +1391,26 @@ Bool_t MakeVertex::Process(Long64_t entry)
double Ex_from_proton = apkin_p.getExc(qqqEfix, theta_recon * 180. / M_PI); double Ex_from_proton = apkin_p.getExc(qqqEfix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(qqqEfixalpha, theta_recon * 180. / M_PI); double Ex_from_alpha = apkin_a.getExc(qqqEfixalpha, theta_recon * 180. / M_PI);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ", 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA0C_label); plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ", 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA_label);
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label); plotter->Fill1D(nA_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA_label);
plotter->Fill1D(nA0C_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label); plotter->Fill1D(nA_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA_label);
plotter->Fill2D(nA0C_label + "_qqq_E_vs_theta_raw_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqevent.Energy1, nA0C_label); plotter->Fill2D(nA_label + "_qqq_E_vs_theta_raw_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqevent.Energy1, nA_label);
plotter->Fill2D(nA0C_label + "_qqq_E_vs_theta_corr_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqEfix, nA0C_label); plotter->Fill2D(nA_label + "_qqq_E_vs_theta_corr_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqEfix, nA_label);
if (vtx_gate != "") if (vtx_gate != "")
{ {
plotter->Fill1D(nA0C_label + "_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA0C_label); plotter->Fill1D(nA_label + "_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA_label);
plotter->Fill1D(nA0C_label + "_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, nA0C_label); plotter->Fill1D(nA_label + "_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, nA_label);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA0C_label); plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA_label);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfixalpha, apSumE * sinTheta, nA0C_label); plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfixalpha, apSumE * sinTheta, nA_label);
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label); plotter->Fill1D(nA_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA_label);
plotter->Fill1D(nA0C_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label); plotter->Fill1D(nA_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA_label);
} }
} }
} }
#endif // DIAG_NA0C_QQQ #endif // DIAG_nA_QQQ
} }
#endif // DIAG_NA0C_SX3 || DIAG_NA0C_QQQ #endif // DIAG_nA_SX3 || DIAG_nA_QQQ
#ifdef DIAG_PC_QQQ #ifdef DIAG_PC_QQQ
for (auto pcevent : PC_Events) for (auto pcevent : PC_Events)

2
anasen_fem/paraview_plotter.py
View File

@ -95,7 +95,7 @@ glyph.ScaleArray = ['POINTS', 'No scale array']
glyph.ScaleFactor = 0.001 glyph.ScaleFactor = 0.001
glyph.GlyphMode = 'Every Nth Point' glyph.GlyphMode = 'Every Nth Point'
glyph.Stride = 24 glyph.Stride = 32
# --- 3. Display the Glyphs --- # --- 3. Display the Glyphs ---
glyph_display = Show(glyph, renderView) glyph_display = Show(glyph, renderView)

10
anasen_fem/run.py
View File

@ -1,10 +1,10 @@
import code import code
import os import os
# val=-178.3 # val=-174.3
val=89.15 val=0
count=11 count=11
while val<178.3+0.1: while val<174.3+0.1:
print(val) print(val)
os.system("python3 wires_gmsh2d_bc.py "+str(val)) os.system("python3 wires_gmsh2d_bc.py "+str(val))
os.system("ElmerGrid 14 2 wires2d.msh -2d") os.system("ElmerGrid 14 2 wires2d.msh -2d")
@ -16,9 +16,9 @@ while val<178.3+0.1:
os.system("cp wires2d/elfield_anasen_t0001.vtu wires2d/vtu_files/elfield_anasen_%02d_%1.4f.vtu"%(count,val)) os.system("cp wires2d/elfield_anasen_t0001.vtu wires2d/vtu_files/elfield_anasen_%02d_%1.4f.vtu"%(count,val))
os.system("cp contour_output.png png/Contour_output_z_%02d_%1.4f.png"%(count,val)) os.system("cp contour_output.png png/Contour_output_z_%02d_%1.4f.png"%(count,val))
os.system("cp Field_output.png png/Field_ouput_z_%02d_%1.4f.png"%(count,val)) os.system("cp Field_output.png png/Field_ouput_z_%02d_%1.4f.png"%(count,val))
val=val+17.83 val=val+17.43
count = count + 1 count = count + 1
break # break
# os.system("tar -cvzf wiress2d/mesh.tar.gz wires2d/mesh_files") # os.system("tar -cvzf wiress2d/mesh.tar.gz wires2d/mesh_files")
# os.system("rm -rf wires2d/mesh_files/*") # os.system("rm -rf wires2d/mesh_files/*")

2
anasen_fem/scalars.dat.names
View File

@ -2,7 +2,7 @@ Metadata for SaveScalars file: ./scalars.dat
Elmer version: 26.2 Elmer version: 26.2
Elmer compilation date: 2026-05-14 Elmer compilation date: 2026-05-14
Solver input file: wires2d.sif Solver input file: wires2d.sif
File started at: 2026/05/18 09:41:59 File started at: 2026/05/19 16:24:56
Variables in columns of matrix: Variables in columns of matrix:
1: res: potential difference 1: res: potential difference

95
anasen_fem/wires_gmsh2d_bc.py
View File

@ -7,7 +7,7 @@ include_needle = True
gmsh.initialize() gmsh.initialize()
gmsh.model.add("adaptive_mesh") gmsh.model.add("adaptive_mesh")
gmsh.option.setNumber('General.NumThreads', 10) gmsh.option.setNumber("General.NumThreads", 10)
# gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfElements", 200000) # gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfElements", 200000)
# gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfNodes", 200000) # gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfNodes", 200000)
# gmsh.option.setNumber("Mesh.Adapt.MaxIter",5) # gmsh.option.setNumber("Mesh.Adapt.MaxIter",5)
@ -17,7 +17,7 @@ gmsh.option.setNumber("Geometry.Tolerance", 4e-2)
# gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) # gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
lc = 0.04 lc = 0.04
#z_loc = -178.3 # z_loc = -174.3
if len(sys.argv) < 2: if len(sys.argv) < 2:
print("Usage: python3 wires_gmsh2d_bc.py <z_locus in mm>") print("Usage: python3 wires_gmsh2d_bc.py <z_locus in mm>")
@ -26,66 +26,71 @@ if len(sys.argv) < 2:
z_loc = float(sys.argv[1]) z_loc = float(sys.argv[1])
wireShift = 4.0 wireShift = 4.0
k = (2 * np.pi / 24.0) k = 2 * np.pi / 24.0
kg = (2 * np.pi / 48.0) kg = k/2.0
# Plane 1 Offsets (-zmax/2)
# Anodes: -6*k (base) - 4*k (correction) = -10*k
offset_a1 = -6 * k - 4 * k
# Cathodes: -6*k (base) + 0.5*k (half-placement)
offset_c1 = -6 * k + (k / 2.0)
# Guard wires: aligned with cathodes
offset_g1 = offset_c1
# Plane 2 Offsets (+zmax/2) with Twist
offset_a2 = offset_a1 + (wireShift * k)
offset_c2 = offset_c1 - (wireShift * k)
offset_g2 = offset_c2
# 1 needle, 24 ic1, 24 ic2, 48 guard wires, 24 anodes, 24 cathodes # 1 needle, 24 ic1, 24 ic2, 48 guard wires, 24 anodes, 24 cathodes
#needle at plane 1 at -zmax/2 # needle at plane 1 at -zmax/2 no rotation
xarr_needle = np.array([0]) xarr_needle = np.array([0])
yarr_needle = np.array([0]) yarr_needle = np.array([0])
#ic1 wires, plane 1 at -zmax/2 # ic1 wires, plane 1 at -zmax/2 no rotation
ki=2*np.pi/24. xarr_i11 = np.array([23 * np.cos(k * i) for i in range(24)])
xarr_i11 = np.array([23 * np.cos(ki * i) for i in range(24)]) yarr_i11 = np.array([23 * np.sin(k * i) for i in range(24)])
yarr_i11 = np.array([23 * np.sin(ki * i) for i in range(24)])
#ic1 wires, plane 1 at -zmax/2 # ic1 wires, plane 1 at -zmax/2 no rotation
xarr_i21 = np.array([23 * np.cos(ki * i + ki/2.0) for i in range(24)]) xarr_i21 = np.array([23 * np.cos(k * i + k / 2.0) for i in range(24)])
yarr_i21 = np.array([23 * np.sin(ki * i + ki/2.0) for i in range(24)]) yarr_i21 = np.array([23 * np.sin(k * i + k / 2.0) for i in range(24)])
# Guard wires (48 total) - aligned with Cathode phasing # --- Coordinate Arrays (Plane 1) ---
offset_g1 = -5 * k - 2 * k - (np.pi / 24.0) # Anodes: -k*i (Left-handed twist direction)
xarrg_1 = np.array([32 * np.cos(kg * i + offset_g1) for i in range(48)])
yarrg_1 = np.array([32 * np.sin(kg * i + offset_g1) for i in range(48)])
# Anodes (24 total) - index increases leftward (-k)
offset_a1 = -5 * k - 5 * k
xarra_1 = np.array([37 * np.cos(-k * i + offset_a1) for i in range(24)]) xarra_1 = np.array([37 * np.cos(-k * i + offset_a1) for i in range(24)])
yarra_1 = np.array([37 * np.sin(-k * i + offset_a1) for i in range(24)]) yarra_1 = np.array([37 * np.sin(-k * i + offset_a1) for i in range(24)])
# Cathodes (24 total) - index increases rightward (+k) # Cathodes: +k*i (Right-handed twist direction)
offset_c1 = -5 * k - 2 * k - (np.pi / 24.0)
xarrc_1 = np.array([42 * np.cos(k * i + offset_c1) for i in range(24)]) xarrc_1 = np.array([42 * np.cos(k * i + offset_c1) for i in range(24)])
yarrc_1 = np.array([42 * np.sin(k * i + offset_c1) for i in range(24)]) yarrc_1 = np.array([42 * np.sin(k * i + offset_c1) for i in range(24)])
# Guard Wires (48 wires, use kg spacing)
xarrg_1 = np.array([32 * np.cos(kg * i + offset_g1) for i in range(48)])
yarrg_1 = np.array([32 * np.sin(kg * i + offset_g1) for i in range(48)])
# needle at plane 2 at zmax/2 # needle at plane 2 at zmax/2
xarr_needle_2 = np.array([0]) xarr_needle_2 = np.array([0])
yarr_needle_2 = np.array([0]) yarr_needle_2 = np.array([0])
# #ic1 wires, plane 2 at zmax/2 # #ic1 wires, plane 2 at zmax/2
xarr_i12 = np.array([23 * np.cos(ki * i) for i in range(24)]) xarr_i12 = np.array([23 * np.cos(k * i) for i in range(24)])
yarr_i12 = np.array([23 * np.sin(ki * i) for i in range(24)]) yarr_i12 = np.array([23 * np.sin(k * i) for i in range(24)])
# #ic2 wires, plane 2 at zmax/2 # #ic2 wires, plane 2 at zmax/2
xarr_i22 = np.array([23 * np.cos(ki * i + ki/2.0) for i in range(24)]) xarr_i22 = np.array([23 * np.cos(k * i + k / 2.0) for i in range(24)])
yarr_i22 = np.array([23 * np.sin(ki * i + ki/2.0) for i in range(24)]) yarr_i22 = np.array([23 * np.sin(k * i + k / 2.0) for i in range(24)])
# Guard wires (48 total) - twists leftward to match cathodes # --- Coordinate Arrays (Plane 2) ---
offset_g2 = offset_g1 - (wireShift * k)
xarrg_2 = np.array([32 * np.cos(kg * i + offset_g2) for i in range(48)])
yarrg_2 = np.array([32 * np.sin(kg * i + offset_g2) for i in range(48)])
# Anodes (24 total) - twists rightward (+shift)
offset_a2 = offset_a1 + (wireShift * k)
xarra_2 = np.array([37 * np.cos(-k * i + offset_a2) for i in range(24)]) xarra_2 = np.array([37 * np.cos(-k * i + offset_a2) for i in range(24)])
yarra_2 = np.array([37 * np.sin(-k * i + offset_a2) for i in range(24)]) yarra_2 = np.array([37 * np.sin(-k * i + offset_a2) for i in range(24)])
# Cathodes (24 total) - twists leftward (-shift)
offset_c2 = offset_c1 - (wireShift * k)
xarrc_2 = np.array([42 * np.cos(k * i + offset_c2) for i in range(24)]) xarrc_2 = np.array([42 * np.cos(k * i + offset_c2) for i in range(24)])
yarrc_2 = np.array([42 * np.sin(k * i + offset_c2) for i in range(24)]) yarrc_2 = np.array([42 * np.sin(k * i + offset_c2) for i in range(24)])
xarrg_2 = np.array([32 * np.cos(kg * i + offset_g2) for i in range(48)])
yarrg_2 = np.array([32 * np.sin(kg * i + offset_g2) for i in range(48)])
direction_needle_x = xarr_needle_2 - xarr_needle direction_needle_x = xarr_needle_2 - xarr_needle
direction_needle_y = yarr_needle_2 - yarr_needle direction_needle_y = yarr_needle_2 - yarr_needle
@ -104,7 +109,7 @@ direction_anodes_y = yarra_2 - yarra_1
direction_cathodes_x = xarrc_2 - xarrc_1 direction_cathodes_x = xarrc_2 - xarrc_1
direction_cathodes_y = yarrc_2 - yarrc_1 direction_cathodes_y = yarrc_2 - yarrc_1
t = (z_loc+178.3)/(2*178.3) #z=-178.3 is 0, z=+178.3 is 1 t = (z_loc + 174.3) / (2 * 174.3) # z=-174.3 is 0, z=+174.3 is 1
xloc_needle = xarr_needle + t * direction_needle_x xloc_needle = xarr_needle + t * direction_needle_x
yloc_needle = yarr_needle + t * direction_needle_y yloc_needle = yarr_needle + t * direction_needle_y
@ -159,8 +164,12 @@ for i, (xa, ya) in enumerate(zip(xloc_a, yloc_a)):
# Place IC wires only if flag is True # Place IC wires only if flag is True
if include_ic_wires: if include_ic_wires:
i1disk = gmsh.model.occ.addDisk(xloc_i1[i], yloc_i1[i], 0, wire_radius, wire_radius) i1disk = gmsh.model.occ.addDisk(
i2disk = gmsh.model.occ.addDisk(xloc_i2[i], yloc_i2[i], 0, wire_radius, wire_radius) xloc_i1[i], yloc_i1[i], 0, wire_radius, wire_radius
)
i2disk = gmsh.model.occ.addDisk(
xloc_i2[i], yloc_i2[i], 0, wire_radius, wire_radius
)
ic1_wires.append(i1disk) ic1_wires.append(i1disk)
ic2_wires.append(i2disk) ic2_wires.append(i2disk)
@ -172,13 +181,17 @@ gmsh.option.setNumber("Geometry.Tolerance", 1e-6)
gmsh.option.setNumber("Geometry.OCCFixDegenerated", 1) gmsh.option.setNumber("Geometry.OCCFixDegenerated", 1)
gmsh.model.occ.synchronize() gmsh.model.occ.synchronize()
# --- Surface Extraction --- # --- Surface Extraction ---
def get_surfs(disks): def get_surfs(disks):
surfs = [] surfs = []
for d in disks: for d in disks:
surfs += [s[1] for s in gmsh.model.getBoundary([(2,d)], oriented=False) if s[0] == 1] surfs += [
s[1] for s in gmsh.model.getBoundary([(2, d)], oriented=False) if s[0] == 1
]
return surfs return surfs
needle_surfs = get_surfs(needle) if include_needle else [] needle_surfs = get_surfs(needle) if include_needle else []
gwire_surfs = get_surfs(guard_wires) gwire_surfs = get_surfs(guard_wires)
awire_surfs = get_surfs(anode_wires) awire_surfs = get_surfs(anode_wires)
@ -187,7 +200,9 @@ i1wire_surfs = get_surfs(ic1_wires) if include_ic_wires else []
i2wire_surfs = get_surfs(ic2_wires) if include_ic_wires else [] i2wire_surfs = get_surfs(ic2_wires) if include_ic_wires else []
all_active_wire_surfs = needle_surfs + gwire_surfs + awire_surfs + cwire_surfs + i1wire_surfs + i2wire_surfs all_active_wire_surfs = (
needle_surfs + gwire_surfs + awire_surfs + cwire_surfs + i1wire_surfs + i2wire_surfs
)
gmsh.model.mesh.embed(1, all_active_wire_surfs, 2, anasen_barrel) gmsh.model.mesh.embed(1, all_active_wire_surfs, 2, anasen_barrel)
f1 = gmsh.model.mesh.field.add("Distance") f1 = gmsh.model.mesh.field.add("Distance")
@ -226,7 +241,7 @@ gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.model.mesh.generate(dim=2) gmsh.model.mesh.generate(dim=2)
# gmsh.model.mesh.refine() # gmsh.model.mesh.refine()
# gmsh.model.mesh.refine() # gmsh.model.mesh.refine()
gmsh.model.mesh.setOrder(3) gmsh.model.mesh.setOrder(2)
gmsh.write("wires2d.msh") gmsh.write("wires2d.msh")
# gmsh.fltk.run() # gmsh.fltk.run()
gmsh.finalize() gmsh.finalize()

67
pccal/anode_gainmatch.C Normal file
View File

@ -0,0 +1,67 @@
void anode_gainmatch(){
TFile *f = new TFile("../results_run16.root");
TH2F *pc_index_h2d = (TH2F*)(f->Get("hRawPC/PC_Index_Vs_Energy"));
std::cout << pc_index_h2d << std::endl;
TCanvas c("c1","c1",0,0,1600,800);
//TCanvas c_g("cg","cg",0,900,400,400);
c.Divide(2,1);
auto c1=c.cd(1);
pc_index_h2d->Draw("COLZ");
pc_index_h2d->GetYaxis()->SetRangeUser(240,5000);
auto c2=c.cd(2);
c2->SetLogy();
TH1F *h_1d=NULL;
int bin_index=1;
std::vector<std::vector<double>> all_peaks;
std::vector<int> found_wire_list;
while(bin_index<=24) {
h_1d=(TH1F*)(pc_index_h2d->ProjectionY("_py",bin_index,bin_index));
auto c1 = c.cd(1);
TBox box(pc_index_h2d->GetXaxis()->GetBinLowEdge(bin_index),0,pc_index_h2d->GetXaxis()->GetBinUpEdge(bin_index),pc_index_h2d->GetYaxis()->GetXmax());
box.SetFillColorAlpha(kYellow+3,0.3);
box.Draw("SAME");
c1->Modified(); c1->Update();
//while(c1->WaitPrimitive());
TSpectrum s;
auto c2 = c.cd(2);
h_1d->Draw();
c2->Modified(); c2->Update();
int npeaks = s.Search(h_1d,8,"",0.02); std::cout << npeaks << std::endl;
if(npeaks>=3) {
std::vector<double> xpeaks(s.GetPositionX(),s.GetPositionX()+npeaks);
std::sort(xpeaks.begin(),xpeaks.end(),std::greater());
found_wire_list.push_back((int)pc_index_h2d->GetXaxis()->GetBinCenter(bin_index));
all_peaks.push_back(xpeaks);
}
while(c2->WaitPrimitive());
bin_index++;
}
c.cd(2)->SetLogy(kFALSE);
gStyle->SetOptFit(1111);
std::ofstream outfile("anode_gm_coeffs.dat");
outfile << found_wire_list.at(0) << " "
<< 1.0 << " "
<< 0.0 << std::endl;
for(int i=0; i<all_peaks.size(); i++){
if(i==1) continue;
TGraph g(all_peaks.at(i).size(), all_peaks.at(i).data(), all_peaks.at(1).data());
auto c2 = c.cd(2);
g.SetMarkerStyle(20);
//g.Print();
g.Draw("AP");
g.Fit("pol1");
outfile << found_wire_list.at(i) << " "
<< ((TF1*)g.FindObject("pol1"))->GetParameter(1) << " "
<< ((TF1*)g.FindObject("pol1"))->GetParameter(0) << std::endl;
c2->Modified();
c2->Update();
while(c2->WaitPrimitive());
}
outfile.close();
f->Close();
return;
}

24
pccal/anode_gm_coeffs.dat Normal file
View File

@ -0,0 +1,24 @@
0 1 0
0 0.937314 -16.871
2 0.965461 -1.54376
3 0.926501 -3.27662
4 0.905634 2.54577
5 0.905634 -11.0387
6 0.853919 6.23079
7 0.945588 -9.54044
8 0.884454 -11.8262
9 0.922501 -3.42538
10 0.903053 9.28069
11 0.914653 9.87642
12 0.965332 13.2526
13 0.923847 -3.41775
14 0.93845 25.9901
15 0.955424 12.324
16 0.95116 4.99595
17 0.910745 2.86648
18 0.941376 4.57217
19 0.871622 932.111
20 1.00624 7.86358
21 0.969834 -45.001
22 0.89304 -31.5635
23 0.933226 4.02193

69
pccal/cathode_gainmatch.C Normal file
View File

@ -0,0 +1,69 @@
void cathode_gainmatch(){
TFile *f = new TFile("../results_run17.root");
TH2F *pc_index_h2d = (TH2F*)(f->Get("hRawPC/PC_Index_Vs_Energy"));
std::cout << pc_index_h2d << std::endl;
TCanvas c("c1","c1",0,0,1600,800);
//TCanvas c_g("cg","cg",0,900,400,400);
c.Divide(2,1);
auto c1=c.cd(1);
pc_index_h2d->Draw("COLZ");
pc_index_h2d->GetYaxis()->SetRangeUser(600,pc_index_h2d->GetYaxis()->GetXmax());
auto c2=c.cd(2);
c2->SetLogy();
TH1F *h_1d=NULL;
int bin_index=25;
std::vector<double> pulser_heights = {0.01,0.05,0.1,0.15,0.2,0.25,0.3,0.5};
std::vector<std::vector<double>> all_peaks;
std::vector<int> found_wire_list;
while(bin_index<=48) {
h_1d=(TH1F*)(pc_index_h2d->ProjectionY("_py",bin_index,bin_index));
auto c1 = c.cd(1);
TBox box(pc_index_h2d->GetXaxis()->GetBinLowEdge(bin_index),0,pc_index_h2d->GetXaxis()->GetBinUpEdge(bin_index),pc_index_h2d->GetYaxis()->GetXmax());
box.SetFillColorAlpha(kYellow+3,0.3);
box.Draw("SAME");
c1->Modified(); c1->Update();
//while(c1->WaitPrimitive());
TSpectrum s;
auto c2 = c.cd(2);
h_1d->Draw();
c2->Modified(); c2->Update();
int npeaks = s.Search(h_1d,20,"",0.1); std::cout << npeaks << std::endl;
if(npeaks==8) {
std::vector<double> xpeaks(s.GetPositionX(),s.GetPositionX()+npeaks);
for(int i=0; i<8; i++) {
std::cout << pc_index_h2d->GetXaxis()->GetBinCenter(bin_index) << " " << xpeaks.at(i) << " " << xpeaks.at(i)/pulser_heights.at(i) << std::endl;
}
std::sort(xpeaks.begin(),xpeaks.end(),std::greater());
found_wire_list.push_back((int)pc_index_h2d->GetXaxis()->GetBinCenter(bin_index));
all_peaks.push_back(xpeaks);
}
while(c2->WaitPrimitive());
bin_index++;
}
c.cd(2)->SetLogy(kFALSE);
gStyle->SetOptFit(1111);
std::ofstream outfile("cathode_gm_coeffs.dat");
outfile << found_wire_list.at(0) << " "
<< 1.0 << " "
<< 0.0 << std::endl;
for(int i=1; i<all_peaks.size(); i++){
TGraph g(all_peaks.at(i).size(), all_peaks.at(i).data(), all_peaks.at(0).data());
auto c2 = c.cd(2);
g.SetMarkerStyle(20);
//g.Print();
g.Draw("AP");
g.Fit("pol1");
outfile << found_wire_list.at(i) << " "
<< ((TF1*)g.FindObject("pol1"))->GetParameter(1) << " "
<< ((TF1*)g.FindObject("pol1"))->GetParameter(0) << std::endl;
c2->Modified();
c2->Update();
while(c2->WaitPrimitive());
}
outfile.close();
f->Close();
return;
}

21
pccal/cathode_gm_coeffs.dat Normal file
View File

@ -0,0 +1,21 @@
24 1 0
25 0.941896 6.16135
26 0.980284 2.86886
27 0.983166 -3.82952
28 0.978704 -2.89713
29 0.964947 2.25786
30 0.94514 0.925074
31 0.977231 1.6493
32 0.919527 5.82742
33 0.972243 2.88061
34 0.928892 7.61384
35 0.947376 -0.644223
36 0.875342 6.066
38 0.970953 6.262
40 0.918408 -3.27891
41 0.913619 4.11288
42 0.954083 2.21261
43 0.993037 5.48924
45 0.926406 -19.719
46 1.00459 5.14574
47 0.942483 5.54183

49
pccal/pc_gm_coeffs.dat Normal file
View File

@ -0,0 +1,49 @@
#Histogram Number Slope Intercept
0 0.937314 -16.871
1 0 0
2 0.965461 -1.54376
3 0.926501 -3.27662
4 0.905634 2.54577
5 0.905634 -11.0387
6 0.853919 6.23079
7 0.945588 -9.54044
8 0.884454 -11.8262
9 0.922501 -3.42538
10 0.903053 9.28069
11 0.914653 9.87642
12 0.965332 13.2526
13 0.923847 -3.41775
14 0.93845 25.9901
15 0.955424 12.324
16 0.95116 4.99595
17 0.910745 2.86648
18 0.941376 4.57217
19 0.871622 932.111
20 1.00624 7.86358
21 0.969834 -45.001
22 0.89304 -31.5635
23 0.933226 4.02193
24 0 0
25 0.941896 6.16135
26 0.980284 2.86886
27 0.983166 -3.82952
28 0.978704 -2.89713
29 0.964947 2.25786
30 0.94514 0.925074
31 0.977231 1.6493
32 0.919527 5.82742
33 0.972243 2.88061
34 0.928892 7.61384
35 0.947376 -0.644223
36 0.875342 6.066
37 0 0
38 0.970953 6.262
39 0 0
40 0.918408 -3.27891
41 0.913619 4.11288
42 0.954083 2.21261
43 0.993037 5.48924
44 0 0
45 0.926406 -19.719
46 1.00459 5.14574
47 0.942483 5.54183

50
pccal/slope_intercept_26Al.dat Normal file
View File

@ -0,0 +1,50 @@
#Histogram Number Slope Intercept
#Histogram Number Slope Intercept
0 0.931015 -1.35431
1 1 -1.87356e-10
2 0.964185 1.49989
3 0.92638 -1.30621
4 0.905569 1.00834
5 0.901182 0.470903
6 0.853932 3.32687
7 0.942785 1.08887
8 0.878904 -0.0107433
9 0.922662 -2.32259
10 0.903343 8.38332
11 0.914227 6.56108
12 0.961008 23.0982
13 0.920976 5.22104
14 0.936584 31.5073
15 0.959044 5.43267
16 0.95263 -0.404053
17 0.90953 4.82833
18 0.940277 10.3629
19 0.86746 -17.8678
20 1.00683 4.76371
21 0.968342 -43.9496
22 0.892882 -32.0742
23 0.933615 1.10704
24 1 -2.89219e-10
25 0.942098 -0.105169
26 0.980862 -0.732032
27 0.982975 -2.22704
28 0.978815 -1.51477
29 0.965245 -2.19515
30 0.945384 -0.892599
31 0.977408 -0.908592
32 0.919546 3.25464
33 0.972194 2.44956
34 0.92852 5.44745
35 0.947098 1.40531
36 0.875491 -1.13145
37 1 0
38 0.970862 2.86019
39 1 0
40 0.91793 -3.80615
41 0.913897 -2.12964
42 0.954014 -0.760604
43 0.993616 -1.40278
44 1 0
45 0.926169 -21.2016
46 1.00577 -2.14281
47 0.943312 -1.26464

85
scratch/sx3z_vs_phiz/scan_offset.C Normal file
View File

@ -0,0 +1,85 @@
#include "testmodel.h"
int quit=0;
void handler(int){quit=0;}
int colors[] = {kSpring+3, kRed, kGreen+3, kBlue+3, kViolet, kOrange, kSpring-7, kAzure-5};
void scan_offset(){
signal(SIGINT,handler);
TCanvas c("c1","c1",0,0,1600,800);
c.Divide(2,1);
TF1 f1("model",model,-200,200,2);
f1.SetNpx(10000);
std::vector<double> pars = {0.0,1.};
f1.SetParameters(pars.data());
f1.SetLineColor(kGreen+2);
f1.SetLineStyle(kLine);
TFile* f=NULL;
std::vector<TFile*> files;
int ctr=0;
for(int i=12; i<=21; i++) {
auto c1=c.cd(1);
c1->SetGrid(1,1);
f = new TFile(Form("../../results_run%d.root",i));
if(i==12) {
//TH2F *h2 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int"));
TH2F *h23 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int_A1C2"));
std::cout << "aaa" << h23 << std::endl;
h23->SetLineColorAlpha(kOrange,0.75);
h23->GetYaxis()->SetRangeUser(-200,200);
h23->Draw("box");
while(gPad->WaitPrimitive());
} else {
//TH2F *h2 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int"));
//TH2F *h2 = (TH2F*)(f->Get("pcz_vs_sx3pczguess_A1C2_strip12"));
TH2F *h2 = (TH2F*)(f->Get("pcz_vs_sx3pczguess_A1C2"));
std::cout << h2 << std::endl;
//TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2"));
if(!h2) continue;
h2->SetTitle(Form("case%d",i));
//h2->Draw("colz same");
h2->SetLineColorAlpha(colors[ctr],0.75);
h2->Draw("box same");
f1.Draw("same");
}
TF1 eqline("x","x",-200,200);
eqline.Draw("SAME");
c1->Modified();
c1->Update();
ctr+=1;
auto c2=c.cd(2);
c2->SetGrid(1,1);
TH2F *h3 = (TH2F*)(f->Get("sx3phi_vs_pcphi1"));
// TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2"));
if(!h3) continue;
h3->SetTitle(Form("case%d",i));
h3->Draw("colz");
eqline.Draw("SAME");
c2->Modified();
c2->Update();
while(gPad->WaitPrimitive());
files.emplace_back(f);
std::cout <<"Test" << std::endl;
if(i==21) {
i=11;
c.Clear();
c.Divide(2,1);
ctr=0;
}
//if(quit) break;
}
for(auto file : files) {
file->Close();
}
}

81
scratch/sx3z_vs_phiz/scan_offset_fix.C Normal file
View File

@ -0,0 +1,81 @@
#include "testmodel.h"
int quit=0;
void handler(int){quit=1;}
int colors[] = {kSpring+3, kRed, kGreen+3, kBlue+3, kViolet, kOrange, kSpring-7, kAzure-5};
void scan_offset_fix(){
signal(SIGINT,handler);
TCanvas c("c1","c1",0,0,1600,800);
c.Divide(2,1);
TF1 f1("model",model,-200,200,2);
f1.SetNpx(10000);
std::vector<double> pars = {0.0,1.};
f1.SetParameters(pars.data());
f1.SetLineColor(kGreen+2);
f1.SetLineStyle(kLine);
TFile* f=NULL;
std::vector<TFile*> files;
int ctr=0;
for(int i=12; i<=21; i++) {
if(i==15) continue;
auto c1=c.cd(1);
c1->SetGrid(1,1);
f = new TFile(Form("../../results_run%d.root",i));
if(i==12) {
//TH2F *h2 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int"));
TH2F *h23 = (TH2F*)(f->Get("pczfix_vs_qqqpczguess_A1C2"));
h23->SetLineColorAlpha(kOrange,0.75);
h23->Draw("box SAME");
} else {
//TH2F *h2 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int"));
//TH2F *h2 = (TH2F*)(f->Get("pcz_vs_sx3pczguess_A1C2_strip12"));
TH2F *h2 = (TH2F*)(f->Get("pczfix_vs_sx3pczguess_A1C2"));
//TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2"));
if(!h2) continue;
h2->SetTitle(Form("case%d",i));
//h2->Draw("colz same");
h2->SetLineColorAlpha(colors[ctr],0.75);
h2->Draw("box same");
//f1.Draw("same");
}
TF1 eqline("x","x",-200,200);
eqline.Draw("SAME");
c1->Modified();
c1->Update();
ctr+=1;
auto c2=c.cd(2);
c2->SetGrid(1,1);
TH2F *h3 = (TH2F*)(f->Get("sx3phi_vs_pcphi1"));
// TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2"));
if(!h3) continue;
h3->SetTitle(Form("case%d",i));
h3->Draw("colz");
eqline.Draw("SAME");
c2->Modified();
c2->Update();
while(gPad->WaitPrimitive());
files.emplace_back(f);
if(i==21) {
i=11;
c.Clear();
c.Divide(2,1);
ctr=0;
}
if(quit) break;
}
for(auto file : files) {
file->Close();
}
}