new vZ recon from known Ex

Armory/Archives/anasen_anode_cathode_hyperboloids.h

@@ -0,0 +1,108 @@
+inline std::tuple<TVector3,TVector3,double> 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);
+}
+
+
+

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);

Armory/anasen_anode_cathode_hyperboloids.h

@@ -1,108 +1,37 @@
 inline std::tuple<TVector3,TVector3,double> find_PC_PathLength(const TVector3& x1, const TVector3& x2) {
     /*
-        Square-aperture detector geometry replacing the original circular one-sheet hyperboloid.
+        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
-        The active surfaces are centered on the z axis and form squares rather than circles:
+        * path length found for a given particle moving along a certain direction from x1 to x2
-        - anodes: 40 mm x 40 mm
+        * Typical arguments here will be x1=r_rhoMin, x2=qqqevent.pos or sx3event.pos
-        - 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.
     */
-
+    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 
-    // 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 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 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;
+		auto C = pow(x1.Perp(),2)/(a*a) - pow(x1.Z(),2)/(c*c) - 1.0;
-        double disc = B*B - 4*A*C;
+		double disc = B*B - 4*A*C;
-        if(disc<0)
+		if(disc<0) 
-            return TVector3(0,0,54321);
+			return TVector3(0,0,54321);
-        else {
+		else {
-            double tsol1 = (-B + TMath::Sqrt(disc))/(2*A);
+			double tsol1 = (-B + TMath::Sqrt(disc))/(2*A);
-            double tsol2 = (-B - TMath::Sqrt(disc))/(2*A);
+			double tsol2 = (-B - TMath::Sqrt(disc))/(2*A);
-            if(tsol1 >= 0  && tsol1 <= t2)
+			if(tsol1 >= 0  && tsol1 <= t2) 
-                return x1+tsol1*dx;
+				return x1+tsol1*dx;
-            else if(tsol2>=0 && tsol2 <= t2)
+			else if(tsol2>=0 && tsol2 <= t2)
-                return x1+tsol2*dx;
+				return x1+tsol2*dx;
-            else
+			else
-                return TVector3(0,0,54321);
+				return TVector3(0,0,54321);
-        }
+		}
-    };
+	};
 
-    TVector3 anode_intersect = onesheet_hyperboloid_intersect(32.0429,301.895);
+	//TODO: Magic numbers here describing waist 'a', and flare 'c' will need better treatment. 
-    TVector3 cathode_intersect = onesheet_hyperboloid_intersect(37.239045,301.895);
+	//Currently, these are derived by fitting the crossover points to R^2/a^2 - z^2/c^2 = 1 for anodes
-    TVector3 gw_intersect = onesheet_hyperboloid_intersect(27.712,301.895);
+	// Cathode a, c values are found by scaling up the anode waist by 43/37, the ratio of the outermost radii
-    if(anode_intersect.Z()!=54321 && cathode_intersect.Z()!=54321 && gw_intersect.Z()!=54321)
+	TVector3 anode_intersect = onesheet_hyperboloid_intersect(32.0429,301.895);
-        return std::tuple(cathode_intersect,gw_intersect,(cathode_intersect-gw_intersect).Mag()*0.1);
+	TVector3 cathode_intersect = onesheet_hyperboloid_intersect(37.239045,301.895);
-    else
+	if(anode_intersect.Z()!=54321 && cathode_intersect.Z()!=54321)
-        return std::tuple(TVector3(0,0,0), TVector3(0,0,0), 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);
-    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);
 }
-
-
-

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");
+
 
 }