diff --git a/MakeVertex.C b/MakeVertex.C index 0db0c50..6b1ebd6 100755 --- a/MakeVertex.C +++ b/MakeVertex.C @@ -1,5 +1,12 @@ #define MakeVertex_cxx +// 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_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, @@ -866,6 +873,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"); } +#ifdef DIAG_1WIRE for (auto sx3event : SX3_Events) { for (int i = 0; i < 24; i++) @@ -940,7 +948,9 @@ 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) { @@ -1225,6 +1235,7 @@ Bool_t MakeVertex::Process(Long64_t entry) }*/ } } // end PC-SX3 coincidence +#endif // DIAG_PC_SX3 /*for(size_t ii=0; ii= 1) { - // 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); + std::string nA0C_label = std::to_string(aClusters.size()) + "A0C"; - // 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; + // 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); - // Loop over SX3_Events directly + auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(allAnodeHits, "ANODE"); + +#ifdef DIAG_NA0C_SX3 for (auto sx3event : SX3_Events) { - - if (sx3event.Time1 - aTime < -150) // Time cut for protons + if (sx3event.Time1 - apTSMaxE < -150) { - // 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); + TVector3 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, sx3event.pos.Phi()); - // 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,33 +1305,19 @@ 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]"; - } - // 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); @@ -1342,59 +1325,37 @@ 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); - // 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"); + 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); if (vtx_gate != "") { - 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"); + 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); } } } +#endif // DIAG_NA0C_SX3 - // Loop over QQQ_Events directly - +#ifdef DIAG_NA0C_QQQ for (auto qqqevent : QQQ_Events) { - if (qqqevent.Time1 - aTime < -150) // Time cut for protons + if (qqqevent.Time1 - apTSMaxE < -150) { - // 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); + TVector3 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, qqqevent.pos.Phi()); - // 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; @@ -1405,38 +1366,22 @@ 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]"; - } - // 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); @@ -1446,95 +1391,28 @@ 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); - // 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"); + 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("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); - } - } + 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) @@ -1543,6 +1421,7 @@ Bool_t MakeVertex::Process(Long64_t entry) plotter->Fill2D("dt_pcA_qqqR_vs_qqqRE", 640, -2000, 2000, 400, 0, 30, qqqevent.Time1 - pcevent.Time1, qqqevent.Energy1); plotter->Fill1D("dt_pcC_qqqW", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2); plotter->Fill2D("phiPC_vs_phiQQQ", 180, -360, 360, 180, -360, 360, qqqevent.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI); + plotter->Fill1D("phiPC-phiQQQ", 180, -180, 180, pcevent.pos.Phi() * 180 / M_PI - qqqevent.pos.Phi() * 180 / M_PI); double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()); /// TMath::Sin((TVector3(51.5,0,128.) - TVector3(0,0,85)).Theta()); TVector3 x2(pcevent.pos); @@ -1595,6 +1474,7 @@ Bool_t MakeVertex::Process(Long64_t entry) plotter->Fill2D("dE2_theta_AnodeQQQR_zoomin", 60, 0, 30, 400, 0, 5000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy1 * sinTheta); plotter->Fill2D("dE2_theta_AnodeQQQR", 90, 0, 90, 400, 0, 20000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy1 * sinTheta); plotter->Fill2D("phiPC_vs_phiQQQ_TimeCut", 180, -360, 360, 180, -360, 360, qqqevent.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI); + plotter->Fill1D("phiPC-phiQQQ_TimeCut", 180, 0, 180, pcevent.pos.Phi() * 180 / M_PI - qqqevent.pos.Phi() * 180 / M_PI); // plotter->Fill2D("E_theta_AnodeQQQR_TC1_PC"+std::to_string(phicut),75,0,90,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1); // plotter->Fill2D("E_theta_zoomin_AnodeQQQR_TC1_PC"+std::to_string(phicut),60,0,30,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1); @@ -1729,6 +1609,7 @@ 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) @@ -2274,4 +2155,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/anasen_fem/run.py b/anasen_fem/run.py index 6a54548..f26f85f 100755 --- a/anasen_fem/run.py +++ b/anasen_fem/run.py @@ -4,13 +4,13 @@ import os # val=-178.3 val=89.15 count=11 -while val<89.15+0.1: +while val<178.3+0.1: print(val) os.system("python3 wires_gmsh2d_bc.py "+str(val)) os.system("ElmerGrid 14 2 wires2d.msh -2d") os.system("ElmerSolver wires2d.sif") os.system("./paraview_plotter.py") - os.system("python3 garfield_sim.py") + # os.system("python3 garfield_sim.py") os.system("cp wires2d.msh wires2d/mesh_files/wires2d%02d_%1.4f.msh"%(count,val)) os.system("cp wires2d.sif wires2d/sif_files/wires2d_%02d_%1.4f.sif"%(count,val)) os.system("cp wires2d/elfield_anasen_t0001.vtu wires2d/vtu_files/elfield_anasen_%02d_%1.4f.vtu"%(count,val)) diff --git a/anasen_fem/scalars.dat.names b/anasen_fem/scalars.dat.names index f0aa591..44cd88c 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/15 17:54:54 +File started at: 2026/05/18 09:41:59 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 ab69087..0c39847 100755 --- a/anasen_fem/wires_gmsh2d_bc.py +++ b/anasen_fem/wires_gmsh2d_bc.py @@ -25,7 +25,9 @@ if len(sys.argv) < 2: z_loc = float(sys.argv[1]) -k=(2*np.pi/24.) +wireShift = 4.0 +k = (2 * np.pi / 24.0) +kg = (2 * np.pi / 48.0) #1 needle, 24 ic1, 24 ic2, 48 guard wires, 24 anodes, 24 cathodes @@ -42,23 +44,20 @@ yarr_i11 = np.array([23 * np.sin(ki * i) for i in range(24)]) 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)]) -#guard wires, plane 1 at -zmax/2 -kg=2*np.pi/48. -offsetg = -4*kg + 2*kg - np.pi/24 #-pi/4 -xarrg_1 = np.array([32*np.cos(kg*i+offsetg) for i in np.arange(0,48)]) -yarrg_1 = np.array([32*np.sin(kg*i+offsetg) for i in np.arange(0,48)]) +# 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, plane 1 at -zmax/2 -k=-2*np.pi/24. -offset = 6*k + 3*k #-pi/2 -xarra_1 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)]) -yarra_1 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)]) +# 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)]) +yarra_1 = np.array([37 * np.sin(-k * i + offset_a1) for i in range(24)]) -#cathodes, plane 1 at -zmax/2 -kc=2*np.pi/48. -offsetc = -4*kc + 2*kc - np.pi/24 #-pi/4 -xarrc_1 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,48)]) -yarrc_1 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,48)]) +# Cathodes (24 total) - index increases rightward (+k) +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)]) +yarrc_1 = np.array([42 * np.sin(k * i + offset_c1) for i in range(24)]) #needle at plane 2 at zmax/2 xarr_needle_2 = np.array([0]) @@ -72,26 +71,20 @@ 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 -# 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)]) +# 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, 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)]) +# 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)]) +yarra_2 = np.array([37 * np.sin(-k * i + offset_a2) for i in range(24)]) -# 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)]) -yarra_2 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,48)]) +# 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)]) direction_needle_x = xarr_needle_2 - xarr_needle direction_needle_y = yarr_needle_2 - yarr_needle @@ -150,9 +143,12 @@ for i, (xn, yn) in enumerate(zip(xloc_needle, yloc_needle)): needle.append(ndisk) #create Guard Wires (48 total) -for i, (xg, yg, xc, yc) in enumerate(zip(xloc_g, yloc_g, xloc_c, yloc_c)): +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) +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)