From 88281cd38680dc6f5f2819ed565ca6399acf65e5 Mon Sep 17 00:00:00 2001 From: tongtongcao Date: Thu, 15 Oct 2026 15:07:21 -0400 Subject: [PATCH 1/2] fix bugs for sign of d0, calculation of residual, and calculation of position and momentum at fixed radius during CVT track finding --- .../java/org/jlab/rec/cvt/fit/HelicalTrackFitter.java | 6 +++--- .../java/org/jlab/rec/cvt/track/TrackSeederXY.java | 8 +++----- .../main/java/org/jlab/rec/cvt/trajectory/Helix.java | 10 ++++------ 3 files changed, 10 insertions(+), 14 deletions(-) diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/fit/HelicalTrackFitter.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/fit/HelicalTrackFitter.java index af77f44273..3a2c7cac70 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/fit/HelicalTrackFitter.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/fit/HelicalTrackFitter.java @@ -152,7 +152,7 @@ public FitStatus fit(List X, List Y, List Z, List X, List Y, List Z, List cand, DetectorType dt) { } private double calcResi(double rho, double d0, double phi0, double xc, double yc) { double r = Math.sqrt(xc*xc+yc*yc); - double par = 1. - ((r * r - d0 * d0) * rho * rho) / (2. * (1. + d0 * Math.abs(rho))); + double par = (2. * (1. + d0 * rho) - (r * r - d0 * d0) * rho * rho) / (2. * Math.abs(1. + d0 * rho)); double newPathLength = Math.abs(Math.acos(par) / rho); int charge = (int) Math.signum(rho); double alpha = -newPathLength * rho; - double x = d0 * charge * Math.sin(phi0) + (charge / Math.abs(rho)) - * (Math.sin(phi0) - Math.cos(alpha) * Math.sin(phi0) - Math.sin(alpha) * Math.cos(phi0)); - double y = -d0 * charge * Math.cos(phi0) - (charge / Math.abs(rho)) - * (Math.cos(phi0) + Math.sin(alpha) * Math.sin(phi0) - Math.cos(alpha) * Math.cos(phi0)); + double x = d0 * Math.sin(phi0) + 1./rho * (Math.sin(phi0) - Math.sin(phi0 + alpha)); + double y = -d0 * Math.cos(phi0) - 1./rho * (Math.cos(phi0) - Math.cos(phi0 + alpha)); double res2 = ((x-xc)*(x-xc)+(y-yc)*(y-yc)); return Math.sqrt(res2); diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java index 107e5e7e2d..435fbbe2fe 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java @@ -273,15 +273,13 @@ public Point3D getPointAtRadius(double r) { } - double par = 1. - ((r * r - d0 * d0) * omega * omega) / (2. * (1. + d0 * Math.abs(omega))); + double par = (2. * (1. - d0 * omega) - (r * r - d0 * d0) * omega * omega) / (2. * Math.abs(1. - d0 * omega)); double newPathLength = Math.abs(Math.acos(par) / omega); double alpha = -newPathLength * omega; - double x = d0 * charge * Math.sin(phi0) + (charge / Math.abs(omega)) - * (Math.sin(phi0) - Math.cos(alpha) * Math.sin(phi0) - Math.sin(alpha) * Math.cos(phi0))+xb; - double y = -d0 * charge * Math.cos(phi0) - (charge / Math.abs(omega)) - * (Math.cos(phi0) + Math.sin(alpha) * Math.sin(phi0) - Math.cos(alpha) * Math.cos(phi0))+yb; + double x = xb - d0 * Math.sin(phi0) + 1./omega * (Math.sin(phi0) - Math.sin(phi0 + alpha)); + double y = yb + d0 * Math.cos(phi0) - 1./omega * (Math.cos(phi0) - Math.cos(phi0 + alpha)); double z = z0 + newPathLength * tandip; return new Point3D(x, y, z); @@ -301,7 +299,7 @@ public Vector3D getTrackDirectionAtRadius(double r) { return new Vector3D(ux, uy, uz); } - double par = 1. - ((r * r - d0 * d0) * omega * omega) / (2. * (1. + d0 * Math.abs(omega))); + double par = (2. * (1. - d0 * omega) - (r * r - d0 * d0) * omega * omega) / (2. * Math.abs(1. - d0 * omega)); double newPathLength = Math.abs(Math.acos(par) / omega); double alpha = newPathLength * omega; From 93962b99598c8a313ca49410db1382327041c855 Mon Sep 17 00:00:00 2001 From: tongtongcao Date: Sat, 17 Oct 2026 19:46:07 -0400 Subject: [PATCH 2/2] update functions to share common codes for future maintenance --- .../org/jlab/rec/cvt/track/TrackSeederXY.java | 19 +++-- .../org/jlab/rec/cvt/trajectory/Helix.java | 73 +++++++++---------- 2 files changed, 46 insertions(+), 46 deletions(-) diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/TrackSeederXY.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/TrackSeederXY.java index 87ef58d66a..b44fa3f100 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/TrackSeederXY.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/TrackSeederXY.java @@ -4,12 +4,14 @@ import java.util.HashMap; import java.util.List; import java.util.Map; +import org.jlab.geom.prim.Point3D; import org.jlab.detector.base.DetectorType; import org.jlab.rec.cvt.Constants; import org.jlab.rec.cvt.bmt.BMTType; import org.jlab.rec.cvt.cross.Cross; import org.jlab.rec.cvt.fit.CircleFitter; import org.jlab.rec.cvt.fit.CircleFitPars; +import org.jlab.rec.cvt.trajectory.Helix; public class TrackSeederXY { @@ -235,17 +237,18 @@ private int countType(List cand, DetectorType dt) { return countsvt; } + + // A method to calculate residual for a cross relative to a circle from fit private double calcResi(double rho, double d0, double phi0, double xc, double yc) { - double r = Math.sqrt(xc*xc+yc*yc); - double par = (2. * (1. + d0 * rho) - (r * r - d0 * d0) * rho * rho) / (2. * Math.abs(1. + d0 * rho)); - double newPathLength = Math.abs(Math.acos(par) / rho); - int charge = (int) Math.signum(rho); + Helix helix = new Helix(-d0, phi0, rho, 0, 0, 0, 0); // d0 from circular fitting is opposite to d0 in Helix + + double r = Math.sqrt(xc*xc+yc*yc); + double newPathLength = helix.getPathLengthAtRadius(r); double alpha = -newPathLength * rho; - - double x = d0 * Math.sin(phi0) + 1./rho * (Math.sin(phi0) - Math.sin(phi0 + alpha)); - double y = -d0 * Math.cos(phi0) - 1./rho * (Math.cos(phi0) - Math.cos(phi0 + alpha)); + + double[] xy = helix.getXYAtRadius(alpha); - double res2 = ((x-xc)*(x-xc)+(y-yc)*(y-yc)); + double res2 = ((xy[0]-xc)*(xy[0]-xc)+(xy[1]-yc)*(xy[1]-yc)); return Math.sqrt(res2); } diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java index 435fbbe2fe..503b40ba81 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java @@ -255,59 +255,56 @@ public int getCharge() { return charge; } - + + // a method to return path length at a fixed radius in xy plane + public double getPathLengthAtRadius(double r){ + double par = (2. * (1. - _dca * _curvature) - (r * r - _dca * _dca) * _curvature * _curvature) / (2. * Math.abs(1. - _dca * _curvature)); + return Math.abs(Math.acos(par) / _curvature); + } + + // a method to return a point (as a vector) at a given radial position. + public double[] getXYAtRadius(double alpha){ + double x = xb - _dca * Math.sin(_phi_at_dca) + 1./_curvature * (Math.sin(_phi_at_dca) - Math.sin(_phi_at_dca + alpha)); + double y = yb + _dca * Math.cos(_phi_at_dca) - 1./_curvature * (Math.cos(_phi_at_dca) - Math.cos(_phi_at_dca + alpha)); + return new double[]{x, y}; + } + + // a method to return a point (as a vector) at a given radial position. public Point3D getPointAtRadius(double r) { - // a method to return a point (as a vector) at a given radial position. - double d0 = _dca; - double omega = _curvature; - double charge = this.getCharge(); - double phi0 = _phi_at_dca; - double tandip = _tandip; - double z0 = _Z0; - if (Math.abs(_curvature)<1E-5) { // R > 100 m, assume it's straight track - double x = -d0 * Math.sin(phi0) + r*Math.cos(phi0)+xb; - double y = d0 * Math.cos(phi0) + r*Math.sin(phi0)+yb; - double z = z0 + r*tandip; + double x = -_dca * Math.sin(_phi_at_dca) + r*Math.cos(_phi_at_dca)+xb; + double y = _dca * Math.cos(_phi_at_dca) + r*Math.sin(_phi_at_dca)+yb; + double z = _Z0 + r*_tandip; return new Point3D(x, y, z); } + double newPathLength = getPathLengthAtRadius(r); - double par = (2. * (1. - d0 * omega) - (r * r - d0 * d0) * omega * omega) / (2. * Math.abs(1. - d0 * omega)); - double newPathLength = Math.abs(Math.acos(par) / omega); - - double alpha = -newPathLength * omega; - - double x = xb - d0 * Math.sin(phi0) + 1./omega * (Math.sin(phi0) - Math.sin(phi0 + alpha)); - double y = yb + d0 * Math.cos(phi0) - 1./omega * (Math.cos(phi0) - Math.cos(phi0 + alpha)); - double z = z0 + newPathLength * tandip; + double alpha = -newPathLength * _curvature; + + double[] xy = getXYAtRadius(alpha); + double z = _Z0 + newPathLength * _tandip; - return new Point3D(x, y, z); + return new Point3D(xy[0], xy[1], z); } - public Vector3D getTrackDirectionAtRadius(double r) { - // a method to return a point (as a vector) at a given radial position. - double d0 = _dca; - double omega = _curvature; - double phi0 = _phi_at_dca; - double tandip = _tandip; - + // a method to return a point (as a vector) at a given radial position. + public Vector3D getTrackDirectionAtRadius(double r) { if (Math.abs(_curvature) <1E-5) { // R > 100 m, assume it's straight track - double ux = Math.cos(phi0); - double uy = Math.sin(phi0); - double uz = tandip; - return new Vector3D(ux, uy, uz); + double ux = Math.cos(_phi_at_dca); + double uy = Math.sin(_phi_at_dca); + double uz = _tandip; + Vector3D trkDir = new Vector3D(ux, uy, uz); + return trkDir.asUnit(); } - double par = (2. * (1. - d0 * omega) - (r * r - d0 * d0) * omega * omega) / (2. * Math.abs(1. - d0 * omega)); - double newPathLength = Math.abs(Math.acos(par) / omega); - - double alpha = newPathLength * omega; + double newPathLength = getPathLengthAtRadius(r); + double alpha = -newPathLength * _curvature; double sintheta = Math.abs(sinTheta()); - double ux = Math.cos(-alpha + phi0) * sintheta; - double uy = Math.sin(-alpha + phi0) * sintheta; + double ux = Math.cos(_phi_at_dca + alpha) * sintheta; + double uy = Math.sin(_phi_at_dca + alpha) * sintheta; double uz = cosTheta(); Vector3D trkDir = new Vector3D(ux, uy, uz);