Skip to content

Correct calculated doca in KF tracking #1108

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
tongtongcao wants to merge 2 commits into
base: development
Choose a base branch
from
Open

Correct calculated doca in KF tracking #1108

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
793e70b
correct calculated doca by alpha
tongtongcao Feb 12, 2026
0209925
Merge branch 'development' into iss1107
baltzell Feb 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/AMeasVecs.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ public double[] dhDoca(int k, StateVec stateVec) {

Surface surf = this.measurements.get(stateVec.k).surface;
Point3D point = new Point3D(stateVec.x, stateVec.y, stateVec.z);
double h = hDoca(point, surf.wireLine[0]);
double h = hDoca(point, surf.wireLine[0], surf.alpha[0]);

double signMeas = 1;
double sign = 1;
Expand All @@ -66,7 +66,7 @@ public double[] dhDoca(int k, StateVec stateVec) {

//USE THE DOUBLE HIT
if(surf.doca[1]!=-99) {
h = hDoca(point, surf.wireLine[1]);
h = hDoca(point, surf.wireLine[1], surf.alpha[1]);

signMeas = Math.signum(surf.doca[1]);
sign = Math.signum(h);
Expand All @@ -78,13 +78,13 @@ public double[] dhDoca(int k, StateVec stateVec) {
}

// Return a signed doca for DC
public double hDoca(Point3D point, Line3D wireLine) {
public double hDoca(Point3D point, Line3D wireLine, double alpha) {

Line3D WL = new Line3D();
WL.copy(wireLine);
WL.copy(WL.distance(point));

return WL.length()*Math.signum(WL.direction().x());
return WL.length()*Math.signum(WL.direction().x()) * Math.cos(Math.toRadians(alpha));
}

public double dh(int k, StateVec stateVec) {
Expand Down
17 changes: 17 additions & 0 deletions common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ public class Surface implements Comparable<Surface> {
public double[] unc = new double[2];
public double[] doca = new double[2];
public Line3D[] wireLine = new Line3D[2];
public double[] alpha = {0., 0.};
public int region;
public int superlayer;
public int nMeas = 1;
Expand Down Expand Up @@ -89,6 +90,22 @@ public Surface(int region, double[] doca, double[] unc, double error, Line3D[] w
materials.add(material_argon);
}

public Surface(int region, double[] doca, double[] unc, double error, Line3D[] wireLine, double[] alpha) {
type = Type.LINEDOCA;
this.region = region;
this.doca = doca;
this.unc = unc;
this.error = error;
this.wireLine = wireLine;
this.alpha = alpha;
Point3D point = new Point3D(0, 0, wireLine[0].origin().z());
this.measPoint = point;
Material material_air = new Material("air", 0, 0, 0, 30400, 0, Units.CM);
Material material_argon = new Material("argon", 0, 0, 0, 14, 0, Units.CM);
materials.add(material_air);
materials.add(material_argon);
}

public Surface(Plane3D plane3d, Point3D refrPoint, Point3D c1, Point3D c2, double accuracy) {
type = Type.PLANEWITHPOINT;
plane = plane3d;
Expand Down
47 changes: 23 additions & 24 deletions ...s/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitter.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -592,7 +592,7 @@ private boolean filter(int k, boolean forward, double annealingFactor) {

double[] K = new double[5];
double V = effectiveVar;
double[] H = mv.H(sVec.x, sVec.y, mVec.surface.measPoint.z(), mVec.surface.wireLine[0]);
double[] H = mv.H(sVec.x, sVec.y, mVec.surface.measPoint.z(), mVec.surface.wireLine[0], mVec.surface.alpha[0]);
Matrix CaInv = this.filterCovMat(H, sVec.CM, V);
if (CaInv != null) {
Matrix5x5.copy(CaInv, cMat);
Expand All @@ -607,7 +607,7 @@ private boolean filter(int k, boolean forward, double annealingFactor) {
}

Point3D point = new Point3D(sVec.x, sVec.y, mVec.surface.measPoint.z());
double h = mv.hDoca(point, mVec.surface.wireLine[0]);
double h = mv.hDoca(point, mVec.surface.wireLine[0], mVec.surface.alpha[0]);

c2 = (effectiveDoca - h) * (effectiveDoca - h) / V;

Expand All @@ -623,7 +623,7 @@ private boolean filter(int k, boolean forward, double annealingFactor) {
+ K[4] * (effectiveDoca - h);

Point3D pointFiltered = new Point3D(x_filt, y_filt, mVec.surface.measPoint.z());
double h0 = mv.hDoca(pointFiltered, mVec.surface.wireLine[0]);
double h0 = mv.hDoca(pointFiltered, mVec.surface.wireLine[0], mVec.surface.alpha[0]);

double residual = effectiveDoca - h0;
updatedWeights_singleHit = daf.calc_updatedWeight_singleHit(residual, annealingFactor);
Expand All @@ -644,7 +644,7 @@ private boolean filter(int k, boolean forward, double annealingFactor) {

double[] K = new double[5];
double V = effectiveVar;
double[] H = mv.H(sVec.x, sVec.y, mVec.surface.measPoint.z(), mVec.surface.wireLine[indexReferenceWire]);
double[] H = mv.H(sVec.x, sVec.y, mVec.surface.measPoint.z(), mVec.surface.wireLine[indexReferenceWire], mVec.surface.alpha[indexReferenceWire]);
Matrix CaInv = this.filterCovMat(H, sVec.CM, V);
if (CaInv != null) {
Matrix5x5.copy(CaInv, cMat);
Expand All @@ -659,7 +659,7 @@ private boolean filter(int k, boolean forward, double annealingFactor) {
}

Point3D point = new Point3D(sVec.x, sVec.y, mVec.surface.measPoint.z());
double h = mv.hDoca(point, mVec.surface.wireLine[indexReferenceWire]);
double h = mv.hDoca(point, mVec.surface.wireLine[indexReferenceWire], mVec.surface.alpha[indexReferenceWire]);

c2 = (effectiveDoca - h) * (effectiveDoca - h) / V;

Expand All @@ -675,8 +675,8 @@ private boolean filter(int k, boolean forward, double annealingFactor) {
+ K[4] * (effectiveDoca - h);

Point3D pointFiltered = new Point3D(x_filt, y_filt, mVec.surface.measPoint.z());
double h0 = mv.hDoca(pointFiltered, mVec.surface.wireLine[0]);
double h1 = mv.hDoca(pointFiltered, mVec.surface.wireLine[1]);
double h0 = mv.hDoca(pointFiltered, mVec.surface.wireLine[0], mVec.surface.alpha[0]);
double h1 = mv.hDoca(pointFiltered, mVec.surface.wireLine[1], mVec.surface.alpha[1]);
double[] residuals = {mVec.surface.doca[0] - h0, mVec.surface.doca[1] - h1};
updatedWeights_doubleHits = daf.calc_updatedWeights_doubleHits(residuals, annealingFactor);
}
Expand Down Expand Up @@ -725,7 +725,7 @@ private boolean filter(int k, boolean forward) {

double[] K = new double[5];
double V = mVec.surface.unc[0] * KFScale;
double[] H = mv.H(sVec.x, sVec.y, mVec.surface.measPoint.z(), mVec.surface.wireLine[0]);
double[] H = mv.H(sVec.x, sVec.y, mVec.surface.measPoint.z(), mVec.surface.wireLine[0], mVec.surface.alpha[0]);
Matrix CaInv = this.filterCovMat(H, sVec.CM, V);
Matrix cMat = new Matrix();
if (CaInv != null) {
Expand All @@ -741,7 +741,7 @@ private boolean filter(int k, boolean forward) {
}

Point3D point = new Point3D(sVec.x, sVec.y, mVec.surface.measPoint.z());
double h = mv.hDoca(point, mVec.surface.wireLine[0]);
double h = mv.hDoca(point, mVec.surface.wireLine[0], mVec.surface.alpha[0]);

double signMeas = 1;
double sign = 1;
Expand Down Expand Up @@ -774,8 +774,7 @@ private boolean filter(int k, boolean forward) {
if (mVec.surface.doca[1] != -99) {
// now filter using the other Hit
V = mVec.surface.unc[1] * KFScale;
H = mv.H(x_filt, y_filt, mVec.surface.measPoint.z(),
mVec.surface.wireLine[1]);
H = mv.H(x_filt, y_filt, mVec.surface.measPoint.z(), mVec.surface.wireLine[1], mVec.surface.alpha[1]);
CaInv = this.filterCovMat(H, cMat, V);
if (CaInv != null) {
for (int i = 0; i < 5; i++) {
Expand All @@ -792,7 +791,7 @@ private boolean filter(int k, boolean forward) {

Point3D point2 = new Point3D(x_filt, y_filt, mVec.surface.measPoint.z());

h = mv.hDoca(point2, mVec.surface.wireLine[1]);
h = mv.hDoca(point2, mVec.surface.wireLine[1], mVec.surface.alpha[1]);

signMeas = Math.signum(mVec.surface.doca[1]);
sign = Math.signum(h);
Expand Down Expand Up @@ -889,7 +888,7 @@ public void calcFinalChisq(int sector, boolean nofilter) {
double V0 = mv.measurements.get(0).surface.unc[0];

Point3D point = new Point3D(svc.x, svc.y, mv.measurements.get(0).surface.measPoint.z());
double h0 = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[0]);
double h0 = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[0], mv.measurements.get(0).surface.alpha[0]);

svc.setProjector(mv.measurements.get(0).surface.wireLine[0].origin().x());
svc.setProjectorDoca(h0);
Expand All @@ -900,7 +899,7 @@ public void calcFinalChisq(int sector, boolean nofilter) {
//USE THE DOUBLE HIT
if (mv.measurements.get(0).surface.doca[1] != -99) {
V0 = mv.measurements.get(0).surface.unc[1];
h0 = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[1]);
h0 = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[1], mv.measurements.get(0).surface.alpha[1]);
res = (mv.measurements.get(0).surface.doca[1] - h0);
chi2 += (mv.measurements.get(0).surface.doca[1] - h0) * (mv.measurements.get(0).surface.doca[1] - h0) / V0;
nRj[mv.measurements.get(0).region - 1] += res * res / mv.measurements.get(0).error;
Expand All @@ -923,7 +922,7 @@ public void calcFinalChisq(int sector, boolean nofilter) {

point = new Point3D(sv.transported(forward).get(k1 + 1).x, sv.transported(forward).get(k1 + 1).y, mv.measurements.get(k1 + 1).surface.measPoint.z());

double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]);
double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0], mv.measurements.get(k1 + 1).surface.alpha[0]);
svc = sv.transported(forward).get(k1 + 1);
path += (forward ? 1 : -1) * svc.deltaPath;
svc.setPathLength(path);
Expand All @@ -936,7 +935,7 @@ public void calcFinalChisq(int sector, boolean nofilter) {
//USE THE DOUBLE HIT
if (mv.measurements.get(k1 + 1).surface.doca[1] != -99) {
V = mv.measurements.get(k1 + 1).surface.unc[1];
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1]);
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1], mv.measurements.get(k1 + 1).surface.alpha[1]);
res = (mv.measurements.get(k1 + 1).surface.doca[1] - h);
chi2 += (mv.measurements.get(k1 + 1).surface.doca[1] - h) * (mv.measurements.get(k1 + 1).surface.doca[1] - h) / V;
nRj[mv.measurements.get(k1 + 1).region - 1] += res * res / V;
Expand Down Expand Up @@ -995,7 +994,7 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) {
double effectiveDoca = daf.get_EffectiveDoca();
double effectiveVar = daf.get_EffectiveVar();

double h = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[0]);
double h = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[0], mv.measurements.get(0).surface.alpha[0]);
double res = (effectiveDoca - h);
chi2 += res*res / effectiveVar;
ndfDAF += daf_weight;
Expand All @@ -1020,20 +1019,20 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) {
double effectiveVar = daf.get_EffectiveVar();
int indexReferenceWire = daf.get_IndexReferenceWire();

double h = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[indexReferenceWire]);
double h = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[indexReferenceWire], mv.measurements.get(0).surface.alpha[indexReferenceWire]);
double res = (effectiveDoca - h);
chi2 += res*res / effectiveVar;
ndfDAF += (daf_weights[0] + daf_weights[1]);

h = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[0]);
h = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[0], mv.measurements.get(0).surface.alpha[0]);
svc.setProjectorDoca(h);
svc.setProjector(mv.measurements.get(0).surface.wireLine[0].origin().x());
svc.setFinalDAFWeight(daf_weights[0]);
svc.setIsDoubleHit(true);
kfStateVecsAlongTrajectory.add(svc);

StateVec svc2 = sv.new StateVec(svc);
h = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[1]);
h = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[1], mv.measurements.get(0).surface.alpha[1]);
svc2.setProjectorDoca(h);
svc2.setProjector(mv.measurements.get(0).surface.wireLine[1].origin().x());
svc2.setFinalDAFWeight(daf_weights[1]);
Expand Down Expand Up @@ -1067,7 +1066,7 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) {
double effectiveDoca = daf.get_EffectiveDoca();
double effectiveVar = daf.get_EffectiveVar();

double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]);
double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0], mv.measurements.get(k1 + 1).surface.alpha[0]);
double res = (effectiveDoca - h);
chi2 += res*res / effectiveVar;
ndfDAF += daf_weight;
Expand All @@ -1092,20 +1091,20 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) {
double effectiveVar = daf.get_EffectiveVar();
int indexReferenceWire = daf.get_IndexReferenceWire();

double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[indexReferenceWire]);
double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[indexReferenceWire], mv.measurements.get(k1 + 1).surface.alpha[indexReferenceWire]);
double res = (effectiveDoca - h);
chi2 += res*res / effectiveVar;
ndfDAF += (daf_weights[0] + daf_weights[1]);

h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]);
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0], mv.measurements.get(k1 + 1).surface.alpha[0]);
svc.setProjectorDoca(h);
svc.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[0].origin().x());
svc.setFinalDAFWeight(daf_weights[0]);
svc.setIsDoubleHit(true);
kfStateVecsAlongTrajectory.add(svc);

StateVec svc2 = sv.new StateVec(svc);
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1]);
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1], mv.measurements.get(k1 + 1).surface.alpha[1]);
svc2.setProjectorDoca(h);
svc2.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[1].origin().x());
svc2.setFinalDAFWeight(daf_weights[1]);
Expand Down
16 changes: 8 additions & 8 deletions ...racking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitterStraight.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -272,7 +272,7 @@ private boolean filter(int k, boolean forward) {

double[] K = new double[5];
double V = mVec.surface.unc[0] * KFScale;
double[] H = mv.H(sVec.x, sVec.y, mVec.surface.measPoint.z(), mVec.surface.wireLine[0]);
double[] H = mv.H(sVec.x, sVec.y, mVec.surface.measPoint.z(), mVec.surface.wireLine[0], mVec.surface.alpha[0]);
Matrix CaInv = this.filterCovMat(H, sVec.CM, V);
Matrix cMat = new Matrix();
if (CaInv != null) {
Expand All @@ -288,7 +288,7 @@ private boolean filter(int k, boolean forward) {
}

Point3D point = new Point3D(sVec.x, sVec.y, mVec.surface.measPoint.z());
double h = mv.hDoca(point, mVec.surface.wireLine[0]);
double h = mv.hDoca(point, mVec.surface.wireLine[0], mVec.surface.alpha[0]);

double signMeas = 1;
double sign = 1;
Expand Down Expand Up @@ -321,7 +321,7 @@ private boolean filter(int k, boolean forward) {
if (mVec.surface.doca[1] != -99) {
// now filter using the other Hit
V = mVec.surface.unc[1] * KFScale;
H = mv.H(x_filt, y_filt, mVec.surface.measPoint.z(), mVec.surface.wireLine[1]);
H = mv.H(x_filt, y_filt, mVec.surface.measPoint.z(), mVec.surface.wireLine[1], mVec.surface.alpha[1]);
CaInv = this.filterCovMat(H, cMat, V);
if (CaInv != null) {
for (int i = 0; i < 5; i++) {
Expand All @@ -338,7 +338,7 @@ private boolean filter(int k, boolean forward) {

Point3D point2 = new Point3D(x_filt, y_filt, mVec.surface.measPoint.z());

h = mv.hDoca(point2, mVec.surface.wireLine[1]);
h = mv.hDoca(point2, mVec.surface.wireLine[1], mVec.surface.alpha[1]);

signMeas = Math.signum(mVec.surface.doca[1]);
sign = Math.signum(h);
Expand Down Expand Up @@ -437,7 +437,7 @@ private void calcFinalChisq(int sector, boolean nofilter) {
double V0 = mv.measurements.get(0).surface.unc[0];

Point3D point = new Point3D(svc.x, svc.y, mv.measurements.get(0).surface.measPoint.z());
double h0 = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[0]);
double h0 = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[0], mv.measurements.get(0).surface.alpha[0]);

svc.setProjector(mv.measurements.get(0).surface.wireLine[0].origin().x());
svc.setProjectorDoca(h0);
Expand All @@ -448,7 +448,7 @@ private void calcFinalChisq(int sector, boolean nofilter) {
//USE THE DOUBLE HIT
if(mv.measurements.get(0).surface.doca[1]!=-99) {
V0 = mv.measurements.get(0).surface.unc[1];
h0 = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[1]);
h0 = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[1], mv.measurements.get(0).surface.alpha[1]);
res = (mv.measurements.get(0).surface.doca[1] - h0);
chi2 += (mv.measurements.get(0).surface.doca[1] - h0) * (mv.measurements.get(0).surface.doca[1] - h0) / V0;
nRj[mv.measurements.get(0).region-1]+=res*res/mv.measurements.get(0).error;
Expand All @@ -472,7 +472,7 @@ private void calcFinalChisq(int sector, boolean nofilter) {

point = new Point3D(sv.transported(forward).get(k1+1).x, sv.transported(forward).get(k1+1).y, mv.measurements.get(k1+1).surface.measPoint.z());

double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]);
double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0], mv.measurements.get(k1 + 1).surface.alpha[0]);
svc = sv.transported(forward).get(k1+1);
path += (forward ? 1 : -1) * svc.deltaPath;
svc.setPathLength(path);
Expand All @@ -485,7 +485,7 @@ private void calcFinalChisq(int sector, boolean nofilter) {
//USE THE DOUBLE HIT
if(mv.measurements.get(k1 + 1).surface.doca[1]!=-99) {
V = mv.measurements.get(k1 + 1).surface.unc[1];
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1]);
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1], mv.measurements.get(k1 + 1).surface.alpha[1]);
res = (mv.measurements.get(k1 + 1).surface.doca[1] - h);
chi2 += (mv.measurements.get(k1 + 1).surface.doca[1] - h) * (mv.measurements.get(k1 + 1).surface.doca[1] - h) / V;
nRj[mv.measurements.get(k1 + 1).region-1]+=res*res/V;
Expand Down
Loading
Loading