diff --git a/common-tools/clas-decay-tools/src/main/java/org/jlab/clas/decay/analysis/Particle.java b/common-tools/clas-decay-tools/src/main/java/org/jlab/clas/decay/analysis/Particle.java index cbea622750..219288a909 100644 --- a/common-tools/clas-decay-tools/src/main/java/org/jlab/clas/decay/analysis/Particle.java +++ b/common-tools/clas-decay-tools/src/main/java/org/jlab/clas/decay/analysis/Particle.java @@ -997,7 +997,7 @@ private double calcUncorrMass(Particle part1, Particle part2) { //momentum not c } else { swim.SetSwimParameters(vxch, vych, vzch, -pxch, -pych, -pzch, -q); - double[] tr1 = swim.SwimToPlaneBoundary((vzvo-buffer),new Vector3D(0,0,1), -1); + double[] tr1 = swim.SwimToPlaneBoundary((vzvo-buffer),new Vector3D(0,0,1)); swim.SetSwimParameters(tr1[0], tr1[1], tr1[2], -tr1[3], -tr1[4], -tr1[5], q); } diff --git a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/ASwim.java b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/ASwim.java new file mode 100644 index 0000000000..75dff76274 --- /dev/null +++ b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/ASwim.java @@ -0,0 +1,52 @@ +package org.jlab.clas.swimtools; + +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Vector3D; + +/** + * + * @author baltzell + */ +public abstract class ASwim extends SwimPars implements ISwim { + + @Override + public double[] SwimToPlaneBoundary(double d_cm, Vector3D n) { + // Normalize the normal - should already be done, but just in case + Vector3D nhat = n.asUnit(); + + // Point on the plane at distance d_cm from origin + Point3D p = new Point3D( + nhat.x() * d_cm, + nhat.y() * d_cm, + nhat.z() * d_cm + ); + return SwimPlane(nhat, p, accuracy); + } + + @Override + public double[] SwimToPlaneLab(double z_cm) { + return SwimPlane(new Vector3D(0,0,1), new Point3D(0,0,z_cm), accuracy); + } + + //@Override + //public double[] SwimToCylinder(double radius) { + // return SwimGenCylinder(new Point3D(0,0,-1), new Point3D(0,0,1), radius, accuracy); + //} + + @Override + public double[] SwimRho(double radius, double accuracy) { + return SwimGenCylinder(new Point3D(0,0,-1), new Point3D(0,0,1), radius, accuracy); + } + + @Override + public double[] SwimToZ(double Z, int dir) { + return SwimPlane(new Vector3D(0,0,dir*1), new Point3D(0,0,Z), accuracy); + } + + @Override + public double[] SwimToBeamLine(double xB, double yB) { + return SwimToLine(new Line3D(xB,yB,-1,xB,yB,1)); + } + +} diff --git a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/AdaptiveSwim.java b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/AdaptiveSwim.java new file mode 100644 index 0000000000..d28ff1487e --- /dev/null +++ b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/AdaptiveSwim.java @@ -0,0 +1,137 @@ +package org.jlab.clas.swimtools; + +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Vector3D; + +import cnuphys.adaptiveSwim.geometry.Line; +import cnuphys.adaptiveSwim.geometry.Point; +import cnuphys.adaptiveSwim.geometry.Vector; +import cnuphys.adaptiveSwim.geometry.Sphere; +import cnuphys.adaptiveSwim.geometry.Cylinder; + +import cnuphys.adaptiveSwim.AdaptiveSwimException; +import cnuphys.adaptiveSwim.AdaptiveSwimResult; +import cnuphys.adaptiveSwim.AdaptiveSwimmer; +import cnuphys.adaptiveSwim.geometry.Plane; + +public class AdaptiveSwim extends ASwim { + + private static double[] convert(AdaptiveSwimResult result, double p) { + + double[] value = null; + + if (result.getStatus() == AdaptiveSwimmer.SWIM_SUCCESS) { + value = new double[8]; + value[0] = result.getUf()[0] * 100; // convert back to cm + value[1] = result.getUf()[1] * 100; // convert back to cm + value[2] = result.getUf()[2] * 100; // convert back to cm + value[3] = result.getUf()[3] * p; // normalized values + value[4] = result.getUf()[4] * p; + value[5] = result.getUf()[5] * p; + value[6] = result.getFinalS() * 100; + value[7] = 0; // Conversion from kG.m to T.cm + } + + return value; + } + + @Override + public double[] SwimRho(double radius, double accuracy) { + + // convert to meters: + radius = radius/100; + + try { + AdaptiveSwimResult result = new AdaptiveSwimResult(false); + PC.AS.swimRho(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, radius, + accuracy/100, _rMax, stepSize, cnuphys.swim.Swimmer.getEps(), result); + return convert(result, _pTot); + } + catch (AdaptiveSwimException e) { + e.printStackTrace(); + } + return null; + } + + @Override + public double[] SwimGenCylinder(Point3D axisPoint1, Point3D axisPoint2, double radius, double accuracy) { + + // convert to meters: + radius = radius/100; + Point a1 = new Point(axisPoint1.x()/100, axisPoint1.y()/100, axisPoint1.z()/100); + Point a2 = new Point(axisPoint2.x()/100, axisPoint2.y()/100, axisPoint2.z()/100); + Cylinder targetCylinder = new Cylinder(new Line(a1,a2), radius); + + try { + AdaptiveSwimResult result = new AdaptiveSwimResult(false); + PC.AS.swimCylinder(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, targetCylinder, + accuracy/100, _rMax, stepSize, cnuphys.swim.Swimmer.getEps(), result); + return convert(result, _pTot); + } + catch (AdaptiveSwimException e) { + e.printStackTrace(); + } + return null; + } + + @Override + public double[] SwimPlane(Vector3D n, Point3D p, double accuracy) { + + // convert to meters: + Vector norm = new Vector(n.asUnit().x(), n.asUnit().y(), n.asUnit().z()); + Point point = new Point(p.x()/100, p.y()/100, p.z()/100); + Plane targetPlane = new Plane(norm, point); + + try { + AdaptiveSwimResult result = new AdaptiveSwimResult(false); + PC.AS.swimPlane(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, targetPlane, + accuracy/100, _rMax, stepSize, cnuphys.swim.Swimmer.getEps(), result); + return convert(result, _pTot); + } + catch (AdaptiveSwimException e) { + e.printStackTrace(); + } + return null; + } + + @Override + public double[] SwimToSphere(double Rad) { + + // convert to meters: + Sphere targetSphere = new Sphere(new Point(0,0,0), Rad/100); + + try { + AdaptiveSwimResult result = new AdaptiveSwimResult(false); + PC.AS.swimSphere(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, targetSphere, + accuracy/100, _rMax, stepSize, cnuphys.swim.Swimmer.getEps(), result); + return convert(result, _pTot); + } + catch (AdaptiveSwimException e) { + e.printStackTrace(); + } + return null; + } + + @Override + public double[] SwimToLine(Line3D l) { + + // convert to meters: + Point a1 = new Point(l.origin().x()/100, l.origin().y()/100, l.origin().z()/100); + Point a2 = new Point(l.end().x()/100, l.end().y()/100, l.end().z()/100); + Line targetLine = new Line(a1, a2); + + try { + AdaptiveSwimResult result = new AdaptiveSwimResult(false); + PC.AS.swimLine(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, targetLine, + accuracy/100, _rMax, stepSize, cnuphys.swim.Swimmer.getEps(), result); + return convert(result, _pTot); + } + catch (AdaptiveSwimException e) { + e.printStackTrace(); + } + return null; + } + + +} diff --git a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/ISwim.java b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/ISwim.java new file mode 100644 index 0000000000..5fd721d968 --- /dev/null +++ b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/ISwim.java @@ -0,0 +1,33 @@ +package org.jlab.clas.swimtools; + +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Vector3D; + +/** + * + * @author baltzell + */ +interface ISwim { + + public double[] SwimToPlaneLab(double z_cm); + + //public double[] SwimToCylinder(double Rad); + + public double[] SwimRho(double radius, double accuracy); + + public double[] SwimGenCylinder(Point3D axisPoint1, Point3D axisPoint2, double radius, double accuracy); + + public double[] SwimPlane(Vector3D n, Point3D p, double accuracy); + + public double[] SwimToSphere(double Rad); + + public double[] SwimToPlaneBoundary(double d_cm, Vector3D n); + + public double[] SwimToBeamLine(double xB, double yB); + + public double[] SwimToLine(Line3D l); + + public double[] SwimToZ(double Z, int dir); + +} diff --git a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/MagFieldsEngine.java b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/MagFieldsEngine.java index ceab50389a..b02d767101 100644 --- a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/MagFieldsEngine.java +++ b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/MagFieldsEngine.java @@ -1,8 +1,9 @@ package org.jlab.clas.swimtools; +import cnuphys.magfield.MagneticFieldInitializationException; import cnuphys.magfield.MagneticFields; +import java.io.FileNotFoundException; import java.util.Arrays; -import java.util.concurrent.atomic.AtomicInteger; import java.util.logging.Level; import java.util.logging.Logger; @@ -14,7 +15,7 @@ public class MagFieldsEngine extends ReconstructionEngine { - public static Logger LOGGER = Logger.getLogger(MagFieldsEngine.class.getName()); + public static final Logger LOGGER = Logger.getLogger(MagFieldsEngine.class.getName()); private String solShift = null; @@ -22,8 +23,6 @@ public MagFieldsEngine() { super("MagFields", "ziegler", "1.0"); } - AtomicInteger Run = new AtomicInteger(0); - /** * Choose one of YAML or ENV values. * @@ -57,19 +56,19 @@ public boolean initializeMagneticFields() { final String mapDir = CLASResources.getResourcePath("etc")+"/data/magfield"; if (torusMap==null) { - LOGGER.log(Level.SEVERE,"["+this.getName()+"] ERROR: torus field is undefined."); + LOGGER.log(Level.SEVERE, "[{0}] ERROR: torus field is undefined.", this.getName()); return false; } if (solenoidMap==null) { - LOGGER.log(Level.SEVERE,"["+this.getName()+"] ERROR: solenoid is undefined."); + LOGGER.log(Level.SEVERE, "[{0}] ERROR: solenoid is undefined.", this.getName()); return false; } try { MagneticFields.getInstance().initializeMagneticFields(mapDir, torusMap, solenoidMap); } - catch (Exception e) { - e.printStackTrace(); + catch (MagneticFieldInitializationException | FileNotFoundException e) { + LOGGER.log(Level.SEVERE,"Magfields error",e); return false; } @@ -77,77 +76,67 @@ public boolean initializeMagneticFields() { solShift = this.getEngineConfigString("magfieldSolenoidShift"); if (solShift != null) { - LOGGER.log(Level.INFO, "[" + this.getName() - + "] run with solenoid z shift in tracking config chosen based on yaml = " + solShift + " cm"); - Swimmer.set_zShift(Float.valueOf(solShift)); + LOGGER.log(Level.INFO, "[{0}] run with solenoid z shift in tracking config chosen based on yaml = {1} cm", new Object[]{this.getName(), solShift}); + Swimmer.set_zShift(Float.parseFloat(solShift)); } else { solShift = System.getenv("COAT_MAGFIELD_SOLENOIDSHIFT"); if (solShift != null) { - LOGGER.log(Level.INFO, "[" + this.getName() - + "] run with solenoid z shift in tracking config chosen based on env = " + solShift + " cm"); - Swimmer.set_zShift(Float.valueOf(solShift)); + LOGGER.log(Level.INFO, "[{0}] run with solenoid z shift in tracking config chosen based on env = {1} cm", new Object[]{this.getName(), solShift}); + Swimmer.set_zShift(Float.parseFloat(solShift)); } } if (solShift == null) { - LOGGER.log(Level.INFO, "[" + this.getName() + "] run with solenoid z shift based on CCDB CD position"); + LOGGER.log(Level.INFO, "[{0}] run with solenoid z shift based on CCDB CD position", this.getName()); // this.solenoidShift = (float) 0; } // torus: String TorX = this.getEngineConfigString("magfieldTorusXShift"); if (TorX != null) { - LOGGER.log(Level.INFO, "[" + this.getName() - + "] run with torus x shift in tracking config chosen based on yaml = " + TorX + " cm"); - Swimmer.setTorXShift(Float.valueOf(TorX)); + LOGGER.log(Level.INFO, "[{0}] run with torus x shift in tracking config chosen based on yaml = {1} cm", new Object[]{this.getName(), TorX}); + Swimmer.setTorXShift(Float.parseFloat(TorX)); } else { TorX = System.getenv("COAT_MAGFIELD_TORUSXSHIFT"); if (TorX != null) { - LOGGER.log(Level.INFO, "[" + this.getName() - + "] run with torus x shift in tracking config chosen based on env = " + TorX + " cm"); - Swimmer.setTorXShift(Float.valueOf(TorX)); + LOGGER.log(Level.INFO, "[{0}] run with torus x shift in tracking config chosen based on env = {1} cm", new Object[]{this.getName(), TorX}); + Swimmer.setTorXShift(Float.parseFloat(TorX)); } } if (TorX == null) { - LOGGER.log(Level.INFO, "[" + this.getName() + "] run with torus x shift in tracking set to 0 cm"); + LOGGER.log(Level.INFO, "[{0}] run with torus x shift in tracking set to 0 cm", this.getName()); // this.solenoidShift = (float) 0; } String TorY = this.getEngineConfigString("magfieldTorusYShift"); if (TorY != null) { - LOGGER.log(Level.INFO, "[" + this.getName() - + "] run with torus y shift in tracking config chosen based on yaml = " + TorY + " cm"); - Swimmer.setTorYShift(Float.valueOf(TorY)); + LOGGER.log(Level.INFO, "[{0}] run with torus y shift in tracking config chosen based on yaml = {1} cm", new Object[]{this.getName(), TorY}); + Swimmer.setTorYShift(Float.parseFloat(TorY)); } else { TorY = System.getenv("COAT_MAGFIELD_TORUSYSHIFT"); if (TorY != null) { - LOGGER.log(Level.INFO, "[" + this.getName() - + "] run with torus y shift in tracking config chosen based on env = " + TorY + " cm"); - Swimmer.setTorYShift(Float.valueOf(TorY)); + LOGGER.log(Level.INFO, "[{0}] run with torus y shift in tracking config chosen based on env = {1} cm", new Object[]{this.getName(), TorY}); + Swimmer.setTorYShift(Float.parseFloat(TorY)); } } if (TorY == null) { - LOGGER.log(Level.INFO, "[" + this.getName() + "] run with torus y shift in tracking set to 0 cm"); - // this.solenoidShift = (float) 0; + LOGGER.log(Level.INFO, "[{0}] run with torus y shift in tracking set to 0 cm", this.getName()); } String TorZ = this.getEngineConfigString("magfieldTorusZShift"); if (TorZ != null) { - LOGGER.log(Level.INFO, "[" + this.getName() - + "] run with torus z shift in tracking config chosen based on yaml = " + TorZ + " cm"); - Swimmer.setTorZShift(Float.valueOf(TorZ)); + LOGGER.log(Level.INFO, "[{0}] run with torus z shift in tracking config chosen based on yaml = {1} cm", new Object[]{this.getName(), TorZ}); + Swimmer.setTorZShift(Float.parseFloat(TorZ)); } else { TorZ = System.getenv("COAT_MAGFIELD_TORUSZSHIFT"); if (TorZ != null) { - LOGGER.log(Level.INFO, "[" + this.getName() - + "] run with torus z shift in tracking config chosen based on env = " + TorZ + " cm"); - Swimmer.setTorZShift(Float.valueOf(TorZ)); + LOGGER.log(Level.INFO, "[{0}] run with torus z shift in tracking config chosen based on env = {1} cm", new Object[]{this.getName(), TorZ}); + Swimmer.setTorZShift(Float.parseFloat(TorZ)); } } if (TorZ == null) { - LOGGER.log(Level.INFO, "[" + this.getName() + "] run with torus z shift in tracking set to 0 cm"); - // this.solenoidShift = (float) 0; + LOGGER.log(Level.INFO, "[{0}] run with torus z shift in tracking set to 0 cm", this.getName()); } return true; @@ -170,9 +159,10 @@ public boolean processDataEvent(DataEvent event) { if (newRun == 0) return true; - if (solShift == null) { // if no shift is set in the yaml file or environment, read from CCDB - // will read target position and assume that is representative of the shift of - // the whole CD + if (solShift == null) { + // if no shift is set in the yaml file or environment, read from CCDB + // will read target position and assume that is representative of the + // shift of the whole CD IndexedTable targetPosition = this.getConstantsManager().getConstants(newRun, "/geometry/shifts/solenoid"); Swimmer.set_zShift((float) targetPosition.getDoubleValue("z", 0, 0, 0)); } diff --git a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/ProbeCollection.java b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/ProbeCollection.java index 043ffec40e..f410927af6 100644 --- a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/ProbeCollection.java +++ b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/ProbeCollection.java @@ -1,8 +1,3 @@ -/* - * To change this license header, choose License Headers in Project Properties. - * To change this template file, choose Tools | Templates - * and open the template in the editor. - */ package org.jlab.clas.swimtools; import cnuphys.magfield.CompositeProbe; diff --git a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Stoppers.java b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Stoppers.java new file mode 100644 index 0000000000..c9deacda41 --- /dev/null +++ b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Stoppers.java @@ -0,0 +1,328 @@ +package org.jlab.clas.swimtools; + +import cnuphys.rk4.IStopper; +import cnuphys.swim.SwimTrajectory; +import java.util.ArrayList; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Vector3D; + +/** + * + * @author ziegler + */ +public class Stoppers { + + public static class CylindricalBoundarySwimStopper implements IStopper { + + private double _finalPathLength = Double.NaN; + private double _Rad; + double max = -1.0; + + /** + * A swim stopper that will stop if the boundary of a plane is crossed + * + * @param Rad the max radial coordinate in meters. + */ + public CylindricalBoundarySwimStopper(double Rad) { + _Rad = Rad; + } + + @Override + public boolean stopIntegration(double t, double[] y) { + double r = Math.sqrt(y[0] * y[0] + y[1] * y[1]) * 100.; + return (r < max || r > _Rad); // stop intergration at closest distance to the cylinder + } + + /** + * Get the final path length in meters + * + * @return the final path length in meters + */ + @Override + public double getFinalT() { + return _finalPathLength; + } + + /** + * Set the final path length in meters + * + * @param finalPathLength + * the final path length in meters + */ + @Override + public void setFinalT(double finalPathLength) { + _finalPathLength = finalPathLength; + } + } + + public static class SphericalBoundarySwimStopper implements IStopper { + + private double _finalPathLength = Double.NaN; + private double _Rad; + + /** + * A swim stopper that will stop if the boundary of a plane is crossed + * + * @param maxR the max radial coordinate in meters. + */ + public SphericalBoundarySwimStopper(double Rad) { + _Rad = Rad; + } + + @Override + public boolean stopIntegration(double t, double[] y) { + double r = Math.sqrt(y[0] * y[0] + y[1] * y[1] + y[2] * y[2]) * 100.; + return r > _Rad; + } + + /** + * Get the final path length in meters + * + * @return the final path length in meters + */ + @Override + public double getFinalT() { + return _finalPathLength; + } + + /** + * Set the final path length in meters + * + * @param finalPathLength + * the final path length in meters + */ + @Override + public void setFinalT(double finalPathLength) { + _finalPathLength = finalPathLength; + } + } + + // added for swimming to outer detectors + public static class PlaneBoundarySwimStopper implements IStopper { + + private double _finalPathLength = Double.NaN; + private double _d; + private Vector3D _n; + private double _dist2plane; + private int _dir; + + /** + * A swim stopper that will stop if the boundary of a plane is crossed + * + * @param d + * @param dir + * @param n + */ + public PlaneBoundarySwimStopper(double d, Vector3D n, int dir) { + _d = d; + _n = n; + _dir = dir; + } + + @Override + public boolean stopIntegration(double t, double[] y) { + double dtrk = y[0] * _n.x() + y[1] * _n.y() + y[2] * _n.z(); + //double accuracy = 20e-6; // 20 microns + return _dir < 0 ? dtrk < _d : dtrk > _d; + } + + @Override + public double getFinalT() { + return _finalPathLength; + } + + /** + * Set the final path length in meters + * + * @param finalPathLength the final path length in meters + */ + @Override + public void setFinalT(double finalPathLength) { + _finalPathLength = finalPathLength; + } + } + + public static class LineSwimStopper implements IStopper { + + private double _finalPathLength = Double.NaN; + private Line3D _l; + private Point3D _p; + double min = 999; + public LineSwimStopper(Line3D l) { + _l =l; + _p = new Point3D(); + } + + @Override + public boolean stopIntegration(double t, double[] y) { + _p.set(y[0]* 100.0, y[1]* 100.0, y[2]* 100.0); + double doca = _l.distance(_p).length(); + if(doca min ); + } + + /** + * Get the final path length in meters + * + * @return the final path length in meters + */ + @Override + public double getFinalT() { + return _finalPathLength; + } + + /** + * Set the final path length in meters + * + * @param finalPathLength + * the final path length in meters + */ + @Override + public void setFinalT(double finalPathLength) { + _finalPathLength = finalPathLength; + } + } + + public static class ZSwimStopper implements IStopper { + + private double _finalPathLength = Double.NaN; + private double _Z; + private int _dir; + + public ZSwimStopper(double Z, int dir) { + _Z = Z; + _dir = dir; + } + + @Override + public boolean stopIntegration(double t, double[] y) { + double z = y[2] * 100.; + return _dir>0 ? z>_Z : z<_Z; + } + + /** + * Get the final path length in meters + * + * @return the final path length in meters + */ + @Override + public double getFinalT() { + return _finalPathLength; + } + + /** + * Set the final path length in meters + * + * @param finalPathLength + * the final path length in meters + */ + @Override + public void setFinalT(double finalPathLength) { + _finalPathLength = finalPathLength; + } + } + + public static class DCASwimStopper implements IStopper { + + public DCASwimStopper(SwimTrajectory swimTraj) { + _swimTraj = swimTraj; + for(int i = 0; i < _swimTraj.size()-1; i++) { + polylines.add(new Line3D(_swimTraj.get(i)[0],_swimTraj.get(i)[1],_swimTraj.get(i)[2], + _swimTraj.get(i+1)[0],_swimTraj.get(i+1)[1],_swimTraj.get(i+1)[2])); + } + } + + private ArrayList polylines = new ArrayList<>(); + private SwimTrajectory _swimTraj; + private double _finalPathLength = Double.NaN; + private double _doca = Double.POSITIVE_INFINITY; + + @Override + public boolean stopIntegration(double t, double[] y) { + + Point3D dcaCand = new Point3D(y[0],y[1],y[2]); + double maxDoca = Double.POSITIVE_INFINITY; + + for(Line3D l : polylines) { + if(l.distance(dcaCand).length() min ; + + } + + /** + * Get the final path length in meters + * + * @return the final path length in meters + */ + @Override + public double getFinalT() { + return _finalPathLength; + } + + /** + * Set the final path length in meters + * + * @param finalPathLength + * the final path length in meters + */ + @Override + public void setFinalT(double finalPathLength) { + _finalPathLength = finalPathLength; + } + } + +} \ No newline at end of file diff --git a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java index 66897dc346..2dcd98ef5e 100644 --- a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java +++ b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java @@ -1,1431 +1,378 @@ -/* - * To change this license header, choose License Headers in Project Properties. - * To change this template file, choose Tools | Templates - * and open the template in the editor. - */ package org.jlab.clas.swimtools; -import cnuphys.adaptiveSwim.AdaptiveSwimException; -import cnuphys.rk4.IStopper; +import org.jlab.geom.prim.Vector3D; +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Line3D; + import cnuphys.rk4.RungeKuttaException; + import cnuphys.swim.SwimTrajectory; import cnuphys.swim.util.Plane; + import cnuphys.swimZ.SwimZException; import cnuphys.swimZ.SwimZResult; import cnuphys.swimZ.SwimZStateVector; -import org.apache.commons.math3.util.FastMath; -import org.jlab.geom.prim.Vector3D; -import org.jlab.geom.prim.Point3D; + import cnuphys.adaptiveSwim.AdaptiveSwimResult; -import cnuphys.adaptiveSwim.AdaptiveSwimmer; -import cnuphys.adaptiveSwim.geometry.Cylinder; -import cnuphys.adaptiveSwim.geometry.Line; -import cnuphys.adaptiveSwim.geometry.Point; -import cnuphys.adaptiveSwim.geometry.Vector; -import java.util.ArrayList; -import java.util.List; -import org.jlab.geom.prim.Line3D; +import cnuphys.rk4.IStopper; +import org.jlab.clas.swimtools.Stoppers.BeamLineSwimStopper; + +import org.jlab.clas.swimtools.Stoppers.CylindricalBoundarySwimStopper; +import org.jlab.clas.swimtools.Stoppers.DCASwimStopper; +import org.jlab.clas.swimtools.Stoppers.LineSwimStopper; +import org.jlab.clas.swimtools.Stoppers.SphericalBoundarySwimStopper; +import org.jlab.clas.swimtools.Stoppers.ZSwimStopper; + /** + * Class for swimming to various surfaces. The input and output units are cm and GeV/c * * @author ziegler */ +public class Swim extends ASwim { -public class Swim { + private static final int STATE_SIZE = 8; + private static final double CM_PER_M = 100.0; + private static final double KGCM_TO_TCM = 0.1; // divide by 10 - private double _x0; - private double _y0; - private double _z0; - private double _phi; - private double _theta; - private double _pTot; - private final double _rMax = 5 + 3; // increase to allow swimming to outer - // detectors - private double _maxPathLength = 9; - private boolean SwimUnPhys = false; //Flag to indicate if track is swimmable - private int _charge; + private SwimTrajectory swimTraj; - final double SWIMZMINMOM = 0.75; // GeV/c - final double MINTRKMOM = 0.05; // GeV/c - double accuracy = 20e-6; // 20 microns - public double stepSize = 5.00 * 1.e-4; // 500 microns - public double distanceBetweenSaves= 100*stepSize; - - private ProbeCollection PC; - - /** - * Class for swimming to various surfaces. The input and output units are cm and GeV/c - */ - public Swim() { - PC = Swimmer.getProbeCollection(Thread.currentThread()); - if (PC == null) { - PC = new ProbeCollection(); - Swimmer.put(Thread.currentThread(), PC); - } - } - - /** - * Set max swimming path length - * - * @param maxPathLength - */ - public void setMaxPathLength(double _maxPathLength) { - this._maxPathLength = _maxPathLength; + public SwimTrajectory getSwimTraj() { + return swimTraj; } - /** - * - * @param direction - * +1 for out -1 for in - * @param x0 - * @param y0 - * @param z0 - * @param thx - * @param thy - * @param p - * @param charge - */ - public void SetSwimParameters(int direction, double x0, double y0, double z0, double thx, double thy, double p, - int charge) { - - // x,y,z in m = swimmer units - _x0 = x0 / 100; - _y0 = y0 / 100; - _z0 = z0 / 100; - this.checkR(_x0, _y0, _z0); - double pz = direction * p / Math.sqrt(thx * thx + thy * thy + 1); - double px = thx * pz; - double py = thy * pz; - _phi = Math.toDegrees(FastMath.atan2(py, px)); - _pTot = Math.sqrt(px * px + py * py + pz * pz); - _theta = Math.toDegrees(Math.acos(pz / _pTot)); - - _charge = direction * charge; + public void setSwimTraj(SwimTrajectory swimTraj) { + this.swimTraj = swimTraj; } - /** - * Sets the parameters used by swimmer based on the input track state vector - * parameters swimming outwards - * - * @param superlayerIdx - * @param layerIdx - * @param x0 - * @param y0 - * @param thx - * @param thy - * @param p - * @param charge - */ - public void SetSwimParameters(int superlayerIdx, int layerIdx, double x0, double y0, double z0, double thx, - double thy, double p, int charge) { - // z at a given DC plane in the tilted coordinate system - // x,y,z in m = swimmer units - _x0 = x0 / 100; - _y0 = y0 / 100; - _z0 = z0 / 100; - this.checkR(_x0, _y0, _z0); - double pz = p / Math.sqrt(thx * thx + thy * thy + 1); - double px = thx * pz; - double py = thy * pz; - _phi = Math.toDegrees(FastMath.atan2(py, px)); - _pTot = Math.sqrt(px * px + py * py + pz * pz); - _theta = Math.toDegrees(Math.acos(pz / _pTot)); - - _charge = charge; + // ------------------------------------------------------------------------ + // Common guards and helpers + // ------------------------------------------------------------------------ + private boolean isInvalid() { + return SwimUnPhys || _pTot < MINTRKMOM; } - /** - * Sets the parameters used by swimmer based on the input track parameters - * - * @param x0 - * @param y0 - * @param z0 - * @param px - * @param py - * @param pz - * @param charge - */ - public void SetSwimParameters(double x0, double y0, double z0, double px, double py, double pz, int charge) { - _x0 = x0 / 100; - _y0 = y0 / 100; - _z0 = z0 / 100; - this.checkR(_x0, _y0, _z0); - _phi = Math.toDegrees(FastMath.atan2(py, px)); - _pTot = Math.sqrt(px * px + py * py + pz * pz); - _theta = Math.toDegrees(Math.acos(pz / _pTot)); - - _charge = charge; - + private double[] newState() { + return new double[STATE_SIZE]; } - /** - * - * @param xcm - * @param ycm - * @param zcm - * @param phiDeg - * @param thetaDeg - * @param p - * @param charge - * @param maxPathLength - */ - public void SetSwimParameters(double xcm, double ycm, double zcm, double phiDeg, double thetaDeg, double p, - int charge, double maxPathLength) { - - _maxPathLength = maxPathLength; - _charge = charge; - _phi = phiDeg; - _theta = thetaDeg; - _pTot = p; - _x0 = xcm / 100; - _y0 = ycm / 100; - _z0 = zcm / 100; - this.checkR(_x0, _y0, _z0); + private void fillFromTrajectory(double[] out, SwimTrajectory traj) { + double[] y = traj.lastElement(); + out[0] = y[0] * CM_PER_M; + out[1] = y[1] * CM_PER_M; + out[2] = y[2] * CM_PER_M; + out[3] = y[3] * _pTot; + out[4] = y[4] * _pTot; + out[5] = y[5] * _pTot; + out[6] = y[6] * CM_PER_M; + out[7] = y[7] * 10.0; // kG·m → T·cm } - /** - * - * @param xcm - * @param ycm - * @param zcm - * @param phiDeg - * @param thetaDeg - * @param p - * @param charge - * @param maxPathLength - * @param Accuracy - * @param StepSize - */ - public void SetSwimParameters(double xcm, double ycm, double zcm, double phiDeg, double thetaDeg, double p, - int charge, double maxPathLength, double Accuracy, double StepSize) { - - _maxPathLength = maxPathLength; - accuracy = Accuracy/100; - stepSize = StepSize/100; - _charge = charge; - _phi = phiDeg; - _theta = thetaDeg; - _pTot = p; - _x0 = xcm / 100; - _y0 = ycm / 100; - _z0 = zcm / 100; - this.checkR(_x0, _y0, _z0); + private void fillFromZResult(double[] out, SwimZResult szr, double bdlKgCm) { + SwimZStateVector last = szr.last(); + double[] p3 = szr.getThreeMomentum(last); + out[0] = last.x; + out[1] = last.y; + out[2] = last.z; + out[3] = p3[0]; + out[4] = p3[1]; + out[5] = p3[2]; + out[6] = szr.getPathLength(); + out[7] = bdlKgCm * KGCM_TO_TCM; } - public double[] SwimToPlaneTiltSecSys(int sector, double z_cm) { - double z = z_cm / 100; // the magfield method uses meters - double[] value = new double[8]; - - if (_pTot < MINTRKMOM || this.SwimUnPhys==true) // fiducial cut - { + private SwimZResult tryZSwimSector(int sector, double z_cm, double[] hdata) { + if (_pTot <= SWIMZMINMOM) return null; + try { + double stepSizeCM = stepSize * CM_PER_M; + SwimZStateVector start = new SwimZStateVector( + _x0 * CM_PER_M, _y0 * CM_PER_M, _z0 * CM_PER_M, + _pTot, _theta, _phi); + return PC.RCF_z.sectorAdaptiveRK(sector, _charge, _pTot, start, z_cm, stepSizeCM, hdata); + } catch (SwimZException e) { return null; } + } - // use a SwimZResult instead of a trajectory (dph) - SwimZResult szr = null; - - SwimTrajectory traj = null; - double hdata[] = new double[3]; - + private SwimZResult tryZSwimLab(double z_cm, double[] hdata) { + if (_pTot <= SWIMZMINMOM) return null; try { + double stepSizeCM = stepSize * CM_PER_M; + SwimZStateVector start = new SwimZStateVector( + _x0 * CM_PER_M, _y0 * CM_PER_M, _z0 * CM_PER_M, + _pTot, _theta, _phi); + return PC.CF_z.adaptiveRK(_charge, _pTot, start, z_cm, stepSizeCM, hdata); + } catch (SwimZException e) { + return null; + } + } - if (_pTot > SWIMZMINMOM) { - - // use the new z swimmer (dph) - // NOTE THE DISTANCE, UNITS FOR swimZ are cm, NOT m like the old - // swimmer (dph) - - double stepSizeCM = stepSize * 100; // convert to cm - - // create the starting SwimZ state vector - SwimZStateVector start = new SwimZStateVector(_x0 * 100, _y0 * 100, _z0 * 100, _pTot, _theta, _phi); - - try { - szr = PC.RCF_z.sectorAdaptiveRK(sector, _charge, _pTot, start, z_cm, stepSizeCM, hdata); - } catch (SwimZException e) { - szr = null; - //System.err.println("[WARNING] Tilted SwimZ Failed for p = " + _pTot); - } - } - - if (szr != null) { - double bdl = szr.sectorGetBDL(sector, PC.RCF_z.getProbe()); - double pathLength = szr.getPathLength(); // already in cm - - SwimZStateVector last = szr.last(); - double p3[] = szr.getThreeMomentum(last); - - value[0] = last.x; // xf in cm - value[1] = last.y; // yz in cm - value[2] = last.z; // zf in cm - value[3] = p3[0]; - value[4] = p3[1]; - value[5] = p3[2]; - value[6] = pathLength; - value[7] = bdl / 10; // convert from kg*cm to T*cm - } else { // use old swimmer. Either low momentum or SwimZ failed. - // (dph) + // ------------------------------------------------------------------------ + // Plane / Z swimmers + // ------------------------------------------------------------------------ - traj = PC.RCF.sectorSwim(sector, _charge, _x0, _y0, _z0, _pTot, _theta, _phi, z, accuracy, _rMax, - _maxPathLength, stepSize, cnuphys.swim.Swimmer.CLAS_Tolerance, hdata); + public double[] SwimToPlaneTiltSecSys(int sector, double z_cm) { + if (isInvalid()) return null; - // traj.computeBDL(sector, rprob); - if(traj==null) - return null; - - traj.sectorComputeBDL(sector, PC.RCP); - // traj.computeBDL(rcompositeField); + double[] hdata = new double[3]; + double[] out = newState(); - double lastY[] = traj.lastElement(); - value[0] = lastY[0] * 100; // convert back to cm - value[1] = lastY[1] * 100; // convert back to cm - value[2] = lastY[2] * 100; // convert back to cm - value[3] = lastY[3] * _pTot; - value[4] = lastY[4] * _pTot; - value[5] = lastY[5] * _pTot; - value[6] = lastY[6] * 100; - value[7] = lastY[7] * 10; - } // use old swimmer - } catch (Exception e) { - e.printStackTrace(); + SwimZResult szr = tryZSwimSector(sector, z_cm, hdata); + if (szr != null) { + double bdl = szr.sectorGetBDL(sector, PC.RCF_z.getProbe()); + fillFromZResult(out, szr, bdl); + return out; } - return value; - } - - public double[] SwimToPlaneTiltSecSysBdlXZPlane(int sector, double z_cm) { - double z = z_cm / 100; // the magfield method uses meters - double[] value = new double[8]; - - if (_pTot < MINTRKMOM || this.SwimUnPhys==true) // fiducial cut - { + try { + double z_m = z_cm / CM_PER_M; + SwimTrajectory traj = PC.RCF.sectorSwim( + sector, _charge, _x0, _y0, _z0, _pTot, _theta, _phi, + z_m, accuracy, _rMax, _maxPathLength, stepSize, + cnuphys.swim.Swimmer.CLAS_Tolerance, hdata); + if (traj == null) return null; + traj.sectorComputeBDL(sector, PC.RCP); + fillFromTrajectory(out, traj); + return out; + } catch (RungeKuttaException e) { return null; } + } - // use a SwimZResult instead of a trajectory (dph) - SwimZResult szr = null; - - SwimTrajectory traj = null; - double hdata[] = new double[3]; - - try { - - if (_pTot > SWIMZMINMOM) { - - // use the new z swimmer (dph) - // NOTE THE DISTANCE, UNITS FOR swimZ are cm, NOT m like the old - // swimmer (dph) - - double stepSizeCM = stepSize * 100; // convert to cm - - // create the starting SwimZ state vector - SwimZStateVector start = new SwimZStateVector(_x0 * 100, _y0 * 100, _z0 * 100, _pTot, _theta, _phi); - - try { - szr = PC.RCF_z.sectorAdaptiveRK(sector, _charge, _pTot, start, z_cm, stepSizeCM, hdata); - } catch (SwimZException e) { - szr = null; - //System.err.println("[WARNING] Tilted SwimZ Failed for p = " + _pTot); - } - } - - if (szr != null) { - double bdl = szr.sectorGetBDLXZPlane(sector, PC.RCF_z.getProbe()); - double pathLength = szr.getPathLength(); // already in cm - - SwimZStateVector last = szr.last(); - double p3[] = szr.getThreeMomentum(last); - - value[0] = last.x; // xf in cm - value[1] = last.y; // yz in cm - value[2] = last.z; // zf in cm - value[3] = p3[0]; - value[4] = p3[1]; - value[5] = p3[2]; - value[6] = pathLength; - value[7] = bdl / 10; // convert from kg*cm to T*cm - } else { // use old swimmer. Either low momentum or SwimZ failed. - // (dph) - - traj = PC.RCF.sectorSwim(sector, _charge, _x0, _y0, _z0, _pTot, _theta, _phi, z, accuracy, _rMax, - _maxPathLength, stepSize, cnuphys.swim.Swimmer.CLAS_Tolerance, hdata); - - // traj.computeBDL(sector, rprob); - if(traj==null) - return null; - - traj.sectorComputeBDL(sector, PC.RCP); - // traj.computeBDL(rcompositeField); - - double lastY[] = traj.lastElement(); - value[0] = lastY[0] * 100; // convert back to cm - value[1] = lastY[1] * 100; // convert back to cm - value[2] = lastY[2] * 100; // convert back to cm - value[3] = lastY[3] * _pTot; - value[4] = lastY[4] * _pTot; - value[5] = lastY[5] * _pTot; - value[6] = lastY[6] * 100; - value[7] = lastY[7] * 10; - } // use old swimmer - } catch (Exception e) { - e.printStackTrace(); - } - return value; + public double[] SwimToPlaneTiltSecSysBdlXZPlane(int sector, double z_cm) { + if (isInvalid()) return null; - } - - /** - * - * @param z_cm - * @return state x,y,z,px,py,pz, pathlength, iBdl at the plane surface - */ - public double[] SwimToPlaneLab(double z_cm) { - double z = z_cm / 100; // the magfield method uses meters - double[] value = new double[8]; + double[] hdata = new double[3]; + double[] out = newState(); - if (_pTot < MINTRKMOM || this.SwimUnPhys==true) // fiducial cut - { - return null; + SwimZResult szr = tryZSwimSector(sector, z_cm, hdata); + if (szr != null) { + double bdl = szr.sectorGetBDLXZPlane(sector, PC.RCF_z.getProbe()); + fillFromZResult(out, szr, bdl); + return out; } - SwimTrajectory traj = null; - double hdata[] = new double[3]; - - // use a SwimZResult instead of a trajectory (dph) - SwimZResult szr = null; try { - - if (_pTot > SWIMZMINMOM) { - - // use the new z swimmer (dph) - // NOTE THE DISTANCE, UNITS FOR swimZ are cm, NOT m like the old - // swimmer (dph) - - double stepSizeCM = stepSize * 100; // convert to cm - - // create the starting SwimZ state vector - SwimZStateVector start = new SwimZStateVector(_x0 * 100, _y0 * 100, _z0 * 100, _pTot, _theta, _phi); - - try { - szr = PC.CF_z.adaptiveRK(_charge, _pTot, start, z_cm, stepSizeCM, hdata); - } catch (SwimZException e) { - szr = null; - //System.err.println("[WARNING] SwimZ Failed for p = " + _pTot); - - } - } - - if (szr != null) { - double bdl = szr.getBDL(PC.CF_z.getProbe()); - double pathLength = szr.getPathLength(); // already in cm - - SwimZStateVector last = szr.last(); - double p3[] = szr.getThreeMomentum(last); - - value[0] = last.x; // xf in cm - value[1] = last.y; // yz in cm - value[2] = last.z; // zf in cm - value[3] = p3[0]; - value[4] = p3[1]; - value[5] = p3[2]; - value[6] = pathLength; - value[7] = bdl / 10; // convert from kg*cm to T*cm - } else { // use old swimmer. Either low momentum or SwimZ failed. - // (dph) - traj = PC.CF.swim(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, z, accuracy, _rMax, _maxPathLength, - stepSize, cnuphys.swim.Swimmer.CLAS_Tolerance, hdata); - if(traj==null) - return null; - traj.computeBDL(PC.CP); - // traj.computeBDL(compositeField); - - double lastY[] = traj.lastElement(); - - value[0] = lastY[0] * 100; // convert back to cm - value[1] = lastY[1] * 100; // convert back to cm - value[2] = lastY[2] * 100; // convert back to cm - value[3] = lastY[3] * _pTot; - value[4] = lastY[4] * _pTot; - value[5] = lastY[5] * _pTot; - value[6] = lastY[6] * 100; - value[7] = lastY[7] * 10; - } // old swimmer - + double z_m = z_cm / CM_PER_M; + SwimTrajectory traj = PC.RCF.sectorSwim( + sector, _charge, _x0, _y0, _z0, _pTot, _theta, _phi, + z_m, accuracy, _rMax, _maxPathLength, stepSize, + cnuphys.swim.Swimmer.CLAS_Tolerance, hdata); + if (traj == null) return null; + traj.sectorComputeBDL(sector, PC.RCP); + fillFromTrajectory(out, traj); + return out; } catch (RungeKuttaException e) { - e.printStackTrace(); + return null; } - return value; - - } - - private void checkR(double _x0, double _y0, double _z0) { - this.SwimUnPhys=false; - if(Math.sqrt(_x0*_x0 + _y0*_y0)>this._rMax || - Math.sqrt(_x0*_x0 + _y0*_y0 + _z0*_z0)>this._maxPathLength) - this.SwimUnPhys=true; } - /** - * Cylindrical stopper - */ - private class CylindricalBoundarySwimStopper implements IStopper { - - private double _finalPathLength = Double.NaN; - private double _Rad; - //boolean cutOff = false; - // check that the track can get to R. Stops at the track radius - //float[] b = new float[3]; - //double x0 = _x0*100; - //double y0 = _y0*100; - //double z0 = _z0*100; - - double max = -1.0; - /** - * A swim stopper that will stop if the boundary of a plane is crossed - * - * @param maxR - * the max radial coordinate in meters. - */ - private CylindricalBoundarySwimStopper(double Rad) { - // DC reconstruction units are cm. Swim units are m. Hence scale by - // 100 - _Rad = Rad; - // check if the track will reach the surface of the cylinder. - //BfieldLab(x0, y0, z0, b); - //double trkR = _pTot*Math.sin(Math.toRadians(_theta))/Math.abs(b[2]*LIGHTVEL); - //double trkXc = x0 + trkR * Math.sin(Math.toRadians(_phi)); - //if(trkR<(Rad+trkXc) && Math.sqrt(x0*x0+y0*y0)<_Rad) { // check only for swimming inside out. - // cutOff=true; - //} - } + @Override + public double[] SwimToPlaneLab(double z_cm) { + if (isInvalid()) return null; - @Override - public boolean stopIntegration(double t, double[] y) { - - double r = Math.sqrt(y[0] * y[0] + y[1] * y[1]) * 100.; -// if(r>max ) -// max = r; -// else System.out.println(r + " " + max + " " + t); - //if(cutOff) { - return (r < max || r > _Rad); // stop intergration at closest distance to the cylinder - //} - //else { - // return (r > _Rad); - //} - } + double[] hdata = new double[3]; + double[] out = newState(); - /** - * Get the final path length in meters - * - * @return the final path length in meters - */ - @Override - public double getFinalT() { - return _finalPathLength; + SwimZResult szr = tryZSwimLab(z_cm, hdata); + if (szr != null) { + double bdl = szr.getBDL(PC.CF_z.getProbe()); + fillFromZResult(out, szr, bdl); + return out; } - /** - * Set the final path length in meters - * - * @param finalPathLength - * the final path length in meters - */ - @Override - public void setFinalT(double finalPathLength) { - _finalPathLength = finalPathLength; + try { + double z_m = z_cm / CM_PER_M; + SwimTrajectory traj = PC.CF.swim( + _charge, _x0, _y0, _z0, _pTot, _theta, _phi, + z_m, accuracy, _rMax, _maxPathLength, stepSize, + cnuphys.swim.Swimmer.CLAS_Tolerance, hdata); + if (traj == null) return null; + traj.computeBDL(PC.CP); + fillFromTrajectory(out, traj); + return out; + } catch (RungeKuttaException e) { + return null; } } - //private final double LIGHTVEL = 0.000299792458 ; - - /** - * - * @param Rad - * @return state x,y,z,px,py,pz, pathlength, iBdl at the surface - */ - public double[] SwimToCylinder(double Rad) { - - double[] value = new double[8]; - if(this.SwimUnPhys) - return null; - - CylindricalBoundarySwimStopper stopper = new CylindricalBoundarySwimStopper(Rad); - - SwimTrajectory st = PC.CF.swim(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, stopper, _maxPathLength, stepSize, - 0.0005); - if(st==null) - return null; - st.computeBDL(PC.CP); - // st.computeBDL(compositeField); - double[] lastY = st.lastElement(); - - value[0] = lastY[0] * 100; // convert back to cm - value[1] = lastY[1] * 100; // convert back to cm - value[2] = lastY[2] * 100; // convert back to cm - value[3] = lastY[3] * _pTot; // normalized values - value[4] = lastY[4] * _pTot; - value[5] = lastY[5] * _pTot; - value[6] = lastY[6] * 100; - value[7] = lastY[7] * 10; // Conversion from kG.m to T.cm - - return value; + // ------------------------------------------------------------------------ + // Geometry-based stoppers + // ------------------------------------------------------------------------ +// @Override + public double[] SwimToCylinder(double radius) { + if (SwimUnPhys) return null; + return swimWithStopper(new CylindricalBoundarySwimStopper(radius)); } - /** - * - * @param radius in cm - * @return state x,y,z,px,py,pz, pathlength, iBdl at the surface - */ - public double[] SwimRho(double radius) { - return SwimRho(radius, accuracy*100); + @Override + public double[] SwimToSphere(double radius) { + if (SwimUnPhys) return null; + return swimWithStopper(new SphericalBoundarySwimStopper(radius)); } - - /** - * - * @param radius in cm - * @param accuracy in cm - * @return state x,y,z,px,py,pz, pathlength, iBdl at the surface - */ - public double[] SwimRho(double radius, double accuracy) { - double[] value = null; + @Override + public double[] SwimToBeamLine(double xB, double yB) { + if (SwimUnPhys) return null; + return swimWithStopper(new BeamLineSwimStopper(xB, yB)); + } - // using adaptive stepsize - if(this.SwimUnPhys) - return null; + @Override + public double[] SwimToLine(Line3D l) { + if (SwimUnPhys) return null; + return swimWithStopper(new LineSwimStopper(l)); + } + private double[] swimWithStopper(IStopper stopper) { try { - - AdaptiveSwimResult result = new AdaptiveSwimResult(false); - - PC.CF.swimRho(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, radius/100, accuracy/100, _rMax, stepSize, cnuphys.swim.Swimmer.CLAS_Tolerance, result); - - if(result.getStatus()==0) { - value = new double[8]; - value[0] = result.getUf()[0] * 100; // convert back to cm - value[1] = result.getUf()[1] * 100; // convert back to cm - value[2] = result.getUf()[2] * 100; // convert back to cm - value[3] = result.getUf()[3] * _pTot; // normalized values - value[4] = result.getUf()[4] * _pTot; - value[5] = result.getUf()[5] * _pTot; - value[6] = result.getFinalS() * 100; - value[7] = 0; // Conversion from kG.m to T.cm - } - - } catch (RungeKuttaException e) { - System.out.println(_charge + " " + _x0 + " " + _y0 + " " + _z0 + " " + _pTot + " " + _theta + " " + _phi); - e.printStackTrace(); + SwimTrajectory traj = PC.CF.swim( + _charge, _x0, _y0, _z0, _pTot, _theta, _phi, + stopper, _maxPathLength, stepSize, 0.0005); + if (traj == null) return null; + traj.computeBDL(PC.CP); + double[] out = newState(); + fillFromTrajectory(out, traj); + return out; + } catch (Exception e) { + e.printStackTrace(); + return null; } - return value; - - } - - /** - * - * @param axisPoint1 in cm - * @param axisPoint2 in cm - * @param radius in cm - * @return swam trajectory to the cylinder - */ - public double[] SwimGenCylinder(Point3D axisPoint1, Point3D axisPoint2, double radius) { - return SwimGenCylinder(axisPoint1, axisPoint2, radius, accuracy*100); } - - /** - * - * @param axisPoint1 in cm - * @param axisPoint2 in cm - * @param radius in cm - * @param accuracy in cm - * @return swam trajectory to the cylinder - */ - public double[] SwimGenCylinder(Point3D axisPoint1, Point3D axisPoint2, double radius, double accuracy) { - double[] value = null; - double[] p1 = new double[3]; - double[] p2 = new double[3]; - p1[0] = axisPoint1.x()/100; - p1[1] = axisPoint1.y()/100; - p1[2] = axisPoint1.z()/100; - p2[0] = axisPoint2.x()/100; - p2[1] = axisPoint2.y()/100; - p2[2] = axisPoint2.z()/100; - - Cylinder targCyl = new Cylinder(p1, p2, radius/100); - // using adaptive stepsize - if(this.SwimUnPhys) - return null; + // ------------------------------------------------------------------------ + // Adaptive swimmers + // ------------------------------------------------------------------------ + @Override + public double[] SwimRho(double radius, double accuracy) { + if (SwimUnPhys) return null; try { - - AdaptiveSwimResult result = new AdaptiveSwimResult(false); - - PC.CF.swimCylinder(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, - p1, p2, radius/100, accuracy/100, _rMax, stepSize, cnuphys.swim.Swimmer.CLAS_Tolerance, result); - - if(result.getStatus()==0) { - value = new double[8]; - value[0] = result.getUf()[0] * 100; // convert back to cm - value[1] = result.getUf()[1] * 100; // convert back to cm - value[2] = result.getUf()[2] * 100; // convert back to cm - value[3] = result.getUf()[3] * _pTot; // normalized values - value[4] = result.getUf()[4] * _pTot; - value[5] = result.getUf()[5] * _pTot; - value[6] = result.getFinalS() * 100; - value[7] = 0; // Conversion from kG.m to T.cm - } - + AdaptiveSwimResult r = new AdaptiveSwimResult(false); + PC.CF.swimRho(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, + radius / CM_PER_M, accuracy / CM_PER_M, + _rMax, stepSize, cnuphys.swim.Swimmer.CLAS_Tolerance, r); + return r.getStatus() == 0 ? adaptiveOut(r) : null; } catch (RungeKuttaException e) { - System.out.println(_charge + " " + _x0 + " " + _y0 + " " + _z0 + " " + _pTot + " " + _theta + " " + _phi); - e.printStackTrace(); + e.printStackTrace(); + return null; } - return value; - } - public double[] SwimPlane(Vector3D n, Point3D p, double accuracy) { - - double[] value = null; - - - // using adaptive stepsize - if(this.SwimUnPhys) - return null; - + @Override + public double[] SwimGenCylinder(Point3D a1, Point3D a2, double radius, double accuracy) { + if (SwimUnPhys) return null; try { - - AdaptiveSwimResult result = new AdaptiveSwimResult(false); - - PC.CF.swimPlane(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, - n.x(),n.y(),n.z(),p.x()/100,p.y()/100,p.z()/100, - accuracy/100, _rMax, stepSize, cnuphys.swim.Swimmer.CLAS_Tolerance, result); - - - if(result.getStatus()==0) { - value = new double[8]; - value[0] = result.getUf()[0] * 100; // convert back to cm - value[1] = result.getUf()[1] * 100; // convert back to cm - value[2] = result.getUf()[2] * 100; // convert back to cm - value[3] = result.getUf()[3] * _pTot; // normalized values - value[4] = result.getUf()[4] * _pTot; - value[5] = result.getUf()[5] * _pTot; - value[6] = result.getFinalS() * 100; - value[7] = 0; // Conversion from kG.m to T.cm - } - + AdaptiveSwimResult r = new AdaptiveSwimResult(false); + double[] p1 = {a1.x() / CM_PER_M, a1.y() / CM_PER_M, a1.z() / CM_PER_M}; + double[] p2 = {a2.x() / CM_PER_M, a2.y() / CM_PER_M, a2.z() / CM_PER_M}; + PC.CF.swimCylinder(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, + p1, p2, radius / CM_PER_M, accuracy / CM_PER_M, + _rMax, stepSize, cnuphys.swim.Swimmer.CLAS_Tolerance, r); + return r.getStatus() == 0 ? adaptiveOut(r) : null; } catch (RungeKuttaException e) { - System.out.println(_charge + " " + _x0 + " " + _y0 + " " + _z0 + " " + _pTot + " " + _theta + " " + _phi); - e.printStackTrace(); - } - return value; - - } - - - private class SphericalBoundarySwimStopper implements IStopper { - - private double _finalPathLength = Double.NaN; - - private double _Rad; - - /** - * A swim stopper that will stop if the boundary of a plane is crossed - * - * @param maxR - * the max radial coordinate in meters. - */ - private SphericalBoundarySwimStopper(double Rad) { - // DC reconstruction units are cm. Swim units are m. Hence scale by - // 100 - _Rad = Rad; - } - - @Override - public boolean stopIntegration(double t, double[] y) { - - double r = Math.sqrt(y[0] * y[0] + y[1] * y[1] + y[2] * y[2]) * 100.; - - return (r > _Rad); - - } - - /** - * Get the final path length in meters - * - * @return the final path length in meters - */ - @Override - public double getFinalT() { - return _finalPathLength; - } - - /** - * Set the final path length in meters - * - * @param finalPathLength - * the final path length in meters - */ - @Override - public void setFinalT(double finalPathLength) { - _finalPathLength = finalPathLength; - } - } - /** - * - * @param Rad - * @return state x,y,z,px,py,pz, pathlength, iBdl at the surface - */ - public double[] SwimToSphere(double Rad) { - - double[] value = new double[8]; - // using adaptive stepsize - if(this.SwimUnPhys==true) - return null; - SphericalBoundarySwimStopper stopper = new SphericalBoundarySwimStopper(Rad); - - SwimTrajectory st = PC.CF.swim(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, stopper, _maxPathLength, stepSize, - 0.0005); - if(st==null) + e.printStackTrace(); return null; - st.computeBDL(PC.CP); - // st.computeBDL(compositeField); - - double[] lastY = st.lastElement(); - - value[0] = lastY[0] * 100; // convert back to cm - value[1] = lastY[1] * 100; // convert back to cm - value[2] = lastY[2] * 100; // convert back to cm - value[3] = lastY[3] * _pTot; // normalized values - value[4] = lastY[4] * _pTot; - value[5] = lastY[5] * _pTot; - value[6] = lastY[6] * 100; - value[7] = lastY[7] * 10; // Conversion from kG.m to T.cm - - return value; - - } - - // added for swimming to outer detectors - private class PlaneBoundarySwimStopper implements IStopper { - - private double _finalPathLength = Double.NaN; - private double _d; - private Vector3D _n; - private double _dist2plane; - private int _dir; - - /** - * A swim stopper that will stop if the boundary of a plane is crossed - * - * @param maxR - * the max radial coordinate in meters. - */ - private PlaneBoundarySwimStopper(double d, Vector3D n, int dir) { - // DC reconstruction units are cm. Swim units are m. Hence scale by - // 100 - _d = d; - _n = n; - _dir = dir; - } - - @Override - public boolean stopIntegration(double t, double[] y) { - double dtrk = y[0] * _n.x() + y[1] * _n.y() + y[2] * _n.z(); - - double accuracy = 20e-6; // 20 microns - // System.out.println(" dist "+dtrk*100+ " state "+y[0]*100+", - // "+y[1]*100+" , "+y[2]*100); - if (_dir < 0) { - return dtrk < _d; - } else { - return dtrk > _d; - } - - } - - @Override - public double getFinalT() { - - return _finalPathLength; - } - - /** - * Set the final path length in meters - * - * @param finalPathLength - * the final path length in meters - */ - @Override - public void setFinalT(double finalPathLength) { - _finalPathLength = finalPathLength; } } - /** - * - * @param d_cm - * @param n - * @param dir - * @return return state x,y,z,px,py,pz, pathlength, iBdl at the plane surface in the lab frame - */ - public double[] SwimToPlaneBoundary(double d_cm, Vector3D n, int dir) { - - double[] value = new double[8]; - if(this.SwimUnPhys) - return null; - double d = d_cm / 100; - - double hdata[] = new double[3]; - // using adaptive stepsize - // the new swim to plane in swimmer - Plane plane = new Plane(n.x(), n.y(), n.z(), d); - SwimTrajectory st; + @Override + public double[] SwimPlane(Vector3D n, Point3D p, double accuracy) { + if (SwimUnPhys) return null; try { - - st = PC.CF.swim(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, plane, accuracy, _maxPathLength, stepSize, - cnuphys.swim.Swimmer.CLAS_Tolerance, hdata); - - st.computeBDL(PC.CP); - - double[] lastY = st.lastElement(); - - value[0] = lastY[0] * 100; // convert back to cm - value[1] = lastY[1] * 100; // convert back to cm - value[2] = lastY[2] * 100; // convert back to cm - value[3] = lastY[3] * _pTot; // normalized values - value[4] = lastY[4] * _pTot; - value[5] = lastY[5] * _pTot; - value[6] = lastY[6] * 100; - value[7] = lastY[7] * 10; // Conversion from kG.m to T.cm - - // System.out.println("\nCOMPARE plane swims DIRECTION = " + - // dir); - // for (int i = 0; i < 8; i++) { - // System.out.print(String.format("%-8.5f ", value[i])); - // } - - + AdaptiveSwimResult r = new AdaptiveSwimResult(false); + PC.CF.swimPlane(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, + n.x(), n.y(), n.z(), + p.x() / CM_PER_M, p.y() / CM_PER_M, p.z() / CM_PER_M, + accuracy / CM_PER_M, _rMax, stepSize, + cnuphys.swim.Swimmer.CLAS_Tolerance, r); + return r.getStatus() == 0 ? adaptiveOut(r) : null; } catch (RungeKuttaException e) { - e.printStackTrace(); - } - return value; - - } - - - - private class BeamLineSwimStopper implements IStopper { - - private double _finalPathLength = Double.NaN; - - private double _xB; - private double _yB; - double min = Double.POSITIVE_INFINITY; - double thetaRad = Math.toRadians(_theta); - double phiRad = Math.toRadians(_phi); - double pz = _pTot * Math.cos(thetaRad); - private BeamLineSwimStopper(double xB, double yB) { - // DC reconstruction units are cm. Swim units are m. Hence scale by - // 100 - _xB = xB; - _yB = yB; - } - - @Override - public boolean stopIntegration(double t, double[] y) { - - double r = Math.sqrt((_xB-y[0]* 100.) * (_xB-y[0]* 100.) + (_yB-y[1]* 100.) * (_yB-y[1]* 100.)); - if(r min ); - - } - - /** - * Get the final path length in meters - * - * @return the final path length in meters - */ - @Override - public double getFinalT() { - return _finalPathLength; - } - - /** - * Set the final path length in meters - * - * @param finalPathLength - * the final path length in meters - */ - @Override - public void setFinalT(double finalPathLength) { - _finalPathLength = finalPathLength; - } - } - - public double[] SwimToBeamLine(double xB, double yB) { - - double[] value = new double[8]; - - if(this.SwimUnPhys==true) - return null; - BeamLineSwimStopper stopper = new BeamLineSwimStopper(xB, yB); - - SwimTrajectory st = PC.CF.swim(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, stopper, _maxPathLength, stepSize, - 0.0005); - if(st==null) + e.printStackTrace(); return null; - st.computeBDL(PC.CP); - // st.computeBDL(compositeField); - - double[] lastY = st.lastElement(); - - value[0] = lastY[0] * 100; // convert back to cm - value[1] = lastY[1] * 100; // convert back to cm - value[2] = lastY[2] * 100; // convert back to cm - value[3] = lastY[3] * _pTot; // normalized values - value[4] = lastY[4] * _pTot; - value[5] = lastY[5] * _pTot; - value[6] = lastY[6] * 100; - value[7] = lastY[7] * 10; // Conversion from kG.m to T.cm - - return value; - - } - - // - - private class LineSwimStopper implements IStopper { - - private double _finalPathLength = Double.NaN; - - private Line3D _l; - private Point3D _p; - double min = 999; - private LineSwimStopper(Line3D l) { - // DC reconstruction units are cm. Swim units are m. Hence scale by - // 100 - _l =l; - _p = new Point3D(); - } - - @Override - public boolean stopIntegration(double t, double[] y) { - _p.set(y[0]* 100.0, y[1]* 100.0, y[2]* 100.0); - double doca = _l.distance(_p).length(); - if(doca min ); - - } - - /** - * Get the final path length in meters - * - * @return the final path length in meters - */ - @Override - public double getFinalT() { - return _finalPathLength; } - - /** - * Set the final path length in meters - * - * @param finalPathLength - * the final path length in meters - */ - @Override - public void setFinalT(double finalPathLength) { - _finalPathLength = finalPathLength; - } - } - - public double[] SwimToLine(Line3D l) { - - double[] value = new double[8]; - - if(this.SwimUnPhys==true) - return null; - LineSwimStopper stopper = new LineSwimStopper(l); - - SwimTrajectory st = PC.CF.swim(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, stopper, _maxPathLength, stepSize, - 0.0005); - if(st==null) - return null; - st.computeBDL(PC.CP); - // st.computeBDL(compositeField); - - double[] lastY = st.lastElement(); - - value[0] = lastY[0] * 100; // convert back to cm - value[1] = lastY[1] * 100; // convert back to cm - value[2] = lastY[2] * 100; // convert back to cm - value[3] = lastY[3] * _pTot; // normalized values - value[4] = lastY[4] * _pTot; - value[5] = lastY[5] * _pTot; - value[6] = lastY[6] * 100; - value[7] = lastY[7] * 10; // Conversion from kG.m to T.cm - - return value; - } - // - - - private void printV(String pfx, double v[]) { - double x = v[0] / 100; - double y = v[1] / 100; - double z = v[2] / 100; - double r = Math.sqrt(x * x + y * y + z * z); - System.out.println(String.format("%s: (%-8.5f, %-8.5f, %-8.5f) R: %-8.5f", pfx, z, y, z, r)); + private double[] adaptiveOut(AdaptiveSwimResult r) { + double[] out = newState(); + out[0] = r.getUf()[0] * CM_PER_M; + out[1] = r.getUf()[1] * CM_PER_M; + out[2] = r.getUf()[2] * CM_PER_M; + out[3] = r.getUf()[3] * _pTot; + out[4] = r.getUf()[4] * _pTot; + out[5] = r.getUf()[5] * _pTot; + out[6] = r.getFinalS() * CM_PER_M; + out[7] = 0.0; + return out; } - /** - * - * @param sector - * @param x_cm - * @param y_cm - * @param z_cm - * @param result B field components in T in the tilted sector system - */ - public void Bfield(int sector, double x_cm, double y_cm, double z_cm, float[] result) { - - PC.RCP.field(sector, (float) x_cm, (float) y_cm, (float) z_cm, result); - // rcompositeField.field((float) x_cm, (float) y_cm, (float) z_cm, - // result); - result[0] = result[0] / 10; - result[1] = result[1] / 10; - result[2] = result[2] / 10; - - } - /** - * - * @param x_cm - * @param y_cm - * @param z_cm - * @param result B field components in T in the lab frame - */ - public void BfieldLab(double x_cm, double y_cm, double z_cm, float[] result) { - - PC.CP.field((float) x_cm, (float) y_cm, (float) z_cm, result); - result[0] = result[0] / 10; - result[1] = result[1] / 10; - result[2] = result[2] / 10; - - } - - - - public double[] AdaptiveSwimPlane(double px, double py, double pz, double nx, double ny, double nz, double accuracy) { -// System.out.println("Don't use yet"); - - double[] value = new double[8]; - - Vector norm = new Vector(nx,ny,nz); - Point point = new Point(px/100,py/100,pz/100); - - cnuphys.adaptiveSwim.geometry.Plane targetPlane = new cnuphys.adaptiveSwim.geometry.Plane(norm, point); - - - // using adaptive stepsize - if(this.SwimUnPhys) - return null; + // ------------------------------------------------------------------------ + // Swim to boundaries in the lab + // ------------------------------------------------------------------------ + @Override + public double[] SwimToPlaneBoundary(double d_cm, Vector3D n) { + if (SwimUnPhys) return null; try { - - AdaptiveSwimResult result = new AdaptiveSwimResult(false); - - PC.AS.swimPlane(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, targetPlane, - accuracy/100, _rMax, stepSize, cnuphys.swim.Swimmer.getEps(), result); - - if(result.getStatus() == AdaptiveSwimmer.SWIM_SUCCESS) { - value[0] = result.getUf()[0] * 100; // convert back to cm - value[1] = result.getUf()[1] * 100; // convert back to cm - value[2] = result.getUf()[2] * 100; // convert back to cm - value[3] = result.getUf()[3] * _pTot; // normalized values - value[4] = result.getUf()[4] * _pTot; - value[5] = result.getUf()[5] * _pTot; - value[6] = result.getFinalS() * 100; - value[7] = 0; // Conversion from kG.m to T.cm - } - else { - return null; - } - - } catch (AdaptiveSwimException e) { - e.printStackTrace(); - } - return value; - - } - - - public double[] AdaptiveSwimCylinder(double a1x, double a1y, double a1z, double a2x, double a2y, double a2z, double radius, double accuracy) { - // System.out.println("Don't use yet"); - double[] value = new double[8]; - - radius = radius/100; - Point a1 = new Point(a1x/100, a1y/100, a1z/100); - Point a2 = new Point(a2x/100, a2y/100, a2z/100); - Line centerLine = new Line(a1, a2); - - cnuphys.adaptiveSwim.geometry.Cylinder targetCylinder = new cnuphys.adaptiveSwim.geometry.Cylinder(centerLine, radius); - - - // using adaptive stepsize - if(this.SwimUnPhys) - return null; - - try { - - AdaptiveSwimResult result = new AdaptiveSwimResult(false); - - PC.AS.swimCylinder(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, targetCylinder, - accuracy/100, _rMax, stepSize, cnuphys.swim.Swimmer.getEps(), result); - - if(result.getStatus() == AdaptiveSwimmer.SWIM_SUCCESS) { - value[0] = result.getUf()[0] * 100; // convert back to cm - value[1] = result.getUf()[1] * 100; // convert back to cm - value[2] = result.getUf()[2] * 100; // convert back to cm - value[3] = result.getUf()[3] * _pTot; // normalized values - value[4] = result.getUf()[4] * _pTot; - value[5] = result.getUf()[5] * _pTot; - value[6] = result.getFinalS() * 100; - value[7] = 0; // Conversion from kG.m to T.cm - } - else { - return null; - } - - } catch (AdaptiveSwimException e) { - e.printStackTrace(); - } - return value; - - } - - public double[] AdaptiveSwimRho(double radius, double accuracy) { - System.out.println("Don't use yet"); - - double[] value = new double[8]; - - radius = radius/100; - // using adaptive stepsize - if(this.SwimUnPhys) + Plane plane = new Plane(n.x(), n.y(), n.z(), d_cm / CM_PER_M); + SwimTrajectory traj = PC.CF.swim( + _charge, _x0, _y0, _z0, _pTot, _theta, _phi, + plane, accuracy, _maxPathLength, stepSize, + cnuphys.swim.Swimmer.CLAS_Tolerance, new double[3]); + if (traj == null) return null; + traj.computeBDL(PC.CP); + double[] out = newState(); + fillFromTrajectory(out, traj); + return out; + } catch (RungeKuttaException e) { + e.printStackTrace(); return null; - - try { - - AdaptiveSwimResult result = new AdaptiveSwimResult(false); - - PC.AS.swimRho(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, radius, - accuracy/100, _rMax, stepSize, cnuphys.swim.Swimmer.getEps(), result); - - if(result.getStatus() == AdaptiveSwimmer.SWIM_SUCCESS) { - value[0] = result.getUf()[0] * 100; // convert back to cm - value[1] = result.getUf()[1] * 100; // convert back to cm - value[2] = result.getUf()[2] * 100; // convert back to cm - value[3] = result.getUf()[3] * _pTot; // normalized values - value[4] = result.getUf()[4] * _pTot; - value[5] = result.getUf()[5] * _pTot; - value[6] = result.getFinalS() * 100; - value[7] = 0; // Conversion from kG.m to T.cm - } - else { - return null; - } - - } catch (AdaptiveSwimException e) { - e.printStackTrace(); } - return value; - } -/** - * - * @param Z - * @return state x,y,z,px,py,pz, pathlength, iBdl at the surface - */ + @Override public double[] SwimToZ(double Z, int dir) { - - double[] value = new double[8]; - //if(this.SwimUnPhys) - // return null; - + if (SwimUnPhys) return null; ZSwimStopper stopper = new ZSwimStopper(Z, dir); - - SwimTrajectory st = PC.CF.swim(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, stopper, _maxPathLength, stepSize, - distanceBetweenSaves); - if(st==null) - return null; - st.computeBDL(PC.CP); - // st.computeBDL(compositeField); - this.setSwimTraj(st); - double[] lastY = st.lastElement(); - - value[0] = lastY[0] * 100; // convert back to cm - value[1] = lastY[1] * 100; // convert back to cm - value[2] = lastY[2] * 100; // convert back to cm - value[3] = lastY[3] * _pTot; // normalized values - value[4] = lastY[4] * _pTot; - value[5] = lastY[5] * _pTot; - value[6] = lastY[6] * 100; - value[7] = lastY[7] * 10; // Conversion from kG.m to T.cm - - return value; - - } - - private SwimTrajectory swimTraj; - - private class ZSwimStopper implements IStopper { - - private double _finalPathLength = Double.NaN; - - private double _Z; - private int _dir; - - private ZSwimStopper(double Z, int dir) { - // The reconstruction units are cm. Swim units are m. Hence scale by - // 100 - _Z = Z; - _dir = dir; - } - - @Override - public boolean stopIntegration(double t, double[] y) { - - double z = y[2] * 100.; - if(_dir>0) { - return (z > _Z); - } else { - return (z<_Z); - } - } - - /** - * Get the final path length in meters - * - * @return the final path length in meters - */ - @Override - public double getFinalT() { - return _finalPathLength; - } - - /** - * Set the final path length in meters - * - * @param finalPathLength - * the final path length in meters - */ - @Override - public void setFinalT(double finalPathLength) { - _finalPathLength = finalPathLength; + try { + SwimTrajectory traj = PC.CF.swim( + _charge, _x0, _y0, _z0, _pTot, _theta, _phi, + stopper, _maxPathLength, stepSize, distanceBetweenSaves); + if (traj == null) return null; + traj.computeBDL(PC.CP); + setSwimTraj(traj); + double[] out = newState(); + fillFromTrajectory(out, traj); + return out; + } catch (Exception e) { + e.printStackTrace(); + return null; } } - /** - * @return the swimTraj - */ - public SwimTrajectory getSwimTraj() { - return swimTraj; - } - - /** - * @param swimTraj the swimTraj to set - */ - public void setSwimTraj(SwimTrajectory swimTraj) { - this.swimTraj = swimTraj; - } - - private class DCASwimStopper implements IStopper { - - public DCASwimStopper(SwimTrajectory swimTraj) { - _swimTraj = swimTraj; - for(int i = 0; i < _swimTraj.size()-1; i++) { - polylines.add(new Line3D(_swimTraj.get(i)[0],_swimTraj.get(i)[1],_swimTraj.get(i)[2], - _swimTraj.get(i+1)[0],_swimTraj.get(i+1)[1],_swimTraj.get(i+1)[2])); - - } - } - - private List polylines = new ArrayList<>(); - private SwimTrajectory _swimTraj; - private double _finalPathLength = Double.NaN; - private double _doca = Double.POSITIVE_INFINITY; - - @Override - public boolean stopIntegration(double t, double[] y) { - - Point3D dcaCand = new Point3D(y[0],y[1],y[2]); - double maxDoca = Double.POSITIVE_INFINITY; - - for(Line3D l : polylines) { - if(l.distance(dcaCand).length()this._rMax || + Math.sqrt(_x0*_x0 + _y0*_y0 + _z0*_z0)>this._maxPathLength) + this.SwimUnPhys=true; + } + +} diff --git a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swimmer.java b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swimmer.java index 7915ce6129..c6344fc593 100644 --- a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swimmer.java +++ b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swimmer.java @@ -1,16 +1,14 @@ package org.jlab.clas.swimtools; -import cnuphys.magfield.MagneticFields; import java.util.HashMap; -import java.util.logging.Level; import java.util.logging.Logger; +import cnuphys.magfield.MagneticFields; /** * * @author ziegler, heddle */ - public class Swimmer { public static Logger LOGGER = Logger.getLogger(Swimmer.class.getName()); diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java index 6f68d554ab..88c4e488b5 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java @@ -526,7 +526,7 @@ public void setTrackPars(Track cand, double theta_n = ((double) (sector - 1)) * Math.toRadians(60.); double x_n = Math.cos(theta_n); double y_n = Math.sin(theta_n); - double[] Vt = dcSwim.SwimToPlaneBoundary(0, new Vector3D(x_n, y_n, 0), -1); + double[] Vt = dcSwim.SwimToPlaneBoundary(0, new Vector3D(x_n, y_n, 0)); if (Vt == null) { return; @@ -672,7 +672,7 @@ public void setTrackPars(Track cand, //double x_n = Math.cos(theta_n); //double y_n = Math.sin(theta_n); //double d = x_n*xB + y_n*yB; - //Vt = dcSwim.SwimToPlaneBoundary(d, new Vector3D(x_n, y_n, 0), -1); + //Vt = dcSwim.SwimToPlaneBoundary(d, new Vector3D(x_n, y_n, 0)); //if(Vt==null) // return; double xOrFix = Vt[0]; diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/Trajectory.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/Trajectory.java index ddb69ebf94..b56c2671f2 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/Trajectory.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/Trajectory.java @@ -218,19 +218,19 @@ public void calcTrajectory(int trackId, Swim dcSwim, Point3D v, Vector3D p, int // swim backward from HTCC to Target if(surface.getDetectorType() == DetectorType.TARGET) { - int dir = -1; + //int dir = -1; float b[] = new float[3]; dcSwim.BfieldLab(v.x(), v.y(), v.z(), b); dcSwim.SetSwimParameters(last.getPoint().x(), last.getPoint().y(), last.getPoint().z(), -last.getMomentum().x(), -last.getMomentum().y(), -last.getMomentum().z(), -q); - double[] tPars = dcSwim.SwimToPlaneBoundary(surface.getD(), surface.getNormal(), dir); + double[] tPars = dcSwim.SwimToPlaneBoundary(surface.getD(), surface.getNormal()); if(tPars==null || surface.distanceFromPlane(tPars[0], tPars[1], tPars[2])>TOLERANCE) return; this.addTrajectoryPoint(tPars[0], tPars[1], tPars[2], -tPars[3], -tPars[4], -tPars[5], last.getPath()-tPars[6], last.getiBdl()-tPars[7], surface); } // swim forward from vertex to FMT and from HTCC to FD planes else { - int dir = 1; + //int dir = 1; double path = 0; double bdl = 0; if(surface.getDetectorType() == DetectorType.FMT) { @@ -242,7 +242,7 @@ public void calcTrajectory(int trackId, Swim dcSwim, Point3D v, Vector3D p, int path = htccPars[6]; bdl = htccPars[7]; } - double[] tPars = dcSwim.SwimToPlaneBoundary(surface.getD(), surface.getNormal(), dir); + double[] tPars = dcSwim.SwimToPlaneBoundary(surface.getD(), surface.getNormal()); if(tPars==null) return; double MAXDIST = -99; if(surface.getDetectorType() == DetectorType.RICH) MAXDIST = -1; diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java index 4b6bb58d2a..f6819ce356 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java @@ -444,7 +444,7 @@ private static double[] getTrajectory(double x, double y, double z, double px, d double d = p.dot(n); if(v.dot(n)