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

modified: MakeVertex.C

Browse Source 
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
This commit is contained in:
Vignesh Sitaraman 2026-05-15 18:01:31 -04:00
parent 17f44bc12e
commit b69dcf39d6
6 changed files with 240 additions and 107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

284
MakeVertex.C
View File

@ -1,13 +1,5 @@
#define MakeVertex_cxx #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] = { 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,
@ -859,7 +851,6 @@ Bool_t MakeVertex::Process(Long64_t entry)
if (QQQ_Events.size() && PC_Events.size()) if (QQQ_Events.size() && PC_Events.size())
plotter->Fill2D("PCEv_vs_QQQEv", 20, 0, 20, 20, 0, 20, 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"); plotter->Fill2D("ac_vs_cc", 20, 0, 20, 20, 0, 20, aClusters.size(), cClusters.size(), "wiremult");
for (auto cluster : aClusters) 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"); 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 (auto sx3event : SX3_Events)
{ {
for (int i = 0; i < 24; i++) for (int i = 0; i < 24; i++)
@ -951,9 +940,7 @@ 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)
{ {
@ -1238,7 +1225,6 @@ 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++) {
@ -1272,28 +1258,42 @@ Bool_t MakeVertex::Process(Long64_t entry)
} }
}*/ }*/
///////////////////nA0C analysis using pseudo-wire (GetPseudoWire + getClosestWirePosAtWirePhi)/////////////////// ///////////////////Single wire analysis for the anodes///////////////////
#if defined(DIAG_NA0C_SX3) || defined(DIAG_NA0C_QQQ) if (aClusters.size() == 1 && cClusters.size() == 0)
if (cClusters.size() == 0 && aClusters.size() >= 1)
{ {
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 // Get the 3D endpoints of the fired twisted anode wire from your geometry class
std::vector<std::tuple<int, double, double>> allAnodeHits; TVector3 a1 = pwinstance.An[aWireID].first;
for (const auto &acluster : aClusters) TVector3 wireVec = pwinstance.An[aWireID].first - pwinstance.An[aWireID].second;
for (const auto &hit : acluster)
allAnodeHits.push_back(hit);
auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(allAnodeHits, "ANODE"); // Loop over SX3_Events directly
#ifdef DIAG_NA0C_SX3
for (auto sx3event : SX3_Events) 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 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,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); 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)
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 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);
@ -1328,37 +1342,59 @@ 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);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_SX3", 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA0C_label); // 5. Fill Diagnostics
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, 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_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label); plotter->Fill1D("1A0C_Ex_from_protons_SX3", 200, -10, 10, Ex_from_proton, "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->Fill1D("1A0C_Ex_from_alphas_SX3", 200, -10, 10, Ex_from_alpha, "1A0C");
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(nA0C_label + "_twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA0C_label); plotter->Fill1D("1A0C_twisted_vertex_recon_SX3" + vtx_gate, 600, -300, 300, vertex_recon, "1A0C");
plotter->Fill1D(nA0C_label + "_twisted_vertex_recon_SX3" + vtx_gate, 600, -300, 300, vertex_recon, nA0C_label); plotter->Fill1D("1A0C_twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), "1A0C");
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->Fill2D("1A0C_dE_Ecorr_Anode_SX3" + vtx_gate, 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(nA0C_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label); 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) 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 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;
@ -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_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)
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 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);
@ -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_proton = apkin_p.getExc(qqqEfix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(qqqEfixalpha, theta_recon * 180. / M_PI); double Ex_from_alpha = apkin_a.getExc(qqqEfixalpha, theta_recon * 180. / M_PI);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ", 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA0C_label); // 5. Fill Diagnostics
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, 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_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label); // plotter->Fill1D("1A0C_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, "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->Fill1D("1A0C_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, "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_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 != "") if (vtx_gate != "")
{ {
plotter->Fill1D(nA0C_label + "_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA0C_label); plotter->Fill1D("1A0C_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), "1A0C");
plotter->Fill1D(nA0C_label + "_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, nA0C_label); plotter->Fill1D("1A0C_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, "1A0C");
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->Fill2D("1A0C_dE_Ecorr_Anode_QQQ" + vtx_gate, 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->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_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label); 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<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)
@ -1610,7 +1729,6 @@ 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)

29
README.md
View File

@ -1,6 +1,6 @@
# ANASEN Analysis # 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 ## 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) │ │ Binary: Mapper (Armory/Mapper.cpp) │
│ Script: ProcessRun.sh <run> <tw> 0 (calls Mapper internally) │ │ Script: ProcessRun.sh <run> <tw> 0 (calls Mapper internally) │
│ Config: mapping.h │ │ Config: mapping.h │
│ Input : eventbuilt ROOT tree │ │ Input : Eventbuilt ROOT tree │
│ Output: Run_NNN_mapped.root │ │ Output: Run_NNN_mapped.root │
│ Translates hardware (digitizer SN, channel) to logical │ │ Translates hardware (digitizer SN, channel) to logical │
│ detector identity (SX3/QQQ/PC, strip/wire number). │ │ detector identity (SX3/QQQ/PC, strip/wire number). │
@ -55,13 +55,22 @@ Raw .fsu files (FSUNSCL digitizer output)
▼ ▼
┌─────────────────────────────────────────────────────────────────┐ ┌─────────────────────────────────────────────────────────────────┐
│ 4. CALIBRATION │ │ 4. CALIBRATION │
│ ├── sx3cal/EXFit.C / EXFit2.C │ │ │
│ │ Fit SX3 front-strip position vs back-strip energy │ | SX3 — two-pass procedure: │
│ │ to extract front/back gain coefficients │ │ Pass 1 — Left/Right matching (sx3cal/LRFit.C) │
│ ├── sx3cal/LRFit.C │ │ │ Start with unity gains: │
│ │ Left-right ratio fit for SX3 position calibration │ │ │ LRFit.C fits the left/right charge ratio │
│ │ Output: sx3cal/{17F,27Al}/ (frontgains.dat, │ │ │ Collate per-detector results into a single rightgains.dat │
│ │ backgains.dat, rightgains.dat per run set) │ | │ │
│ 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 │ │ ├── GainMatchQQQ.C │
│ │ QQQ ring/wedge gain matching │ │ │ QQQ ring/wedge gain matching │
│ │ Output: qqq_GainMatch.dat │ │ │ Output: qqq_GainMatch.dat │

2
anasen_fem/run.py
View File

@ -4,7 +4,7 @@ import os
# val=-178.3 # val=-178.3
val=89.15 val=89.15
count=11 count=11
while val<178.3+0.1: while val<89.15+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")

6
anasen_fem/scalars.dat.names
View File

@ -1,8 +1,8 @@
Metadata for SaveScalars file: ./scalars.dat Metadata for SaveScalars file: ./scalars.dat
Elmer version: 26.1 Elmer version: 26.2
Elmer compilation date: 2026-03-15 Elmer compilation date: 2026-05-14
Solver input file: wires2d.sif 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: Variables in columns of matrix:
1: res: potential difference 1: res: potential difference

20
anasen_fem/wires_gmsh2d_bc.py
View File

@ -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)]) 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, plane 2 at +zmax/2
offsetg = offsetg-3*kg # 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)]) 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)]) yarrg_2 = np.array([32*np.sin(kg*i+offsetg) for i in np.arange(0,48)])
#anodes, plane 2 at +zmax/2 # anodes, plane 2 at +zmax/2
offset = offset-3*k # 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)]) 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)]) yarra_2 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
#cathodes, plane2 at +zmax/2 # cathodes, plane 2 at +zmax/2
offsetc = offsetc-3*kc # 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)]) 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_x = xarr_needle_2 - xarr_needle
direction_needle_y = yarr_needle_2 - yarr_needle direction_needle_y = yarr_needle_2 - yarr_needle

4
run_sx3.sh
View File

@ -20,7 +20,7 @@ fi
export DATASET="27Al" export DATASET="27Al"
export flip180="0" 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; #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=230.0;
export timecut_low=400.0; export timecut_low=400.0;
#export timecut_high=400.0; #export timecut_high=400.0;
@ -37,7 +37,7 @@ fi
#protons+gas, 27Al #protons+gas, 27Al
#export flip180="1" #export flip180="1"
#export flip180="0" #export flip180="0"
if [[ 1 -eq 0 ]]; then if [[ 1 -eq 1 ]]; then
export flipa=0 export flipa=0
export anode_offset=0 export anode_offset=0
export source_vertex=-200.0; #put the 'source' on the entrance window export source_vertex=-200.0; #put the 'source' on the entrance window