diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitter.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitter.java index c723fb2a3..cfdce11c3 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitter.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitter.java @@ -708,7 +708,9 @@ private boolean filter(int k, boolean forward, double annealingFactor) { return false; } } - + + // Since no vertex inforamtion, the starting point for path length is the final point at the last layer. + // After vertex information is obtained, transition for the starting point from the final point to vertex will be taken. private boolean filter(int k, boolean forward) { StateVec sVec = sv.transported(forward).get(k); org.jlab.clas.tracking.kalmanfilter.AMeasVecs.MeasVec mVec = mv.measurements.get(k); @@ -857,7 +859,9 @@ public Matrix filterCovMat(double[] H, Matrix Ci, double V) { private void calcFinalChisq(int sector) { calcFinalChisq(sector, false); } - + + // Since no vertex inforamtion, the starting point for path length is the final point at the last layer. + // After vertex information is obtained, transition for the starting point from the final point to vertex will be taken. private void calcFinalChisq(int sector, boolean nofilter) { int k = svzLength - 1; this.chi2 = 0; @@ -880,9 +884,9 @@ private void calcFinalChisq(int sector, boolean nofilter) { sv.transport(sector, k, 0, sVec, mv, this.getSwimmer(), forward); StateVec svc = sv.transported(forward).get(0); - path += svc.deltaPath; + path += (forward ? 1 : -1) * svc.deltaPath; svc.setPathLength(path); - + double V0 = mv.measurements.get(0).surface.unc[0]; Point3D point = new Point3D(svc.x, svc.y, mv.measurements.get(0).surface.measPoint.z()); @@ -922,7 +926,7 @@ private void calcFinalChisq(int sector, boolean nofilter) { double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]); svc = sv.transported(forward).get(k1 + 1); - path += svc.deltaPath; + path += (forward ? 1 : -1) * svc.deltaPath; svc.setPathLength(path); svc.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[0].origin().x()); svc.setProjectorDoca(h); @@ -975,7 +979,7 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) { sv.transport(sector, k, 0, sVec, mv, this.getSwimmer(), forward); StateVec svc = sv.transported(forward).get(0); - path += svc.deltaPath; + path += (forward ? 1 : -1) * svc.deltaPath; svc.setPathLength(path); Point3D point = new Point3D(svc.x, svc.y, mv.measurements.get(0).surface.measPoint.z()); @@ -1047,7 +1051,7 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) { } svc = sv.transported(forward).get(k1 + 1); - path += svc.deltaPath; + path += (forward ? 1 : -1) * svc.deltaPath; svc.setPathLength(path); 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()); @@ -1116,6 +1120,10 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) { public Matrix propagateToVtx(int sector, double Zf) { return sv.transport(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer()); } + + public double getDeltaPathToVtx(int sector, double Zf) { + return sv.getDeltaPath(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer()); + } //Todo: apply the common funciton to replace current function above @Override diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitterStraight.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitterStraight.java index 232bcaa7d..fb563fcc8 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitterStraight.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitterStraight.java @@ -406,6 +406,8 @@ private void calcFinalChisq(int sector) { calcFinalChisq(sector, false); } + // Since no vertex inforamtion, the starting point for path length is the final point at the last layer. + // After vertex information is obtained, transition for the starting point from the final point to vertex will be taken. private void calcFinalChisq(int sector, boolean nofilter) { int k = svzLength - 1; this.chi2 = 0; @@ -426,11 +428,11 @@ private void calcFinalChisq(int sector, boolean nofilter) { kfStateVecsAlongTrajectory = new ArrayList<>(); if (sVec != null && sVec.CM != null) { - boolean forward = false; + boolean forward = false; sv.transport(sector, k, 0, sVec, mv, this.getSwimmer(), forward); StateVec svc = sv.transported(forward).get(0); - path += svc.deltaPath; + path += (forward ? 1 : -1) * svc.deltaPath; svc.setPathLength(path); double V0 = mv.measurements.get(0).surface.unc[0]; @@ -473,7 +475,7 @@ private void calcFinalChisq(int sector, boolean nofilter) { double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]); svc = sv.transported(forward).get(k1+1); - path += svc.deltaPath; + path += (forward ? 1 : -1) * svc.deltaPath; svc.setPathLength(path); svc.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[0].origin().x()); svc.setProjectorDoca(h); @@ -504,6 +506,10 @@ public Matrix propagateToVtx(int sector, double Zf) { return sv.transport(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer()); } + public double getDeltaPathToVtx(int sector, double Zf) { + return sv.getDeltaPath(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer()); + } + @Override public void runFitter(AStateVecs sv, AMeasVecs mv) { throw new UnsupportedOperationException("Not supported yet."); diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java index f69dd0bfc..269e966c2 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java @@ -53,7 +53,7 @@ public void initFromHB(StateVec initSV, double beta) { this.trackTrajS.clear(); this.trackTrajT.put(0, new StateVec(initSV)); } - + /** * * @param sector @@ -150,6 +150,104 @@ public Matrix transport(int sector, int i, double Zf, StateVec iVec, AMeasVecs m return fVec.CM; + } + + /** + * + * @param sector + * @param i initial state vector index + * @param Zf + * @param iVec state vector at the initial index + * @param mv measurements + */ + public double getDeltaPath(int sector, int i, double Zf, StateVec iVec, AMeasVecs mv, Swim swimmer) { // s = signed step-size + + double stepSize = 1.0; + StateVec fVec = new StateVec(0); + fVec.x = iVec.x; + fVec.y = iVec.y; + fVec.z = iVec.z; + fVec.tx = iVec.tx; + fVec.ty = iVec.ty; + fVec.Q = iVec.Q; + fVec.B = iVec.B; + Matrix5x5.copy(iVec.CM, fVec.CM); + + double s = 0; + double zInit = mv.measurements.get(i).surface.measPoint.z(); + double BatMeas = iVec.B; + + double z = zInit; + + while (Math.signum(Zf - zInit) * z < Math.signum(Zf - zInit) * Zf) { + z = fVec.z; + if (z == Zf) { + break; + } + + double x = fVec.x; + double y = fVec.y; + double tx = fVec.tx; + double ty = fVec.ty; + double Q = fVec.Q; + double dPath = fVec.deltaPath; + Matrix cMat = new Matrix(); + Matrix5x5.copy(fVec.CM, cMat); + s = Math.signum(Zf - zInit) * stepSize; + + // LOGGER.log(Level.FINE, " from "+(float)Z[i]+" to "+(float)Z[f]+" at "+(float)z+" By is "+bf[1]+" B is "+Math.sqrt(bf[0]*bf[0]+bf[1]*bf[1]+bf[2]*bf[2])/Bmax+" stepSize is "+s); + if (Math.signum(Zf - zInit) * (z + s) > Math.signum(Zf - zInit) * Zf) { + s = Math.signum(Zf - zInit) * Math.abs(Zf - z); + } + + //rk.RK4transport(sector, Q, x, y, z, tx, ty, s, swimmer, cMat, fVec, dPath); + rk.RK4transport(sector, s, swimmer, cMat, fVec); + + // Q process noise matrix estimate + double p = Math.abs(1. / iVec.Q); + + double X0 = this.getX0(mv.measurements.get(i).surface, z, Z); + double t_ov_X0 = Math.abs(s) / X0;//path length in radiation length units = t/X0 [true path length/ X0] ; Ar radiation length = 14 cm + + double beta = this.beta; + if (beta > 1.0 || beta <= 0) { + beta = 1.0; + } + + double sctRMS = 0; + + ////// Todo: Modify multi-scattering or remove it; After update, some parameters, like iteration termintion chonditions, may need to be updated. + // Speed of light should be 1 + // From one measurement site to another, F and Q should be calculated separaetely with multiple steps; and then C' = FTCF + Q + if (Math.abs(s) > 0) { + sctRMS = ((0.0136) / (beta * PhysicsConstants.speedOfLight() * p)) * Math.sqrt(t_ov_X0) + * (1 + 0.038 * Math.log(t_ov_X0)); + } + + double cov_txtx = (1 + tx * tx) * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; + double cov_tyty = (1 + ty * ty) * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; + double cov_txty = tx * ty * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; + + fMS.set( + 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, + 0, 0, cov_txtx, cov_txty, 0, + 0, 0, cov_txty, cov_tyty, 0, + 0, 0, 0, 0, 0 + ); + + Matrix5x5.copy(fVec.CM, copyMatrix); + Matrix5x5.add(copyMatrix, fMS, fVec.CM); + + if (Math.abs(fVec.B - BatMeas) < 0.0001) { + stepSize *= 2; + } + + BatMeas = fVec.B; + } + + return fVec.deltaPath; + } /** 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 6d0256005..6e8ca4cb7 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 @@ -907,8 +907,6 @@ private List findStraightTracks(CrossList crossList, DCGeant4Factory DcDe continue; } - List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef); - StateVec fn = new StateVec(); if (!kFZRef.setFitFailed && kFZRef.finalStateVec != null) { fn.set(kFZRef.finalStateVec.x, kFZRef.finalStateVec.y, kFZRef.finalStateVec.tx, kFZRef.finalStateVec.ty); @@ -927,6 +925,11 @@ private List findStraightTracks(CrossList crossList, DCGeant4Factory DcDe cand.set_FitConvergenceStatus(kFZRef.ConvStatus); cand.set_Id(cands.size() + 1); cand.set_CovMat(kFZRef.finalStateVec.CM); + + Point3D VTCS = cand.get(cand.size()-1).getCoordsInTiltedSector(cand.get_Vtx0().x(), cand.get_Vtx0().y(), cand.get_Vtx0().z()); + double deltaPathToVtx = kFZRef.getDeltaPathToVtx(cand.get(cand.size()-1).get_Sector(), VTCS.z()); + + List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef, deltaPathToVtx); cand.setStateVecs(kfStateVecsAlongTrajectory); // add candidate to list of tracks @@ -1071,7 +1074,6 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete kFZRef.init(measSurfaces, initSV); kFZRef.runFitter(); - List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef); if (kFZRef.finalStateVec == null) { continue; @@ -1093,15 +1095,21 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete cand.set_FitNDF(kFZRef.NDF); cand.set_FitConvergenceStatus(kFZRef.ConvStatus); - cand.set_CovMat(kFZRef.finalStateVec.CM); - cand.setStateVecs(kfStateVecsAlongTrajectory); - + cand.set_CovMat(kFZRef.finalStateVec.CM); + cand.setFinalStateVec(fitStateVec); cand.set_Id(cands.size() + 1); this.setTrackPars(cand, traj, trjFind, fitStateVec, fitStateVec.getZ(), DcDetector, dcSwim); + + Point3D VTCS = cand.get(cand.size()-1).getCoordsInTiltedSector(cand.get_Vtx0().x(), cand.get_Vtx0().y(), cand.get_Vtx0().z()); + double deltaPathToVtx = kFZRef.getDeltaPathToVtx(cand.get(cand.size()-1).get_Sector(), VTCS.z()); + + List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef, deltaPathToVtx); + cand.setStateVecs(kfStateVecsAlongTrajectory); + // add candidate to list of tracks if (cand.fit_Successful = true) { cands.add(cand); @@ -1117,7 +1125,7 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete return cands; } - public List setKFStateVecsAlongTrajectory(KFitter kFZRef) { + public List setKFStateVecsAlongTrajectory(KFitter kFZRef, double deltaPathToVtx) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -1125,7 +1133,7 @@ public List setKFStateVecsAlongTrajectory(K org.jlab.rec.dc.trajectory.StateVec sv = new org.jlab.rec.dc.trajectory.StateVec(svc.x, svc.y, svc.tx, svc.ty); sv.setZ(svc.z); sv.setB(svc.B); - sv.setPathLength(svc.getPathLength()); + sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex sv.setProjector(svc.getProjector()); sv.setProjectorDoca(svc.getProjectorDoca()); kfStateVecsAlongTrajectory.add(sv); @@ -1134,7 +1142,7 @@ public List setKFStateVecsAlongTrajectory(K return kfStateVecsAlongTrajectory; } - public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef) { + public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef, double deltaPathToVtx) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -1142,7 +1150,7 @@ public List setKFStateVecsAlongTrajectory(K org.jlab.rec.dc.trajectory.StateVec sv = new org.jlab.rec.dc.trajectory.StateVec(svc.x, svc.y, svc.tx, svc.ty); sv.setZ(svc.z); sv.setB(svc.B); - sv.setPathLength(svc.getPathLength()); + sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex sv.setProjector(svc.getProjector()); sv.setProjectorDoca(svc.getProjectorDoca()); kfStateVecsAlongTrajectory.add(sv); diff --git a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java index 2bd9d1e9f..851dc8a16 100644 --- a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java +++ b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java @@ -235,7 +235,6 @@ public boolean processDataEvent(DataEvent event) { getInitState(TrackArray1, measSurfaces.get(0).measPoint.z(), initSV, kFZRef, dcSwim, new float[3]); kFZRef.initFromHB(measSurfaces, initSV, TrackArray1.get(0).get(0).get(0).get_Beta()); kFZRef.runFitter(); - List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef); StateVec fn = new StateVec(); if (kFZRef.setFitFailed==false && kFZRef.finalStateVec!=null) { @@ -253,15 +252,18 @@ public boolean processDataEvent(DataEvent event) { TrackArray1.set_FitChi2(kFZRef.chi2); TrackArray1.set_FitNDF(kFZRef.NDF); - TrackArray1.setStateVecs(kfStateVecsAlongTrajectory); TrackArray1.set_FitConvergenceStatus(kFZRef.ConvStatus); if (TrackArray1.get_Vtx0().toVector3D().mag() > 500) { continue; } - + // get CovMat at vertex - Point3D VTCS = crosses.get(0).getCoordsInSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z()); + Point3D VTCS = crosses.get(0).getCoordsInTiltedSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z()); TrackArray1.set_CovMat(kFZRef.propagateToVtx(crosses.get(0).get_Sector(), VTCS.z())); + + double deltaPathToVtx = kFZRef.getDeltaPathToVtx(TrackArray1.get(TrackArray1.size()-1).get_Sector(), VTCS.z()); + List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef, deltaPathToVtx); + TrackArray1.setStateVecs(kfStateVecsAlongTrajectory); if (TrackArray1.isGood()) { trkcands.add(TrackArray1); @@ -277,8 +279,6 @@ public boolean processDataEvent(DataEvent event) { kFZRef.initFromHB(measSurfaces, initSV, TrackArray1.get(0).get(0).get(0).get_Beta(), useDAF); kFZRef.runFitter(useDAF); - List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef); - StateVec fn = new StateVec(); if (kFZRef.setFitFailed==false && kFZRef.finalStateVec!=null) { // set the state vector at the last measurement site @@ -296,15 +296,18 @@ public boolean processDataEvent(DataEvent event) { TrackArray1.set_FitChi2(kFZRef.chi2); TrackArray1.set_FitNDF(kFZRef.NDF); TrackArray1.set_NDFDAF(kFZRef.getNDFDAF()); - TrackArray1.setStateVecs(kfStateVecsAlongTrajectory); TrackArray1.set_FitConvergenceStatus(kFZRef.ConvStatus); if (TrackArray1.get_Vtx0().toVector3D().mag() > 500) { continue; } // get CovMat at vertex - Point3D VTCS = crosses.get(0).getCoordsInSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z()); + Point3D VTCS = crosses.get(0).getCoordsInTiltedSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z()); TrackArray1.set_CovMat(kFZRef.propagateToVtx(crosses.get(0).get_Sector(), VTCS.z())); + + double deltaPathToVtx = kFZRef.getDeltaPathToVtx(TrackArray1.get(TrackArray1.size()-1).get_Sector(), VTCS.z()); + List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef, deltaPathToVtx); + TrackArray1.setStateVecs(kfStateVecsAlongTrajectory); if (TrackArray1.isGood()) { trkcands.add(TrackArray1); @@ -361,7 +364,7 @@ public boolean processDataEvent(DataEvent event) { return true; } - public List setKFStateVecsAlongTrajectory(KFitter kFZRef) { + public List setKFStateVecsAlongTrajectory(KFitter kFZRef, double deltaPathToVtx) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -369,7 +372,7 @@ public List setKFStateVecsAlongTrajectory(K org.jlab.rec.dc.trajectory.StateVec sv = new org.jlab.rec.dc.trajectory.StateVec(svc.x, svc.y, svc.tx, svc.ty); sv.setZ(svc.z); sv.setB(svc.B); - sv.setPathLength(svc.getPathLength()); + sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex sv.setProjector(svc.getProjector()); sv.setProjectorDoca(svc.getProjectorDoca()); sv.setDAFWeight(svc.getFinalDAFWeight()); @@ -380,7 +383,7 @@ public List setKFStateVecsAlongTrajectory(K return kfStateVecsAlongTrajectory; } - public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef) { + public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef, double deltaPathToVtx) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -388,7 +391,7 @@ public List setKFStateVecsAlongTrajectory(K org.jlab.rec.dc.trajectory.StateVec sv = new org.jlab.rec.dc.trajectory.StateVec(svc.x, svc.y, svc.tx, svc.ty); sv.setZ(svc.z); sv.setB(svc.B); - sv.setPathLength(svc.getPathLength()); + sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex sv.setProjector(svc.getProjector()); sv.setProjectorDoca(svc.getProjectorDoca()); kfStateVecsAlongTrajectory.add(sv);