From bb05baf89d72a881aab312b08b7bd1ff48545a6c Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Tue, 19 May 2026 17:57:54 -0400 Subject: [PATCH] modified: Armory/ClassPW.h changes made to explain the andoe and cathode shifts in a more intuitive geometric manner 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 --- Armory/ClassPW.h | 6 +- MakeVertex.C | 66 ++++----- anasen_fem/paraview_plotter.py | 2 +- anasen_fem/run.py | 10 +- anasen_fem/scalars.dat.names | 2 +- anasen_fem/wires_gmsh2d_bc.py | 181 +++++++++++++------------ pccal/anode_gainmatch.C | 67 +++++++++ pccal/anode_gm_coeffs.dat | 24 ++++ pccal/cathode_gainmatch.C | 69 ++++++++++ pccal/cathode_gm_coeffs.dat | 21 +++ pccal/pc_gm_coeffs.dat | 49 +++++++ pccal/slope_intercept_26Al.dat | 50 +++++++ scratch/sx3z_vs_phiz/scan_offset.C | 85 ++++++++++++ scratch/sx3z_vs_phiz/scan_offset_fix.C | 81 +++++++++++ 14 files changed, 587 insertions(+), 126 deletions(-) create mode 100644 pccal/anode_gainmatch.C create mode 100644 pccal/anode_gm_coeffs.dat create mode 100644 pccal/cathode_gainmatch.C create mode 100644 pccal/cathode_gm_coeffs.dat create mode 100644 pccal/pc_gm_coeffs.dat create mode 100644 pccal/slope_intercept_26Al.dat create mode 100644 scratch/sx3z_vs_phiz/scan_offset.C create mode 100644 scratch/sx3z_vs_phiz/scan_offset_fix.C diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h index 2ed454b..7af4985 100755 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -153,9 +153,9 @@ inline void PW::ConstructGeo() std::pair p1; // anode std::pair q1; // cathode - double k = TMath::TwoPi() / 24.; // 48 solder thru holes, wires in every other one - double offset_a1 = -5 * k - 5 * k; - double offset_c1 = -5 * k - 2 * k - TMath::TwoPi() / 48; // correct for a half-turn + double k = TMath::TwoPi() / 24.; // 48 solder thru holes, wires in every other one + 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 = -6 * k + k / 2.0; // correct for a half-turn // std::cerr << "Here!" << std::endl; // #include "../scratch/testing.h" double offset_a2 = offset_a1 + wireShift * k; diff --git a/MakeVertex.C b/MakeVertex.C index 6b1ebd6..2e5662d 100755 --- a/MakeVertex.C +++ b/MakeVertex.C @@ -3,8 +3,8 @@ // Comment out any line below to disable that diagnostic section #define DIAG_1WIRE // raw per-wire dE vs Si (no PC required) #define DIAG_PC_SX3 // PC-SX3 coincidence analysis -// #define DIAG_NA0C_SX3 // nA0C (n>=1) pseudo-wire vertex with SX3 -// #define DIAG_NA0C_QQQ // nA0C (n>=1) pseudo-wire vertex with QQQ +#define DIAG_NA_SX3 // nA (n>=1) pseudo-wire vertex with SX3 +#define DIAG_NA_QQQ // nA (n>=1) pseudo-wire vertex with QQQ #define DIAG_PC_QQQ // PC-QQQ coincidence analysis 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 (cClusters.size() == 0 && aClusters.size() >= 1) +#if defined(DIAG_NA_SX3) || defined(DIAG_NA_QQQ) + 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 std::vector> allAnodeHits; @@ -1284,7 +1284,7 @@ Bool_t MakeVertex::Process(Long64_t entry) auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(allAnodeHits, "ANODE"); -#ifdef DIAG_NA0C_SX3 +#ifdef DIAG_NA_SX3 for (auto sx3event : SX3_Events) { 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_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->Fill1D(nA0C_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label); - plotter->Fill1D(nA0C_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_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(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 + "_dE_Ecorr_Anode_SX3", 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA_label); + plotter->Fill1D(nA_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA_label); + plotter->Fill1D(nA_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA_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(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 != "") { - plotter->Fill1D(nA0C_label + "_twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA0C_label); - plotter->Fill1D(nA0C_label + "_twisted_vertex_recon_SX3" + vtx_gate, 600, -300, 300, vertex_recon, nA0C_label); - plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_SX3" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA0C_label); - plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_SX3_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efixalpha, apSumE * sinTheta, nA0C_label); - plotter->Fill1D(nA0C_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label); - plotter->Fill1D(nA0C_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label); + plotter->Fill1D(nA_label + "_twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA_label); + plotter->Fill1D(nA_label + "_twisted_vertex_recon_SX3" + vtx_gate, 600, -300, 300, vertex_recon, nA_label); + plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_SX3" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA_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(nA_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA_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) { 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_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->Fill1D(nA0C_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label); - plotter->Fill1D(nA0C_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_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(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 + "_dE_Ecorr_Anode_QQQ", 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA_label); + plotter->Fill1D(nA_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA_label); + plotter->Fill1D(nA_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA_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(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 != "") { - plotter->Fill1D(nA0C_label + "_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA0C_label); - plotter->Fill1D(nA0C_label + "_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, nA0C_label); - plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA0C_label); - plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfixalpha, apSumE * sinTheta, nA0C_label); - plotter->Fill1D(nA0C_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label); - plotter->Fill1D(nA0C_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label); + plotter->Fill1D(nA_label + "_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA_label); + plotter->Fill1D(nA_label + "_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, nA_label); + plotter->Fill2D(nA_label + "_dE_Ecorr_Anode_QQQ" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA_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(nA_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA_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 for (auto pcevent : PC_Events) diff --git a/anasen_fem/paraview_plotter.py b/anasen_fem/paraview_plotter.py index e2762b4..d492ced 100755 --- a/anasen_fem/paraview_plotter.py +++ b/anasen_fem/paraview_plotter.py @@ -95,7 +95,7 @@ glyph.ScaleArray = ['POINTS', 'No scale array'] glyph.ScaleFactor = 0.001 glyph.GlyphMode = 'Every Nth Point' -glyph.Stride = 24 +glyph.Stride = 32 # --- 3. Display the Glyphs --- glyph_display = Show(glyph, renderView) diff --git a/anasen_fem/run.py b/anasen_fem/run.py index f26f85f..1e886a6 100755 --- a/anasen_fem/run.py +++ b/anasen_fem/run.py @@ -1,10 +1,10 @@ import code import os -# val=-178.3 -val=89.15 +# val=-174.3 +val=0 count=11 -while val<178.3+0.1: +while val<174.3+0.1: print(val) os.system("python3 wires_gmsh2d_bc.py "+str(val)) 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 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)) - val=val+17.83 + val=val+17.43 count = count + 1 - break + # break # os.system("tar -cvzf wiress2d/mesh.tar.gz wires2d/mesh_files") # os.system("rm -rf wires2d/mesh_files/*") diff --git a/anasen_fem/scalars.dat.names b/anasen_fem/scalars.dat.names index 44cd88c..a1c6859 100755 --- a/anasen_fem/scalars.dat.names +++ b/anasen_fem/scalars.dat.names @@ -2,7 +2,7 @@ Metadata for SaveScalars file: ./scalars.dat Elmer version: 26.2 Elmer compilation date: 2026-05-14 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: 1: res: potential difference diff --git a/anasen_fem/wires_gmsh2d_bc.py b/anasen_fem/wires_gmsh2d_bc.py index 0c39847..be212a9 100755 --- a/anasen_fem/wires_gmsh2d_bc.py +++ b/anasen_fem/wires_gmsh2d_bc.py @@ -1,23 +1,23 @@ import numpy as np -import gmsh,sys +import gmsh, sys # --- Configuration Flags --- -include_ic_wires = True +include_ic_wires = True include_needle = True gmsh.initialize() gmsh.model.add("adaptive_mesh") -gmsh.option.setNumber('General.NumThreads', 10) -#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfElements", 200000) -#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfNodes", 200000) -#gmsh.option.setNumber("Mesh.Adapt.MaxIter",5) -#gmsh.option.setNumber("Mesh.MeshSizeMin", 5e-3) -#gmsh.option.setNumber("Mesh.MeshSizeMax", 10.0) +gmsh.option.setNumber("General.NumThreads", 10) +# gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfElements", 200000) +# gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfNodes", 200000) +# gmsh.option.setNumber("Mesh.Adapt.MaxIter",5) +# gmsh.option.setNumber("Mesh.MeshSizeMin", 5e-3) +# gmsh.option.setNumber("Mesh.MeshSizeMax", 10.0) gmsh.option.setNumber("Geometry.Tolerance", 4e-2) -#gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) +# gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) lc = 0.04 -#z_loc = -178.3 +# z_loc = -174.3 if len(sys.argv) < 2: print("Usage: python3 wires_gmsh2d_bc.py ") @@ -26,66 +26,71 @@ if len(sys.argv) < 2: z_loc = float(sys.argv[1]) wireShift = 4.0 -k = (2 * np.pi / 24.0) -kg = (2 * np.pi / 48.0) +k = 2 * np.pi / 24.0 +kg = k/2.0 -#1 needle, 24 ic1, 24 ic2, 48 guard wires, 24 anodes, 24 cathodes +# 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 -#needle at plane 1 at -zmax/2 +# 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 + +# needle at plane 1 at -zmax/2 no rotation xarr_needle = np.array([0]) yarr_needle = np.array([0]) -#ic1 wires, plane 1 at -zmax/2 -ki=2*np.pi/24. -xarr_i11 = np.array([23 * np.cos(ki * 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 no rotation +xarr_i11 = np.array([23 * np.cos(k * i) for i in range(24)]) +yarr_i11 = np.array([23 * np.sin(k * i) for i in range(24)]) -#ic1 wires, plane 1 at -zmax/2 -xarr_i21 = np.array([23 * np.cos(ki * i + ki/2.0) for i in range(24)]) -yarr_i21 = np.array([23 * np.sin(ki * i + ki/2.0) for i in range(24)]) +# ic1 wires, plane 1 at -zmax/2 no rotation +xarr_i21 = np.array([23 * np.cos(k * i + k / 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 -offset_g1 = -5 * k - 2 * k - (np.pi / 24.0) -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 +# --- Coordinate Arrays (Plane 1) --- +# Anodes: -k*i (Left-handed twist direction) 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)]) -# Cathodes (24 total) - index increases rightward (+k) -offset_c1 = -5 * k - 2 * k - (np.pi / 24.0) +# Cathodes: +k*i (Right-handed twist direction) 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)]) -#needle at plane 2 at zmax/2 +# 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 xarr_needle_2 = np.array([0]) yarr_needle_2 = np.array([0]) # #ic1 wires, plane 2 at zmax/2 -xarr_i12 = np.array([23 * np.cos(ki * i) for i in range(24)]) -yarr_i12 = np.array([23 * np.sin(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(k * i) for i in range(24)]) # #ic2 wires, plane 2 at zmax/2 -xarr_i22 = np.array([23 * np.cos(ki * i + ki/2.0) for i in range(24)]) -yarr_i22 = np.array([23 * np.sin(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(k * i + k / 2.0) for i in range(24)]) -# Guard wires (48 total) - twists leftward to match cathodes -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) +# --- Coordinate Arrays (Plane 2) --- 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)]) -# 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)]) 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_y = yarr_needle_2 - yarr_needle @@ -104,81 +109,89 @@ direction_anodes_y = yarra_2 - yarra_1 direction_cathodes_x = xarrc_2 - xarrc_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 -yloc_needle = yarr_needle + t*direction_needle_y -xloc_i1 = xarr_i11 + t*direction_ic1_x -yloc_i1 = yarr_i11 + t*direction_ic1_y -xloc_i2 = xarr_i21 + t*direction_ic2_x -yloc_i2 = yarr_i21 + t*direction_ic2_y -xloc_g = xarrg_1 + t*direction_guard_x -yloc_g = yarrg_1 + t*direction_guard_y -xloc_a = xarra_1 + t*direction_anodes_x -yloc_a = yarra_1 + t*direction_anodes_y -xloc_c = xarrc_1 + t*direction_cathodes_x -yloc_c = yarrc_1 + t*direction_cathodes_y +xloc_needle = xarr_needle + t * direction_needle_x +yloc_needle = yarr_needle + t * direction_needle_y +xloc_i1 = xarr_i11 + t * direction_ic1_x +yloc_i1 = yarr_i11 + t * direction_ic1_y +xloc_i2 = xarr_i21 + t * direction_ic2_x +yloc_i2 = yarr_i21 + t * direction_ic2_y +xloc_g = xarrg_1 + t * direction_guard_x +yloc_g = yarrg_1 + t * direction_guard_y +xloc_a = xarra_1 + t * direction_anodes_x +yloc_a = yarra_1 + t * direction_anodes_y +xloc_c = xarrc_1 + t * direction_cathodes_x +yloc_c = yarrc_1 + t * direction_cathodes_y # wire_radius_a = 0.018 #mm # wire_radius_c = 0.0762 #mm # wire_radius_g = 0.0762 #mm -wire_radius = 0.254 #mm -needle=[] +wire_radius = 0.254 # mm +needle = [] ic1_wires = [] ic2_wires = [] guard_wires = [] anode_wires = [] cathode_wires = [] -iw1_tags = [(3,i) for i in range(24)] -iw2_tags = [(3,i+24) for i in range(24)] -gw_tags = [(3,i+48) for i in range(48)] -aw_tags = [(3,i) for i in range(24)] -cw_tags = [(3,i+24) for i in range(24)] +iw1_tags = [(3, i) for i in range(24)] +iw2_tags = [(3, i + 24) for i in range(24)] +gw_tags = [(3, i + 48) for i in range(48)] +aw_tags = [(3, i) for i in range(24)] +cw_tags = [(3, i + 24) for i in range(24)] -#for i,[xa,ya,xc,yc] in enumerate(zip(xarra_1,yarra_1,xarrc_1,yarrc_1)): -#create Hot Needle (1 total) +# for i,[xa,ya,xc,yc] in enumerate(zip(xarra_1,yarra_1,xarrc_1,yarrc_1)): +# create Hot Needle (1 total) for i, (xn, yn) in enumerate(zip(xloc_needle, yloc_needle)): if include_needle: ndisk = gmsh.model.occ.addDisk(xn, yn, 0, wire_radius, wire_radius) - needle.append(ndisk) + needle.append(ndisk) -#create Guard Wires (48 total) +# create Guard Wires (48 total) for xg, yg in zip(xloc_g, yloc_g): gdisk = gmsh.model.occ.addDisk(xg, yg, 0, wire_radius, wire_radius) guard_wires.append(gdisk) -#create Cathode Wires (24 total) +# create Cathode Wires (24 total) for xc, yc in zip(xloc_c, yloc_c): cdisk = gmsh.model.occ.addDisk(xc, yc, 0, wire_radius, wire_radius) cathode_wires.append(cdisk) -#create IC Anode and Cathode Wires (24 total each) +# create IC Anode and Cathode Wires (24 total each) for i, (xa, ya) in enumerate(zip(xloc_a, yloc_a)): adisk = gmsh.model.occ.addDisk(xa, ya, 0, wire_radius, wire_radius) anode_wires.append(adisk) - + # Place IC wires only if flag is True if include_ic_wires: - i1disk = gmsh.model.occ.addDisk(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) + i1disk = gmsh.model.occ.addDisk( + 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) ic2_wires.append(i2disk) -anasen_barrel = gmsh.model.occ.addDisk(0,0,0,500,500) -#gmsh.model.occ.synchronize() -#gmsh.model.mesh.embed(1,anode_wires+cathode_wires,2,anasen_barrel) +anasen_barrel = gmsh.model.occ.addDisk(0, 0, 0, 500, 500) +# gmsh.model.occ.synchronize() +# gmsh.model.mesh.embed(1,anode_wires+cathode_wires,2,anasen_barrel) gmsh.option.setNumber("Geometry.Tolerance", 1e-6) gmsh.option.setNumber("Geometry.OCCFixDegenerated", 1) gmsh.model.occ.synchronize() + # --- Surface Extraction --- def get_surfs(disks): surfs = [] 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 + needle_surfs = get_surfs(needle) if include_needle else [] gwire_surfs = get_surfs(guard_wires) awire_surfs = get_surfs(anode_wires) @@ -187,17 +200,19 @@ i1wire_surfs = get_surfs(ic1_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) - + f1 = gmsh.model.mesh.field.add("Distance") gmsh.model.mesh.field.setNumbers(f1, "CurvesList", all_active_wire_surfs) f2 = gmsh.model.mesh.field.add("Threshold") gmsh.model.mesh.field.setNumber(f2, "InField", f1) -gmsh.model.mesh.field.setNumber(f2, "SizeMin", 0.05) # Fine mesh near wires +gmsh.model.mesh.field.setNumber(f2, "SizeMin", 0.05) # Fine mesh near wires gmsh.model.mesh.field.setNumber(f2, "SizeMax", 5.0) # Large mesh in empty space -gmsh.model.mesh.field.setNumber(f2, "DistMin", 0.5) # Apply SizeMin within 1mm +gmsh.model.mesh.field.setNumber(f2, "DistMin", 0.5) # Apply SizeMin within 1mm gmsh.model.mesh.field.setNumber(f2, "DistMax", 15.0) # Transition to SizeMax by 20mm gmsh.model.mesh.field.setAsBackgroundMesh(f2) @@ -226,7 +241,7 @@ gmsh.option.setNumber("Mesh.Algorithm", 6) gmsh.model.mesh.generate(dim=2) # gmsh.model.mesh.refine() # gmsh.model.mesh.refine() -gmsh.model.mesh.setOrder(3) +gmsh.model.mesh.setOrder(2) gmsh.write("wires2d.msh") -#gmsh.fltk.run() +# gmsh.fltk.run() gmsh.finalize() diff --git a/pccal/anode_gainmatch.C b/pccal/anode_gainmatch.C new file mode 100644 index 0000000..fa267be --- /dev/null +++ b/pccal/anode_gainmatch.C @@ -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> all_peaks; + std::vector 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 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; iGetParameter(1) << " " + << ((TF1*)g.FindObject("pol1"))->GetParameter(0) << std::endl; + c2->Modified(); + c2->Update(); + while(c2->WaitPrimitive()); + } + outfile.close(); + f->Close(); + return; +} diff --git a/pccal/anode_gm_coeffs.dat b/pccal/anode_gm_coeffs.dat new file mode 100644 index 0000000..a241cfe --- /dev/null +++ b/pccal/anode_gm_coeffs.dat @@ -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 diff --git a/pccal/cathode_gainmatch.C b/pccal/cathode_gainmatch.C new file mode 100644 index 0000000..d52e32a --- /dev/null +++ b/pccal/cathode_gainmatch.C @@ -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 pulser_heights = {0.01,0.05,0.1,0.15,0.2,0.25,0.3,0.5}; + std::vector> all_peaks; + std::vector 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 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; iGetParameter(1) << " " + << ((TF1*)g.FindObject("pol1"))->GetParameter(0) << std::endl; + c2->Modified(); + c2->Update(); + while(c2->WaitPrimitive()); + } + outfile.close(); + f->Close(); + return; +} diff --git a/pccal/cathode_gm_coeffs.dat b/pccal/cathode_gm_coeffs.dat new file mode 100644 index 0000000..b1aaa77 --- /dev/null +++ b/pccal/cathode_gm_coeffs.dat @@ -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 diff --git a/pccal/pc_gm_coeffs.dat b/pccal/pc_gm_coeffs.dat new file mode 100644 index 0000000..92d002f --- /dev/null +++ b/pccal/pc_gm_coeffs.dat @@ -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 diff --git a/pccal/slope_intercept_26Al.dat b/pccal/slope_intercept_26Al.dat new file mode 100644 index 0000000..875aaec --- /dev/null +++ b/pccal/slope_intercept_26Al.dat @@ -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 diff --git a/scratch/sx3z_vs_phiz/scan_offset.C b/scratch/sx3z_vs_phiz/scan_offset.C new file mode 100644 index 0000000..f6ea1b0 --- /dev/null +++ b/scratch/sx3z_vs_phiz/scan_offset.C @@ -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 pars = {0.0,1.}; + f1.SetParameters(pars.data()); + f1.SetLineColor(kGreen+2); + f1.SetLineStyle(kLine); + + + + + TFile* f=NULL; + std::vector 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(); + } +} diff --git a/scratch/sx3z_vs_phiz/scan_offset_fix.C b/scratch/sx3z_vs_phiz/scan_offset_fix.C new file mode 100644 index 0000000..20dde89 --- /dev/null +++ b/scratch/sx3z_vs_phiz/scan_offset_fix.C @@ -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 pars = {0.0,1.}; + f1.SetParameters(pars.data()); + f1.SetLineColor(kGreen+2); + f1.SetLineStyle(kLine); + + + + + TFile* f=NULL; + std::vector 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(); + } +}