true if we crossed the boundary, in which case we should
+ * terminate and interpolate to the intersection.
+ */
+ public abstract boolean crossedBoundary(double newS, double[] newU);
+
+}
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12CylinderListener.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12CylinderListener.java
new file mode 100644
index 0000000000..0ac98004bd
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12CylinderListener.java
@@ -0,0 +1,57 @@
+package cnuphys.CLAS12Swim;
+
+import cnuphys.CLAS12Swim.geometry.Cylinder;
+
+/**
+ * A listener for swimming to the surface of a fixed infinite cylinder
+ */
+
+public class CLAS12CylinderListener extends CLAS12BoundaryListener {
+
+ // the target cylinder
+ private Cylinder _targetCylinder;
+
+ // starting inside or outside
+ private boolean _inside;
+
+ /**
+ * Create a CLAS12 boundary target cylinder listener, for swimming to a fixed
+ * infinite cylinder
+ *
+ * @param ivals the initial values of the swim
+ * @param targetCylinder the target infinite cylinder
+ * @param accuracy the desired accuracy (cm)
+ * @param sMax the final or max path length (cm)
+ */
+ public CLAS12CylinderListener(CLAS12Values ivals, Cylinder targetCylinder, double accuracy, double sMax) {
+ super(ivals, accuracy, sMax);
+ _targetCylinder = targetCylinder;
+ _inside = _targetCylinder.isInside(ivals.x, ivals.y, ivals.z);
+ _canMakeStraightLine = false;
+ }
+
+ @Override
+ public boolean accuracyReached(double newS, double[] newU) {
+ double dist = _targetCylinder.distance(newU[0], newU[1], newU[2]);
+ return dist < _accuracy;
+ }
+
+ @Override
+ public boolean crossedBoundary(double newS, double[] newU) {
+ boolean newInside = _targetCylinder.isInside(newU[0], newU[1], newU[2]);
+ return newInside != _inside;
+ }
+
+ /**
+ * Get the absolute distance to the target (boundary) in cm.
+ *
+ * @param newS the new path length
+ * @param newU the new state vector
+ * @return the distance to the target (boundary) in cm.
+ */
+ @Override
+ public double distanceToTarget(double newS, double[] newU) {
+ return _targetCylinder.distance(newU[0], newU[1], newU[2]);
+ }
+
+}
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12DOCAListener.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12DOCAListener.java
new file mode 100644
index 0000000000..7cf8b264ee
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12DOCAListener.java
@@ -0,0 +1,94 @@
+package cnuphys.CLAS12Swim;
+
+/**
+ * This is an abstract class to be extended by classes that to a distance of
+ * closest approach. The assumption is that the first doca is the only one. i.e.
+ * we are not dealing with low energy particles looping about.
+ */
+
+public abstract class CLAS12DOCAListener extends CLAS12Listener {
+
+ // the requested accuracy in cm
+ protected final double _accuracy;
+
+ // current doca
+ protected double _currentDOCA = Double.POSITIVE_INFINITY;
+
+ /**
+ * Create a CLAS12 boundary crossing listener
+ *
+ * @param ivals the initial values of the swim
+ * @param accuracy the accuracy (cm)
+ * @param sMax the final or max path length (cm)
+ */
+ public CLAS12DOCAListener(CLAS12Values ivals, double accuracy, double sMax) {
+ super(ivals, sMax);
+ _accuracy = accuracy;
+ }
+
+ /**
+ * Reset the current DOCA to infinity
+ */
+ @Override
+ public void reset() {
+ super.reset();
+ _currentDOCA = Double.POSITIVE_INFINITY;
+ }
+
+ /**
+ * Get the requested accuracy (on on difference in successive docas)in cm.
+ *
+ * @return the requested accuracy
+ */
+ public double getAccuracy() {
+ return _accuracy;
+ }
+
+ /**
+ * Get the current estimate of the doca
+ *
+ * @return the current doca
+ */
+ public double getCurrentDOCA() {
+ return _currentDOCA;
+ }
+
+ /**
+ * Called when a new step is taken in the ODE solving process.
+ *
+ * @param newS The new path length after the step.
+ * @param newU The new state vector after the step.
+ * @return A boolean indicating whether to continue (true) or stop (false) the
+ * integration.
+ */
+ @Override
+ public boolean newStep(double newS, double[] newU) {
+ accept(newS, newU);
+
+ double doca = doca(newS, newU);
+
+ if (doca > _currentDOCA) { // getting farther
+ _status = CLAS12Swimmer.SWIM_SUCCESS;
+ return false;
+ }
+
+ // have we reached the max path length?
+ if (newS >= _sMax) {
+ _status = CLAS12Swimmer.SWIM_TARGET_MISSED;
+ return false;
+ }
+
+ _currentDOCA = doca;
+ return true;
+ }
+
+ /**
+ * Get the absolute distance to the target (boundary) in cm.
+ *
+ * @param newS the new path length
+ * @param newU the new state vector
+ * @return the distance to the target (boundary) in cm.
+ */
+ public abstract double doca(double newS, double[] newU);
+
+}
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Listener.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Listener.java
new file mode 100644
index 0000000000..00874d54c8
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Listener.java
@@ -0,0 +1,254 @@
+package cnuphys.CLAS12Swim;
+
+/**
+ * Base ODE step listener used by {@link CLAS12Swimmer} to observe integration progress
+ * and optionally terminate the swim.
+ * + * This “basic” listener never requests termination on its own; the integrator will + * stop when the configured maximum path length {@code sMax} is reached. Specialized + * subclasses (e.g., z, rho, plane, cylinder) override termination checks to stop when + * a target condition is met within a requested accuracy. + *
+ * + *true if the listener can make a straight line to the
+ * target for a neutral particle.
+ */
+ public boolean canMakeStraightLine() {
+ return _canMakeStraightLine;
+ }
+
+ /**
+ * Get the final (target) or maximum path length in cm
+ *
+ * @return the final (target) or maximum path length in cm
+ */
+ public double getSMax() {
+ return _sMax;
+ }
+
+ /*
+ * Basic initialization and reset
+ */
+ public void reset() {
+ _status = CLAS12Swimmer.SWIM_SWIMMING;
+ _trajectory.clear();
+ _trajectory.add(0., _initialVaues.getU());
+ }
+
+ /**
+ * Called when a new step is taken in the ODE solving process.
+ *
+ * @param newS The new path length after the step.
+ * @param newU The new state vector after the step.
+ * @return A boolean indicating whether to continue (true) or stop (false) the
+ * integration.
+ */
+ @Override
+ public boolean newStep(double newS, double[] newU) {
+
+ accept(newS, newU);
+
+ // if we are done, set the status
+ if (Math.abs(newS - _sMax) < TINY) {
+ _status = CLAS12Swimmer.SWIM_SUCCESS;
+ }
+
+ // base always continues, the solve with integrate to sMax and stop
+ return true;
+ }
+
+ /**
+ * Accept the next step.
+ *
+ * @param newS The new path length after the step.
+ * @param newU The new state vector after the step.
+ */
+ protected void accept(double newS, double[] newU) {
+ _trajectory.add(newS, newU);
+ }
+
+ /**
+ * Get the trajectory
+ *
+ * @return the trajectory
+ */
+ public CLAS12Trajectory getTrajectory() {
+ return _trajectory;
+ }
+
+ /**
+ * Get the initial values
+ *
+ * @return the initial values
+ */
+ public CLAS12Values getIvals() {
+ return _initialVaues;
+ }
+
+ /**
+ * Get the current state vector
+ *
+ * @return the current state vector
+ */
+ public double[] getU() {
+ return _trajectory.get(_trajectory.size() - 1);
+ }
+
+ /**
+ * Get the state vector at the given index
+ *
+ * @param index the index
+ * @return the state vector
+ */
+ public double[] getU(int index) {
+ return _trajectory.get(index);
+ }
+
+ /**
+ * Get the number of integration steps
+ *
+ * @return the number of integration steps
+ */
+ public int getNumStep() {
+ return _trajectory.size();
+ }
+
+ /**
+ * Get the current path length
+ *
+ * @return the current path length in cm
+ */
+ public double getS() {
+ return _trajectory.getS(_trajectory.size() - 1);
+ }
+
+ /**
+ * Get the path length at the given index
+ *
+ * @param index the index
+ * @return the path length in cm
+ */
+ public double getS(int index) {
+ return _trajectory.getS(index);
+ }
+
+ /**
+ * Get the status of the swim. The values are the CLAS12Swimmer constants:
+ * SWIM_SUCCESS or SWIM_TARGET_MISSED.
+ *
+ * @return the status
+ */
+ public int getStatus() {
+ return _status;
+ }
+
+ /**
+ * Set the status
+ *
+ * @param status the status. The values are the CLAS12Swimmer constants:
+ * SWIM_SUCCESS or SWIM_TARGET_MISSED.
+ * @see CLAS12Swimmer
+ */
+ public void setStatus(int status) {
+ _status = status;
+ }
+
+ /**
+ * Get the status of the swim as a string
+ *
+ * @return the status of the swim as a string
+ */
+ public String statusString() {
+ String s = CLAS12Swimmer.resultNames.get(_status);
+ if (s == null) {
+ s = "Unknown (" + _status + ")";
+ }
+ return s;
+ }
+
+ /**
+ * Add a second point creating a straight line. This is only used when
+ * "swimming" neutral particles. This can be overridden to stop the straight
+ * line at a target.
+ */
+ public void straightLine() {
+
+ double xo = _initialVaues.x;
+ double yo = _initialVaues.y;
+ double zo = _initialVaues.z;
+ double theta = _initialVaues.theta;
+ double phi = _initialVaues.phi;
+ double sf = _sMax;
+
+ double sintheta = Math.sin(Math.toRadians(theta));
+ double costheta = Math.cos(Math.toRadians(theta));
+ double sinphi = Math.sin(Math.toRadians(phi));
+ double cosphi = Math.cos(Math.toRadians(phi));
+
+ double xf = xo + sf * sintheta * cosphi;
+ double yf = yo + sf * sintheta * sinphi;
+ double zf = zo + sf * costheta;
+
+ _trajectory.addPoint(xf, yf, zf, theta, phi, sf);
+ _status = CLAS12Swimmer.SWIM_SUCCESS;
+
+ }
+
+}
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12PlaneListener.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12PlaneListener.java
new file mode 100644
index 0000000000..4810c69f6c
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12PlaneListener.java
@@ -0,0 +1,67 @@
+package cnuphys.CLAS12Swim;
+
+import java.io.File;
+import java.io.FileNotFoundException;
+
+import cnuphys.CLAS12Swim.geometry.Plane;
+import cnuphys.magfield.MagneticFieldInitializationException;
+import cnuphys.magfield.MagneticFields;
+import cnuphys.magfield.MagneticFields.FieldType;
+
+/**
+ * A listener for swimming to the surface of a fixed infinite plane
+ */
+
+public class CLAS12PlaneListener extends CLAS12BoundaryListener {
+
+ // the target plane
+ private Plane _targetPlane;
+
+ // the starting sign. When this changes we have crossed.
+ private double _startSign;
+
+ /**
+ * Create a CLAS12 boundary target plane listener, for swimming to a fixed
+ * infinite plane
+ *
+ * @param ivals the initial values of the swim
+ * @param targetPlane the target infinite plane
+ * @param accuracy the desired accuracy (cm)
+ * @param sMax the final or max path length (cm)
+ */
+ public CLAS12PlaneListener(CLAS12Values ivals, Plane targetPlane, double accuracy, double sMax) {
+ super(ivals, accuracy, sMax);
+ _targetPlane = targetPlane;
+ _startSign = _targetPlane.sign(ivals.x, ivals.y, ivals.z);
+ _canMakeStraightLine = false;
+ }
+
+ @Override
+ public boolean crossedBoundary(double newS, double[] newU) {
+ int sign = _targetPlane.sign(newU[0], newU[1], newU[2]);
+
+ if (sign != _startSign) {
+ return true;
+ }
+ return false;
+ }
+
+ @Override
+ public boolean accuracyReached(double newS, double[] newU) {
+ double distance = _targetPlane.distance(newU[0], newU[1], newU[2]);
+ return distance < _accuracy;
+ }
+
+ /**
+ * Get the absolute distance to the target (boundary) in cm.
+ *
+ * @param newS the new path length
+ * @param newU the new state vector
+ * @return the distance to the target (boundary) in cm.
+ */
+ @Override
+ public double distanceToTarget(double newS, double[] newU) {
+ return _targetPlane.distance(newU[0], newU[1], newU[2]);
+ }
+
+}
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12RhoListener.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12RhoListener.java
new file mode 100644
index 0000000000..b524fb35cf
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12RhoListener.java
@@ -0,0 +1,133 @@
+package cnuphys.CLAS12Swim;
+
+import cnuphys.magfield.FastMath;
+
+/**
+ * A listener for swimming to a fixed cylindrical radius (rho).
+ */
+public class CLAS12RhoListener extends CLAS12BoundaryListener {
+
+ // the target rho (cm)
+ private double _rhoTarget;
+
+ // the starting sign. When this changes we have crossed.
+ private double _startSign;
+
+ /**
+ * Create a CLAS12 boundary target Z listener, for swimming to a fixed z
+ *
+ * @param ivals the initial values of the swim
+ * @param rhoTarget the target rho (cylindrical r) (cm)
+ * @param accuracy the desired accuracy (cm)
+ * @param sMax the final or max path length (cm)
+ */
+ public CLAS12RhoListener(CLAS12Values ivals, double rhoTarget, double accuracy, double sMax) {
+ super(ivals, accuracy, sMax);
+ _rhoTarget = rhoTarget;
+
+ double x = ivals.x;
+ double y = ivals.y;
+ _startSign = sign(Math.hypot(x, y));
+ }
+
+ @Override
+ public boolean crossedBoundary(double newS, double[] newU) {
+ int sign = sign(rho(newU));
+
+ if (sign != _startSign) {
+ return true;
+ }
+ return false;
+ }
+
+ // the rho (cylindrical r) of the state vector in cm
+ private double rho(double u[]) {
+ double x = u[0];
+ double y = u[1];
+ return FastMath.hypot(x, y);
+ }
+
+ @Override
+ public boolean accuracyReached(double newS, double[] newU) {
+ double dRho = Math.abs(rho(newU) - _rhoTarget);
+ return dRho < _accuracy;
+ }
+
+ // left or right of the target rho?
+ private int sign(double rho) {
+ return (rho < _rhoTarget) ? -1 : 1;
+ }
+
+ /**
+ * Get the absolute distance to the target (boundary) in cm.
+ *
+ * @param newS the new path length
+ * @param newU the new state vector
+ * @return the distance to the target (boundary) in cm.
+ */
+ @Override
+ public double distanceToTarget(double newS, double[] newU) {
+ return Math.abs(rho(newU) - _rhoTarget);
+ }
+
+ /**
+ * Add a second point creating a straight line to the target rho
+ */
+ @Override
+ public void straightLine() {
+
+ double u[] = _trajectory.get(_trajectory.size() - 1);
+ double s = _trajectory.getS(_trajectory.size() - 1);
+
+ double u2[] = findPoint(u[0], u[1], u[2], u[3], u[4], u[5], _rhoTarget);
+
+ double dx = u2[0] - u[0];
+ double dy = u2[1] - u[1];
+ double dz = u2[2] - u[2];
+ double ds = Math.sqrt(dx * dx + dy * dy + dz * dz);
+
+ _trajectory.add(s + ds, u2);
+ _status = CLAS12Swimmer.SWIM_SUCCESS;
+
+ }
+
+ private double[] findPoint(double x0, double y0, double z0, double tx, double ty, double tz, double rTarget) {
+ // Calculate the coefficients of the quadratic equation
+ double a = tx * tx + ty * ty;
+ double b = 2 * (x0 * tx + y0 * ty);
+ double c = x0 * x0 + y0 * y0 - rTarget * rTarget;
+
+ // Solve the quadratic equation
+ double discriminant = b * b - 4 * a * c;
+ if (discriminant < 0) {
+ return null; // No real solutions, rTarget cannot be reached
+ }
+
+ // Find the two possible values of s
+ double t1 = (-b + Math.sqrt(discriminant)) / (2 * a);
+ double t2 = (-b - Math.sqrt(discriminant)) / (2 * a);
+
+ // Choose the appropriate t (the one with the smaller positive value)
+
+ if (t1 < 0 && t2 < 0) {
+ return null;
+ }
+
+ double t;
+ if (t1 < 0) {
+ t = t2;
+ } else if (t2 < 0) {
+ t = t1;
+ } else {
+ t = Math.min(t1, t2);
+ }
+
+ // Calculate the resulting point
+ double x = x0 + tx * t;
+ double y = y0 + ty * t;
+ double z = z0 + tz * t;
+
+ return new double[] { x, y, z, tx, ty, tz };
+ }
+
+}
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12SphereListener.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12SphereListener.java
new file mode 100644
index 0000000000..eb939adcb9
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12SphereListener.java
@@ -0,0 +1,57 @@
+package cnuphys.CLAS12Swim;
+
+import cnuphys.CLAS12Swim.geometry.Sphere;
+
+/**
+ * A listener for swimming to the surface of a fixed sphere
+ */
+
+public class CLAS12SphereListener extends CLAS12BoundaryListener {
+
+ // the target sphere
+ private Sphere _targetSphere;
+
+ // starting inside or outside
+ private boolean _inside;
+
+ /**
+ * Create a CLAS12 boundary target sphere listener, for swimming to a fixed
+ * sphere
+ *
+ * @param ivals the initial values of the swim
+ * @param targetSphere the target infinite sphere
+ * @param accuracy the desired accuracy (cm)
+ * @param sMax the final or max path length (cm)
+ */
+ public CLAS12SphereListener(CLAS12Values ivals, Sphere targetSphere, double accuracy, double sMax) {
+ super(ivals, accuracy, sMax);
+ _targetSphere = targetSphere;
+ _inside = targetSphere.isInside(ivals.x, ivals.y, ivals.z);
+ _canMakeStraightLine = false;
+ }
+
+ @Override
+ public boolean accuracyReached(double newS, double[] newU) {
+ double dist = _targetSphere.distance(newU[0], newU[1], newU[2]);
+ return dist < _accuracy;
+ }
+
+ @Override
+ public boolean crossedBoundary(double newS, double[] newU) {
+ boolean newInside = _targetSphere.isInside(newU[0], newU[1], newU[2]);
+ return newInside != _inside;
+ }
+
+ /**
+ * Get the absolute distance to the target (boundary) in cm.
+ *
+ * @param newS the new path length
+ * @param newU the new state vector
+ * @return the distance to the target (boundary) in cm.
+ */
+ @Override
+ public double distanceToTarget(double newS, double[] newU) {
+ return _targetSphere.distance(newU[0], newU[1], newU[2]);
+ }
+
+}
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12SwimResult.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12SwimResult.java
new file mode 100644
index 0000000000..63fc1ecfca
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12SwimResult.java
@@ -0,0 +1,209 @@
+package cnuphys.CLAS12Swim;
+
+/**
+ * Result container returned by all {@code CLAS12Swimmer} swim operations.
+ * + * A {@code CLAS12SwimResult} encapsulates: + *
+ * This class is a passive data container and is not thread-safe for mutation. + *
+ */ + +public class CLAS12SwimResult { + + private CLAS12Listener _listener; + + public CLAS12SwimResult(CLAS12Listener listener) { + _listener = listener; + } + + /** + * Get the trajectory + * + * @return the trajectory + */ + public CLAS12Trajectory getTrajectory() { + return _listener.getTrajectory(); + } + + /** + * Get the initial values of the swim + * + * @return the initial values + */ + public CLAS12Values getInitialValues() { + return _listener.getIvals(); + } + + /** + * Get the final values of the swim + * + * @return the final values + */ + public CLAS12Values getFinalValues() { + double u[] = _listener.getU(); + int q = _listener.getIvals().q; + double p = _listener.getIvals().p; + return new CLAS12Values(q, p, u); + } + + /** + * Get the path length in cm + * + * @return the path length in cm + */ + public double getPathLength() { + return _listener.getS(); + } + + /** + * Get the termination status code for the swim. + *+ * The value is one of the {@link CLAS12Swimmer} status constants, e.g. + * {@link CLAS12Swimmer#SWIM_SUCCESS}, {@link CLAS12Swimmer#SWIM_TARGET_MISSED}, + * or {@link CLAS12Swimmer#BELOW_MIN_MOMENTUM}. + *
+ * + * @return integer status code indicating how the swim terminated + */ + public int getStatus() { + return _listener.getStatus(); + } + + /** + * Get a copy of the final state vector. + *+ * The state vector {@code u} has length 6 and is interpreted as: + *
+ * The returned array is a defensive copy and may be modified by the caller. + *
+ * + * @return a copy of the final state vector + */ + public double[] getFinalU() { + double[] u = _listener.getU(); + return (u == null) ? null : u.clone(); + } + + /** + * Determine whether the swim terminated successfully. + *+ * A successful termination indicates that the swimmer reached the requested target condition + * (such as a surface, target {@code z}, target {@code ρ}, distance of closest approach, + * or maximum path length) within the specified accuracy, and without encountering an internal + * failure condition. + *
+ *+ * If this method returns {@code false}, the final state stored in this result still represents + * the particle state at the point where the swim terminated (for example, due to exceeding + * {@code sMax} or failing to converge). + *
+ * + * @return {@code true} if the swim terminated successfully; {@code false} otherwise + */ + public boolean isSuccess() { + return getStatus() == CLAS12Swimmer.SWIM_SUCCESS; + } + + /** + * Get the final rho in cm + * + * @return the final rho in cm + */ + public double getFinalRho() { + return Math.hypot(_listener.getU()[0], _listener.getU()[1]); + } + + /** + * Get the status of the swim as a string + * + * @return the status of the swim as a string + */ + public String statusString() { + int status = getStatus(); + String s = CLAS12Swimmer.resultNames.get(status); + if (s == null) { + s = "Unknown (" + status + ")"; + } + return s; + } + + /** + * Get the number of integration steps + * + * @return the number of integration steps + */ + public int getNStep() { + return _listener.getNumStep(); + } + + /** + * Get a summary of the results of the swim + */ + @Override + public String toString() { + StringBuffer sb = new StringBuffer(2000); + CLAS12Values ivalues = getInitialValues(); + CLAS12Values fvalues = getFinalValues(); + + double norm = ivalues.p / fvalues.p; // should be 1.0 + + sb.append("Swim results:\n"); + sb.append("Status: " + statusString() + "\n"); + sb.append("Initial values:\n"); + sb.append("charge = " + ivalues.q + "\n"); + sb.append(String.format("vertex = (%.4f, %.4f, %.4f) cm\n", ivalues.x, ivalues.y, ivalues.z)); + + sb.append(String.format("momentum = %.4f GeV/c\n", ivalues.p)); + sb.append(String.format("theta = %.4f deg\n", ivalues.theta)); + sb.append(String.format("phi = %.4f deg\n", ivalues.phi)); + sb.append("--------\nFinal values:\n"); + sb.append(String.format("location = (%.4f, %.4f, %.4f) cm\n", fvalues.x, fvalues.y, fvalues.z)); + sb.append(String.format("momentum = %.4f GeV/c\n", fvalues.p)); + sb.append(String.format("norm = %.4f (should be 1)\n", norm)); + sb.append(String.format("theta = %.4f deg\n", fvalues.theta)); + sb.append(String.format("phi = %.4f deg\n", fvalues.phi)); + sb.append(String.format("rho = %.4f cm\n", Math.hypot(fvalues.x, fvalues.y))); + sb.append(String.format("path length = %.4f cm\n", getPathLength())); + sb.append(String.format("number of steps = %d\n", getNStep())); + + return sb.toString(); + } + +} diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Swimmer.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Swimmer.java new file mode 100644 index 0000000000..97e86ced83 --- /dev/null +++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Swimmer.java @@ -0,0 +1,1105 @@ +package cnuphys.CLAS12Swim; + +import java.util.Hashtable; + +import cnuphys.CLAS12Swim.geometry.Cylinder; +import cnuphys.CLAS12Swim.geometry.Plane; +import cnuphys.CLAS12Swim.geometry.Sphere; +import cnuphys.magfield.FieldProbe; +import cnuphys.magfield.RotatedCompositeProbe; + +import org.apache.commons.math3.ode.FirstOrderDifferentialEquations; +import org.apache.commons.math3.ode.events.EventHandler; +import org.apache.commons.math3.ode.sampling.StepHandler; +import org.apache.commons.math3.ode.sampling.StepInterpolator; +import org.apache.commons.math3.ode.nonstiff.ClassicalRungeKuttaIntegrator; +import org.apache.commons.math3.ode.nonstiff.DormandPrince54Integrator; + +/** + * The CLAS12 swimmer implementation, based on the Apache Commons Math 3.6.1 + * ODE solvers. + * + *Units: positions in cm, momentum in GeV/c, angles in degrees, path length in cm.
+ */ + @Override + public CLAS12SwimResult swimZ(int q, + double xo, double yo, double zo, + double p, double theta, double phi, + double zTarget, double accuracy, + double sMax, double h, double tolerance) { + + final CLAS12Values ivals = new CLAS12Values(q, xo, yo, zo, p, theta, phi); + final CLAS12ZListener listener = new CLAS12ZListener(ivals, zTarget, accuracy, sMax); + + // Momentum guard + if (p < minMomentum) { + listener.setStatus(CLAS12Swimmer.BELOW_MIN_MOMENTUM); + return new CLAS12SwimResult(listener); + } + + // Neutral shortcut + if (q == 0 && listener.canMakeStraightLine()) { + listener.straightLine(); + return new CLAS12SwimResult(listener); + } + + // Initial state y = [x,y,z, tx,ty,tz] + final double[] y = ivals.getU().clone(); + final SwimEquations ode = new SwimEquations(q, p, probe); + + final double targetMiss = accuracy; // cm + final double successTol = accuracy; // cm + + // Commons Math requires per-component tolerances. Interpret the provided "tolerance" knob + // as a position absolute tolerance in cm. + final double absPos = Math.max(1.0e-12, tolerance); // cm + final double absDir = 1.0e-10; // dimensionless + final double rel = 1.0e-12; + + final double[] absTol = new double[] { absPos, absPos, absPos, absDir, absDir, absDir }; + final double[] relTol = new double[] { rel, rel, rel, rel, rel, rel }; + + // Adaptive integrator: allow step-size growth up to maxStepSize (do NOT cap at h) + final DormandPrince54Integrator integrator = + new DormandPrince54Integrator( + Math.max(minStepSize, 1e-12), + Math.max(maxStepSize, minStepSize), + absTol, + relTol + ); + + // Use h only as an initial step-size guess + final double h0 = Math.max(minStepSize, Math.min(Math.abs(h), maxStepSize)); + integrator.setInitialStepSize(h0); + + // Record the trajectory at each accepted step + integrator.addStepHandler(new StepHandler() { + @Override + public void init(double s0, double[] y0, double sEnd) { + // no-op + } + + @Override + public void handleStep(StepInterpolator interpolator, boolean isLast) { + final double s = interpolator.getCurrentTime(); + final double[] state = interpolator.getInterpolatedState().clone(); + listener.accept(s, state); + } + }); + + // Event: stop at zTarget + final HitFlag hit = new HitFlag(); + + final EventHandler zEvent = new EventHandler() { + @Override + public void init(double s0, double[] y0, double sEnd) { } + + @Override + public double g(double s, double[] y) { + return y[2] - zTarget; + } + + @Override + public Action eventOccurred(double s, double[] y, boolean increasing) { + hit.hit = true; + return Action.STOP; + } + + @Override + public void resetState(double s, double[] y) { } + }; + + final double maxCheckInterval = Math.max(0.5, h0); + final double eventConv = Math.max(1.0e-12, targetMiss); + integrator.addEventHandler(zEvent, maxCheckInterval, eventConv, 200); + + // Integrate + double sFinal; + try { + sFinal = integrator.integrate(ode, 0.0, y, sMax, y); + } catch (Exception ex) { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + return new CLAS12SwimResult(listener); + } + + // Ensure final point captured even if the last step handler didn't run as expected + listener.accept(sFinal, y.clone()); + + // Status + if (hit.hit && Math.abs(listener.getU()[2] - zTarget) <= successTol) { + listener.setStatus(CLAS12Swimmer.SWIM_SUCCESS); + } else { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + } + + return new CLAS12SwimResult(listener); + } + + + // ------------------------------------------------------------------------- + // The rest of the interface: skeleton placeholders + // ------------------------------------------------------------------------- + + @Override + public CLAS12SwimResult swim(int q, double xo, double yo, double zo, double p, double theta, double phi, + double sMax, double h, double tolerance) { + + final CLAS12Values ivals = new CLAS12Values(q, xo, yo, zo, p, theta, phi); + final CLAS12Listener listener = new CLAS12Listener(ivals, sMax); + + // Momentum guard + if (p < minMomentum) { + listener.setStatus(CLAS12Swimmer.BELOW_MIN_MOMENTUM); + return new CLAS12SwimResult(listener); + } + + // Neutral shortcut: if permitted, do an exact straight-line propagation to sMax + if (q == 0 && listener.canMakeStraightLine()) { + listener.straightLine(); + return new CLAS12SwimResult(listener); + } + + // Initial state y0 = [x,y,z, tx,ty,tz] + final double[] y0 = ivals.getU(); + final double[] y = y0.clone(); + + final FirstOrderDifferentialEquations ode = new SwimEquations(q, p, probe); + + // Per-component tolerances + final double absPos = Math.max(1.0e-12, tolerance); // cm + final double absDir = 1.0e-10; // dimensionless + final double rel = 1.0e-12; + + final double[] absTol = new double[] { absPos, absPos, absPos, absDir, absDir, absDir }; + final double[] relTol = new double[] { rel, rel, rel, rel, rel, rel }; + + // Adaptive integrator: allow step-size growth up to maxStepSize; use h only as initial guess. + final DormandPrince54Integrator integrator = + new DormandPrince54Integrator( + Math.max(minStepSize, 1e-12), + Math.max(maxStepSize, minStepSize), + absTol, + relTol + ); + + final double h0 = Math.max(minStepSize, Math.min(Math.abs(h), maxStepSize)); + integrator.setInitialStepSize(h0); + + // Record trajectory at accepted steps + integrator.addStepHandler(new StepHandler() { + @Override + public void init(double t0, double[] y0, double t) { + // listener.reset() already added the initial point + } + + @Override + public void handleStep(StepInterpolator interpolator, boolean isLast) { + final double s = interpolator.getCurrentTime(); + final double[] state = interpolator.getInterpolatedState().clone(); + listener.accept(s, state); + } + }); + + try { + integrator.integrate(ode, 0.0, y, sMax, y); + } catch (Exception ex) { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + return new CLAS12SwimResult(listener); + } + + // For the basic swim, reaching sMax is considered success. + listener.setStatus(CLAS12Swimmer.SWIM_SUCCESS); + return new CLAS12SwimResult(listener); + } + + @Override + public CLAS12SwimResult swimFixed(int q, double xo, double yo, double zo, double p, double theta, double phi, + double sMax, double h) { + if (h <= 0.0) { + throw new IllegalArgumentException("Fixed step size must be positive"); + } + + final CLAS12Values ivals = new CLAS12Values(q, xo, yo, zo, p, theta, phi); + final CLAS12Listener listener = new CLAS12Listener(ivals, sMax); + + if (p < minMomentum) { + listener.setStatus(CLAS12Swimmer.BELOW_MIN_MOMENTUM); + return new CLAS12SwimResult(listener); + } + + if (q == 0 && listener.canMakeStraightLine()) { + listener.straightLine(); + return new CLAS12SwimResult(listener); + } + + final double[] y = ivals.getU().clone(); + final FirstOrderDifferentialEquations ode = new SwimEquations(q, p, probe); + final ClassicalRungeKuttaIntegrator integrator = new ClassicalRungeKuttaIntegrator(h); + + integrator.addStepHandler(new StepHandler() { + @Override + public void init(double s0, double[] y0, double sEnd) { + // listener reset already recorded the initial point + } + + @Override + public void handleStep(StepInterpolator interpolator, boolean isLast) { + listener.accept(interpolator.getCurrentTime(), interpolator.getInterpolatedState().clone()); + } + }); + + try { + integrator.integrate(ode, 0.0, y, sMax, y); + listener.setStatus(CLAS12Swimmer.SWIM_SUCCESS); + } catch (Exception ex) { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + } + + return new CLAS12SwimResult(listener); + } + + @Override + public CLAS12SwimResult swimCylinder(int q, double xo, double yo, double zo, double p, double theta, double phi, + double[] p1, double[] p2, double r, double accuracy, double sMax, double h, + double tolerance) { + Cylinder targetCylinder = new Cylinder(p1, p2, r); + return swimCylinder(q, xo, yo, zo, p, theta, phi, targetCylinder, accuracy, sMax, h, tolerance); + } + + @Override + public CLAS12SwimResult swimCylinder(int q, double xo, double yo, double zo, double p, double theta, double phi, + Cylinder targetCylinder, double accuracy, double sMax, double h, + double tolerance) { + + // If the target cylinder is centered on the z axis, this is exactly a rho swim. + if (targetCylinder.centeredOnZ()) { + return swimRho(q, xo, yo, zo, p, theta, phi, targetCylinder.radius, accuracy, sMax, h, tolerance); + } + + final CLAS12Values ivals = new CLAS12Values(q, xo, yo, zo, p, theta, phi); + final CLAS12CylinderListener listener = + new CLAS12CylinderListener(ivals, targetCylinder, accuracy, sMax); + + // Momentum guard + if (p < minMomentum) { + listener.setStatus(CLAS12Swimmer.BELOW_MIN_MOMENTUM); + return new CLAS12SwimResult(listener); + } + + // Neutral shortcut is intentionally not used here unless the cylinder reduces to rho. + // CLAS12CylinderListener disables exact straight-line handling. + + final double[] y = ivals.getU().clone(); + final SwimEquations ode = new SwimEquations(q, p, probe); + + final double targetMiss = accuracy; + final double successTol = accuracy; + + final double absPos = Math.max(1.0e-12, tolerance); + final double absDir = 1.0e-10; + final double rel = 1.0e-12; + + final double[] absTol = { absPos, absPos, absPos, absDir, absDir, absDir }; + final double[] relTol = { rel, rel, rel, rel, rel, rel }; + + final DormandPrince54Integrator integrator = + new DormandPrince54Integrator( + Math.max(minStepSize, 1.0e-12), + Math.max(maxStepSize, minStepSize), + absTol, + relTol + ); + + final double h0 = Math.max(minStepSize, Math.min(Math.abs(h), maxStepSize)); + integrator.setInitialStepSize(h0); + + integrator.addStepHandler(new StepHandler() { + @Override + public void init(double s0, double[] y0, double sEnd) { } + + @Override + public void handleStep(StepInterpolator interpolator, boolean isLast) { + final double s = interpolator.getCurrentTime(); + final double[] state = interpolator.getInterpolatedState().clone(); + listener.accept(s, state); + } + }); + + // Cylinder event: signed distance to surface = 0. + // Keep checks reasonably frequent to reduce any chance of stepping over the cylinder. + final HitFlag hit = new HitFlag(); + + final EventHandler cylinderEvent = new EventHandler() { + @Override + public void init(double s0, double[] y0, double sEnd) { } + + @Override + public double g(double s, double[] y) { + return targetCylinder.signedDistance(y[0], y[1], y[2]); + } + + @Override + public Action eventOccurred(double s, double[] y, boolean increasing) { + hit.hit = true; + return Action.STOP; + } + + @Override + public void resetState(double s, double[] y) { } + }; + + final double maxCheckInterval = Math.min(Math.max(0.5, h0), 5.0); + final double eventConv = Math.max(1.0e-12, targetMiss); + + integrator.addEventHandler(cylinderEvent, maxCheckInterval, eventConv, 200); + + double sFinal; + try { + sFinal = integrator.integrate(ode, 0.0, y, sMax, y); + } catch (Exception ex) { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + return new CLAS12SwimResult(listener); + } + + // Ensure final point is recorded + listener.accept(sFinal, y.clone()); + + // Status + double dist = targetCylinder.distance(listener.getU()[0], + listener.getU()[1], + listener.getU()[2]); + + if (hit.hit && dist <= successTol) { + listener.setStatus(CLAS12Swimmer.SWIM_SUCCESS); + } else { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + } + + return new CLAS12SwimResult(listener); + } + + @Override + public CLAS12SwimResult swimSphere(int q, double xo, double yo, double zo, double p, double theta, double phi, + double[] center, double r, double accuracy, double sMax, double h, + double tolerance) { + return swimSphere(q, xo, yo, zo, p, theta, phi, new Sphere(center, r), accuracy, sMax, h, tolerance); + } + + @Override + public CLAS12SwimResult swimSphere(int q, double xo, double yo, double zo, double p, double theta, double phi, + Sphere targetSphere, double accuracy, double sMax, double h, + double tolerance) { + final CLAS12Values ivals = new CLAS12Values(q, xo, yo, zo, p, theta, phi); + final CLAS12SphereListener listener = new CLAS12SphereListener(ivals, targetSphere, accuracy, sMax); + + if (p < minMomentum) { + listener.setStatus(BELOW_MIN_MOMENTUM); + return new CLAS12SwimResult(listener); + } + + final double[] y = ivals.getU().clone(); + final SwimEquations ode = new SwimEquations(q, p, probe); + final DormandPrince54Integrator integrator = createAdaptiveIntegrator(h, tolerance); + final HitFlag hit = new HitFlag(); + + integrator.addStepHandler(recordSteps(listener)); + integrator.addEventHandler(new EventHandler() { + @Override public void init(double s0, double[] y0, double sEnd) { } + @Override public double g(double s, double[] state) { + return targetSphere.signedDistance(state[0], state[1], state[2]); + } + @Override public Action eventOccurred(double s, double[] state, boolean increasing) { + hit.hit = true; + return Action.STOP; + } + @Override public void resetState(double s, double[] state) { } + }, Math.max(0.5, initialStep(h)), Math.max(1.0e-12, accuracy), 200); + + try { + double sFinal = integrator.integrate(ode, 0.0, y, sMax, y); + listener.accept(sFinal, y.clone()); + } catch (Exception ex) { + listener.setStatus(SWIM_TARGET_MISSED); + return new CLAS12SwimResult(listener); + } + + listener.setStatus(hit.hit && targetSphere.distance(y[0], y[1], y[2]) <= accuracy + ? SWIM_SUCCESS : SWIM_TARGET_MISSED); + return new CLAS12SwimResult(listener); + } + + @Override + public CLAS12SwimResult swimPlane(int q, double xo, double yo, double zo, double p, double theta, double phi, + double nx, double ny, double nz, double px, double py, double pz, + double accuracy, double sMax, double h, double tolerance) { + return swimPlane(q, xo, yo, zo, p, theta, phi, new Plane(nx, ny, nz, px, py, pz), accuracy, sMax, h, + tolerance); + } + + @Override + public CLAS12SwimResult swimPlane(int q, double xo, double yo, double zo, double p, double theta, double phi, + double[] norm, double[] point, double accuracy, double sMax, double h, + double tolerance) { + return swimPlane(q, xo, yo, zo, p, theta, phi, new Plane(norm, point), accuracy, sMax, h, tolerance); + } + + /** + * Swim to a fixed plane, stopping when the trajectory intersects the plane + * or when the path length {@code sMax} is reached. + * + *Units: positions in cm, momentum in GeV/c, angles in degrees.
+ */ + @Override + public CLAS12SwimResult swimPlane(int q, + double xo, double yo, double zo, + double p, double theta, double phi, + Plane plane, + double accuracy, + double sMax, + double h, + double tolerance) { + + final CLAS12Values ivals = new CLAS12Values(q, xo, yo, zo, p, theta, phi); + final CLAS12PlaneListener listener = + new CLAS12PlaneListener(ivals, plane, accuracy, sMax); + + // Momentum guard + if (p < minMomentum) { + listener.setStatus(CLAS12Swimmer.BELOW_MIN_MOMENTUM); + return new CLAS12SwimResult(listener); + } + + // Neutral shortcut + if (q == 0 && listener.canMakeStraightLine()) { + listener.straightLine(); + return new CLAS12SwimResult(listener); + } + + // Initial state y = [x,y,z, tx,ty,tz] + final double[] y = ivals.getU().clone(); + final SwimEquations ode = new SwimEquations(q, p, probe); + + final double targetMiss = accuracy; + final double successTol = accuracy; + final double absPos = Math.max(1.0e-12, tolerance); + final double absDir = 1.0e-10; + final double rel = 1.0e-12; + + final double[] absTol = { + absPos, absPos, absPos, + absDir, absDir, absDir + }; + final double[] relTol = { + rel, rel, rel, rel, rel, rel + }; + + final DormandPrince54Integrator integrator = + new DormandPrince54Integrator( + Math.max(minStepSize, 1.0e-12), + Math.max(maxStepSize, minStepSize), + absTol, + relTol + ); + + // Initial step-size guess + final double h0 = Math.max(minStepSize, Math.min(Math.abs(h), maxStepSize)); + integrator.setInitialStepSize(h0); + + // Record every accepted step (matches swimZ behavior) + integrator.addStepHandler(new StepHandler() { + @Override + public void init(double s0, double[] y0, double sEnd) { } + + @Override + public void handleStep(StepInterpolator interpolator, boolean isLast) { + final double s = interpolator.getCurrentTime(); + final double[] state = interpolator.getInterpolatedState().clone(); + listener.accept(s, state); + } + }); + + // ------------------------------------------------------------ + // Plane event: signed distance = 0 + // ------------------------------------------------------------ + final HitFlag hit = new HitFlag(); + + final EventHandler planeEvent = new EventHandler() { + + @Override + public void init(double s0, double[] y0, double sEnd) { } + + @Override + public double g(double s, double[] y) { + return plane.signedDistance(y[0], y[1], y[2]); + } + + @Override + public Action eventOccurred(double s, double[] y, boolean increasing) { + hit.hit = true; + return Action.STOP; + } + + @Override + public void resetState(double s, double[] y) { } + }; + + final double maxCheckInterval = Math.max(0.5, h0); + final double eventConv = Math.max(1.0e-12, targetMiss); + + integrator.addEventHandler(planeEvent, maxCheckInterval, eventConv, 200); + + // ------------------------------------------------------------ + // Integrate + // ------------------------------------------------------------ + double sFinal; + try { + sFinal = integrator.integrate(ode, 0.0, y, sMax, y); + } catch (Exception ex) { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + return new CLAS12SwimResult(listener); + } + + // Ensure final point is recorded + listener.accept(sFinal, y.clone()); + + // ------------------------------------------------------------ + // Status + // ------------------------------------------------------------ + double dist = plane.distance(listener.getU()[0], + listener.getU()[1], + listener.getU()[2]); + + if (hit.hit && dist <= successTol) { + listener.setStatus(CLAS12Swimmer.SWIM_SUCCESS); + } else { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + } + + return new CLAS12SwimResult(listener); + } + + /** + * Swim to a fixed target z (cm) in a given CLAS12 sector using the RotatedComposite field. + *+ * This operation is only valid when the active probe is a + * {@link RotatedCompositeProbe}. If not, it prints an error and returns {@code null}. + *
+ */ + @Override + public CLAS12SwimResult sectorSwimZ(int sector, int q, + double xo, double yo, double zo, + double p, double theta, double phi, + double zTarget, double accuracy, + double sMax, double h, double tolerance) { + + // Must use the rotated field. + if (!(probe instanceof RotatedCompositeProbe)) { + System.err.println("sectorSwimZ only valid with RotatedCompositeProbe."); + return null; + } + + final CLAS12Values ivals = new CLAS12Values(q, xo, yo, zo, p, theta, phi); + final CLAS12ZListener listener = new CLAS12ZListener(ivals, zTarget, accuracy, sMax); + + if (p < minMomentum) { + listener.setStatus(CLAS12Swimmer.BELOW_MIN_MOMENTUM); + return new CLAS12SwimResult(listener); + } + + if (q == 0 && listener.canMakeStraightLine()) { + listener.straightLine(); + return new CLAS12SwimResult(listener); + } + + final double[] y = ivals.getU().clone(); + + final SectorSwimEquations ode = + new SectorSwimEquations(sector, q, p, probe); + + final double targetMiss = accuracy; + final double successTol = accuracy; + final double absPos = Math.max(1.0e-12, tolerance); + final double absDir = 1.0e-10; + final double rel = 1.0e-12; + + final double[] absTol = { absPos, absPos, absPos, absDir, absDir, absDir }; + final double[] relTol = { rel, rel, rel, rel, rel, rel }; + + final DormandPrince54Integrator integrator = + new DormandPrince54Integrator( + Math.max(minStepSize, 1.0e-12), + Math.max(maxStepSize, minStepSize), + absTol, + relTol + ); + + final double h0 = Math.max(minStepSize, Math.min(Math.abs(h), maxStepSize)); + integrator.setInitialStepSize(h0); + + integrator.addStepHandler(new StepHandler() { + @Override + public void init(double s0, double[] y0, double sEnd) { } + + @Override + public void handleStep(StepInterpolator interpolator, boolean isLast) { + double s = interpolator.getCurrentTime(); + double[] state = interpolator.getInterpolatedState().clone(); + listener.accept(s, state); + } + }); + + final HitFlag hit = new HitFlag(); + + EventHandler zEvent = new EventHandler() { + + @Override + public void init(double s0, double[] y0, double sEnd) { } + + @Override + public double g(double s, double[] y) { + return y[2] - zTarget; + } + + @Override + public Action eventOccurred(double s, double[] y, boolean increasing) { + hit.hit = true; + return Action.STOP; + } + + @Override + public void resetState(double s, double[] y) { } + }; + + double maxCheckInterval = Math.max(0.5, h0); + double eventConv = Math.max(1.0e-12, targetMiss); + + integrator.addEventHandler(zEvent, maxCheckInterval, eventConv, 200); + + double sFinal; + + try { + sFinal = integrator.integrate(ode, 0.0, y, sMax, y); + } catch (Exception ex) { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + return new CLAS12SwimResult(listener); + } + + listener.accept(sFinal, y.clone()); + + if (hit.hit && Math.abs(listener.getU()[2] - zTarget) <= successTol) { + listener.setStatus(CLAS12Swimmer.SWIM_SUCCESS); + } else { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + } + + return new CLAS12SwimResult(listener); + } + + + @Override + public CLAS12SwimResult swimRho(int q, double xo, double yo, double zo, double p, double theta, double phi, + double rhoTarget, double accuracy, double sMax, double h, double tolerance) { + + final CLAS12Values ivals = new CLAS12Values(q, xo, yo, zo, p, theta, phi); + final CLAS12RhoListener listener = new CLAS12RhoListener(ivals, rhoTarget, accuracy, sMax); + + // Momentum guard + if (p < minMomentum) { + listener.setStatus(CLAS12Swimmer.BELOW_MIN_MOMENTUM); + return new CLAS12SwimResult(listener); + } + + // Neutral shortcut: listener can compute exact straight-line intersection with rho target + if (q == 0 && listener.canMakeStraightLine()) { + listener.straightLine(); + return new CLAS12SwimResult(listener); + } + + // Initial state y0 = [x,y,z, tx,ty,tz] + final double[] y0 = ivals.getU(); + final double[] y = y0.clone(); + + final FirstOrderDifferentialEquations ode = new SwimEquations(q, p, probe); + + final double targetMiss = accuracy; + final double absPos = Math.max(1e-12, tolerance); + final double absDir = 1e-10; + final double rel = 1e-12; + + final double[] absTol = new double[] { absPos, absPos, absPos, absDir, absDir, absDir }; + final double[] relTol = new double[] { rel, rel, rel, rel, rel, rel }; + + final DormandPrince54Integrator integrator = + new DormandPrince54Integrator( + Math.max(minStepSize, 1e-12), + Math.max(maxStepSize, minStepSize), + absTol, + relTol + ); + + final double h0 = Math.max(minStepSize, Math.min(Math.abs(h), maxStepSize)); + integrator.setInitialStepSize(h0); + + // Record trajectory at accepted steps + integrator.addStepHandler(new StepHandler() { + @Override + public void init(double t0, double[] y0, double t) { + // listener.reset() already added the initial point + } + + @Override + public void handleStep(StepInterpolator interpolator, boolean isLast) { + final double s = interpolator.getCurrentTime(); + final double[] state = interpolator.getInterpolatedState().clone(); + listener.accept(s, state); + } + }); + + // Event: rho(s) - rhoTarget = 0, where rho = sqrt(x^2 + y^2) + final HitFlag hit = new HitFlag(); + + final EventHandler rhoEventHandler = new EventHandler() { + @Override + public void init(double t0, double[] y0, double t) { + // nothing + } + + @Override + public double g(double s, double[] y) { + final double rho = Math.hypot(y[0], y[1]); + return rho - rhoTarget; + } + + @Override + public Action eventOccurred(double s, double[] y, boolean increasing) { + hit.hit = true; + return Action.STOP; + } + + @Override + public void resetState(double s, double[] y) { + // no reset + } + }; + + integrator.addEventHandler( + rhoEventHandler, + Math.max(0.5, h0), + Math.max(1e-12, targetMiss), + 200 + ); + + try { + integrator.integrate(ode, 0.0, y, sMax, y); + } catch (Exception ex) { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + return new CLAS12SwimResult(listener); + } + + // Decide final status + final double rhoFinal = Math.hypot(listener.getU()[0], listener.getU()[1]); + if (hit.hit && Math.abs(rhoFinal - rhoTarget) <= targetMiss) { + listener.setStatus(CLAS12Swimmer.SWIM_SUCCESS); + } else { + listener.setStatus(CLAS12Swimmer.SWIM_TARGET_MISSED); + } + + return new CLAS12SwimResult(listener); + } + + @Override + public CLAS12SwimResult swimZLine(int q, double xo, double yo, double zo, double p, double theta, double phi, + double xb, double yb, double accuracy, double sMax, double h, double tolerance) { + CLAS12Values ivals = new CLAS12Values(q, xo, yo, zo, p, theta, phi); + return swimToClosestApproach(ivals, new CLAS12ZLineListener(ivals, xb, yb, accuracy, sMax), + xb, yb, accuracy, sMax, h, tolerance); + } + + @Override + public CLAS12SwimResult swimBeamline(int q, double xo, double yo, double zo, double p, double theta, double phi, + double accuracy, double sMax, double h, double tolerance) { + CLAS12Values ivals = new CLAS12Values(q, xo, yo, zo, p, theta, phi); + return swimToClosestApproach(ivals, new CLAS12BeamlineListener(ivals, accuracy, sMax), + 0.0, 0.0, accuracy, sMax, h, tolerance); + } + + // ------------------------------------------------------------------------- + // Internal helpers + // ------------------------------------------------------------------------- + + private CLAS12SwimResult swimToClosestApproach(CLAS12Values ivals, CLAS12DOCAListener listener, + double xb, double yb, double accuracy, double sMax, + double h, double tolerance) { + if (ivals.p < minMomentum) { + listener.setStatus(BELOW_MIN_MOMENTUM); + return new CLAS12SwimResult(listener); + } + + final double[] y = ivals.getU().clone(); + if ((y[0] - xb) * y[3] + (y[1] - yb) * y[4] >= 0.0) { + listener.setStatus(SWIM_SUCCESS); + return new CLAS12SwimResult(listener); + } + + final DormandPrince54Integrator integrator = createAdaptiveIntegrator(h, tolerance); + final HitFlag hit = new HitFlag(); + integrator.addStepHandler(recordSteps(listener)); + integrator.addEventHandler(new EventHandler() { + @Override public void init(double s0, double[] y0, double sEnd) { } + @Override public double g(double s, double[] state) { + return (state[0] - xb) * state[3] + (state[1] - yb) * state[4]; + } + @Override public Action eventOccurred(double s, double[] state, boolean increasing) { + hit.hit = true; + return Action.STOP; + } + @Override public void resetState(double s, double[] state) { } + }, Math.max(0.5, initialStep(h)), Math.max(1.0e-12, accuracy), 200); + + try { + double sFinal = integrator.integrate(new SwimEquations(ivals.q, ivals.p, probe), 0.0, y, sMax, y); + listener.accept(sFinal, y.clone()); + } catch (Exception ex) { + listener.setStatus(SWIM_TARGET_MISSED); + return new CLAS12SwimResult(listener); + } + + listener.setStatus(hit.hit ? SWIM_SUCCESS : SWIM_TARGET_MISSED); + return new CLAS12SwimResult(listener); + } + + private DormandPrince54Integrator createAdaptiveIntegrator(double h, double tolerance) { + double absPos = Math.max(1.0e-12, tolerance); + double[] absTol = {absPos, absPos, absPos, 1.0e-10, 1.0e-10, 1.0e-10}; + double[] relTol = {1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12}; + DormandPrince54Integrator integrator = new DormandPrince54Integrator( + Math.max(minStepSize, 1.0e-12), Math.max(maxStepSize, minStepSize), absTol, relTol); + integrator.setInitialStepSize(initialStep(h)); + return integrator; + } + + private double initialStep(double h) { + return Math.max(minStepSize, Math.min(Math.abs(h), maxStepSize)); + } + + private static StepHandler recordSteps(final CLAS12Listener listener) { + return new StepHandler() { + @Override public void init(double s0, double[] y0, double sEnd) { } + @Override public void handleStep(StepInterpolator interpolator, boolean isLast) { + listener.accept(interpolator.getCurrentTime(), interpolator.getInterpolatedState().clone()); + } + }; + } + + /** + * ODE system for Cartesian CLAS12 swimming. + * Independent variable is path length s in cm. + */ + private static final class SwimEquations implements FirstOrderDifferentialEquations { + + private final FieldProbe probe; + private final double alpha; // 1/(kG*cm) + private final float[] b = new float[3]; + + SwimEquations(int q, double p, FieldProbe probe) { + this.probe = probe; + // Curvature scale for field in kG and distance in cm. + // alpha = 1.0e-14 * q * C / p (units: 1/(kG*cm)) + this.alpha = 1.0e-14 * q * CLAS12Swimmer.C / p; + } + + @Override + public int getDimension() { + return 6; + } + + @Override + public void computeDerivatives(double s, double[] y, double[] yDot) { + double Bx = 0.0, By = 0.0, Bz = 0.0; + + if (probe != null) { + probe.field((float) y[0], (float) y[1], (float) y[2], b); + Bx = b[0]; + By = b[1]; + Bz = b[2]; + } + + // dr/ds = t + yDot[0] = y[3]; + yDot[1] = y[4]; + yDot[2] = y[5]; + + // dt/ds = alpha * (t x B) + yDot[3] = alpha * (y[4] * Bz - y[5] * By); + yDot[4] = alpha * (y[5] * Bx - y[3] * Bz); + yDot[5] = alpha * (y[3] * By - y[4] * Bx); + } + } + + /** + * Sector-aware ODE system for sector-dependent swimming with a {@link RotatedCompositeProbe}. + *+ * The independent variable is the path length {@code s} in cm. + *
+ * + *+ * This implementation tries to call a sector-aware method on the probe via reflection: + * {@code field(int sector, float x, float y, float z, float[] b)}. + * If not found (or invocation fails), it falls back to {@code probe.field(x,y,z,b)}. + *
+ */ + private static final class SectorSwimEquations implements FirstOrderDifferentialEquations { + + private final int sector; + private final FieldProbe probe; + private final double alpha; // 1/(kG*cm) + private final float[] b = new float[3]; + + + // Cached reflective call (lazy init) + private transient java.lang.reflect.Method sectorFieldMethod; + private transient boolean searched = false; + + SectorSwimEquations(int sector, int q, double p, FieldProbe probe) { + this.sector = sector; + this.probe = probe; + this.alpha = 1.0e-14 * q * CLAS12Swimmer.C / p; + } + + @Override + public int getDimension() { + return 6; + } + + @Override + public void computeDerivatives(double s, double[] y, double[] yDot) { + + double Bx = 0.0, By = 0.0, Bz = 0.0; + + if (probe != null) { + if (!searched) { + searched = true; + sectorFieldMethod = findSectorFieldMethod(probe.getClass()); + } + + boolean ok = false; + + if (sectorFieldMethod != null) { + try { + // signature: (int, float, float, float, float[]) + sectorFieldMethod.invoke(probe, sector, (float) y[0], (float) y[1], (float) y[2], b); + ok = true; + } catch (Throwable t) { + // Disable and fall back for remainder of this swim + sectorFieldMethod = null; + } + } + + if (!ok) { + probe.field((float) y[0], (float) y[1], (float) y[2], b); + } + + + Bx = b[0]; + By = b[1]; + Bz = b[2]; + } + + // dr/ds = t + yDot[0] = y[3]; + yDot[1] = y[4]; + yDot[2] = y[5]; + + // dt/ds = alpha * (t x B) + yDot[3] = alpha * (y[4] * Bz - y[5] * By); + yDot[4] = alpha * (y[5] * Bx - y[3] * Bz); + yDot[5] = alpha * (y[3] * By - y[4] * Bx); + } + + private static java.lang.reflect.Method findSectorFieldMethod(Class cls) { + try { + return cls.getMethod("field", int.class, float.class, float.class, float.class, float[].class); + } catch (NoSuchMethodException e) { + return null; + } + } + } + + + private static final class HitFlag { + boolean hit = false; + } +} diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Trajectory.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Trajectory.java new file mode 100644 index 0000000000..56c7de172c --- /dev/null +++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Trajectory.java @@ -0,0 +1,191 @@ +package cnuphys.CLAS12Swim; + +import java.util.Arrays; + +import cnuphys.magfield.FieldProbe; +import cnuphys.magfield.RotatedCompositeProbe; +import cnuphys.swim.SwimTrajectory; + +@SuppressWarnings("serial") +public class CLAS12Trajectory extends SwimTrajectory { + + private double[] _s = new double[200]; + private int _sSize = 0; + + private double _bdlValue = Double.NaN; + + public CLAS12Trajectory(CLAS12Values initialValues) { + super(initialValues.toGeneratedParticleRecord(), 200); + } + + public void add(double s, double[] u) { + addS(s); + super.add(u); + _bdlValue = Double.NaN; + } + + public void replaceLastPoint(double s, double[] u) { + if (_sSize > 0) { + int index = _sSize - 1; + removePoint(index); + add(s, u); + } + } + + public void removeLastPoint() { + if (_sSize > 0) { + removePoint(_sSize - 1); + _bdlValue = Double.NaN; + } + } + + public void removePoint(int index) { + if (index >= 0 && index < _sSize) { + System.arraycopy(_s, index + 1, _s, index, _sSize - index - 1); + _sSize--; + remove(index); + _bdlValue = Double.NaN; + } + } + + public double getS(int index) { + return _s[index]; + } + + public int getSSize() { + return _sSize; + } + + public String sizeReport() { + return String.format("State vector size: %d Pathlength size: %d", size(), _sSize); + } + + @Override + public void clear() { + super.clear(); + _sSize = 0; + _bdlValue = Double.NaN; + } + + @Override + public boolean add(double[] u) { + throw new UnsupportedOperationException("Use add(s, u) instead."); + } + + @Override + public boolean add(double[] u, double s) { + throw new UnsupportedOperationException("Use add(s, u) instead."); + } + + @Override + public void add(double xo, double yo, double zo, double p, double theta, double phi) { + throw new UnsupportedOperationException("Use addPoint instead."); + } + + public void addPoint(double x, double y, double z, double theta, double phi, double s) { + double thetaRad = Math.toRadians(theta); + double phiRad = Math.toRadians(phi); + double sinTheta = Math.sin(thetaRad); + + double[] u = new double[6]; + u[0] = x; + u[1] = y; + u[2] = z; + u[3] = sinTheta * Math.cos(phiRad); + u[4] = sinTheta * Math.sin(phiRad); + u[5] = Math.cos(thetaRad); + + add(s, u); + } + + @Override + public double getR(int index) { + if ((index < 0) || (index >= size())) { + return Double.NaN; + } + + double[] v = get(index); + return Math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + } + + @Override + public double getComputedBDL() { + if (Double.isNaN(_bdlValue)) { + computeBDL(FieldProbe.factory()); + } + return _bdlValue; + } + + @Override + public void computeBDL(FieldProbe probe) { + if (!(probe instanceof RotatedCompositeProbe) && Double.isNaN(_bdlValue) && size() >= 2) { + _bdlValue = 0.0; + int n = size(); + double[] dr = new double[3]; + float[] b = new float[3]; + double[] bxdl = new double[3]; + + for (int i = 0; i < n - 1; i++) { + double[] p0 = get(i); + double[] p1 = get(i + 1); + + dr[0] = p1[0] - p0[0]; + dr[1] = p1[1] - p0[1]; + dr[2] = p1[2] - p0[2]; + + float xavg = (float) ((p0[0] + p1[0]) * 0.5); + float yavg = (float) ((p0[1] + p1[1]) * 0.5); + float zavg = (float) ((p0[2] + p1[2]) * 0.5); + + probe.field(xavg, yavg, zavg, b); + cross(b, dr, bxdl); + _bdlValue += vecmag(bxdl); + } + } + } + + @Override + public void sectorComputeBDL(int sector, RotatedCompositeProbe probe) { + if (Double.isNaN(_bdlValue) && size() >= 2) { + _bdlValue = 0.0; + int n = size(); + double[] dr = new double[3]; + float[] b = new float[3]; + double[] bxdl = new double[3]; + + for (int i = 0; i < n - 1; i++) { + double[] p0 = get(i); + double[] p1 = get(i + 1); + + dr[0] = p1[0] - p0[0]; + dr[1] = p1[1] - p0[1]; + dr[2] = p1[2] - p0[2]; + + float xavg = (float) ((p0[0] + p1[0]) * 0.5); + float yavg = (float) ((p0[1] + p1[1]) * 0.5); + float zavg = (float) ((p0[2] + p1[2]) * 0.5); + + probe.field(sector, xavg, yavg, zavg, b); + cross(b, dr, bxdl); + _bdlValue += vecmag(bxdl); + } + } + } + + private void addS(double s) { + if (_sSize >= _s.length) { + _s = Arrays.copyOf(_s, _s.length * 2); + } + _s[_sSize++] = s; + } + + private static void cross(float[] a, double[] b, double[] out) { + out[0] = a[1] * b[2] - a[2] * b[1]; + out[1] = a[2] * b[0] - a[0] * b[2]; + out[2] = a[0] * b[1] - a[1] * b[0]; + } + + private static double vecmag(double[] a) { + return Math.sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]); + } +} diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Values.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Values.java new file mode 100644 index 0000000000..a7713bef39 --- /dev/null +++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12Values.java @@ -0,0 +1,155 @@ +package cnuphys.CLAS12Swim; + +import cnuphys.lund.GeneratedParticleRecord; +import cnuphys.magfield.FastMath; + +/** + * A class to hold the initial or final values for a swim. + */ +public class CLAS12Values { + + /** The integer charge */ + public final int q; + + /** The coordinate x in cm */ + public final double x; + + /** The y coordinate in cm */ + public final double y; + + /** The z coordinate of in cm */ + public final double z; + + /** The momentum in GeV/c */ + public final double p; + + /** The DIRECTIONAL polar angle in degrees, i.e. theta component of p */ + public final double theta; + + /** The azimuthal angle in degrees */ + public final double phi; + + public final double tx; + public final double ty; + public final double tz; + + /** + * Store the initial conditions of a swim + * + * @param q The integer charge + * @param xo The x coordinate of the vertex in cm + * @param yo The y coordinate of the vertex in cm + * @param zo The z coordinate of the vertex in cm + * @param p The momentum in GeV/c + * @param theta The DITECTIONAL polar angle in degrees + * @param phi The DIRECTIONAL azimuthal angle in degrees + */ + public CLAS12Values(int q, double xo, double yo, double zo, double p, double theta, double phi) { + this.q = q; + this.x = xo; + this.y = yo; + this.z = zo; + this.p = p; + this.theta = theta; + this.phi = phi; + double thetaRad = Math.toRadians(theta); + double phiRad = Math.toRadians(phi); + double sinTheta = Math.sin(thetaRad); + tx = sinTheta * Math.cos(phiRad); + ty = sinTheta * Math.sin(phiRad); + tz = Math.cos(thetaRad); + } + + /** + * Get the POSITIONAL values from a state vector. The state vector is the vector + * of that is the dependent variable in the integration. The anlges are + * positional, not the directional angles for the momentum. + * + * @param q the integer charge. Must be supplied, not part of the state vector. + * It shouldn't change, but we assume this is the original momentum, so + * we mutliply by the state vector norm of the t components, which + * should be 1 since we have magnetic field only. + * @param p the momentum in GeV/c + * @param u the state vector + */ + public CLAS12Values(int q, double p, double[] u) { + this.q = q; + x = u[0]; + y = u[1]; + z = u[2]; + + tx = u[3]; + ty = u[4]; + tz = u[5]; + + // norm should be 1 + double norm = Math.sqrt(u[3] * u[3] + u[4] * u[4] + u[5] * u[5]); + + this.p = norm * p; + + // directional theta and phi + theta = FastMath.acos2Deg(u[5]); + phi = FastMath.atan2Deg(u[4], u[3]); + } + + /** + * Get the values as a state vector used in integration + * + * @return the values as a state vector + */ + public double[] getU() { + double uo[] = new double[6]; + + double thetaRad = Math.toRadians(theta); + double phiRad = Math.toRadians(phi); + double sinTheta = Math.sin(thetaRad); + + double tx = sinTheta * Math.cos(phiRad); // px/p + double ty = sinTheta * Math.sin(phiRad); // py/p + double tz = Math.cos(thetaRad); // pz/p + + // set uf to the starting state vector + uo[0] = x; + uo[1] = y; + uo[2] = z; + uo[3] = tx; + uo[4] = ty; + uo[5] = tz; + return uo; + } + + /** + * Copy constructor + * + * @param src the source initial values + */ + public CLAS12Values(CLAS12Values src) { + this(src.q, src.x, src.y, src.z, src.p, src.theta, src.phi); + } + + @Override + public String toString() { + return String.format("Q: %d\n", q) + String.format("xo: %10.7e cm\n", x) + String.format("yo: %10.7e cm\n", y) + + String.format("zo: %10.7e cm\n", z) + String.format("p: %10.7e GeV/c\n", p) + + String.format("theta: %10.7f deg\n", theta) + String.format("phi: %10.7f deg", phi); + } + + /** + * Convert to a GeneratedParticleRecord for backwards compatibility + * + * @return a GeneratedParticleRecord corresponding to this data + */ + public GeneratedParticleRecord toGeneratedParticleRecord() { + return new GeneratedParticleRecord(q, x, y, z, p, theta, phi); + } + + /** + * A raw string for output, just numbers no units + * + * @return a raw string for output + */ + public String toStringRaw() { + return String.format("%-7.4f %-7.4f %-7.4f %-6.3f %-6.3f %-6.3f", x, y, z, p, theta, phi); + } + +} diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12ZLineListener.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12ZLineListener.java new file mode 100644 index 0000000000..914183e653 --- /dev/null +++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12ZLineListener.java @@ -0,0 +1,33 @@ +package cnuphys.CLAS12Swim; + +public class CLAS12ZLineListener extends CLAS12DOCAListener { + + private double _xb; //x offset in cm + private double _yb; //y offset in cm + + /** + * Create a CLAS12 swim to an "offest beamline" listener. The offset + * beamline is a line parallel to the z-axis, but offset in the x and y + * directions by _xb and _yb. + * + * @param ivals the initial values of the swim + * @param xb the x offset (cm) + * @param yb the y offset (cm) + * @param accuracy the accuracy (cm) (on on difference in successive docas) + * @param sMax the final or max path length (cm) + */ + public CLAS12ZLineListener(CLAS12Values ivals, double xb, double yb, double accuracy, double sMax) { + super(ivals, accuracy, sMax); + _xb = xb; + _yb = yb; + } + + + @Override + public double doca(double newS, double[] newU) { + double dx = newU[0] - _xb; + double dy = newU[1] - _yb; + return Math.hypot(dx, dy); + } + +} diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12ZListener.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12ZListener.java new file mode 100644 index 0000000000..5565ee5c6a --- /dev/null +++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/CLAS12ZListener.java @@ -0,0 +1,118 @@ +package cnuphys.CLAS12Swim; + +/** + * A listener for swimming to a fixed value of z + */ +public class CLAS12ZListener extends CLAS12BoundaryListener { + + // the target z (cm) + private double _zTarget; + + // the starting sign. When this changes we have crossed. + private double _startSign; + + /** + * Create a CLAS12 boundary target Z listener, for swimming to a fixed z + * + * @param ivals the initial values of the swim + * @param zTarget the target z (cm) + * @param accuracy the desired accuracy (cm) + * @param sMax the final or max path length (cm) + */ + public CLAS12ZListener(CLAS12Values ivals, double zTarget, double accuracy, double sMax) { + super(ivals, accuracy, sMax); + _zTarget = zTarget; + _startSign = sign(ivals.z); + } + + @Override + public boolean crossedBoundary(double newS, double[] newU) { + double newZ = newU[2]; + int sign = sign(newZ); + + if (sign != _startSign) { + return true; + } + return false; + } + + @Override + public boolean accuracyReached(double newS, double[] newU) { + double dZ = Math.abs(newU[2] - _zTarget); + return dZ < _accuracy; + } + + // left or right of the target Z? + private int sign(double z) { + return (z < _zTarget) ? -1 : 1; + } + + /** + * Get the absolute distance to the target (boundary) in cm. + * + * @param newS the new path length + * @param newU the new state vector + * @return the distance to the target (boundary) in cm. + */ + @Override + public double distanceToTarget(double newS, double[] newU) { + return Math.abs(newU[2] - _zTarget); + } + + /** + * Add a second point creating a straight line to the target z + */ + @Override + public void straightLine() { + + double u[] = _trajectory.get(_trajectory.size() - 1); + double s = _trajectory.getS(_trajectory.size() - 1); + + double u2[] = findPoint(u[0], u[1], u[2], u[3], u[4], u[5], _zTarget); + + double dx = u2[0] - u[0]; + double dy = u2[1] - u[1]; + double dz = u2[2] - u[2]; + double ds = Math.sqrt(dx * dx + dy * dy + dz * dz); + + _trajectory.add(s + ds, u2); + _status = CLAS12Swimmer.SWIM_SUCCESS; + + } + + /** + * Finds the point along the line of velocity where the z coordinate reaches + * zTarget. + * + * @param x0 Starting x coordinate + * @param y0 Starting y coordinate + * @param z0 Starting z coordinate + * @param tx x component of the unit direction vector + * @param ty y component of the unit direction vector + * @param tz z component of the unit direction vector + * @param zTarget The target z coordinate to reach + * @return The point [x, y, z] where the z coordinate reaches zTarget, or null + * if it never reaches. + */ + private double[] findPoint(double x0, double y0, double z0, double tx, double ty, double tz, double zTarget) { + // Check if the line is parallel to the z-plane (tz = 0) + if (tz == 0) { + if (z0 == zTarget) { + // The entire line is on the plane where z = zTarget + return new double[] { x0, y0, zTarget }; + } else { + // The line will never reach zTarget + return null; + } + } + + // Calculate the parameter (s) at which z coordinate reaches zTarget + double t = (zTarget - z0) / tz; + + // Calculate the x and y coordinates at this point + double x = x0 + tx * t; + double y = y0 + ty * t; + + return new double[] { x, y, zTarget, tx, ty, tz }; + } +} diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/ICLAS12Swimmer.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/ICLAS12Swimmer.java new file mode 100644 index 0000000000..fb7593ad5b --- /dev/null +++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/ICLAS12Swimmer.java @@ -0,0 +1,369 @@ +package cnuphys.CLAS12Swim; + +import cnuphys.CLAS12Swim.geometry.Cylinder; +import cnuphys.CLAS12Swim.geometry.Plane; +import cnuphys.CLAS12Swim.geometry.Sphere; +import cnuphys.magfield.FieldProbe; + +/** + * Public API for the CLAS12 charged-particle swimmer. + *+ * A “swim” numerically propagates a charged particle through the magnetic field from an + * initial vertex and direction/momentum until a termination condition is reached + * (path length, target surface, target z, target rho/beamline, etc.). + *
+ * + *+ *
+ * If the swim starts inside the cylinder, it will terminate immediately (subject to the + * implementation’s handling of that case in {@link CLAS12SwimResult}). + *
+ * + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param p1 first point on the cylinder centerline: {x,y,z} in cm + * @param p2 second point on the cylinder centerline: {x,y,z} in cm + * @param r cylinder radius in cm + * @param accuracy desired accuracy in cm for reaching the surface + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult swimCylinder(int q, double xo, double yo, double zo, double p, double theta, double phi, + double p1[], double p2[], double r, double accuracy, double sMax, double h, + double tolerance); + + /** + * Swim a particle to the surface of a target cylinder. + * The cylinder is specified by a {@link Cylinder} object (typically treated as infinite in length). + * + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param targetCylinder target cylinder + * @param accuracy desired accuracy in cm for reaching the surface + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult swimCylinder(int q, double xo, double yo, double zo, double p, double theta, double phi, + Cylinder targetCylinder, double accuracy, double sMax, double h, + double tolerance); + + /** + * Swim a particle to the surface of a target sphere. + * The sphere is defined by a center point {@code center} as a {@code double[3]} in {x,y,z} cm, + * and radius {@code r} in cm. + *+ * If the swim starts inside the sphere, it will terminate immediately (subject to the + * implementation’s handling of that case in {@link CLAS12SwimResult}). + *
+ * + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param center sphere center: {x,y,z} in cm + * @param r sphere radius in cm + * @param accuracy desired accuracy in cm for reaching the surface + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult swimSphere(int q, double xo, double yo, double zo, double p, double theta, double phi, + double center[], double r, double accuracy, double sMax, double h, + double tolerance); + + /** + * Swim a particle to the surface of a target sphere. + * The sphere is specified by a {@link Sphere} object. + * + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param targetSphere target sphere + * @param accuracy desired accuracy in cm for reaching the surface + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult swimSphere(int q, double xo, double yo, double zo, double p, double theta, double phi, + Sphere targetSphere, double accuracy, double sMax, double h, + double tolerance); + + /** + * Swim a particle until it intersects a target plane or until {@code sMax} is reached. + * The plane is defined by the components of a normal vector and the components of a point on the plane. + * + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param nx plane normal x-component + * @param ny plane normal y-component + * @param nz plane normal z-component + * @param px x-component of a point on the plane (cm) + * @param py y-component of a point on the plane (cm) + * @param pz z-component of a point on the plane (cm) + * @param accuracy desired accuracy in cm for reaching the plane + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult swimPlane(int q, double xo, double yo, double zo, double p, double theta, double phi, + double nx, double ny, double nz, double px, double py, double pz, + double accuracy, double sMax, double h, double tolerance); + + /** + * Swim a particle until it intersects a target plane or until {@code sMax} is reached. + * The plane is defined by a normal vector and a point on the plane. + * + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param norm plane normal vector {nx, ny, nz} + * @param point a point on the plane {px, py, pz} in cm + * @param accuracy desired accuracy in cm for reaching the plane + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult swimPlane(int q, double xo, double yo, double zo, double p, double theta, double phi, + double norm[], double point[], double accuracy, double sMax, double h, + double tolerance); + + /** + * Swim a particle until it intersects a target plane or until {@code sMax} is reached. + * The plane is specified by a {@link Plane} object. + * + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param targetPlane target plane + * @param accuracy desired accuracy in cm for reaching the plane + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult swimPlane(int q, double xo, double yo, double zo, double p, double theta, double phi, + Plane targetPlane, double accuracy, double sMax, double h, + double tolerance); + + /** + * Swim to a target {@code z} (cm) in a sector coordinate system. + *+ * Important: this is only valid if the underlying field/probe is a rotated composite + * field implementation (your {@code CLAS12Swimmer} uses {@link cnuphys.magfield.RotatedCompositeProbe} + * internally for sector coordinate transforms). + *
+ * The swim is terminated when the particle reaches {@code zTarget} or if {@code sMax} is reached. + * + * @param sector sector number in [1..6] + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param zTarget target z position in cm + * @param accuracy desired accuracy in cm + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult sectorSwimZ(int sector, int q, double xo, double yo, double zo, double p, double theta, + double phi, double zTarget, double accuracy, double sMax, double h, + double tolerance); + + /** + * Swim to a target {@code z} (cm). + * The swim is terminated when the particle reaches {@code zTarget} or if {@code sMax} is reached. + * + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param zTarget target z position in cm + * @param accuracy desired accuracy in cm + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult swimZ(int q, double xo, double yo, double zo, double p, double theta, double phi, + double zTarget, double accuracy, double sMax, double h, double tolerance); + + /** + * Swim to a target cylindrical radius {@code rho} (cm), i.e. to the surface of an infinite cylinder + * about the z-axis. + * The swim is terminated when the particle reaches {@code rhoTarget} or if {@code sMax} is reached. + * + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param rhoTarget target rho (radius) in cm + * @param accuracy desired accuracy in cm + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult swimRho(int q, double xo, double yo, double zo, double p, double theta, double phi, + double rhoTarget, double accuracy, double sMax, double h, double tolerance); + + /** + * Swim to an "offset beamline" listener. + * The offset beamline is a line parallel to the z-axis, offset in x and y by {@code xb} and {@code yb}. + * The goal is to swim to the distance of closest approach (DOCA) to this offset line. + * Swim terminates when successive DOCA estimates differ by less than {@code accuracy}. + * + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param xb beamline x offset in cm + * @param yb beamline y offset in cm + * @param accuracy desired DOCA convergence accuracy in cm + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult swimZLine(int q, double xo, double yo, double zo, double p, double theta, double phi, + double xb, double yb, double accuracy, double sMax, double h, double tolerance); + + /** + * Swim to the beamline (defined by {@code rho = 0}), i.e. find the distance of closest approach (DOCA) + * to the z-axis. + * Swim terminates when successive DOCA estimates differ by less than {@code accuracy}. + * + * @param q particle charge in integer units of e + * @param xo initial x vertex position in cm + * @param yo initial y vertex position in cm + * @param zo initial z vertex position in cm + * @param p initial momentum magnitude in GeV/c + * @param theta initial polar angle in degrees + * @param phi initial azimuthal angle in degrees + * @param accuracy desired DOCA convergence accuracy in cm + * @param sMax maximum path length in cm + * @param h initial step size in cm + * @param tolerance desired tolerance; the integrator adapts step size to meet this tolerance + * @return the result of the swim + */ + CLAS12SwimResult swimBeamline(int q, double xo, double yo, double zo, double p, double theta, double phi, + double accuracy, double sMax, double h, double tolerance); +} diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/ODEStepListener.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/ODEStepListener.java new file mode 100644 index 0000000000..6cbc0d0c0d --- /dev/null +++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/ODEStepListener.java @@ -0,0 +1,16 @@ +package cnuphys.CLAS12Swim; + +/** + * Interface for listening to steps taken by an ODE solver. + */ +public interface ODEStepListener { + /** + * Called when a new step is taken in the ODE solving process. + * + * @param newT The new independent variable after the step. + * @param newY The new state vector after the step. + * @return A boolean indicating whether to continue (true) or stop (false) the + * integration. + */ + boolean newStep(double newT, double[] newY); +} diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Cylinder.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Cylinder.java new file mode 100644 index 0000000000..0f89c9d215 --- /dev/null +++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Cylinder.java @@ -0,0 +1,109 @@ +package cnuphys.CLAS12Swim.geometry; + + +/** + * An INFINITE cylinder is defined by a centerline and a radius + * @author heddle + * + */ +public class Cylinder { + + //the centerline + private Line _centerLine; + + //the radius + public double radius; + + /** + * Create a cylinder + * @param centerLine the center line + * @param radius the radius + */ + public Cylinder(Line centerLine, double radius) { + _centerLine = new Line(centerLine); + this.radius = radius; + } + + /** + * Create a cylinder + * @param p1 one point of center line as an xyz array + * @param p2 another point of center line as an xyz array + * @param radius + */ + public Cylinder(double[] p1, double[] p2, double radius) { + this(new Line(p1, p2), radius); + } + + + /** + * Get the shortest distance between the surface of this infinite cylinder and a point. + * If the value is negative, we are inside the cylinder. + * @param p a point + * @return the perpendicular distance + */ + public double signedDistance(Point p) { + double lineDist = _centerLine.distance(p); + return lineDist - radius; + } + + /** + * Set the path length of the swim + * @deprecated Use {@link Cylinder#signedDistance} instead. + * @param p a point + * @return the perpendicular distance + */ + @Deprecated + public double distance(Point p) { + double lineDist = _centerLine.distance(p); + return lineDist - radius; + } + + /** + * Get the shortest distance between the surface of this infinite cylinder and a point. + * If the value is negative, we are inside the cylinder. + * @param x the x coordinate + * @param y the y coordinate + * @param z the z coordinate + * @return the perpendicular distance + */ + public double signedDistance(double x, double y, double z) { + Point p = new Point(x, y, z); + return signedDistance(p); + } + + /** + * Get the shortest absolute distance between the surface of this infinite cylinder and a point. + * @param x the x coordinate + * @param y the y coordinate + * @param z the z coordinate + * @return the perpendicular distance + */ + public double distance(double x, double y, double z) { + Point p = new Point(x, y, z); + return Math.abs(signedDistance(p)); + } + + /** + * Is the point inside the cylinder? + * @param x the x coordinate + * @param y the y coordinate + * @param z the z coordinate + * @returntrue if the point is inside the cylinder.
+ */
+ public boolean isInside(double x, double y, double z) {
+ return signedDistance(x, y, z) < 0;
+ }
+
+ /**
+ * Is the cylinder centered on the z axis?
+ * @return true if the cylinder is centered on the z axis.
+ */
+ public boolean centeredOnZ() {
+ double x0 = _centerLine.getP0().x;
+ double y0 = _centerLine.getP0().y;
+ double x1 = _centerLine.getP1().x;
+ double y1 = _centerLine.getP1().y;
+ return (x0 == 0) && (y0 == 0) && (x1 == 0) && (y1 == 0);
+ }
+
+}
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Line.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Line.java
new file mode 100644
index 0000000000..ebb96780c8
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Line.java
@@ -0,0 +1,157 @@
+package cnuphys.CLAS12Swim.geometry;
+
+/**
+ * 3D line of the form p(t) = po + t*dp where p(t) is a point on the line, po is
+ * one point and dp = (p1-po) where p1 is another point. If this is an infinite
+ * line, the t = [-infinity, infinity]. If this is a directed segment, t = [0,
+ * 1]
+ *
+ * @author heddle
+ *
+ */
+public class Line {
+
+ private Point _po; // "start" point on the line
+ private Vector _dp; // segment from start to end
+ private double _dpLen; // length of start to end segment
+
+ /**
+ * Create a line from two points on the line. If this is a directed line
+ * segment, the line will go from po to p1.
+ *
+ * @param po one point
+ * @param p1 the other point
+ */
+ public Line(Point po, Point p1) {
+ _po = new Point(po);
+ _dp = new Vector(Point.difference(p1, po));
+ _dpLen = _dp.length();
+ }
+
+ /**
+ * Create a line from two ponts passed as arrays
+ * @param p1 one point as an xyz array
+ * @param p2 another point as an xyz array
+ */
+ public Line(double[] p1, double[] p2) {
+ this(new Point(p1), new Point(p2));
+ }
+
+ /**
+ * Copy constructor
+ * @param line the line to copy
+ */
+ public Line(Line line) {
+ _po = new Point(line._po);
+ _dp = new Vector(line._dp);
+ _dpLen = line._dpLen;
+ }
+
+ /**
+ * Create a line from the origing in the direction of a vector
+ * @param v the vector
+ */
+ public Line(Vector v) {
+ this(new Point(0,0,0), new Point(v.x, v.y, v.z));
+ }
+
+ /**
+ * Get the po "start" point. This is just an arbitrary point on an infinite
+ * line, but the starting point if this is a directed line segment
+ *
+ * @return the "starting" point.
+ */
+ public Point getP0() {
+ return _po;
+ }
+
+ /**
+ * Get the p1-po "dP" segment
+ *
+ * @return dP = p1 - po
+ */
+ public Vector getDelP() {
+ return _dp;
+ }
+
+ /**
+ * Get the p1 "end" point. This is just an arbitrary point on an infinite line,
+ * but the end point if this is a directed line segment
+ *
+ * @return the "end" point.
+ */
+ public Point getP1() {
+ return new Point(_po.x + _dp.x, _po.y + _dp.y, _po.z + _dp.z);
+ }
+
+ /**
+ * Get a point on the line
+ *
+ * @param t the t parameter. If this is a directed line segment, t should be
+ * restricted to [0, 1]
+ * @return a point on the line
+ */
+ public Point getP(double t) {
+ Point p = new Point();
+ getP(t, p);
+ return p;
+ }
+
+ /**
+ * Get a point on the line (in place)
+ *
+ * @param t the t parameter. If this is a directed line segment, t should be
+ * restricted to [0, 1]
+ * @param p upon return, a point on the line
+ */
+ public void getP(double t, Point p) {
+ p.x = _po.x + t * _dp.x;
+ p.y = _po.y + t * _dp.y;
+ p.z = _po.z + t * _dp.z;
+ }
+
+ /**
+ * Get the shortest distance between this line (as an infinite line) and a point
+ *
+ * @param p a point
+ * @return the perpendicular distance
+ */
+ public double distance(Point p) {
+ Vector ap = new Vector(Point.difference(p, _po));
+ Vector c = Vector.cross(ap, _dp);
+ return c.length() / _dpLen;
+ }
+
+ /**
+ * Find the point on the line closest to the given point
+ * @param p the given point
+ * @return the point on the line closest to the given point
+ */
+ public Point closestPointOnLine(Point p) {
+ Point pointVec = p.subtract(_po);
+ double t = pointVec.dot(_dp) / _dp.dot(_dp);
+ return _po.add(_dp.scale(t));
+ }
+
+
+
+ /**
+ * Get a String representation
+ *
+ * @return a String representation of the Line
+ */
+ @Override
+ public String toString() {
+ return "Line from " + getP0() + " to " + getP1();
+ }
+
+ /**
+ * Get the center of the line
+ *
+ * @return the center of the line
+ */
+ public Point getCenter() {
+ return getP(0.5);
+ }
+
+}
\ No newline at end of file
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Plane.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Plane.java
new file mode 100644
index 0000000000..623ec5873b
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Plane.java
@@ -0,0 +1,325 @@
+package cnuphys.CLAS12Swim.geometry;
+
+/**
+ * A plane is defined by the equation (r - ro).norm = 0 Where r is an arbitrary
+ * point on the plane, ro is a given point on the plane and norm is the normal
+ * to the plane
+ *
+ * @author heddle
+ *
+ */
+public class Plane {
+
+ /** Effectively zero */
+ private static final double TINY = 1.0e-20;
+
+ // for the form ax + by + cz = d;
+ public final double a;
+ public final double b;
+ public final double c;
+ public final double d;
+
+ private double _denom = Double.NaN;
+
+ /**
+ * Create a plane from a normal vector and a point on the plane
+ *
+ * @param norm the normal vector
+ * @param p0 a point in the plane
+ * @return the plane that contains p and its normal is norm
+ */
+ public Plane(Vector anorm, Point p0) {
+ // lets make it a unit vector
+ Vector norm = anorm.unitVector();
+ a = norm.x; // A
+ b = norm.y; // B
+ c = norm.z; // C
+ d = a * p0.x + b * p0.y + c * p0.z; // D
+ }
+
+ /**
+ * Create a plane from the coefficients of the equation ax + by + cz = d
+ *
+ * @param a the a coefficient
+ * @param b the b coefficient
+ * @param c the c coefficient
+ * @param d the d coefficient
+ */
+ public Plane(double a, double b, double c, double d) {
+ this.a = a;
+ this.b = b;
+ this.c = c;
+ this.d = d;
+ }
+
+ /**
+ * Create a plane from the normal vector in an array of doubles and a point in
+ * the plane in an array, both (x, y, z)
+ *
+ * @param norm the normal
+ * @param point the point in the plane
+ */
+ public Plane(double norm[], double point[]) {
+
+ this(new Vector(norm[0], norm[1], norm[2]), new Point(point[0], point[1], point[2]));
+ }
+
+ /**
+ * Create a plane from a normal vector and a point on the plane
+ *
+ * @param nx x component of normal vector
+ * @param ny y component of normal vector
+ * @param nz z component of normal vector
+ * @param px x component of point on plane
+ * @param py y component of point on plane
+ * @param pz z component of point on plane
+ */
+ public Plane(double nx, double ny, double nz, double px, double py, double pz) {
+
+ this(new Vector(nx, ny, nz), new Point(px, py, pz));
+ }
+
+ /**
+ * Create a line from two points and then get the intersection with the plane
+ *
+ * @param p1 one point
+ * @param p2 another point
+ * @param p will hold the intersection, NaNs if no intersection
+ * @return the t parameter. If NaN it means the line is parallel to the plane.
+ * If t [0,1] then the segment intersects the plane. If t outside [0, 1]
+ * the infinite line intersects the plane, but not the segment
+ */
+ public double interpolate(Point p1, Point p2, Point p) {
+ Line line = new Line(p1, p2);
+ return lineIntersection(line, p);
+ }
+
+ /**
+ * Distance from a point to the plane
+ *
+ * @param x the x coordinate
+ * @param y the y coordinate
+ * @param z the z coordinate
+ * @return the signed distance (indicates which side you are on where norm
+ * defines positive side)
+ */
+ public double distance(double x, double y, double z) {
+ return Math.abs(signedDistance(x, y, z));
+ }
+
+ /**
+ * Signed distance from a point to the plane
+ *
+ * @param x the x coordinate
+ * @param y the y coordinate
+ * @param z the z coordinate
+ * @return the signed distance (indicates which side you are on where norm
+ * defines positive side)
+ */
+ public double signedDistance(double x, double y, double z) {
+ if (Double.isNaN(_denom)) {
+ _denom = Math.sqrt(a * a + b * b + c * c);
+ }
+ return (a * x + b * y + c * z - d) / _denom;
+ }
+
+ /**
+ * Compute the intersection of an infinite line with the plane
+ *
+ * @param line the line
+ * @param intersection will hold the point of intersection
+ * @return the t parameter. If NaN it means the line is parallel to the plane.
+ * If t [0,1] then the segment intersects the plane. If t outside [0, 1]
+ * the infinite line intersects the plane, but not the segment
+ */
+ public double lineIntersection(Line line, Point intersection) {
+ // Direction vector of the line
+ Vector lineDir = line.getDelP();
+
+ Point p0 = line.getP0();
+
+ // Check if the line is parallel to the plane
+ double dotProduct = a * lineDir.x + b * lineDir.y + c * lineDir.z;
+ if (Math.abs(dotProduct) < TINY) {
+ System.err.println("The line is parallel to the plane in Plane.findLinePlaneIntersection.");
+ return Double.NaN;
+ }
+
+ // Parameter t in the parametric line equation
+ double t = (d - a * p0.x - b * p0.y - c * p0.z) / dotProduct;
+
+ line.getP(t, intersection);
+ return t;
+ }
+
+ /**
+ * Get whether the point is to the left, right or (exactly) on the plane
+ *
+ * @param x the x coordinate
+ * @param y the y coordinate
+ * @param z the z coordinate
+ * @return +1 if to the left, -1 if to the right, 0 if on the plane
+ */
+ public int sign(double x, double y, double z) {
+ double result = a * x + b * y + c * z;
+
+ if (result > d) {
+ return +1;
+ } else if (result < d) {
+ return -1;
+ } else {
+ return 0;
+ }
+ }
+
+ /**
+ * Create a plane of constant azimuthal angle phi
+ *
+ * @param phi the azimuthal angle in degrees
+ * @return the plane of constant phi
+ */
+ public static Plane constantPhiPlane(double phi) {
+ phi = Math.toRadians(phi);
+
+ double cphi = Math.cos(phi);
+ double sphi = Math.sin(phi);
+
+ // point in the plane
+ Point p = new Point(cphi, sphi, 0);
+
+ // normal
+ Vector norm = new Vector(sphi, -cphi, 0);
+
+ return new Plane(norm, p);
+ }
+
+ @Override
+ public String toString() {
+ String pstr = String.format("abcd = [%10.6G, %10.6G, %10.6G, %10.6G]", a, b, c, d);
+ return pstr;
+ }
+
+ // is the value essentially 0?
+ private boolean tiny(double v) {
+ return Math.abs(v) < TINY;
+ }
+
+ /**
+ * Find some coordinates suitable for drawing the plane as a Quad in 3D
+ *
+ * @param scale an arbitrary big number, a couple times bigger than the drawing
+ * extent
+ * @return the jogl coordinates for drawing a Quad
+ */
+ public float[] planeQuadCoordinates(float scale) {
+
+ int[] i1 = { -1, -1, 1, 1 };
+ int[] i2 = { -1, 1, 1, -1 };
+
+ if (tiny(a) && tiny(b) && tiny(c)) {
+ return null;
+ }
+
+ float[] coords = new float[12];
+
+ if (tiny(b) && tiny(c)) { // constant x plane
+ float fx = (float) (d / a);
+ for (int k = 0; k < 4; k++) {
+ int j = 3 * k;
+
+ float y = scale * i1[k];
+ float z = scale * i2[k];
+
+ coords[j] = fx;
+ coords[j + 1] = y;
+ coords[j + 2] = z;
+ }
+
+ } else if (tiny(a) && tiny(c)) { // constant y plane
+ float fy = (float) (d / b);
+ for (int k = 0; k < 4; k++) {
+ int j = 3 * k;
+
+ float x = scale * i1[k];
+ float z = scale * i2[k];
+
+ coords[j] = x;
+ coords[j + 1] = fy;
+ coords[j + 2] = z;
+ }
+ } else if (tiny(a) && tiny(b)) { // constant z plane
+ float fz = (float) (d / c);
+ for (int k = 0; k < 4; k++) {
+ int j = 3 * k;
+
+ float x = scale * i1[k];
+ float y = scale * i2[k];
+
+ coords[j] = x;
+ coords[j + 1] = y;
+ coords[j + 2] = fz;
+ }
+ }
+
+ else if (tiny(a)) {
+ for (int k = 0; k < 4; k++) {
+ int j = 3 * k;
+
+ float x = scale * i1[k];
+ float y = scale * i2[k];
+ float z = (float) ((d - b * y) / c);
+
+ coords[j] = x;
+ coords[j + 1] = y;
+ coords[j + 2] = z;
+ }
+ }
+
+ else if (tiny(b)) {
+ for (int k = 0; k < 4; k++) {
+ int j = 3 * k;
+
+ float x = scale * i1[k];
+ float y = scale * i2[k];
+ float z = (float) ((d - a * x) / c);
+
+ coords[j] = x;
+ coords[j + 1] = y;
+ coords[j + 2] = z;
+ }
+
+ }
+
+ else if (tiny(c)) {
+ for (int k = 0; k < 4; k++) {
+ int j = 3 * k;
+
+ float x = scale * i1[k];
+ float z = scale * i2[k];
+ float y = (float) ((d - a * x) / b);
+
+ coords[j] = x;
+ coords[j + 1] = y;
+ coords[j + 2] = z;
+ }
+
+ }
+
+ else { // general case, no small constants
+ for (int k = 0; k < 4; k++) {
+ int j = 3 * k;
+
+ float x = scale * i1[k];
+ float y = scale * i2[k];
+ float z = (float) ((d - a * x - b * y) / c);
+
+ coords[j] = x;
+ coords[j + 1] = y;
+ coords[j + 2] = z;
+ }
+ }
+
+ return coords;
+ }
+
+}
\ No newline at end of file
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Point.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Point.java
new file mode 100644
index 0000000000..00bad9c3f7
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Point.java
@@ -0,0 +1,186 @@
+package cnuphys.CLAS12Swim.geometry;
+
+public class Point {
+
+ /** x component */
+ public double x;
+ /** y component */
+ public double y;
+ /** z component */
+ public double z;
+
+ /**
+ * Create a point at the origin
+ */
+ public Point() {
+ this(0, 0, 0);
+ }
+
+ /**
+ * Copy constructor
+ *
+ * @param p the point to copy
+ */
+ public Point(Point p) {
+ this(p.x, p.y, p.z);
+ }
+
+ /**
+ * Create a point
+ *
+ * @param x the x coordinate
+ * @param y the y coordinate
+ * @param z the z coordinate
+ */
+ public Point(double x, double y, double z) {
+ this.x = x;
+ this.y = y;
+ this.z = z;
+ }
+
+ /**
+ * Create a point from an xyz array
+ * @param p the point as an xyz array
+ */
+ public Point(double[] p) {
+ this(p[0], p[1], p[2]);
+ }
+
+ /**
+ * Set the components of the point (vector)
+ *
+ * @param x the x component
+ * @param y the y component
+ * @param z the z component
+ */
+ public void set(double x, double y, double z) {
+ this.x = x;
+ this.y = y;
+ this.z = z;
+ }
+
+ /**
+ * Set the components of the point (vector)
+ * @param p the point to use to set (e.g., copy)
+ */
+ public void set(Point p) {
+ this.x = p.x;
+ this.y = p.y;
+ this.z = p.z;
+ }
+
+ /**
+ * Get the difference between two points
+ *
+ * @param a one point
+ * @param b the other point
+ * @return the difference between two points a - b
+ */
+ public static Point difference(Point a, Point b) {
+ return new Point(a.x - b.x, a.y - b.y, a.z - b.z);
+ }
+
+ /**
+ * Method to subtract another point from this point
+ *
+ * @param other the point to subtract
+ * @return the difference between this point and the other point
+ */
+ public Point subtract(Point other) {
+ return new Point(x - other.x, y - other.y, z - other.z);
+ }
+
+ /**
+ * Get the in-place difference between two points
+ *
+ * @param a one point
+ * @param b the other point
+ * @param c upon return the difference between two points a - b
+ */
+ public static void difference(Point a, Point b, Point c) {
+ c.set(a.x - b.x, a.y - b.y, a.z - b.z);
+ }
+
+ /**
+ * The dot product of this "vector" with another vector
+ *
+ * @param v the other vector or point
+ * @return the dot product
+ */
+ public double dot(Point v) {
+ return x * v.x + y * v.y + z * v.z;
+ }
+
+ /**
+ * The dot product of two vectors or points
+ *
+ * @param a one vector or point
+ * @param b the other vector or point
+ * @return the dot product
+ */
+ public static double dot(Point a, Point b) {
+ return a.dot(b);
+ }
+
+ /**
+ * Get a string representation of the Point
+ *
+ * @return a String representation
+ */
+ @Override
+ public String toString() {
+ return String.format("(%-10.6f, %-10.6f, %-10.6f)", x, y, z);
+ }
+
+ /**
+ * Compute the distance to another point
+ * @param x the x coordinate of the other point
+ * @param y the y coordinate of the other point
+ * @param z the z coordinate of the other point
+ * @return the distance between the points
+ */
+ public double distance(double x, double y, double z) {
+ double dx = x - this.x;
+ double dy = y - this.y;
+ double dz = z - this.z;
+ return Math.sqrt(dx*dx + dy*dy + dz*dz);
+ }
+
+ /**
+ * Compute the distance to another point
+ * @param p the other point
+ * @return the distance between the points
+ */
+ public double distance(Point p) {
+ return distance(p.x, p.y, p.z);
+ }
+
+ /** Method to calculate distance between two points
+ *
+ * @param p1 one point
+ * @param p2 the other point
+ * @return the distance between the two points
+ */
+ public static double distance(Point p1, Point p2) {
+ return p1.distance(p2);
+ }
+
+ /**
+ * Add another point to this point (i.e., vector addition)
+ * @param other the other point
+ * @return the sum of the two points
+ */
+ public Point add(Point other) {
+ return new Point(x + other.x, y + other.y, z + other.z);
+ }
+
+ /**
+ * Scale this point by a scalar
+ * @param scalar the scalar multiplier
+ * @return the scaled point
+ */
+ public Point scale(double scalar) {
+ return new Point(x * scalar, y * scalar, z * scalar);
+ }
+
+}
\ No newline at end of file
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Sphere.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Sphere.java
new file mode 100644
index 0000000000..c7ecff95ab
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Sphere.java
@@ -0,0 +1,155 @@
+package cnuphys.CLAS12Swim.geometry;
+
+/**
+ * A sphere centered at an arbitrary point
+ *
+ * @author heddle
+ *
+ */
+public class Sphere {
+
+ // the center
+ private Point _center;
+
+ // the radius
+ private double _radius;
+
+ /**
+ * Create a sphere
+ *
+ * @param center the center of the sphere
+ * @param radius the radius of the sphere
+ */
+ public Sphere(Point center, double radius) {
+ _center = new Point(center);
+ _radius = radius;
+ }
+
+ /**
+ * Create a sphere
+ *
+ * @param center the center of the sphere as an xyz array
+ * @param radius the radius of the sphere
+ */
+ public Sphere(double[] center, double radius) {
+ this(new Point(center[0], center[1], center[2]), radius);
+ }
+
+ /**
+ * Create a sphere centered on the origin
+ *
+ * @param radius the radius of the sphere
+ */
+ public Sphere(double radius) {
+ this(new Point(0, 0, 0), radius);
+ }
+
+ /**
+ * Get the radius of the sphere
+ *
+ * @return the radius of the sphere
+ */
+ public double getRadius() {
+ return _radius;
+ }
+
+ /**
+ * Get the shortest distance between the surface of this sphere and a point. If
+ * the value is negative, we are inside the sphere.
+ *
+ * @param p a point
+ * @return the distance to the sphere
+ */
+ public double signedDistance(Point p) {
+ double centDist = _center.distance(p);
+ return centDist - _radius;
+ }
+
+ /**
+ * Get the shortest distance between the surface of this sphere and a point. If
+ * the value is negative, we are inside the sphere.
+ *
+ * @param x the x coordinate
+ * @param y the y coordinate
+ * @param z the z coordinate
+ * @return the distance to the sphere
+ */
+ public double signedDistance(double x, double y, double z) {
+ Point p = new Point(x, y, z);
+ return signedDistance(p);
+ }
+
+ /**
+ * Get the shortest absolute distance between the surface of this infinite cylinder and a point.
+ * @param x the x coordinate
+ * @param y the y coordinate
+ * @param z the z coordinate
+ * @return the perpendicular distance
+ */
+ public double distance(double x, double y, double z) {
+ Point p = new Point(x, y, z);
+ return Math.abs(signedDistance(p));
+ }
+
+ /**
+ * Is the point inside the sphere?
+ * @param x the x coordinate
+ * @param y the y coordinate
+ * @param z the z coordinate
+ * @return true if the point is inside the sphere.
+ */
+ public boolean isInside(double x, double y, double z) {
+ return signedDistance(x, y, z) < 0;
+ }
+
+ /**
+ * Check whether a segment intersects the sphere
+ *
+ * @param x1 x coordinate of one end of segment
+ * @param y1 y coordinate of one end of segment
+ * @param z1 z coordinate of one end of segment
+ * @param x2 x coordinate of other end of segment
+ * @param y2 y coordinate of other end of segment
+ * @param z2 z coordinate of other end of segment
+ * @return true if the segment intersects the sphere
+ */
+ public boolean segmentIntersects(double x1, double y1, double z1, double x2, double y2, double z2) {
+ return (distToSegment(0, 0, 0, x1, y1, z1, x2, y2, z2) < _radius);
+ }
+
+ /**
+ * The closest distance of a line segment to a point
+ *
+ * @param px x coordinate of point
+ * @param py y coordinate of point
+ * @param pz z coordinate of point
+ * @param x1 x coordinate of one end of segment
+ * @param y1 y coordinate of one end of segment
+ * @param z1 z coordinate of one end of segment
+ * @param x2 x coordinate of other end of segment
+ * @param y2 y coordinate of other end of segment
+ * @param z2 z coordinate of other end of segment
+ * @return the closest distance of the segment to point p
+ */
+ private double distToSegment(double px, double py, double pz, double x1, double y1, double z1, double x2, double y2,
+ double z2) {
+
+ double line_dist = distSq(x1, y1, z1, x2, y2, z2);
+ if (line_dist == 0) {
+ return distSq(px, py, pz, x1, y1, z1);
+ }
+ double t = ((px - x1) * (x2 - x1) + (py - y1) * (y2 - y1) + (pz - z1) * (z2 - z1)) / line_dist;
+ t = Math.max(0, Math.min(1, t));
+ return Math.sqrt(distSq(px, py, pz, x1 + t * (x2 - x1), y1 + t * (y2 - y1), z1 + t * (z2 - z1)));
+ }
+
+ // the square of the distance between two points
+ private double distSq(double x1, double y1, double z1, double x2, double y2, double z2) {
+ double dx = x2 - x1;
+ double dy = y2 - y1;
+ double dz = z2 - z1;
+ return dx * dx + dy * dy + dz * dz;
+
+ }
+
+}
diff --git a/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Vector.java b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Vector.java
new file mode 100644
index 0000000000..8c00dfbd4d
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/main/java/cnuphys/CLAS12Swim/geometry/Vector.java
@@ -0,0 +1,142 @@
+package cnuphys.CLAS12Swim.geometry;
+
+/**
+ * Ordinary 3D vector
+ *
+ * @author heddle
+ *
+ */
+public class Vector extends Point {
+
+ /** Effectively zero */
+ private static final double TINY = 1.0e-20;
+
+
+ /**
+ * Create a new vector with a zero components
+ */
+ public Vector() {
+ }
+
+ /**
+ * Create a Vector from a point
+ *
+ * @param p the point
+ */
+ public Vector(Point p) {
+ this(p.x, p.y, p.z);
+ }
+
+ /**
+ * Create a vector
+ *
+ * @param x the x coordinate
+ * @param y the y coordinate
+ * @param z the z coordinate
+ */
+ public Vector(double x, double y, double z) {
+ super(x, y, z);
+ }
+
+ /**
+ * The square of the length of the vector
+ *
+ * @return the square of the length of the vector
+ */
+ public double lengthSquared() {
+ return x * x + y * y + z * z;
+ }
+
+ /**
+ * The length of the vector
+ *
+ * @return the length of the vector
+ */
+ public double length() {
+ return Math.sqrt(lengthSquared());
+ }
+
+ /**
+ * The cross product of two vectors
+ *
+ * @param a one vector
+ * @param b other vector
+ * @return c = a x b
+ */
+ public static Vector cross(Vector a, Vector b) {
+ Vector c = new Vector();
+ cross(a, b, c);
+ return c;
+ }
+
+ /**
+ * The in-place cross product of two vectors
+ *
+ * @param a one vector
+ * @param b other vector
+ * @param c on return c = a x b
+ */
+ public static void cross(Vector a, Vector b, Vector c) {
+ c.x = a.y * b.z - a.z * b.y;
+ c.y = a.z * b.x - a.x * b.z;
+ c.z = a.x * b.y - a.y * b.x;
+ }
+
+
+ /**
+ * Get a unit vector in the same direction as this
+ *
+ * @return a unit vector
+ */
+ public Vector unitVector() {
+ double len = length();
+ if (len < TINY) {
+ return null;
+ }
+
+ return new Vector(x / len, y / len, z / len);
+ }
+
+ /**
+ * Multiplies each element of a vector by a scalar.
+ *
+ * @param vector The vector to be multiplied.
+ * @param scalar The scalar value for multiplication.
+ * @return The resulting vector after multiplication.
+ */
+ public static double[] scalarMultiply(double[] vector, double scalar) {
+ double[] result = new double[vector.length];
+ for (int i = 0; i < vector.length; i++) {
+ result[i] = vector[i] * scalar;
+ }
+ return result;
+ }
+
+ /**
+ * Adds multiple vectors together element-wise.
+ *
+ * @param vectors An array of vectors to be added.
+ * @return The resulting vector after addition.
+ */
+ public static double[] addVectors(double[]... vectors) {
+ if (vectors.length == 0) {
+ throw new IllegalArgumentException("At least one vector is required for addition.");
+ }
+
+ int length = vectors[0].length;
+ for (double[] vector : vectors) {
+ if (vector.length != length) {
+ throw new IllegalArgumentException("All vectors must be of the same length.");
+ }
+ }
+
+ double[] result = new double[length];
+ for (double[] vector : vectors) {
+ for (int i = 0; i < length; i++) {
+ result[i] += vector[i];
+ }
+ }
+ return result;
+ }
+
+}
\ No newline at end of file
diff --git a/common-tools/cnuphys/clas12-swimmer/src/test/java/cnuphys/CLAS12Swim/CLAS12SwimmerTest.java b/common-tools/cnuphys/clas12-swimmer/src/test/java/cnuphys/CLAS12Swim/CLAS12SwimmerTest.java
new file mode 100644
index 0000000000..d5076a3dd0
--- /dev/null
+++ b/common-tools/cnuphys/clas12-swimmer/src/test/java/cnuphys/CLAS12Swim/CLAS12SwimmerTest.java
@@ -0,0 +1,121 @@
+package cnuphys.CLAS12Swim;
+
+import static org.junit.Assert.assertEquals;
+import static org.junit.Assert.assertFalse;
+import static org.junit.Assert.assertTrue;
+
+import org.junit.Test;
+
+import cnuphys.magfield.ZeroProbe;
+import cnuphys.CLAS12Swim.geometry.Sphere;
+
+public class CLAS12SwimmerTest {
+
+ private static final double POSITION_TOLERANCE = 1.0e-7;
+
+ @Test
+ public void chargedParticleFollowsStraightLineInZeroField() {
+ CLAS12Swimmer swimmer = new CLAS12Swimmer(new ZeroProbe());
+ double pathLength = 100.0;
+ double theta = 60.0;
+ double phi = 30.0;
+
+ CLAS12SwimResult result = swimmer.swim(1, 1.0, 2.0, 3.0, 1.0, theta, phi,
+ pathLength, 0.01, 1.0e-9);
+
+ assertTrue(result.statusString(), result.isSuccess());
+ assertEquals(pathLength, result.getPathLength(), POSITION_TOLERANCE);
+
+ double sinTheta = Math.sin(Math.toRadians(theta));
+ double[] finalState = result.getFinalU();
+ assertEquals(1.0 + pathLength * sinTheta * Math.cos(Math.toRadians(phi)), finalState[0],
+ POSITION_TOLERANCE);
+ assertEquals(2.0 + pathLength * sinTheta * Math.sin(Math.toRadians(phi)), finalState[1],
+ POSITION_TOLERANCE);
+ assertEquals(3.0 + pathLength * Math.cos(Math.toRadians(theta)), finalState[2], POSITION_TOLERANCE);
+ }
+
+ @Test
+ public void neutralParticleReachesTargetZExactly() {
+ CLAS12Swimmer swimmer = new CLAS12Swimmer(new ZeroProbe());
+
+ CLAS12SwimResult result = swimmer.swimZ(0, 0.0, 0.0, 0.0, 1.0, 60.0, 0.0,
+ 50.0, 1.0e-8, 200.0, 0.01, 1.0e-9);
+
+ assertTrue(result.statusString(), result.isSuccess());
+ assertEquals(50.0, result.getFinalU()[2], POSITION_TOLERANCE);
+ assertEquals(100.0, result.getPathLength(), POSITION_TOLERANCE);
+ }
+
+ @Test
+ public void momentumBelowThresholdIsRejected() {
+ CLAS12Swimmer swimmer = new CLAS12Swimmer(new ZeroProbe());
+
+ CLAS12SwimResult result = swimmer.swim(1, 0.0, 0.0, 0.0, 1.0e-6, 45.0, 0.0,
+ 100.0, 0.01, 1.0e-9);
+
+ assertFalse(result.isSuccess());
+ assertEquals(CLAS12Swimmer.BELOW_MIN_MOMENTUM, result.getStatus());
+ assertEquals("BELOW_MIN_MOMENTUM", result.statusString());
+ }
+
+ @Test
+ public void fixedStepSwimUsesRequestedStepSize() {
+ CLAS12Swimmer swimmer = new CLAS12Swimmer(new ZeroProbe());
+
+ CLAS12SwimResult result = swimmer.swimFixed(1, 0.0, 0.0, 0.0, 1.0, 90.0, 0.0,
+ 10.0, 2.0);
+
+ assertTrue(result.statusString(), result.isSuccess());
+ assertEquals(10.0, result.getPathLength(), POSITION_TOLERANCE);
+ assertEquals(10.0, result.getFinalU()[0], POSITION_TOLERANCE);
+ assertEquals(6, result.getNStep());
+ }
+
+ @Test
+ public void planeArrayOverloadDelegatesToPlaneSwim() {
+ CLAS12Swimmer swimmer = new CLAS12Swimmer(new ZeroProbe());
+
+ CLAS12SwimResult result = swimmer.swimPlane(1, 0.0, 0.0, 0.0, 1.0, 60.0, 0.0,
+ new double[] {0.0, 0.0, 1.0}, new double[] {0.0, 0.0, 25.0},
+ 1.0e-7, 100.0, 0.01, 1.0e-9);
+
+ assertTrue(result.statusString(), result.isSuccess());
+ assertEquals(25.0, result.getFinalU()[2], POSITION_TOLERANCE);
+ }
+
+ @Test
+ public void sphereSwimStopsAtSurface() {
+ CLAS12Swimmer swimmer = new CLAS12Swimmer(new ZeroProbe());
+ CLAS12SwimResult result = swimmer.swimSphere(1, 0.0, 0.0, 0.0, 1.0, 90.0, 0.0,
+ new Sphere(25.0), 1.0e-7, 100.0, 0.01, 1.0e-9);
+
+ assertTrue(result.statusString(), result.isSuccess());
+ assertEquals(25.0, result.getFinalU()[0], POSITION_TOLERANCE);
+ assertEquals(25.0, result.getPathLength(), POSITION_TOLERANCE);
+ }
+
+ @Test
+ public void zLineSwimFindsClosestApproach() {
+ CLAS12Swimmer swimmer = new CLAS12Swimmer(new ZeroProbe());
+ CLAS12SwimResult result = swimmer.swimZLine(1, 12.0, 7.0, 0.0, 1.0, 90.0, 180.0,
+ 2.0, 3.0, 1.0e-7, 100.0, 0.01, 1.0e-9);
+
+ assertTrue(result.statusString(), result.isSuccess());
+ assertEquals(2.0, result.getFinalU()[0], POSITION_TOLERANCE);
+ assertEquals(7.0, result.getFinalU()[1], POSITION_TOLERANCE);
+ assertEquals(10.0, result.getPathLength(), POSITION_TOLERANCE);
+ }
+
+ @Test
+ public void beamlineSwimFindsClosestApproach() {
+ CLAS12Swimmer swimmer = new CLAS12Swimmer(new ZeroProbe());
+ CLAS12SwimResult result = swimmer.swimBeamline(1, 10.0, 5.0, 0.0, 1.0, 90.0, 180.0,
+ 1.0e-7, 100.0, 0.01, 1.0e-9);
+
+ assertTrue(result.statusString(), result.isSuccess());
+ assertEquals(0.0, result.getFinalU()[0], POSITION_TOLERANCE);
+ assertEquals(5.0, result.getFinalU()[1], POSITION_TOLERANCE);
+ assertEquals(10.0, result.getPathLength(), POSITION_TOLERANCE);
+ }
+}
diff --git a/common-tools/cnuphys/pom.xml b/common-tools/cnuphys/pom.xml
index e3a0dc1eb1..34ce8d1d5c 100644
--- a/common-tools/cnuphys/pom.xml
+++ b/common-tools/cnuphys/pom.xml
@@ -18,6 +18,7 @@