Skip to content

Commit 793e70b

Browse files
tongtongcaotongtongcao
committed
correct calculated doca by alpha
1 parent 09b42a5 commit 793e70b

6 files changed

Lines changed: 64 additions & 44 deletions

File tree

‎common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/AMeasVecs.java‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ public double[] dhDoca(int k, StateVec stateVec) {
5050

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

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

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

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

8080
// Return a signed doca for DC
81-
public double hDoca(Point3D point, Line3D wireLine) {
81+
public double hDoca(Point3D point, Line3D wireLine, double alpha) {
8282

8383
Line3D WL = new Line3D();
8484
WL.copy(wireLine);
8585
WL.copy(WL.distance(point));
8686

87-
return WL.length()*Math.signum(WL.direction().x());
87+
return WL.length()*Math.signum(WL.direction().x()) * Math.cos(Math.toRadians(alpha));
8888
}
8989

9090
public double dh(int k, StateVec stateVec) {

‎common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,7 @@ public class Surface implements Comparable<Surface> {
4646
public double[] unc = new double[2];
4747
public double[] doca = new double[2];
4848
public Line3D[] wireLine = new Line3D[2];
49+
public double[] alpha = {0., 0.};
4950
public int region;
5051
public int superlayer;
5152
public int nMeas = 1;
@@ -89,6 +90,22 @@ public Surface(int region, double[] doca, double[] unc, double error, Line3D[] w
8990
materials.add(material_argon);
9091
}
9192

93+
public Surface(int region, double[] doca, double[] unc, double error, Line3D[] wireLine, double[] alpha) {
94+
type = Type.LINEDOCA;
95+
this.region = region;
96+
this.doca = doca;
97+
this.unc = unc;
98+
this.error = error;
99+
this.wireLine = wireLine;
100+
this.alpha = alpha;
101+
Point3D point = new Point3D(0, 0, wireLine[0].origin().z());
102+
this.measPoint = point;
103+
Material material_air = new Material("air", 0, 0, 0, 30400, 0, Units.CM);
104+
Material material_argon = new Material("argon", 0, 0, 0, 14, 0, Units.CM);
105+
materials.add(material_air);
106+
materials.add(material_argon);
107+
}
108+
92109
public Surface(Plane3D plane3d, Point3D refrPoint, Point3D c1, Point3D c2, double accuracy) {
93110
type = Type.PLANEWITHPOINT;
94111
plane = plane3d;

‎common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitter.java‎

Lines changed: 23 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -592,7 +592,7 @@ private boolean filter(int k, boolean forward, double annealingFactor) {
592592

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

795-
h = mv.hDoca(point2, mVec.surface.wireLine[1]);
794+
h = mv.hDoca(point2, mVec.surface.wireLine[1], mVec.surface.alpha[1]);
796795

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

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

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

924923
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());
925924

926-
double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]);
925+
double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0], mv.measurements.get(k1 + 1).surface.alpha[0]);
927926
svc = sv.transported(forward).get(k1 + 1);
928927
path += (forward ? 1 : -1) * svc.deltaPath;
929928
svc.setPathLength(path);
@@ -936,7 +935,7 @@ public void calcFinalChisq(int sector, boolean nofilter) {
936935
//USE THE DOUBLE HIT
937936
if (mv.measurements.get(k1 + 1).surface.doca[1] != -99) {
938937
V = mv.measurements.get(k1 + 1).surface.unc[1];
939-
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1]);
938+
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1], mv.measurements.get(k1 + 1).surface.alpha[1]);
940939
res = (mv.measurements.get(k1 + 1).surface.doca[1] - h);
941940
chi2 += (mv.measurements.get(k1 + 1).surface.doca[1] - h) * (mv.measurements.get(k1 + 1).surface.doca[1] - h) / V;
942941
nRj[mv.measurements.get(k1 + 1).region - 1] += res * res / V;
@@ -995,7 +994,7 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) {
995994
double effectiveDoca = daf.get_EffectiveDoca();
996995
double effectiveVar = daf.get_EffectiveVar();
997996

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

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

1028-
h = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[0]);
1027+
h = mv.hDoca(point, mv.measurements.get(0).surface.wireLine[0], mv.measurements.get(0).surface.alpha[0]);
10291028
svc.setProjectorDoca(h);
10301029
svc.setProjector(mv.measurements.get(0).surface.wireLine[0].origin().x());
10311030
svc.setFinalDAFWeight(daf_weights[0]);
10321031
svc.setIsDoubleHit(true);
10331032
kfStateVecsAlongTrajectory.add(svc);
10341033

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

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

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

1100-
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]);
1099+
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0], mv.measurements.get(k1 + 1).surface.alpha[0]);
11011100
svc.setProjectorDoca(h);
11021101
svc.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[0].origin().x());
11031102
svc.setFinalDAFWeight(daf_weights[0]);
11041103
svc.setIsDoubleHit(true);
11051104
kfStateVecsAlongTrajectory.add(svc);
11061105

11071106
StateVec svc2 = sv.new StateVec(svc);
1108-
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1]);
1107+
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1], mv.measurements.get(k1 + 1).surface.alpha[1]);
11091108
svc2.setProjectorDoca(h);
11101109
svc2.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[1].origin().x());
11111110
svc2.setFinalDAFWeight(daf_weights[1]);

‎common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitterStraight.java‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -272,7 +272,7 @@ private boolean filter(int k, boolean forward) {
272272

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

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

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

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

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

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

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

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

473473
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());
474474

475-
double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]);
475+
double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0], mv.measurements.get(k1 + 1).surface.alpha[0]);
476476
svc = sv.transported(forward).get(k1+1);
477477
path += (forward ? 1 : -1) * svc.deltaPath;
478478
svc.setPathLength(path);
@@ -485,7 +485,7 @@ private void calcFinalChisq(int sector, boolean nofilter) {
485485
//USE THE DOUBLE HIT
486486
if(mv.measurements.get(k1 + 1).surface.doca[1]!=-99) {
487487
V = mv.measurements.get(k1 + 1).surface.unc[1];
488-
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1]);
488+
h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[1], mv.measurements.get(k1 + 1).surface.alpha[1]);
489489
res = (mv.measurements.get(k1 + 1).surface.doca[1] - h);
490490
chi2 += (mv.measurements.get(k1 + 1).surface.doca[1] - h) * (mv.measurements.get(k1 + 1).surface.doca[1] - h) / V;
491491
nRj[mv.measurements.get(k1 + 1).region-1]+=res*res/V;

0 commit comments

Comments
 (0)