From b69dcf39d6899c95ea981b64b4e9e4cfc00c4c1b Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Fri, 15 May 2026 18:01:31 -0400 Subject: [PATCH] modified: MakeVertex.C modified: README.md modified: anasen_fem/run.py modified: anasen_fem/scalars.dat.names modified: anasen_fem/wires_gmsh2d_bc.py modified: run_sx3.sh --- MakeVertex.C | 286 ++++++++++++++++++++++++---------- README.md | 29 ++-- anasen_fem/run.py | 2 +- anasen_fem/scalars.dat.names | 6 +- anasen_fem/wires_gmsh2d_bc.py | 20 ++- run_sx3.sh | 4 +- 6 files changed, 240 insertions(+), 107 deletions(-) diff --git a/MakeVertex.C b/MakeVertex.C index 923f32d..0db0c50 100755 --- a/MakeVertex.C +++ b/MakeVertex.C @@ -1,13 +1,5 @@ #define MakeVertex_cxx -// Comment out any line below to disable that diagnostic section -// #define DIAG_WIREMULT //] anode/cathode cluster multiplicity plots -#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_PC_QQQ // PC-QQQ coincidence analysis - Int_t colors[40] = { kBlack, kRed, kGreen, kBlue, kYellow, kMagenta, kCyan, kOrange, kSpring, kTeal, kAzure, kViolet, kPink, kGray, kWhite, @@ -859,7 +851,6 @@ Bool_t MakeVertex::Process(Long64_t entry) if (QQQ_Events.size() && PC_Events.size()) plotter->Fill2D("PCEv_vs_QQQEv", 20, 0, 20, 20, 0, 20, QQQ_Events.size(), PC_Events.size()); -#ifdef DIAG_WIREMULT plotter->Fill2D("ac_vs_cc", 20, 0, 20, 20, 0, 20, aClusters.size(), cClusters.size(), "wiremult"); for (auto cluster : aClusters) { @@ -874,9 +865,7 @@ Bool_t MakeVertex::Process(Long64_t entry) { plotter->Fill2D("ac_vs_cc_ign0", 20, 0, 20, 20, 0, 20, aClusters.size(), cClusters.size(), "wiremult"); } -#endif // DIAG_WIREMULT -#ifdef DIAG_1WIRE for (auto sx3event : SX3_Events) { for (int i = 0; i < 24; i++) @@ -951,9 +940,7 @@ Bool_t MakeVertex::Process(Long64_t entry) } } // for 'i' loop } -#endif // DIAG_1WIRE -#ifdef DIAG_PC_SX3 bool PCSX3PhiCut = false; for (auto pcevent : PC_Events) { @@ -1238,7 +1225,6 @@ Bool_t MakeVertex::Process(Long64_t entry) }*/ } } // end PC-SX3 coincidence -#endif // DIAG_PC_SX3 /*for(size_t ii=0; ii= 1) + if (aClusters.size() == 1 && cClusters.size() == 0) { - std::string nA0C_label = std::to_string(aClusters.size()) + "A0C"; + // Extract the primary anode hit properties + auto anodeHit = aClusters.front().front(); + int aWireID = std::get<0>(anodeHit); + double aEnergy = std::get<1>(anodeHit); + double aTime = std::get<2>(anodeHit); - // Flatten all anode clusters into a combined hit list for the pseudo-wire - std::vector> allAnodeHits; - for (const auto &acluster : aClusters) - for (const auto &hit : acluster) - allAnodeHits.push_back(hit); + // Get the 3D endpoints of the fired twisted anode wire from your geometry class + TVector3 a1 = pwinstance.An[aWireID].first; + TVector3 wireVec = pwinstance.An[aWireID].first - pwinstance.An[aWireID].second; - auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(allAnodeHits, "ANODE"); - -#ifdef DIAG_NA0C_SX3 + // Loop over SX3_Events directly for (auto sx3event : SX3_Events) { - if (sx3event.Time1 - apTSMaxE < -150) - { - TVector3 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, sx3event.pos.Phi()); + if (sx3event.Time1 - aTime < -150) // Time cut for protons + { + // 1. Define the plane of the track (Z-axis to SX3 hit) + TVector3 planeNormal(-TMath::Sin(sx3event.pos.Phi()), TMath::Cos(sx3event.pos.Phi()), 0.0); + + // 2. Find intersection of the twisted wire with this track plane + double dot_wireVec = wireVec.Dot(planeNormal); + + // Prevent divide-by-zero if wire is perfectly parallel to the track plane + if (TMath::Abs(dot_wireVec) < 1e-6) + continue; + + double t_intersect = -(a1.Dot(planeNormal)) / dot_wireVec; + + // Calculate the exact 3D coordinate on the wire that matches the SX3 phi + TVector3 pcz_intersect = a1 + t_intersect * wireVec; + + // 3. Reconstruct Vertex Z double deltaRho = sx3event.pos.Perp() - pcz_intersect.Perp(); double deltaZ = sx3event.pos.Z() - pcz_intersect.Z(); double vertex_recon = sx3event.pos.Z() - sx3event.pos.Perp() * (deltaZ / deltaRho); @@ -1308,19 +1308,33 @@ Bool_t MakeVertex::Process(Long64_t entry) Kinematics apkin_a(26.981538, 4.002603, 4.002603, 26.981538, beam_energy_at_vertex); std::string vtx_gate = ""; - if (vertex_recon >= -176.0 && vertex_recon < -100.0) - vtx_gate = "_Z[-176_to_-100]"; - else if (vertex_recon >= -100.0 && vertex_recon < -50.0) - vtx_gate = "_Z[-100_to_-50]"; - else if (vertex_recon >= -50.0 && vertex_recon < 0.0) - vtx_gate = "_Z[-50_to_0]"; - else if (vertex_recon >= 0.0 && vertex_recon < 50.0) - vtx_gate = "_Z[0_to_50]"; - else if (vertex_recon >= 50.0 && vertex_recon < 100.0) - vtx_gate = "_Z[50_to_100]"; - else if (vertex_recon >= 100.0 && vertex_recon < 176.0) - vtx_gate = "_Z[100_to_176]"; + if (vertex_recon >= -176.0 && vertex_recon < -100.0) + { + vtx_gate = "_Z[-176_to_-100]"; + } + else if (vertex_recon >= -100.0 && vertex_recon < -50.0) + { + vtx_gate = "_Z[-100_to_-50]"; + } + else if (vertex_recon >= -50.0 && vertex_recon < 0.0) + { + vtx_gate = "_Z[-50_to_0]"; + } + else if (vertex_recon >= 0.0 && vertex_recon < 50.0) + { + vtx_gate = "_Z[0_to_50]"; + } + else if (vertex_recon >= 50.0 && vertex_recon < 100.0) + { + vtx_gate = "_Z[50_to_100]"; + } + else if (vertex_recon >= 100.0 && vertex_recon < 176.0) + { + vtx_gate = "_Z[100_to_176]"; + } + + // 4. Energy Loss Correction in Silicon double path_length = (sx3event.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1; double sx3Efix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1) - path_length); double sx3Efixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(sx3event.Energy1) - path_length); @@ -1328,37 +1342,59 @@ Bool_t MakeVertex::Process(Long64_t entry) double theta_recon = (sx3event.pos - TVector3(0, 0, vertex_recon)).Theta(); double sinTheta = TMath::Sin(theta_recon); + // Now these functions will use the correct, event-specific beam energy! 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); + // 5. Fill Diagnostics + plotter->Fill2D("1A0C_dE_Ecorr_Anode_SX3", 400, 0, 30, 800, 0, 40000, sx3Efix, aEnergy * sinTheta, "1A0C"); + plotter->Fill1D("1A0C_Ex_from_protons_SX3", 200, -10, 10, Ex_from_proton, "1A0C"); + plotter->Fill1D("1A0C_Ex_from_alphas_SX3", 200, -10, 10, Ex_from_alpha, "1A0C"); + + plotter->Fill2D("1A0C_sx3_E_vs_theta_raw_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3event.Energy1, "1A0C"); + plotter->Fill2D("1A0C_sx3_E_vs_theta_corr_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3Efix, "1A0C"); 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("1A0C_twisted_vertex_recon_SX3" + vtx_gate, 600, -300, 300, vertex_recon, "1A0C"); + plotter->Fill1D("1A0C_twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), "1A0C"); + + plotter->Fill2D("1A0C_dE_Ecorr_Anode_SX3" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efix, aEnergy * sinTheta, "1A0C"); + + plotter->Fill1D("1A0C_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, "1A0C"); + plotter->Fill1D("1A0C_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, "1A0C"); + + // Track where on the wire the hit occurred (0 to 1 is inside the physical PC) + // plotter->Fill1D("1A0C_wire_t_parameter" + vtx_gate, 200, -0.5, 1.5, t_intersect, "1A0C"); } } } -#endif // DIAG_NA0C_SX3 -#ifdef DIAG_NA0C_QQQ + // Loop over QQQ_Events directly + for (auto qqqevent : QQQ_Events) { - if (qqqevent.Time1 - apTSMaxE < -150) + if (qqqevent.Time1 - aTime < -150) // Time cut for protons { - TVector3 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, qqqevent.pos.Phi()); + // 1. Define the plane of the track (Z-axis to SX3 hit) + TVector3 planeNormal(-TMath::Sin(qqqevent.pos.Phi()), TMath::Cos(qqqevent.pos.Phi()), 0.0); + // 2. Find intersection of the twisted wire with this track plane + double dot_wireVec = wireVec.Dot(planeNormal); + + // Prevent divide-by-zero if wire is perfectly parallel to the track plane + if (TMath::Abs(dot_wireVec) < 1e-6) + continue; + + double t_intersect_QQQ = -(a1.Dot(planeNormal)) / dot_wireVec; + + // Calculate the exact 3D coordinate on the wire that matches the SX3 phi + TVector3 pcz_intersect = a1 + t_intersect_QQQ * wireVec; + + // 3. Reconstruct Vertex Z double deltaRho = qqqevent.pos.Perp() - pcz_intersect.Perp(); double deltaZ = qqqevent.pos.Z() - pcz_intersect.Z(); + double vertex_recon = qqqevent.pos.Z() - qqqevent.pos.Perp() * (deltaZ / deltaRho); double z_entrance = -274.3; @@ -1369,22 +1405,38 @@ Bool_t MakeVertex::Process(Long64_t entry) Kinematics apkin_p(26.981538, 4.002603, 1.007825, 29.973770, beam_energy_at_vertex); Kinematics apkin_a(26.981538, 4.002603, 4.002603, 26.981538, beam_energy_at_vertex); + // ========================================================== std::string vtx_gate = ""; - if (vertex_recon >= -176.0 && vertex_recon < -100.0) - vtx_gate = "_Z[-176_to_-100]"; - else if (vertex_recon >= -100.0 && vertex_recon < -50.0) - vtx_gate = "_Z[-100_to_-50]"; - else if (vertex_recon >= -50.0 && vertex_recon < 0.0) - vtx_gate = "_Z[-50_to_0]"; - else if (vertex_recon >= 0.0 && vertex_recon < 50.0) - vtx_gate = "_Z[0_to_50]"; - else if (vertex_recon >= 50.0 && vertex_recon < 100.0) - vtx_gate = "_Z[50_to_100]"; - else if (vertex_recon >= 100.0 && vertex_recon < 176.0) - vtx_gate = "_Z[100_to_176]"; + if (vertex_recon >= -176.0 && vertex_recon < -100.0) + { + vtx_gate = "_Z[-176_to_-100]"; + } + else if (vertex_recon >= -100.0 && vertex_recon < -50.0) + { + vtx_gate = "_Z[-100_to_-50]"; + } + else if (vertex_recon >= -50.0 && vertex_recon < 0.0) + { + vtx_gate = "_Z[-50_to_0]"; + } + else if (vertex_recon >= 0.0 && vertex_recon < 50.0) + { + vtx_gate = "_Z[0_to_50]"; + } + else if (vertex_recon >= 50.0 && vertex_recon < 100.0) + { + vtx_gate = "_Z[50_to_100]"; + } + else if (vertex_recon >= 100.0 && vertex_recon < 176.0) + { + vtx_gate = "_Z[100_to_176]"; + } + + // 4. Energy Loss Correction in Silicon double path_length = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1; + double qqqEfix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1) - path_length); double qqqEfixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy2) - path_length); @@ -1394,28 +1446,95 @@ 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); + // 5. Fill Diagnostics + // plotter->Fill2D("1A0C_dE_Ecorr_Anode_QQQ", 400, 0, 30, 800, 0, 40000, qqqEfix, aEnergy * sinTheta, "1A0C"); + // plotter->Fill1D("1A0C_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, "1A0C"); + // plotter->Fill1D("1A0C_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, "1A0C"); + // plotter->Fill2D("1A0C_qqq_E_vs_theta_raw_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqevent.Energy1, "1A0C"); + // plotter->Fill2D("1A0C_qqq_E_vs_theta_corr_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqEfix, "1A0C"); 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("1A0C_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), "1A0C"); + plotter->Fill1D("1A0C_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, "1A0C"); + + plotter->Fill2D("1A0C_dE_Ecorr_Anode_QQQ" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfix, aEnergy * sinTheta, "1A0C"); + plotter->Fill2D("1A0C_dE_Ecorr_Anode_QQQ_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfixalpha, aEnergy * sinTheta, "1A0C"); + plotter->Fill1D("1A0C_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, "1A0C"); + plotter->Fill1D("1A0C_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, "1A0C"); + + // Track where on the wire the hit occurred (0 to 1 is inside the physical PC) + plotter->Fill1D("1A0C_wire_t_parameter_QQQ" + vtx_gate, 200, -0.5, 1.5, t_intersect_QQQ, "1A0C"); + } + } + + ///////////////////nA0C analysis using pseudo-wire (GetPseudoWire + getClosestWirePosAtWirePhi)/////////////////// + if (cClusters.size() == 0 && aClusters.size() > 0) + { + std::string nA0C_label = std::to_string(aClusters.size()) + "A0C"; + // Flatten all anode clusters into a combined hit list for the pseudo-wire + std::vector> allAnodeHits; + for (const auto &acluster : aClusters) + for (const auto &hit : acluster) + allAnodeHits.push_back(hit); + auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(allAnodeHits, "ANODE"); + for (auto qqqevent : QQQ_Events) + { + if (qqqevent.Time1 - apTSMaxE < -150) + { + // Use pseudo-wire to find the PC hit position at the QQQ phi + TVector3 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, qqqevent.pos.Phi()); + // Reconstruct vertex Z + double deltaRho = qqqevent.pos.Perp() - pcz_intersect.Perp(); + double deltaZ = qqqevent.pos.Z() - pcz_intersect.Z(); + double vertex_recon = qqqevent.pos.Z() - qqqevent.pos.Perp() * (deltaZ / deltaRho); + double z_entrance = -274.3; + double beam_path_length = TMath::Abs(vertex_recon - z_entrance) * 0.1; + double initial_beam_energy = 72.0; + double beam_energy_at_vertex = cm_to_MeV_beam->Eval( + MeV_to_cm_beam->Eval(initial_beam_energy) - beam_path_length); + Kinematics apkin_p(26.981538, 4.002603, 1.007825, 29.973770, beam_energy_at_vertex); + Kinematics apkin_a(26.981538, 4.002603, 4.002603, 26.981538, beam_energy_at_vertex); + std::string vtx_gate = ""; + if (vertex_recon >= -176.0 && vertex_recon < -100.0) + vtx_gate = "_Z[-176_to_-100]"; + else if (vertex_recon >= -100.0 && vertex_recon < -50.0) + vtx_gate = "_Z[-100_to_-50]"; + else if (vertex_recon >= -50.0 && vertex_recon < 0.0) + vtx_gate = "_Z[-50_to_0]"; + else if (vertex_recon >= 0.0 && vertex_recon < 50.0) + vtx_gate = "_Z[0_to_50]"; + else if (vertex_recon >= 50.0 && vertex_recon < 100.0) + vtx_gate = "_Z[50_to_100]"; + else if (vertex_recon >= 100.0 && vertex_recon < 176.0) + vtx_gate = "_Z[100_to_176]"; + double path_length = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1; + double qqqEfix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1) - path_length); + double qqqEfixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy2) - path_length); + double theta_recon = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Theta(); + double sinTheta = TMath::Sin(theta_recon); + 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); + 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); + } + } } } } -#endif // DIAG_NA0C_QQQ } -#endif // DIAG_NA0C_SX3 || DIAG_NA0C_QQQ -#ifdef DIAG_PC_QQQ for (auto pcevent : PC_Events) { for (auto qqqevent : QQQ_Events) @@ -1610,7 +1729,6 @@ Bool_t MakeVertex::Process(Long64_t entry) } } } // end PC QQQ coincidence -#endif // DIAG_PC_QQQ // HALFTIME! Can stop here in future versions // return kTRUE; if (anodeHits.size() >= 1 && cathodeHits.size() >= 1) @@ -2156,4 +2274,4 @@ void protonAlphaHistograms(HistPlotter *plotter, std::vector QQQ_Events, } // end QQQ_Events for loop, end sidetrack a(p,p) return; -} \ No newline at end of file +} diff --git a/README.md b/README.md index fceeb2a..c274f6d 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # ANASEN Analysis -Analysis code for the **Array for Nuclear Astrophysics and Structure with Exotic Nuclei (ANASEN)** detector at FSU. Processes raw FSUNSCL data through event building, channel mapping, calibration, and physics-level vertex reconstruction for transfer reaction experiments (e.g. ²⁷Al(α,p) and ¹⁷F(α,p)). +Analysis code for the **Array for Nuclear Astrophysics and Structure with Exotic Nuclei (ANASEN)** detector at FSU. Processes raw .fsu data through event building, channel mapping, calibration, and physics-level vertex reconstruction for transfer reaction experiments. --- @@ -19,7 +19,7 @@ The PC uses 24 twisted anode wires and 24 cathode wires. Wire geometry, crossove ## Full Analysis Chain ``` -Raw .fsu files (FSUNSCL digitizer output) +Raw .fsu files (FSU digitizer output) │ ▼ ┌─────────────────────────────────────────────────────────────────┐ @@ -46,7 +46,7 @@ Raw .fsu files (FSUNSCL digitizer output) │ Binary: Mapper (Armory/Mapper.cpp) │ │ Script: ProcessRun.sh 0 (calls Mapper internally) │ │ Config: mapping.h │ -│ Input : eventbuilt ROOT tree │ +│ Input : Eventbuilt ROOT tree │ │ Output: Run_NNN_mapped.root │ │ Translates hardware (digitizer SN, channel) to logical │ │ detector identity (SX3/QQQ/PC, strip/wire number). │ @@ -55,13 +55,22 @@ Raw .fsu files (FSUNSCL digitizer output) ▼ ┌─────────────────────────────────────────────────────────────────┐ │ 4. CALIBRATION │ -│ ├── sx3cal/EXFit.C / EXFit2.C │ -│ │ Fit SX3 front-strip position vs back-strip energy │ -│ │ to extract front/back gain coefficients │ -│ ├── sx3cal/LRFit.C │ -│ │ Left-right ratio fit for SX3 position calibration │ -│ │ Output: sx3cal/{17F,27Al}/ (frontgains.dat, │ -│ │ backgains.dat, rightgains.dat per run set) │ +│ │ +| SX3 — two-pass procedure: │ +│ Pass 1 — Left/Right matching (sx3cal/LRFit.C) │ +│ │ Start with unity gains: │ +│ │ LRFit.C fits the left/right charge ratio │ +│ │ Collate per-detector results into a single rightgains.dat │ +| │ │ +│ Pass 2 — Back/Front gain matching (sx3cal/EXFit.C) │ +│ │ Run on data that is unity-gain sorted but L/R matched │ +│ │ EXFit.C : │ +│ │ 1) gain-matches the back strips (backgains.dat) │ +│ │ 2) corrects dynamic range non-linearity in the fronts │ +│ │ (frontgains.dat) │ +│ │ Run for every detector, collate into master backgains.dat │ +│ │ and frontgains.dat. │ +│ │ | │ ├── GainMatchQQQ.C │ │ │ QQQ ring/wedge gain matching │ │ │ Output: qqq_GainMatch.dat │ diff --git a/anasen_fem/run.py b/anasen_fem/run.py index c3fbb09..6a54548 100755 --- a/anasen_fem/run.py +++ b/anasen_fem/run.py @@ -4,7 +4,7 @@ import os # val=-178.3 val=89.15 count=11 -while val<178.3+0.1: +while val<89.15+0.1: print(val) os.system("python3 wires_gmsh2d_bc.py "+str(val)) os.system("ElmerGrid 14 2 wires2d.msh -2d") diff --git a/anasen_fem/scalars.dat.names b/anasen_fem/scalars.dat.names index 68b1a00..f0aa591 100755 --- a/anasen_fem/scalars.dat.names +++ b/anasen_fem/scalars.dat.names @@ -1,8 +1,8 @@ Metadata for SaveScalars file: ./scalars.dat -Elmer version: 26.1 -Elmer compilation date: 2026-03-15 +Elmer version: 26.2 +Elmer compilation date: 2026-05-14 Solver input file: wires2d.sif -File started at: 2026/04/27 17:44:16 +File started at: 2026/05/15 17:54:54 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 9896518..ab69087 100755 --- a/anasen_fem/wires_gmsh2d_bc.py +++ b/anasen_fem/wires_gmsh2d_bc.py @@ -72,20 +72,26 @@ yarr_i12 = np.array([23 * np.sin(ki * i) for i in range(24)]) 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)]) -#guard wires, plane 2 at +zmax/2 -offsetg = offsetg-3*kg +# guard wires, plane 2 at +zmax/2 +# Old 3-wire shift: offsetg = offsetg - 3*kg +# For a 4-wire shift (relative to the 24-wire geometry, 4 anodes = 8 guard positions): +offsetg = offsetg - 8 * kg xarrg_2 = np.array([32*np.cos(kg*i+offsetg) for i in np.arange(0,48)]) yarrg_2 = np.array([32*np.sin(kg*i+offsetg) for i in np.arange(0,48)]) -#anodes, plane 2 at +zmax/2 -offset = offset-3*k +# anodes, plane 2 at +zmax/2 +# Old 3-wire shift: offset = offset - 3*k +# For a 4-wire shift: +offset = offset - 4 * k xarra_2 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)]) yarra_2 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)]) -#cathodes, plane2 at +zmax/2 -offsetc = offsetc-3*kc +# cathodes, plane 2 at +zmax/2 +# Old 3-wire shift: offsetc = offsetc - 3*kc +# For a 4-wire shift (matching guard wire rotation): +offsetc = offsetc - 4 * kc xarrc_2 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,48)]) -yarrc_2 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,48)]) +yarra_2 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,48)]) direction_needle_x = xarr_needle_2 - xarr_needle direction_needle_y = yarr_needle_2 - yarr_needle diff --git a/run_sx3.sh b/run_sx3.sh index b2fb43f..51409d7 100755 --- a/run_sx3.sh +++ b/run_sx3.sh @@ -20,7 +20,7 @@ fi export DATASET="27Al" export flip180="0" #root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run09.root; -if [[ 1 -eq 1 ]]; then +if [[ 1 -eq 0 ]]; then #export timecut_low=230.0; export timecut_low=400.0; #export timecut_high=400.0; @@ -37,7 +37,7 @@ fi #protons+gas, 27Al #export flip180="1" #export flip180="0" -if [[ 1 -eq 0 ]]; then +if [[ 1 -eq 1 ]]; then export flipa=0 export anode_offset=0 export source_vertex=-200.0; #put the 'source' on the entrance window