From 34b9fd5f3d3202dc43c3bc4e6912b0aa92810d1a Mon Sep 17 00:00:00 2001 From: James Szalkie Date: Mon, 21 Sep 2026 11:33:59 -0400 Subject: [PATCH] new vZ recon from known Ex --- .../anasen_anode_cathode_hyperboloids.h | 108 ++++++++++++++ Armory/anasenMS.cpp | 14 +- Armory/anasen_anode_cathode_hyperboloids.h | 133 ++++-------------- Armory/testMacro.C | 7 + 4 files changed, 157 insertions(+), 105 deletions(-) create mode 100644 Armory/Archives/anasen_anode_cathode_hyperboloids.h diff --git a/Armory/Archives/anasen_anode_cathode_hyperboloids.h b/Armory/Archives/anasen_anode_cathode_hyperboloids.h new file mode 100644 index 0000000..96c1c60 --- /dev/null +++ b/Armory/Archives/anasen_anode_cathode_hyperboloids.h @@ -0,0 +1,108 @@ +inline std::tuple find_PC_PathLength(const TVector3& x1, const TVector3& x2) { + /* + Square-aperture detector geometry replacing the original circular one-sheet hyperboloid. + The active surfaces are centered on the z axis and form squares rather than circles: + - anodes: 40 mm x 40 mm + - cathodes: 42 mm x 42 mm + The beam is kept on-axis, so the square openings enclose the beam while still defining + the anode/cathode boundaries in the x-y plane. + */ + + /* + // Original circular hyperboloid implementation retained for reference. + // This section is intentionally commented out to avoid using the circular geometry. + TVector3 dx = x2-x1; // direction vector + double t2 = 1.0; //The value of 't' at the destination point, by definition: t=(z(t)-z0)/dz + auto onesheet_hyperboloid_intersect = [&](double a, double c) { + auto A = pow(dx.Perp(),2)/(a*a) - pow(dx.Z(),2)/(c*c); + auto B = 2*(dx.X()*x1.X()+dx.Y()*x1.Y())/(a*a) - 2*(dx.Z()*x1.Z())/(c*c); + auto C = pow(x1.Perp(),2)/(a*a) - pow(x1.Z(),2)/(c*c) - 1.0; + double disc = B*B - 4*A*C; + if(disc<0) + return TVector3(0,0,54321); + else { + double tsol1 = (-B + TMath::Sqrt(disc))/(2*A); + double tsol2 = (-B - TMath::Sqrt(disc))/(2*A); + if(tsol1 >= 0 && tsol1 <= t2) + return x1+tsol1*dx; + else if(tsol2>=0 && tsol2 <= t2) + return x1+tsol2*dx; + else + return TVector3(0,0,54321); + } + }; + + TVector3 anode_intersect = onesheet_hyperboloid_intersect(32.0429,301.895); + TVector3 cathode_intersect = onesheet_hyperboloid_intersect(37.239045,301.895); + TVector3 gw_intersect = onesheet_hyperboloid_intersect(27.712,301.895); + if(anode_intersect.Z()!=54321 && cathode_intersect.Z()!=54321 && gw_intersect.Z()!=54321) + return std::tuple(cathode_intersect,gw_intersect,(cathode_intersect-gw_intersect).Mag()*0.1); + else + return std::tuple(TVector3(0,0,0), TVector3(0,0,0), 54321); + */ + + const double anode_half_width = 20.0; // 40 mm square side length + const double cathode_half_width = 21.0; // 42 mm square side length + const double guard_half_width = 18.0; // square guard opening centered on beam axis, still encloses beam + const double z_start = -100.0; // square detector opening begins here along the beam axis + const double invalid_z = 54321.0; + const double eps = 1.0e-12; + + TVector3 dx = x2 - x1; + + auto square_intersection = [&](const TVector3& start, const TVector3& direction, double half_width) { + double best_t = 1.0; + bool found = false; + + auto consider_t = [&](double t, bool x_hit, bool y_hit) { + if (t < -eps || t > 1.0 + eps) return; + + TVector3 point = start + t * direction; + double abs_x = std::abs(point.X()); + double abs_y = std::abs(point.Y()); + + bool inside_square = (abs_x <= half_width + eps) && (abs_y <= half_width + eps); + bool on_boundary = (x_hit && std::abs(abs_x - half_width) <= eps) || (y_hit && std::abs(abs_y - half_width) <= eps); + bool beyond_start = point.Z() >= z_start - eps; + + if (inside_square && on_boundary && beyond_start && t >= 0.0 && t <= 1.0) { + if (!found || t < best_t) { + best_t = t; + found = true; + } + } + }; + + if (std::abs(direction.X()) > eps) { + double t_plus_x = (half_width - start.X()) / direction.X(); + double t_minus_x = (-half_width - start.X()) / direction.X(); + consider_t(t_plus_x, true, false); + consider_t(t_minus_x, true, false); + } + + if (std::abs(direction.Y()) > eps) { + double t_plus_y = (half_width - start.Y()) / direction.Y(); + double t_minus_y = (-half_width - start.Y()) / direction.Y(); + consider_t(t_plus_y, false, true); + consider_t(t_minus_y, false, true); + } + + if (!found) { + return TVector3(0, 0, invalid_z); + } + + return start + best_t * direction; + }; + + TVector3 anode_intersect = square_intersection(x1, dx, anode_half_width); + TVector3 cathode_intersect = square_intersection(x1, dx, cathode_half_width); + TVector3 gw_intersect = square_intersection(x1, dx, guard_half_width); + + if (anode_intersect.Z() != invalid_z && cathode_intersect.Z() != invalid_z && gw_intersect.Z() != invalid_z) + return std::tuple(cathode_intersect, gw_intersect, (cathode_intersect - gw_intersect).Mag() * 0.1); + else + return std::tuple(TVector3(0, 0, 0), TVector3(0, 0, 0), invalid_z); +} + + + diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index e13526a..c305af0 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -162,7 +162,7 @@ int main(int argc, char **argv){ transfer.Seta(4, 2); // 4He target transfer.Setb(1, 1); // outgoing proton from the primary transfer transfer.SetB(30, 14); // 30Si* heavy product - double beamE = 100; //56.1 + double beamE = 56.1; //56.1 const ReactionConfig reactionConfig = transfer.GetRectionConfig(); const double beamA = reactionConfig.beamA; // mass number of 14N beam @@ -415,11 +415,12 @@ int main(int argc, char **argv){ tree1->Branch("qqqUp", &qqqUp, "qqqUp/I"); tree1->Branch("qqqBk", &qqqBk, "qqqBk/I"); - double EBeam_Kin_gs=NAN, EBeam_Kin_2_2=NAN, EBeam_Kin_3_4=NAN, Ex_recon=NAN; + double EBeam_Kin_gs=NAN, EBeam_Kin_2_2=NAN, EBeam_Kin_3_4=NAN, Ex_recon=NAN, EbeamRecon=NAN; tree1->Branch("EBeam_Kin", &EBeam_Kin_gs, "EBeam_Kin/D"); tree1->Branch("EBeam_Kin_2.2", &EBeam_Kin_2_2, "EBeam_Kin_2.2/D"); tree1->Branch("EBeam_Kin_3.4", &EBeam_Kin_3_4, "EBeam_Kin_3.4/D"); tree1->Branch("Ex_recon", &Ex_recon, "Ex_recon/D"); + tree1->Branch("EbeamRecon", &EbeamRecon, "EbeamRecon/D"); // reconstructed angles from PW track fit, method 1 and 2 double reTheta, rePhi; @@ -436,10 +437,11 @@ int main(int argc, char **argv){ tree1->Branch("sigma_a", &sigma_a, "sigma_a/D"); tree1->Branch("theta_recon", &theta_recon, "theta_recon/D"); - double vX_recon, vY_recon, vZ_recon; + double vX_recon, vY_recon, vZ_recon, vZ_from_beam; tree1->Branch("vX_recon", &vX_recon, "vX_recon/D"); tree1->Branch("vY_recon", &vY_recon, "vY_recon/D"); tree1->Branch("vZ_recon", &vZ_recon, "vZ_recon/D"); + tree1->Branch("vZ_from_beam", &vZ_from_beam, "vZ_from_beam/D"); // reconstructed vertex Z from PW fit double z0; @@ -499,6 +501,7 @@ int main(int argc, char **argv){ vY_recon = TMath::QuietNaN(); vZ_recon = TMath::QuietNaN(); Ex_recon = TMath::QuietNaN(); + EbeamRecon = TMath::QuietNaN(); sx3ID = -1; sx3Up = -1; sx3Dn = -1; @@ -735,15 +738,20 @@ int main(int argc, char **argv){ distance_sx3); Ex_recon = apkin_27Al.getExc(originalEnergy, theta_recon); + EbeamRecon = apkin_27Al.getEbeam_givenQ(originalEnergy, Ex, theta_recon); + vZ_from_beam = elossBeamInverse->Eval(EbeamRecon) - 450; /*Checklist: anode smudge, sx3 smudge, beam position off axis, beam angle*/ + } else { theta_recon = TMath::QuietNaN(); vX_recon = TMath::QuietNaN(); vY_recon = TMath::QuietNaN(); vZ_recon = TMath::QuietNaN(); Ex_recon = TMath::QuietNaN(); + EbeamRecon = TMath::QuietNaN(); originalEnergy = TMath::QuietNaN(); + vZ_from_beam = TMath::QuietNaN(); } //EBeam_Kin_gs = apkin_27Al.getEbeam_givenQ(Esx3, 0.0, thetab); //EBeam_Kin_2_2 = apkin_27Al.getEbeam_givenQ(Esx3, 2.2, thetab); diff --git a/Armory/anasen_anode_cathode_hyperboloids.h b/Armory/anasen_anode_cathode_hyperboloids.h index 96c1c60..08843d6 100644 --- a/Armory/anasen_anode_cathode_hyperboloids.h +++ b/Armory/anasen_anode_cathode_hyperboloids.h @@ -1,108 +1,37 @@ inline std::tuple find_PC_PathLength(const TVector3& x1, const TVector3& x2) { /* - Square-aperture detector geometry replacing the original circular one-sheet hyperboloid. - The active surfaces are centered on the z axis and form squares rather than circles: - - anodes: 40 mm x 40 mm - - cathodes: 42 mm x 42 mm - The beam is kept on-axis, so the square openings enclose the beam while still defining - the anode/cathode boundaries in the x-y plane. + Function that finds the path length between anode and cathode surfaces, both one-sheet hyperboloids of form (x*x+y*y)/(a*a) - (z*z)/(c*c) = 1 + * path length found for a given particle moving along a certain direction from x1 to x2 + * Typical arguments here will be x1=r_rhoMin, x2=qqqevent.pos or sx3event.pos */ - - /* - // Original circular hyperboloid implementation retained for reference. - // This section is intentionally commented out to avoid using the circular geometry. - TVector3 dx = x2-x1; // direction vector - double t2 = 1.0; //The value of 't' at the destination point, by definition: t=(z(t)-z0)/dz + TVector3 dx = x2-x1;// direction vector + double t2 = 1.0; //The value of 't' at the destination point, by definition: t=(z(t)-z0)/dz auto onesheet_hyperboloid_intersect = [&](double a, double c) { - auto A = pow(dx.Perp(),2)/(a*a) - pow(dx.Z(),2)/(c*c); - auto B = 2*(dx.X()*x1.X()+dx.Y()*x1.Y())/(a*a) - 2*(dx.Z()*x1.Z())/(c*c); - auto C = pow(x1.Perp(),2)/(a*a) - pow(x1.Z(),2)/(c*c) - 1.0; - double disc = B*B - 4*A*C; - if(disc<0) - return TVector3(0,0,54321); - else { - double tsol1 = (-B + TMath::Sqrt(disc))/(2*A); - double tsol2 = (-B - TMath::Sqrt(disc))/(2*A); - if(tsol1 >= 0 && tsol1 <= t2) - return x1+tsol1*dx; - else if(tsol2>=0 && tsol2 <= t2) - return x1+tsol2*dx; - else - return TVector3(0,0,54321); - } - }; + auto A = pow(dx.Perp(),2)/(a*a) - pow(dx.Z(),2)/(c*c); + auto B = 2*(dx.X()*x1.X()+dx.Y()*x1.Y())/(a*a) - 2*(dx.Z()*x1.Z())/(c*c); + auto C = pow(x1.Perp(),2)/(a*a) - pow(x1.Z(),2)/(c*c) - 1.0; + double disc = B*B - 4*A*C; + if(disc<0) + return TVector3(0,0,54321); + else { + double tsol1 = (-B + TMath::Sqrt(disc))/(2*A); + double tsol2 = (-B - TMath::Sqrt(disc))/(2*A); + if(tsol1 >= 0 && tsol1 <= t2) + return x1+tsol1*dx; + else if(tsol2>=0 && tsol2 <= t2) + return x1+tsol2*dx; + else + return TVector3(0,0,54321); + } + }; - TVector3 anode_intersect = onesheet_hyperboloid_intersect(32.0429,301.895); - TVector3 cathode_intersect = onesheet_hyperboloid_intersect(37.239045,301.895); - TVector3 gw_intersect = onesheet_hyperboloid_intersect(27.712,301.895); - if(anode_intersect.Z()!=54321 && cathode_intersect.Z()!=54321 && gw_intersect.Z()!=54321) - return std::tuple(cathode_intersect,gw_intersect,(cathode_intersect-gw_intersect).Mag()*0.1); - else - return std::tuple(TVector3(0,0,0), TVector3(0,0,0), 54321); - */ - - const double anode_half_width = 20.0; // 40 mm square side length - const double cathode_half_width = 21.0; // 42 mm square side length - const double guard_half_width = 18.0; // square guard opening centered on beam axis, still encloses beam - const double z_start = -100.0; // square detector opening begins here along the beam axis - const double invalid_z = 54321.0; - const double eps = 1.0e-12; - - TVector3 dx = x2 - x1; - - auto square_intersection = [&](const TVector3& start, const TVector3& direction, double half_width) { - double best_t = 1.0; - bool found = false; - - auto consider_t = [&](double t, bool x_hit, bool y_hit) { - if (t < -eps || t > 1.0 + eps) return; - - TVector3 point = start + t * direction; - double abs_x = std::abs(point.X()); - double abs_y = std::abs(point.Y()); - - bool inside_square = (abs_x <= half_width + eps) && (abs_y <= half_width + eps); - bool on_boundary = (x_hit && std::abs(abs_x - half_width) <= eps) || (y_hit && std::abs(abs_y - half_width) <= eps); - bool beyond_start = point.Z() >= z_start - eps; - - if (inside_square && on_boundary && beyond_start && t >= 0.0 && t <= 1.0) { - if (!found || t < best_t) { - best_t = t; - found = true; - } - } - }; - - if (std::abs(direction.X()) > eps) { - double t_plus_x = (half_width - start.X()) / direction.X(); - double t_minus_x = (-half_width - start.X()) / direction.X(); - consider_t(t_plus_x, true, false); - consider_t(t_minus_x, true, false); - } - - if (std::abs(direction.Y()) > eps) { - double t_plus_y = (half_width - start.Y()) / direction.Y(); - double t_minus_y = (-half_width - start.Y()) / direction.Y(); - consider_t(t_plus_y, false, true); - consider_t(t_minus_y, false, true); - } - - if (!found) { - return TVector3(0, 0, invalid_z); - } - - return start + best_t * direction; - }; - - TVector3 anode_intersect = square_intersection(x1, dx, anode_half_width); - TVector3 cathode_intersect = square_intersection(x1, dx, cathode_half_width); - TVector3 gw_intersect = square_intersection(x1, dx, guard_half_width); - - if (anode_intersect.Z() != invalid_z && cathode_intersect.Z() != invalid_z && gw_intersect.Z() != invalid_z) - return std::tuple(cathode_intersect, gw_intersect, (cathode_intersect - gw_intersect).Mag() * 0.1); - else - return std::tuple(TVector3(0, 0, 0), TVector3(0, 0, 0), invalid_z); + //TODO: Magic numbers here describing waist 'a', and flare 'c' will need better treatment. + //Currently, these are derived by fitting the crossover points to R^2/a^2 - z^2/c^2 = 1 for anodes + // Cathode a, c values are found by scaling up the anode waist by 43/37, the ratio of the outermost radii + TVector3 anode_intersect = onesheet_hyperboloid_intersect(32.0429,301.895); + TVector3 cathode_intersect = onesheet_hyperboloid_intersect(37.239045,301.895); + if(anode_intersect.Z()!=54321 && cathode_intersect.Z()!=54321) + return std::tuple(cathode_intersect,anode_intersect,(cathode_intersect-anode_intersect).Mag()*0.1); + else + return std::tuple(TVector3(0,0,0), TVector3(0,0,0), 54321); } - - - diff --git a/Armory/testMacro.C b/Armory/testMacro.C index 211b073..b7edee9 100644 --- a/Armory/testMacro.C +++ b/Armory/testMacro.C @@ -42,5 +42,12 @@ tree1->Draw("theta_recon:thetab", "sx3ID >=0 && Esx3 > 0", ""); c2->SaveAs("theta_recon_vs_thetab.png"); + new TCanvas("c3", "Canvas", 900, 600); + tree1->Draw("vZ_recon:vZ", "sx3ID >=0 && !TMath::IsNaN(vZ_recon)", ""); + TF1 *line = new TF1("line", "x", gPad->GetUxmin(), gPad->GetUxmax()); + line->SetLineColor(kRed); line->Draw("SAME"); + + c3->SaveAs("vZ_recon_vs_vZ.png"); + }