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

modified: MakeVertex.C

Browse Source 
modified:   anasen_fem/run.py
	modified:   anasen_fem/scalars.dat.names
	modified:   anasen_fem/wires_gmsh2d_bc.py
This commit is contained in:
Vignesh Sitaraman 2026-05-19 10:28:41 -04:00
parent ce0bedf704
commit d8593b2f55
4 changed files with 92 additions and 215 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

235
MakeVertex.C
View File

@ -1,5 +1,12 @@
#define MakeVertex_cxx #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] = { Int_t colors[40] = {
kBlack, kRed, kGreen, kBlue, kYellow, kMagenta, kCyan, kOrange, kBlack, kRed, kGreen, kBlue, kYellow, kMagenta, kCyan, kOrange,
kSpring, kTeal, kAzure, kViolet, kPink, kGray, kWhite, 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"); 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 (auto sx3event : SX3_Events)
{ {
for (int i = 0; i < 24; i++) for (int i = 0; i < 24; i++)
@ -940,7 +948,9 @@ Bool_t MakeVertex::Process(Long64_t entry)
} }
} // for 'i' loop } // for 'i' loop
} }
#endif // DIAG_1WIRE
#ifdef DIAG_PC_SX3
bool PCSX3PhiCut = false; bool PCSX3PhiCut = false;
for (auto pcevent : PC_Events) for (auto pcevent : PC_Events)
{ {
@ -1225,6 +1235,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
}*/ }*/
} }
} // end PC-SX3 coincidence } // end PC-SX3 coincidence
#endif // DIAG_PC_SX3
/*for(size_t ii=0; ii<QQQ_Events.size(); ii++) { /*for(size_t ii=0; ii<QQQ_Events.size(); ii++) {
for(size_t jj=ii+1; jj<QQQ_Events.size(); jj++) { for(size_t jj=ii+1; jj<QQQ_Events.size(); jj++) {
@ -1258,42 +1269,28 @@ Bool_t MakeVertex::Process(Long64_t entry)
} }
}*/ }*/
///////////////////Single wire analysis for the anodes/////////////////// ///////////////////nA0C analysis using pseudo-wire (GetPseudoWire + getClosestWirePosAtWirePhi)///////////////////
if (aClusters.size() == 1 && cClusters.size() == 0) #if defined(DIAG_NA0C_SX3) || defined(DIAG_NA0C_QQQ)
if (cClusters.size() == 0 && aClusters.size() >= 1)
{ {
// Extract the primary anode hit properties std::string nA0C_label = std::to_string(aClusters.size()) + "A0C";
auto anodeHit = aClusters.front().front();
int aWireID = std::get<0>(anodeHit);
double aEnergy = std::get<1>(anodeHit);
double aTime = std::get<2>(anodeHit);
// Get the 3D endpoints of the fired twisted anode wire from your geometry class // Flatten all anode clusters into a combined hit list for the pseudo-wire
TVector3 a1 = pwinstance.An[aWireID].first; std::vector<std::tuple<int, double, double>> allAnodeHits;
TVector3 wireVec = pwinstance.An[aWireID].first - pwinstance.An[aWireID].second; 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) for (auto sx3event : SX3_Events)
{ {
if (sx3event.Time1 - apTSMaxE < -150)
if (sx3event.Time1 - aTime < -150) // Time cut for protons
{ {
// 1. Define the plane of the track (Z-axis to SX3 hit) TVector3 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, sx3event.pos.Phi());
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 deltaRho = sx3event.pos.Perp() - pcz_intersect.Perp();
double deltaZ = sx3event.pos.Z() - pcz_intersect.Z(); double deltaZ = sx3event.pos.Z() - pcz_intersect.Z();
double vertex_recon = sx3event.pos.Z() - sx3event.pos.Perp() * (deltaZ / deltaRho); 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); Kinematics apkin_a(26.981538, 4.002603, 4.002603, 26.981538, beam_energy_at_vertex);
std::string vtx_gate = ""; std::string vtx_gate = "";
if (vertex_recon >= -176.0 && vertex_recon < -100.0) if (vertex_recon >= -176.0 && vertex_recon < -100.0)
{
vtx_gate = "_Z[-176_to_-100]"; vtx_gate = "_Z[-176_to_-100]";
}
else if (vertex_recon >= -100.0 && vertex_recon < -50.0) else if (vertex_recon >= -100.0 && vertex_recon < -50.0)
{
vtx_gate = "_Z[-100_to_-50]"; vtx_gate = "_Z[-100_to_-50]";
}
else if (vertex_recon >= -50.0 && vertex_recon < 0.0) else if (vertex_recon >= -50.0 && vertex_recon < 0.0)
{
vtx_gate = "_Z[-50_to_0]"; vtx_gate = "_Z[-50_to_0]";
}
else if (vertex_recon >= 0.0 && vertex_recon < 50.0) else if (vertex_recon >= 0.0 && vertex_recon < 50.0)
{
vtx_gate = "_Z[0_to_50]"; vtx_gate = "_Z[0_to_50]";
}
else if (vertex_recon >= 50.0 && vertex_recon < 100.0) else if (vertex_recon >= 50.0 && vertex_recon < 100.0)
{
vtx_gate = "_Z[50_to_100]"; vtx_gate = "_Z[50_to_100]";
}
else if (vertex_recon >= 100.0 && vertex_recon < 176.0) else if (vertex_recon >= 100.0 && vertex_recon < 176.0)
{
vtx_gate = "_Z[100_to_176]"; 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 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 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); 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 theta_recon = (sx3event.pos - TVector3(0, 0, vertex_recon)).Theta();
double sinTheta = TMath::Sin(theta_recon); 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_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);
// 5. Fill Diagnostics plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_SX3", 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA0C_label);
plotter->Fill2D("1A0C_dE_Ecorr_Anode_SX3", 400, 0, 30, 800, 0, 40000, sx3Efix, aEnergy * sinTheta, "1A0C"); plotter->Fill1D(nA0C_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label);
plotter->Fill1D("1A0C_Ex_from_protons_SX3", 200, -10, 10, Ex_from_proton, "1A0C"); plotter->Fill1D(nA0C_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label);
plotter->Fill1D("1A0C_Ex_from_alphas_SX3", 200, -10, 10, Ex_from_alpha, "1A0C"); 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("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 != "") if (vtx_gate != "")
{ {
plotter->Fill1D("1A0C_twisted_vertex_recon_SX3" + vtx_gate, 600, -300, 300, vertex_recon, "1A0C"); plotter->Fill1D(nA0C_label + "_twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA0C_label);
plotter->Fill1D("1A0C_twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), "1A0C"); 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("1A0C_dE_Ecorr_Anode_SX3" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efix, aEnergy * sinTheta, "1A0C"); 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("1A0C_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, "1A0C"); plotter->Fill1D(nA0C_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label);
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
// Loop over QQQ_Events directly #ifdef DIAG_NA0C_QQQ
for (auto qqqevent : QQQ_Events) 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 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, qqqevent.pos.Phi());
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 deltaRho = qqqevent.pos.Perp() - pcz_intersect.Perp();
double deltaZ = qqqevent.pos.Z() - pcz_intersect.Z(); double deltaZ = qqqevent.pos.Z() - pcz_intersect.Z();
double vertex_recon = qqqevent.pos.Z() - qqqevent.pos.Perp() * (deltaZ / deltaRho); double vertex_recon = qqqevent.pos.Z() - qqqevent.pos.Perp() * (deltaZ / deltaRho);
double z_entrance = -274.3; 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_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); Kinematics apkin_a(26.981538, 4.002603, 4.002603, 26.981538, beam_energy_at_vertex);
// ==========================================================
std::string vtx_gate = ""; std::string vtx_gate = "";
if (vertex_recon >= -176.0 && vertex_recon < -100.0) if (vertex_recon >= -176.0 && vertex_recon < -100.0)
{
vtx_gate = "_Z[-176_to_-100]"; vtx_gate = "_Z[-176_to_-100]";
}
else if (vertex_recon >= -100.0 && vertex_recon < -50.0) else if (vertex_recon >= -100.0 && vertex_recon < -50.0)
{
vtx_gate = "_Z[-100_to_-50]"; vtx_gate = "_Z[-100_to_-50]";
}
else if (vertex_recon >= -50.0 && vertex_recon < 0.0) else if (vertex_recon >= -50.0 && vertex_recon < 0.0)
{
vtx_gate = "_Z[-50_to_0]"; vtx_gate = "_Z[-50_to_0]";
}
else if (vertex_recon >= 0.0 && vertex_recon < 50.0) else if (vertex_recon >= 0.0 && vertex_recon < 50.0)
{
vtx_gate = "_Z[0_to_50]"; vtx_gate = "_Z[0_to_50]";
}
else if (vertex_recon >= 50.0 && vertex_recon < 100.0) else if (vertex_recon >= 50.0 && vertex_recon < 100.0)
{
vtx_gate = "_Z[50_to_100]"; vtx_gate = "_Z[50_to_100]";
}
else if (vertex_recon >= 100.0 && vertex_recon < 176.0) else if (vertex_recon >= 100.0 && vertex_recon < 176.0)
{
vtx_gate = "_Z[100_to_176]"; 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 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 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 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_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);
// 5. Fill Diagnostics plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ", 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA0C_label);
// plotter->Fill2D("1A0C_dE_Ecorr_Anode_QQQ", 400, 0, 30, 800, 0, 40000, qqqEfix, aEnergy * sinTheta, "1A0C"); plotter->Fill1D(nA0C_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label);
// plotter->Fill1D("1A0C_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, "1A0C"); plotter->Fill1D(nA0C_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label);
// plotter->Fill1D("1A0C_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, "1A0C"); 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("1A0C_qqq_E_vs_theta_raw_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqevent.Energy1, "1A0C"); 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("1A0C_qqq_E_vs_theta_corr_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqEfix, "1A0C");
if (vtx_gate != "") if (vtx_gate != "")
{ {
plotter->Fill1D("1A0C_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), "1A0C"); plotter->Fill1D(nA0C_label + "_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA0C_label);
plotter->Fill1D("1A0C_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, "1A0C"); 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("1A0C_dE_Ecorr_Anode_QQQ" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfix, aEnergy * sinTheta, "1A0C"); plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfixalpha, apSumE * sinTheta, nA0C_label);
plotter->Fill2D("1A0C_dE_Ecorr_Anode_QQQ_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfixalpha, aEnergy * sinTheta, "1A0C"); plotter->Fill1D(nA0C_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label);
plotter->Fill1D("1A0C_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, "1A0C"); plotter->Fill1D(nA0C_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label);
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<std::tuple<int, double, double>> 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 pcevent : PC_Events)
{ {
for (auto qqqevent : QQQ_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->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->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->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()); 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); 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_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("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->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_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); // 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 } // end PC QQQ coincidence
#endif // DIAG_PC_QQQ
// HALFTIME! Can stop here in future versions // HALFTIME! Can stop here in future versions
// return kTRUE; // return kTRUE;
if (anodeHits.size() >= 1 && cathodeHits.size() >= 1) if (anodeHits.size() >= 1 && cathodeHits.size() >= 1)
@ -2274,4 +2155,4 @@ void protonAlphaHistograms(HistPlotter *plotter, std::vector<Event> QQQ_Events,
} // end QQQ_Events for loop, end sidetrack a(p,p) } // end QQQ_Events for loop, end sidetrack a(p,p)
return; return;
} }

4
anasen_fem/run.py
View File

@ -4,13 +4,13 @@ import os
# val=-178.3 # val=-178.3
val=89.15 val=89.15
count=11 count=11
while val<89.15+0.1: while val<178.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")
os.system("ElmerSolver wires2d.sif") os.system("ElmerSolver wires2d.sif")
os.system("./paraview_plotter.py") 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.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.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)) os.system("cp wires2d/elfield_anasen_t0001.vtu wires2d/vtu_files/elfield_anasen_%02d_%1.4f.vtu"%(count,val))

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/15 17:54:54 File started at: 2026/05/18 09:41:59
Variables in columns of matrix: Variables in columns of matrix:
1: res: potential difference 1: res: potential difference

66
anasen_fem/wires_gmsh2d_bc.py
View File

@ -25,7 +25,9 @@ if len(sys.argv) < 2:
z_loc = float(sys.argv[1]) 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 #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)]) 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)]) yarr_i21 = np.array([23 * np.sin(ki * i + ki/2.0) for i in range(24)])
#guard wires, plane 1 at -zmax/2 # Guard wires (48 total) - aligned with Cathode phasing
kg=2*np.pi/48. offset_g1 = -5 * k - 2 * k - (np.pi / 24.0)
offsetg = -4*kg + 2*kg - np.pi/24 #-pi/4 xarrg_1 = np.array([32 * np.cos(kg * i + offset_g1) for i in range(48)])
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 + offset_g1) for i in range(48)])
yarrg_1 = np.array([32*np.sin(kg*i+offsetg) for i in np.arange(0,48)])
#anodes, plane 1 at -zmax/2 # Anodes (24 total) - index increases leftward (-k)
k=-2*np.pi/24. offset_a1 = -5 * k - 5 * k
offset = 6*k + 3*k #-pi/2 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) for i in np.arange(0,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) for i in np.arange(0,24)])
#cathodes, plane 1 at -zmax/2 # Cathodes (24 total) - index increases rightward (+k)
kc=2*np.pi/48. offset_c1 = -5 * k - 2 * k - (np.pi / 24.0)
offsetc = -4*kc + 2*kc - np.pi/24 #-pi/4 xarrc_1 = np.array([42 * np.cos(k * i + offset_c1) for i in range(24)])
xarrc_1 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,48)]) yarrc_1 = np.array([42 * np.sin(k * i + offset_c1) for i in range(24)])
yarrc_1 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,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])
@ -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)]) 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)]) yarr_i22 = np.array([23 * np.sin(ki * i + ki/2.0) for i in range(24)])
# guard wires, plane 2 at +zmax/2 # Guard wires (48 total) - twists leftward to match cathodes
# Old 3-wire shift: offsetg = offsetg - 3*kg offset_g2 = offset_g1 - (wireShift * k)
# For a 4-wire shift (relative to the 24-wire geometry, 4 anodes = 8 guard positions): xarrg_2 = np.array([32 * np.cos(kg * i + offset_g2) for i in range(48)])
offsetg = offsetg - 8 * kg yarrg_2 = np.array([32 * np.sin(kg * i + offset_g2) for i in range(48)])
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 # Anodes (24 total) - twists rightward (+shift)
# Old 3-wire shift: offset = offset - 3*k offset_a2 = offset_a1 + (wireShift * k)
# For a 4-wire shift: xarra_2 = np.array([37 * np.cos(-k * i + offset_a2) for i in range(24)])
offset = offset - 4 * k yarra_2 = np.array([37 * np.sin(-k * i + offset_a2) for i in range(24)])
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, plane 2 at +zmax/2 # Cathodes (24 total) - twists leftward (-shift)
# Old 3-wire shift: offsetc = offsetc - 3*kc offset_c2 = offset_c1 - (wireShift * k)
# For a 4-wire shift (matching guard wire rotation): xarrc_2 = np.array([42 * np.cos(k * i + offset_c2) for i in range(24)])
offsetc = offsetc - 4 * kc yarrc_2 = np.array([42 * np.sin(k * i + offset_c2) for i in range(24)])
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)])
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
@ -150,9 +143,12 @@ for i, (xn, yn) in enumerate(zip(xloc_needle, yloc_needle)):
needle.append(ndisk) needle.append(ndisk)
#create Guard Wires (48 total) #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) gdisk = gmsh.model.occ.addDisk(xg, yg, 0, wire_radius, wire_radius)
guard_wires.append(gdisk) 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) cdisk = gmsh.model.occ.addDisk(xc, yc, 0, wire_radius, wire_radius)
cathode_wires.append(cdisk) cathode_wires.append(cdisk)