Skip to content
Open
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ public abstract class AKFitter {

// return variables
public boolean setFitFailed = false;
public double chi2;
public double chi2 = Double.POSITIVE_INFINITY;
public int NDF;
public int NDF0;
public int numIter;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ public class Surface implements Comparable<Surface> {
public double swimAccuracy;
public boolean passive = false;
public double hemisphere = 1;
public Point3D rayInterc = new Point3D(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);

public double[] unc = new double[2];
public double[] doca = new double[2];
Expand Down Expand Up @@ -421,11 +422,18 @@ public void setTransformation(Transformation3D transform) {

@Override
public int compareTo(Surface o) {
if (this.index > o.index) {
return 1;
if(this.rayInterc.y()==Double.POSITIVE_INFINITY) {
if (this.index > o.index) {
return 1;
} else {
return -1;
}
} else {
return -1;
}
if (this.rayInterc.y() < o.rayInterc.y()) {
return 1;
} else {
return -1;
}
}
}

}
Original file line number Diff line number Diff line change
@@ -1,13 +1,14 @@
package org.jlab.clas.tracking.kalmanfilter.helical;

import java.util.ArrayList;
import java.util.List;
import org.jlab.clas.swimtools.Swim;
import org.jlab.clas.tracking.kalmanfilter.AMeasVecs;
import org.jlab.clas.tracking.kalmanfilter.AMeasVecs.MeasVec;
import org.jlab.clas.tracking.kalmanfilter.AStateVecs;
import org.jlab.clas.tracking.kalmanfilter.Surface;
import org.jlab.clas.tracking.kalmanfilter.Units;
import org.jlab.clas.tracking.trackrep.Helix;
import org.jlab.clas.tracking.utilities.MatrixOps;
import org.jlab.geom.prim.Line3D;
import org.jlab.geom.prim.Point3D;
import org.jlab.geom.prim.Vector3D;
Expand Down Expand Up @@ -42,9 +43,8 @@ public boolean getStateVecPosAtMeasSite(StateVec sv, AMeasVecs.MeasVec mv, Swim

Point3D st = new Point3D(sv.x, sv.y, sv.z);
Vector3D stu = new Vector3D(sv.px,sv.py,sv.pz).asUnit();

Line3D toPln = new Line3D(st, stu);
if(mv.surface.plane!=null) {
Line3D toPln = new Line3D(st, stu);
Point3D inters = new Point3D();
int ints = mv.surface.plane.intersection(toPln, inters);
sv.x = inters.x();
Expand All @@ -53,22 +53,16 @@ public boolean getStateVecPosAtMeasSite(StateVec sv, AMeasVecs.MeasVec mv, Swim
sv.path = inters.distance(st);
}
else if(mv.surface.cylinder!=null) {
mv.surface.toLocal().apply(st);
mv.surface.toLocal().apply(stu);
double r = mv.surface.cylinder.baseArc().radius();
double delta = Math.sqrt((st.x()*stu.x()+st.y()*stu.y())*(st.x()*stu.x()+st.y()*stu.y())-(-r*r+st.x()*st.x()+st.y()*st.y())*(stu.x()*stu.x()+stu.y()*stu.y()));
double l = (-(st.x()*stu.x()+st.y()*stu.y())+delta)/(stu.x()*stu.x()+stu.y()*stu.y());
if(Math.signum(st.y()+l*stu.y())!=mv.hemisphere) {
l = (-(st.x()*stu.x()+st.y()*stu.y())-delta)/(stu.x()*stu.x()+stu.y()*stu.y());
}
Point3D inters = new Point3D(st.x()+l*stu.x(),st.y()+l*stu.y(),st.z()+l*stu.z());
mv.surface.toGlobal().apply(inters);
// RDV: should switch to use clas-geometry intersection method, not done now to alwys return a value
sv.x = inters.x();
sv.y = inters.y();
sv.z = inters.z();
sv.path = inters.distance(st);
}
List<Point3D> inters = new ArrayList();
int ints = mv.surface.cylinder.intersectionRay(toPln, inters);
if(ints<1) return false;

sv.x = inters.get(0).x() ;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is the first intersection always the right one?

sv.y = inters.get(0).y() ;
sv.z = inters.get(0).z() ;
sv.path = inters.get(0).distance(st);

}
} else {
if(swim==null) { // applicable only to planes parallel to the z -axis
Helix helix = sv.getHelix(xref, yref);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -54,22 +54,25 @@ public void runFitter(AStateVecs sv, AMeasVecs mv) {
double newchisq = this.calc_chi2(sv, mv);
// if curvature is 0, fit failed
if(Double.isNaN(newchisq) || sv.smoothed().get(0)==null) {
this.setFitFailed = true;
this.setFitFailed = true;
break;
}
else if(newchisq < this.chi2) {
this.chi2 = newchisq;
finalStateVec = sv.new StateVec(sv.smoothed().get(0));
this.setTrajectory(sv, mv);
setFitFailed = false;
//System.out.println(finalStateVec.toString());
//setFitFailed = false; //allow for new iteration to be worse but save the track from the best iteration.
// There are instances where the next iteration is worse at the 1 per 10000 level. This would loose the track...
}
// stop if chi2 got worse
else {
break;
}
}
if(!this.setFitFailed) {
finalStateVec = sv.new StateVec(sv.smoothed().get(0));
}
//if(!this.setFitFailed) {
// finalStateVec = sv.new StateVec(sv.smoothed().get(0));
//}
}

@Override
Expand Down
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
package org.jlab.clas.tracking.kalmanfilter.straight;

import java.util.ArrayList;
import java.util.List;
import org.jlab.clas.swimtools.Swim;
import org.jlab.clas.tracking.kalmanfilter.AMeasVecs;
import org.jlab.clas.tracking.kalmanfilter.AMeasVecs.MeasVec;
Expand All @@ -25,7 +27,7 @@ public boolean getStateVecPosAtMeasSite(StateVec vec, MeasVec mv, Swim swim) {

Point3D ref = new Point3D(vec.x0,vec.y0,vec.z0);
Vector3D u = new Vector3D(vec.tx, 1, vec.tz).asUnit();

Line3D toPln = new Line3D(ref,u);
if(mv.k==0) {
vec.x = vec.x0-vec.y0*vec.tx;
vec.y = 0;
Expand All @@ -34,7 +36,6 @@ public boolean getStateVecPosAtMeasSite(StateVec vec, MeasVec mv, Swim swim) {
}
else if(mv.hemisphere!=0) {
if(mv.surface.plane!=null) {
Line3D toPln = new Line3D(ref,u);
Point3D inters = new Point3D();
int ints = mv.surface.plane.intersection(toPln, inters);
vec.x = inters.x() ;
Expand All @@ -43,22 +44,29 @@ else if(mv.hemisphere!=0) {
vec.path = inters.distance(ref);
}
if(mv.surface.cylinder!=null) {
mv.surface.toLocal().apply(ref);
mv.surface.toLocal().apply(u);
double r = 0.5*(mv.surface.cylinder.baseArc().radius()+mv.surface.cylinder.highArc().radius());
double delta = Math.sqrt((ref.x()*u.x()+ref.y()*u.y())*(ref.x()*u.x()+ref.y()*u.y())
-(-r*r+ref.x()*ref.x()+ref.y()*ref.y())*(u.x()*u.x()+u.y()*u.y()));
double l = (-(ref.x()*u.x()+ref.y()*u.y())+delta)/(u.x()*u.x()+u.y()*u.y());
if(Math.signum(ref.y()+l*u.y())!=mv.hemisphere) {
l = (-(ref.x()*u.x()+ref.y()*u.y())-delta)/(u.x()*u.x()+u.y()*u.y());
List<Point3D> inters = new ArrayList();
int ints = mv.surface.cylinder.intersection(toPln, inters);
if(ints<1) {
if(!mv.surface.passive) {
return false;
} else {
mv.skip=true;
return true;
}

Point3D cylInt = new Point3D(ref.x()+l*u.x(),ref.y()+l*u.y(),ref.z()+l*u.z());
mv.surface.toGlobal().apply(cylInt);
vec.x = cylInt.x();
vec.y = cylInt.y();
vec.z = cylInt.z();
vec.path = cylInt.distance(ref);
} else {
int index =0;
if(ints>1) {//pick the intersection that is closest to the global fit ray intersection
double dh0 = Math.abs(inters.get(0).y()-mv.surface.rayInterc.y());
double dh1 = Math.abs(inters.get(1).y()-mv.surface.rayInterc.y());
if(dh1<dh0)
index=1; //if the second intersection is closest use its index
}
vec.x = inters.get(index).x() ;
vec.y = inters.get(index).y() ;
vec.z = inters.get(index).z() ;
vec.path = inters.get(index).distance(ref);
}

}
return true;
}
Expand Down Expand Up @@ -172,7 +180,6 @@ public void init(double x0, double z0, double tx, double tz, Units units, double
this.trackTrajB.clear();
this.trackTrajS.clear();
this.trackTrajT.put(0, new StateVec(initSV));

}

@Override
Expand Down
Loading
Loading